[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: A problem with range objects and floating point numbers
From: |
Daniel J Sebald |
Subject: |
Re: A problem with range objects and floating point numbers |
Date: |
Fri, 26 Sep 2014 12:53:48 -0500 |
User-agent: |
Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv:1.9.2.24) Gecko/20111108 Fedora/3.1.16-1.fc14 Thunderbird/3.1.16 |
On 09/26/2014 02:43 AM, s.jo wrote:
Daniel Sebald wrote
[snip]
To be factorable doesn't necessarily require that the floating point
step size be an exact number system equivalent, i.e., that .1 equal
exactly 1/10. It just requires the machine arithmetic to be as expected
when producing an operation result, i.e.,
octave-cli:1> rem(-2,.1) == 0
ans = 1
octave-cli:2> rem(0,.1) == 0
ans = 1
octave-cli:3> [-20:1:0]*.1
ans =
Columns 1 through 5:
-2.00000 -1.90000 -1.80000 -1.70000 -1.60000
Columns 6 through 10:
-1.50000 -1.40000 -1.30000 -1.20000 -1.10000
Columns 11 through 15:
-1.00000 -0.90000 -0.80000 -0.70000 -0.60000
Columns 16 through 20:
-0.50000 -0.40000 -0.30000 -0.20000 -0.10000
Column 21:
0.00000
Jo, what do you get for the above scripts?
Dan
Dan and others,
Sorry for late reply.
Your integerized stepping seems to be a decimation process.
You used step-size 1. Any bigger integers such as 10 or 10^17 are allowed in
the decimation process.
By decimation, I take it you mean somehow converting the number to true
decimal representation, similar to the decimal arithmetic toolbox Oliver
alluded to. I don't think that's the case.
Comparing the results of colon operators with floating-point numbers,
Matlab seems to use such decimation process, but Octave does not.
I agree that this decimation feature of colon operator is far from any IEEE
standard.
Also I point out that the results of linspace function and colon operator in
Matlab are different.
This may imply decimation processes are different.
By inspection, the linspace result is much close to decimal number.
I guess that the decimation process of Matlab's colon operator is simpler
for speed.
No, I doubt that. The linspace and colon operator are probably two
slightly different algorithms. For linspace the spacing is exactly an
integer, by definition. That's not necessarily so for the colon
operator where the spacing can result in a fraction of the step size at
the last interval.
Looking at the code a bit, I see that the range (colon) operator does
attempt to use integer multiplication as much as possible, e.g., from
Range::matrix_value:
cache.resize (1, rng_nelem);
double b = rng_base;
double increment = rng_inc;
if (rng_nelem > 0)
{
// The first element must always be *exactly* the base.
// E.g, -0 would otherwise become +0 in the loop (-0 +
0*increment).
cache(0) = b;
for (octave_idx_type i = 1; i < rng_nelem; i++)
cache(i) = b + i * increment;
}
So, the accumulation of addition errors doesn't seem to be an issue here.
So far, I learned that integer range can be generated 'safely' with colon
operator.
Regarding to floating-point number range, we may want to use the function
"linspace" for the best decimation.
I don't know who is responsible for the colon operator source code.
But I think that if there is a document for comparing colon operator with
linspace function,
that will be helpful.
There is another question: I found that (-20:1:1)*.1 and [-20:1:1]*0.1 have
different results.
The range with [...] is better decimated than the case of (...).
Interesting. You may have it right there. And in fact, it appears the
conversation in the bug report
https://savannah.gnu.org/bugs/?42589
turns to this issue. John describes how a "Range" class retains the
description of the range rather than converting it to a matrix class.
That actually may be a good idea because it could possibly allow
retaining double precision if needed somehow (e.g., machines where the
single precision is markedly different from "double"). Whereas
converting the range to a matrix might make the result single precision
depending on compiler settings or whatnot. I'm just speaking in the
abstract without any example.
So, in your example, the difference is (-20:1:1) is a Range description
(base, increment, limit), where [-20:1:1] is a matrix_value
implementation of that range (all the actual 21 values).
Let's try to discern which of these might be the more accurate result:
t_1 = (-20:1:1)*.1;
t_2 = [-20:1:1]*.1;
t_f = [-2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 ...
-1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1];
t_d = double(zeros(1,21));
t_d(1) = -2; t_d(2) = -1.9; t_d(3) = -1.8; t_d(4) = -1.7; t_d(5) = -1.6;
t_d(6) = -1.5; t_d(7) = -1.4; t_d(8) = -1.3; t_d(9) = -1.2; t_d(10) = -1.1;
t_d(11) = -1; t_d(12) = -0.9; t_d(13) = -0.8; t_d(14) = -0.7; t_d(15) = -0.6;
t_d(16) = -0.5; t_d(17) = -0.4; t_d(18) = -0.3; t_d(19) = -0.2; t_d(20) = -0.1;
t_d(21) = 0; t_d(22) = 0.1;
max(abs(t_d - t_f))
ans = 0
norm(abs(t_d - t_f))
ans = 0
max(abs(t_d - t_1))
ans = 2.2204e-16
norm(abs(t_d - t_1))
ans = 4.3088e-16
max(abs(t_d - t_2))
ans = 2.2204e-16
norm(abs(t_d - t_2))
ans = 4.7429e-16
max(abs(t_1 - t_2))
ans = 2.2204e-16
[abs(t_d-t_1)' abs(t_d-t_2)']
ans =
0.0000e+00 0.0000e+00
0.0000e+00 2.2204e-16
0.0000e+00 0.0000e+00
0.0000e+00 2.2204e-16
0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00
0.0000e+00 2.2204e-16
2.2204e-16 0.0000e+00
0.0000e+00 2.2204e-16
0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00
1.1102e-16 0.0000e+00
2.2204e-16 0.0000e+00
0.0000e+00 1.1102e-16
1.1102e-16 1.1102e-16
0.0000e+00 0.0000e+00
1.1102e-16 0.0000e+00
1.6653e-16 5.5511e-17
5.5511e-17 0.0000e+00
1.3878e-16 0.0000e+00
0.0000e+00 0.0000e+00
0.0000e+00 0.0000e+00
From the norm(abs(t_d - t_f)) = 0 result, I'd say my machine/compiler
is using double precision by default. Also, there is a difference
between (A:1:B)*d and [A:1:B]*d but no implementation seems better than
the other. If anything (A:1:B)*d is marginally better than [A:1:B]*d,
but statistically speaking the norm difference of 4.3405e-17 is probably
not conclusive. However, if I were to guess, I'd say (A:1:B)*d retains
accuracy better because it doesn't cast to a IEEE float representation
and then multiply. The computations may be all done with more bits in
the ALU. But we're talking the very last bit here, I believe.
More interestingly, if I transpose the range array such as (-20:1:1)'*.1,
then I have the same result in [-20:1:1]'*.1.
The transpose likely forces the Range description into a matrix_value
implementation.
Is this a bug? Did someone say that there is a patch for this?
I am useing Ocatve 3.8.2 on cygwin.
OK, Cygwin. So there's where there could be some differences in
compiler settings that make your computations single precision and
therefore showing larger errors in the arithmetic than what I'm seeing.
Try the script lines above on your machine and see what you get. You
probably should have the latest code with the Range patch applied. In
any case, though, if you can identify a single/double precision issue,
it might suggest that the Octave code isn't retaining precision...or it
could simply be that on your system Matlab is using double precision at
compilation whereas Octave is not.
Dan
- Re: A problem with range objects and floating point numbers, (continued)
- Re: A problem with range objects and floating point numbers, Oliver Heimlich, 2014/09/24
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/24
- Re: A problem with range objects and floating point numbers, Oliver Heimlich, 2014/09/25
- Re: A problem with range objects and floating point numbers, Markus, 2014/09/25
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/25
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/25
Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/24
- Re: A problem with range objects and floating point numbers, s.jo, 2014/09/26
- Re: A problem with range objects and floating point numbers,
Daniel J Sebald <=
- Re: A problem with range objects and floating point numbers, s.jo, 2014/09/28
- Re: A problem with range objects and floating point numbers, s.jo, 2014/09/29
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/29
- Re: A problem with range objects and floating point numbers, s.jo, 2014/09/29
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/29
- Re: A problem with range objects and floating point numbers, s.jo, 2014/09/29
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/29
- Re: A problem with range objects and floating point numbers, John W. Eaton, 2014/09/29
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/29
- Re: A problem with range objects and floating point numbers, Daniel J Sebald, 2014/09/29