reyssat wrote:
> Hi,
>
> I have some questions and remarks about precision of floats and bigfloats.
>
> I construct
>
> f(i):=block([u], u:sqrt(10^(2*i)+1)/10^i, print(bfloat(u)),print(float(u)));
>
> then :
> (%i3) f(2);
> 1.000049998750063b0
> 1.000049998750062
>
> **** this is ok, just slightly strange that float gives a better result
> than bfloat (the right answer is 1.000049998750062496..., so should be
> rounded to ...0062 and not ...0063).
>
I think bfloats try to be accurate, but there's no guarantee that it is
accurate to the last bit, especially for the special functions.
> By the way, I cannot find in the documentation (at least in chapter 10,
> floating point) the definition and precision of a float.
>
The precision of a float depends on the double-float type of the
underlying lisp. But I think everyone uses IEEE double-floats, so the
precision is well defined.
>
> (%i5) f(200);
> 1.0b0
> Maxima encountered a Lisp error:
>
> Error in PROGN [or a callee]: Can't print a non-number.
>
> Automatically continuing.
> To reenable the Lisp debugger set *debugger-hook* to nil.
>
> **** bfloat is ok. But float cannot calculate this since some
> intermediate result is too big (once again, is the limit for the
> exponent given in the documentation ?), although the result itself is a
>
No I don't think the limit is in the documentation. But the largest
IEEE number is 2^1024 minus one bit. There is some package that has
this information, but I can't remember the name of the package.
> perfectly nice number. Would it be possible to improve the float
> algorithm to avoid this problem ?
>
Don't think that's really possible in general. It's up to you to
arrange the computations to reduce overflow. But in this particular
case, we could probably make maxima a little smarter when evaluating
float(sqrt(10^2*i+1)). We could compute x=isqrt(10^(2*i)+1) and
evaluate sqrt(float((10^(2*i)+1)/x)) instead. That would allow larger
numbers to be used, but eventually there would be overflow.
>
> (%i6) f(1000);
> 1.0b0
> 0.0
>
> **** this time, no error message, just a stupid result from float ! is
> it a bug, or just a feature ? bfloat is ok.
>
This is clearly a bug. I don't get this with cmucl, which complains
about 10<many 0's>01 is not a double-float.
I notice that maxima is pretty slow in computing sqrt(10^2000+1). I
guess it's taking a long time looking for any squares that could be
factored out of the sqrt.
Ray