egallesio / STklos

STklos Scheme
http://stklos.net
GNU General Public License v2.0
69 stars 17 forks source link

Make `sqrt` work on (some) bignums #541

Closed jpellegrini closed 1 year ago

jpellegrini commented 1 year ago

While it's not as advantageous as for log, it would also be interesring to have sqrt work on bignums (so long as the result fits a double, of course). This patch does that.

1. create a rational_epsilon and move computation of DBL_TRUE_MIN

In order to use Newton's method to get an approximation (and not the exact result), we need an epsilon. We compute it as a rational less than the least representable number. For that, we need to move the computation of DBL_TRUE_MIN from (scheme flonum) into number.c.

3. Change my_sqrt_exact to approximate when possible

First, we use GMP's mpz_sqrtrem to check if the square root is exact (this seems better than checking back the result, as was being done before).

If not, then we use Newton's method, applying STklos C arithmetic (sub2, add2, div2, mul2, STk_abs...)

4. Include some tests

Some new tests for square roots of bignums are added.

This one, for example, would result in +inf.0, but now returns an approximation:

(sqrt (* (expt 2 700) (expt 2 700) 3))  =>  9.1108226361989e+210
jpellegrini commented 1 year ago

This page: https://github.com/schemedoc/surveys/blob/master/surveys/log-requires-floating-point.md has a survey on the results of (sqrt (+ (expt 2 2030) 111111111111111111111111117)) in several Scheme implementations. This is what this PR would change for STklos.