Enhancements to bignums (including a big speedup to `bignum2double` and `exact->inexact`)

[jpellegrini]

Hi @egallesio !

One more PR... :) The second one seems interesting, because bignum2double is used a lot in the code...

1. Use mpz_sgn instead of mpz_cmp_si when possible

At least for positive? and negative? there seems to be some speedup - and it seems simpler.

2. Speed up bignum2double

This patch adds only 20 lines of code, and a lot more of comments... :)

Currently, bignum2double does not use the function mpz_get_d since it gives an unspecified value when converting a number which is +inf or -inf.

I have verified that it also brings rounding problems: if we use mpz_get_d, the resulting double will be rounded towards zero, but R7RS requires that we pick the closest number.

What this patch does

1. Treat infinities separately

If the result does not fit a double, we return $\pm$inf.0 We don't use the result of mpz_get_d because it does not guarantee that an inf will be returned in this case.

Checks if the number fits into a double, comparing it to $\pm$DBL_MAX, and directly returns $\pm$inf.0 if the number doesn't fit.

2. Use mpz_get_d and find the closest integer representable as double

A very large integer may not be representable as a float.

mpz_get_d always "rounds towards zero" -- that is, it returns the closest float to n that is between 0 and n, but R7RS requires 'inexact' to return the closest number. So we need to adapt.

Suppose there are two representable integers around n, but n itself is not representable. Call those integers 'below' and 'above':

           0                        below       n  above
        ---|--------------------------|---------|----|-----
                                      v
                                   returned
                                    by GMP

As the figure shows, even if there is an integer 'above' that is closer to n, the number 'below' (closer to 0) will be returned. Note that the whole picture could be reflected around zero if n is negative, so "above" actually means "farthest from zero" but not "largest":

           above  n       below                        0
        -----|----|---------|--------------------------|---
                                                       |
                                                neg <--+--> pos

We of course can be sure that there is no integer between 'below' and 'above'.

So we do the following:

1. Get the next representable double with nextafter (the one starting from below, but away from zero). We call this one 'above'

2. Convert both below and above back into bignums (!), as 'zbelow' and 'zabove'

3. Measure (using the GMP) the distances ABS(zabove - n) and ABS(n - zbelow) and if n is closer to zabove, we return ceil(above) or floor(above), depending on the sign. If it's closer to below, we return below.

The speedup

Average of 10 runs is shown.

1. bignum does not fit a double

In this case the time for the currnet STklos implementation is proportional to the size of the number. We take $3^{431} \cdot 2^{560} \cdot 5 \cdot 11$ as an example.

(let ((b  (* (expt 3 431) 5 11 (expt 2 560) ))
      (a 0))
  (time
    (dotimes (i 10_000_000)
      (set! a (inexact b))))
  a)
current code: 34145.194 ms
this patch:    1574.554 ms

2. bignum does fit a double

Both the STklos implementation and the proposed patch do depend on the size of the number, but there is a great speedup.

;; For "not very large bignums"
(let ((b  (* (expt 3 31) 5 11 (expt 2 60) ))
      (a 0))
  (time
    (dotimes (i 10_000_000)
      (set! a (inexact b))))
  a)
current code: 5720.891 ms
this patch:   3867.123 ms

;; For "really large bignums", still representable as doubles:
(let ((b  (* (expt 3 31) 5 11 (expt 2 760) ))
      (a 0))
  (time
   (dotimes (i 10_000_000)
     (set! a (inexact b))))
  a)
current code: 29688.393 ms
this patch:    7489.038 ms

(From almost 30s to 7.5s)

Tests

Both the numbers test in STklos and the tests for SRFI 144 detect rounding problems, and indeed, previous attempts for preparing this patch were wrong, and it was only because of those tests that I realized that.

[jpellegrini]

By the way... Guile also does have an optimized version of exact->inexact, but they use the low-level GMP functions that work on the internal GMP representation for numbers (the functions that deal with number "limbs") -- I thought this would be a bit too complex and didn't go that way.

[egallesio]

Another great contribution. Thanks, @jpellegrini. PR is merged.