Menu
▾
▴
#4659 mod gives incorrect results for bignums
As I understand, mod works for floats:
(%i4) mod(10,3.1);
(%o4) 0.6999999999999993
However,
(%i1) X: 10^100;
(%o1) 100000000000000000000000000000000000000000000000000000000000000000000000\
00000000000000000000000000000
(%i2) mod(X,2.1);
(%o2) 0.0
(%i3) mod(X*10,21);
(%o3) 19
So, the first result, from mod(X,2.1) should have been 1.9.
This is clearly because of the IEEE representation of Y. Perhaps working with the mantissa/exponent in order to get the proper result would work?
- Extract mantissa m, exponent e from Y
- Divide X by 2^e
- convert mantissa to p/q, p = round(m·2^k), q = 2^k, with k being enough bits to represent m
- t = mod(X' , q) * mod(p, q)
- result = (r / 2^e) + (t / q)
Not sure, but seems like this would work. Maybe check if the result is >= Y, and subtract Y once, if necessary (due to rounding)
Floating point arithmetic is "contagious". That is, mod(10^100,2.1) is equivalent to mod(float(10^100),2.1). But neither 10^100 nor 2.1 is represented exactly as a floating point number:
rationalize(float(10^100))-10^100 => ~10^83
rationalize(2.1)-21/10 => ~10^-16.
To get an exact answer, use mod(10^100,21/10) => 19/10.
Note that rationalize converts a floating point number to the exact rational it represents, e.g., rationalize(0.1) => 3602879701896397/36028797018963968 and rationalize(0.1)-1/10 => 1/180143985094819840 = ~10^-16.
Yes - but there probably are ways to get more precise results (that is what I was trying with that reasoning up there). I'll see if I get a MWE for a possible implementation
I suspect that by the time you've come up with some algorithm, you've probably duplicated most of the work of what Stavros is suggesting by using rationalize.
I personally think the answer should be 0. Why? Because
floor(float(10^100)/2.1)*2.1is actually equal tofloat(10^100). There just aren't enough bits in a float. However, you could use bfloats for everything with a large enough value for fpprec.Last edit: Raymond Toy 6 days ago
I had forgotten about rationalize.
What I think is that Maxima could be "gentle" enough to detect a bignum in the first argument, then do rationalize on the second, get the result, and turn it back into a float. You'd get 19.0, which is a mathematically correct answer (although not verifiable through IEEE arithmetic, as you mentioned).
Last edit: Jeronimo Pellegrini 6 days ago
For some reason, sourceforge is unresponsive although other sites are fine
(I'm on airplane wifi).
Anyway, I think what Pellegrini is looking for is
float(mod(10^100,rat(2.1))).
Having mod automagically use bfloat or rat arithmetic for bignum args
would suggest that we do the same thing for all float arithmetic, which
we're clearly not going to do.
We could introduce an "exact float" type, say 123.45x01, which would denote
the rational 1234/10. But what is sqrt(123.45x01) or sin(3.14x0)?
Hello!
It's not a specific problem I'm working on. I was actually using Maxima to verify some math doing in another system, and realized that it wasn't returning the (mathematically) correct answer when the second arg is float. So I reported it.
But since you (Stavros) mentioned a wider view of how float math is handled, then I agree! It should either work all the time, or not.
Sorry about the noise then.
(I am still having trouble posting to the Sourceforge bug)
How is Maxima supposed to know that 2.1 is intended to represent the
rational 21/10? After all, its actual value is 21/10+1/(2^515), and it may
have been the result of some long calculation, not from the input string
"2.1". The exact, mathematically correct value of mod(10^100,2.1) is
actually float(mod(10^100,rationalize(2.1))) = 1.43... By using
rationalize*, we're using the exact rational value of 2.1.
Large numbers aren't even necessary to show this. For example, mod(37,3.7)
=> 3.6999999999999957 although of course mod(37,37/10) = 0. That is
because 3.7 is actually 37/10+1/(2^50*5).
Maxima does have a way of finding close rational approximations to
floating point numbers. So for example rat(3.7) => 37/10. See ?
ratepsilon for details.
Floating point arithmetic is "contagious". That is, mod(10^100,2.1) is
equivalent to mod(float(10^100),2.1). But neither 10^100 nor 2.1 is
represented exactly as a floating point number:
rationalize(float(10^100))-10^100 => ~10^83
rationalize(2.1)-21/10 => ~10^-16.
To get an exact answer, use mod(10^100,21/10) => 19/10.
Note that rationalize converts a floating point number to the exact
rational it represents, e.g., rationalize(0.1)
=> 3602879701896397/36028797018963968 and rationalize(0.1)-1/10 =>
1/180143985094819840 = ~10^-16.
On Sat, Jan 10, 2026, 09:23 Jeronimo Pellegrini via Maxima-bugs maxima-bugs@lists.sourceforge.net wrote:
Related
Bugs:
#4659