Better use of GMP in `expt` (fixes crash)

jpellegrini commented 1 year ago

1. Include missing calls to `mpz_init`

In exact_exponent_expt, some calls were missing, and STklos would crash on this, for example:

(expt (expt 2 100) -2)

This patch also includes a test for that.

If an exponent is reeealy too large:

(expt fx-greatest 100000)

STklos will report a SIGABRT from the GMP.

Anyway, when the exponent is too large, several Scheme implementations either crash (Chez, Guile, Sagittarius), or report an error like "our of memory" or "exponent too large" (Gosh, MIT). Kawa and Bigloo return an imprecise infinity -- so the behavior of STklos should be OK.

2. Don't create and initialize a new `mpz_t` number unnecessarily

Use mpz_ui_pow_ui instead of creating and setting the value of a new mpz number to then call mpz_pow_ui.

Since that function only accpets unsigned integers, we check the sign, pass the number as unsigned, then if the base was negative and the exponent was odd, we multiply by -1.
Call mpz_init(res) for the tc_bignum case.

The first change seems to make exponentiation with fixnum base slightly faster:

;; fixnum result:
(let ((a 0))
  (time
   (repeat 10000000
     (set! a (expt (expt 2 5) (expt 2 3)))))
  a)
;; 6545.442 => 5446.367

;; bignum result:
(let ((a 0))
  (time
   (repeat 10000000
     (set! a (expt (expt 2 10) (expt 2 5)))))
  a)
;; 8859.299 => 7478.944

(Average of 20 runs)

The sign change didn't make much of a difference in timings (That is, timings seem similar for iterated computation of (expt -32 9) and (expt 32 9).

3. Call `mpz_clear` before returning from `expt`

This actually makes it faster:

(let ((a (* (expt 2 60)))
      (b 30)
      (c #f))
  (time
   (repeat 5_000_000
           (set! c (expt a b))))
  c)

4660.171 ms => 2900.074 ms

jpellegrini commented 1 year ago

Force-pushed to fix conflicts due to merging of other PRs. Still passes all tests! :)

egallesio commented 1 year ago

Great, great, great.

Thanks @jpellegrini for this PR.

egallesio / STklos