Make `expt` more precise

[jpellegrini]

If the exponent is a floating point number that happens to represent an integer precisely (that is, with fractinal part equal to exactly zero), turn it into exact, perform the operation, then turn it back into inexact.

Before this patch:

stklos> (real-precision 50)
stklos> (expt 2 3.0)
7.9999999999999982236431605997495353221893310546875
stklos> (exact (expt 2 3.0))
4503599627370495/562949953421312

After this patch:

stklos> (real-precision 50)
stklos> (expt 2 3.0)
8.0
stklos> (exact (expt 2 3.0))
8

[jpellegrini]

Actually, this would force the use of the GMP exponentiation without need. I think the reason why $2^{(3.0)}$ is imprecise is because of the fallthrough from the tc_real case to the tc_complex case.

    case tc_real:     if (zerop(y)) return double2real(1.0);
                      if (zerop(x)) return (x==MAKE_INT(0)) ? x : double2real(0.0);
                      /* FALLTHROUGH */
    case tc_complex:  if (zerop(x)) {
                  ...
                      } else return my_exp(mul2(my_log(x),y));

So expt is being calculated using one multiplication, one log and one exp... But pow is more precise, and would probably be better when x is real and wen the result fits a double:

#include <stdio.h>
#include <math.h>

int main() {
  double base = 2.0;
  double expo = 3.0;
  printf("base = %.60f\n", base);
  printf("expo = %.60f\n", expo);
  printf("exp(log(expo)*base) = %.60f\n", exp(log(base) * expo));
  printf("pow(base,expo)      = %.60f\n", pow(base,expo));
}

The result:

./precision 
base = 2.000000000000000000000000000000000000000000000000000000000000
expo = 3.000000000000000000000000000000000000000000000000000000000000
exp(log(expo)*base) = 7.999999999999998223643160599749535322189331054687500000000000
pow(base,expo)      = 8.000000000000000000000000000000000000000000000000000000000000

I'll fix this PR with a better solution.

[jpellegrini]

I think it's good now...

@egallesio I'll be busier than usual, so I'll stop sending new PRs for a while, probably. But I'll fix anything that shows up in the already submitted ones.

[jpellegrini]

One observation: I tried to add SRFI 194 before, and I found two problems: the first was that math operations were slow -- and I think my recent PRs somehow help with that; the other problem was that some tests wouldn't pass, apparently due to lack of precision. Maybe after all these PRs are merged (if you accept them), we can add that SRFI to STklos. :)

[jpellegrini]

This fixes the real case, but for rationals there is still a problem:

stklos> (expt 4 3/2)
7.9999999999999982236431605997495353221893310546875

Will fix in a different PR, in the future.

[jpellegrini]

Idea for the rationals (future, not this PR):

To compute $x^{N/D}$, there are two paths:

1. compute $y = x^{1/D}$, then do the usual exponentiation algorithm to get $y^N$
2. compute $y = x^N$, then compute $y^{1/D}$

I didn't think much about it (not enough time right now), but it may be simple to decide which one will better at preserving precision.

Then an algorithm would be needed for computing $y^{1/D}$, or $x^{1/D}$ -- Newton's method for n-th root would be the natural candidate (?)

There are other little subtleties (how and when do we detect that we need to answer "zero" or "inf" due to under/overflow in the inexact result?)

[egallesio]

That's great! Merged. Thanks, @jpellegrini .