janet-lang / janet

A dynamic language and bytecode vm
https://janet-lang.org
MIT License
3.45k stars 223 forks source link

Make floating point `rem` more consistent #1229

Closed primo-ppcg closed 1 year ago

primo-ppcg commented 1 year ago

One thing I've noticed is that % with a non-integer divisor is not consistent in terms of modular arithmetic:

(def d 3.2)
(for n -99 100
  (assert (= n (+ (* d (math/trunc (/ n d))) (% n d)))))

This will fail for any number of divisors, in particular when the binary representation is not exact.

It also has some interesting artifacts:

(pp (seq [n :range [0 10]] (% n 0.8)))
# ->
@[0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2]

0.8 evenly divides both 4 and 8, but the remainder produced is 0.8 (or rather, a value minisculely smaller than 0.8), and not zero. I also suspect that this may differ across platforms, depending on the implementation of fmod, although I'm not absolutely sure of that.

It's true that floating point arithmetic is not meant to be exact, but that doesn't prevent us from defining the remainder in a way that is consistent, both in terms of modular arithmetic, and also across platforms:

diff --git a/src/core/vm.c b/src/core/vm.c
@@ -743,4 +743,5 @@
             double x1 = janet_unwrap_number(op1);
             double x2 = janet_unwrap_number(op2);
-            stack[A] = janet_wrap_number(fmod(x1, x2));
+            double intres = x2 * trunc(x1 / x2);
+            stack[A] = janet_wrap_number(x1 - intres);
             vm_pcnext();

Which parallels the current implementation of mod. I've implemented this on my local fork for the purpose of testing. Aside from resolving both issues above, this also appears to be slightly faster, at least on my system:


(use spork/test)

(timeit-loop [:timeout 2]
  (for n 0 1000 (% n 3.1415926)))
# ->
master  2.000s, 15.75µs/body
local   2.000s, 11.12µs/body
bakpakin commented 1 year ago

0.8 evenly divides both 4 and 8, but the remainder produced is 0.8 (or rather, a value minisculely smaller than 0.8), and not zero. I also suspect that this may differ across platforms, depending on the implementation of fmod, although I'm not absolutely sure of that.

Remainder and Mod are usually defined differently - remainder specifically is not defined in terms of modular arithmetic. Use mod instead of % for that.

bakpakin commented 1 year ago

Also, 0.8 is not representable exactly in ieee-754, so yes there is going to be some inexactness there.

primo-ppcg commented 1 year ago

remainder specifically is not defined in terms of modular arithmetic

Yes, I think you're right. I've just been mentioning issues as I encounter them. Sometimes it's user error 😉