In an email interchange with Mike Hansen, I was reminded of this bug...
> (%i5) ceiling(factorial(50)/exp(1));
> (%o5) 11188719610782480421414879249141773426630319613740326700720324608
>
> The correct result should be
> 11188719610782480504630258070757734324011354208865721592720336801 .
> What is the best way to go about fixing this issue?
from dec,2007.
Here's a fix. Use this program, for ceiling computed using bigfloat.
ceilbf(x,q):=block([?fpprec:
q+ceiling(log(abs(x))/log(2))],ceiling(bfloat(x)));
where q is enough precision so that ceiling can work to the nearest integer.
I'm guessing that the number of bits in the answer + 8 more for slop should
be
a good deal better. That is, q>=8.
some code should perhaps be copied out of ceiling to make sure that x
evaluates to a number.
Also note that ?fpprec is like fpprec, except in bits, generously rounded
up, instead of decimal.
setting fpprec:1000 sets ?fpprec to 3324.
the program could be written
ceilbf10(x,q):=block[fpprec:q+ceil(log(abs(x))/log(10))] ...
alternatively, where q is extra decimal digits.
perhaps ceiling should be described as a heuristic using double-precision,
and certainly not to be
trusted outside [-10^16,10^16].
ceilbf depends on Maxima being able to compute the result to the nearest
integer using arithmetic that carries
q extra bits/digits.
You could try ceilbf twice, with q=8 and q=16; see if you get the same
result :)
RJF