`math:sqrt` and `math:log` could have wider domains

[jpellegrini]

Currently, math:log and math:sqrt use the C functions log and sqrt. This is fine, but in the presence of bignums (as Erlang has), it becomes possible to wider their domains:

1. $10^{2480}$ cannot be used as an argument to math:log, only because it is not representable as an IEEE double float. But it is representable as an Erlang bignum, and the result is an IEEE double float:
   $log(10^{2480})$ is close to $5710.411030625233$
   So it would perhaps be possible to make math:log accept larger numbers (a lot of them, since log grows, well, logarithmically :)

2. $2^{2030}+111111111111111111111111117$ does not fit a double either. But similarly, it is an Erlang integer, and its square root does fit an IEEE double float (it is close to $3.511119404027961 \times 10^{305}$).

I'm not saying, of course, that getting this to work is "easy". Anyway, here's a survey on how many different Scheme implementations deal with it in different ways (but no implementation details here):

https://github.com/schemedoc/surveys/blob/master/surveys/log-requires-floating-point.md

And I can explain how STklos Scheme gets it done (because I have implemented these two enhancements):

1. STkloslog of bignums is very simple, but it uses the GMP (I know Erlang doesn't):

static double my_bignum_log(SCM z) {
  /* The GMP function mpz_get_d_2exp is similar to the C library
     function frexp: it returns a float D between 0.5 and 1.0,
     and sets an exponent E. Then,
     N = D * 2^E
     is the (possibly truncated) value of the original number.
     Then calculating the log is simple:

     log(N) = log(D*2^E)
            = log(D) + log(2^E)
            = log(D) + log(2) * E                               */
  long   expo;
  double d = mpz_get_d_2exp(&expo, BIGNUM_VAL(z));
  return log(d) + log(2) * (double) expo;
}

But since Erlang has its own bignum implementation, maybe it won't be hard to get this working (?)

2. STklos sqrt tries a couple of methods (like identifying perfect squares etc), and when all fail, it uses Newton's method. It is slow (because it is computed using bignum rationals), but we get some more answers for the function.

    /* Ok, we tried everything. There's only the slow path now! */
    SCM x = x0;
    SCM x_new = x0;
    SCM err = x0;

    /* Approximate the square root... Essentially, Newton's method,
       but coded using STklos' internal sub2, div2, mul2, abs
       functions.  */
    while(STk_numgt2(err, rational_epsilon) > 0 &&
          isfinite(REAL_VAL(exact2inexact(x_new)))) {
        x_new = sub2(x, div2(sub2(mul2(x, x), z),
                             mul2(x, MAKE_INT(2))));
        err = STk_abs(sub2(x_new, x));
        x = x_new;
    }
    /* Return inexact, because if we got here, the square of this
       result will not equal to z (it's an approximation, so it would
       be strange to give an "exact" result that is not "exactly" the
       result of the operation). But for floating-point, it is
       acceptable to offer an approximation.  */
    return exact2inexact(x_new);

Here mul2, div2, sub2 are generic 2-argument functions, and x is a bignum (so the computation is fully rational until the last moment, when exact2inexact converts it to a double).

Also,

    SCM x = x0;

Here, x0 is computed earlier using mpz_sqrt, so it is the truncated square root.

For reference, the math code for STklos (the C part) is here:

https://github.com/egallesio/STklos/blob/master/src/number.c

By no means I believe Erlang "needs" this -- I see it as a low-priority enhancement; and I also wouldn't say that my implementation for STklos makes sense for Erlang... Just wanted to share, in case the Erlang implementors find it to be useful.

[RaimoNiskanen]

I think that it could be nice to extend math:log/1, math:log2/1 and math:log10/1 since that would extend the input range from "all numbers except too large integers" to "all numbers", which is a clear simplification.

It is probably not very hard to do this with help from Erlang's bignum implementation (erts/emulator/beam/big.[ch]).

But extending math:sqrt/1's input range from "all numbers expect too large integers" to "all numbers except twice as large integers than before" seems less useful.

And this is, as you say, a low priority enhancement...

[jpellegrini]

But extending math:sqrt/1's input range from "all numbers expect too large integers" to "all numbers except twice as large integers than before" seems less useful.

Yes, I agree. I proposed it because then the library would be more consistent, and implementing Newton's method is easy.

Besides those, there are

1. the trigonometric functions... But those are really hard to get working for large arguments (range reduction for sin/cos seems easy, but I wouldn't know how to do that without a hugely precise approximation of pi). Also, math:atan already returns the best IEEE possible approximation for pi/2 when given huge numbers. :)
2. math:pow. One could say that math:pow(2,100) should return an integer (which is possible in Erlang). But since its type signature says the result is float, its behavior seems ok (signals an error).

[RaimoNiskanen]

Agreed. There is also math:fmod(X, Y) that could be augmented to do float(X rem Y) if X and Y are integers, but that would leave out X big integer and Y float that also has a valid floating point result for any big integer X.

[michalmuskala]

One thing to keep in mind is that the compiler assumes certain types for arguments and results of the math functions. Therefore, to keep the compatibility of allowing older-compiler modules to run on the next 2 OTP releases, this needs to be careful not to change the domains of those functions in a way that would violate compiler's assumptions.

[jpellegrini]

this needs to be careful not to change the domains of those functions in a way that would violate compiler's assumptions

Yes. What I proposed was not to change type signatures at all. It's just that when passing a huge bignum as "t_number", it may trigger an error unecessarily, and that restriction could be relaxed.

• ❓ cos, sin, tan - could be done, but it's hard. Range reduction would make this work, but to apply range reduction properly, one needs a high precision approximation of pi (I have implemented this in C, using GMP bigfloats, for example)
• ✅ acos, asin - not needed, since Erlang does not support complexes
• ✅ atan, atan2 - hmm, it already returns the IEEE value for $\pi/2$ when the angle is too large, and accepts bignums
• ✅ sinh, cosh - no change possible: if x is already a huge number, it becomes even larger.
• ❓ acosh, asinh- these could have their domain extended. For example, $asinh(10^{500})$ is $\sim 1151.98$, which fits a double float. But I'm not sure how to implement this one.
• ✅ atanh - not relevant, since Erlang does not support complexes.
• ✅ tanh - already works for huge numbers
• ✅ erf, erfc - already work for huge numbers
• ✅ exp - no change possible: if x is already a huge number, $e^x$ is even larger, and won't fit a float
• ❗ log, log2, log10 - this is the nicest part. These can have their domain extended hugely.
• ❗ sqrt - it is possible to widen the domain a little bit, but not as much as in the case of logs.
• ✅ pow - here @michalmuskala 's comment is relevant. It must return a float, so no change can be done. And there is a separate issue for a new integer exponentiation operator (ERL-834: math:pow limitations #3978)
• ✅ ceil, floor - already works for all bignums
• ❓ fmod - as @RaimoNiskanen said, it may be extended
• ✅ pi, tau - obvisouly not possible to change, as they take no argument :)

So, without changing types:

• the logs - possible. Huge benefit.
• sqrt - possible. Small effect.
• acosh, asinh - possible. I've never tried.
• fmod - possible. Not sure how to deal with the issue @RaimoNiskanen raised.

But anyway... As I said, not a big priority thing. :)

[jpellegrini]

[ this comment is not that important; see later about rationalizing Y ]

fmod - possible. Not sure how to deal with the issue @RaimoNiskanen raised.

IEEE floats complicate things a little here.

In Maxima:

(%i4) mod(10,2.1);
(%o4)                         1.5999999999999996
(%i5) mod(100,21);
(%o5)                                 16
(%i6) mod(100,21.0);
(%o6)                                16.0

Yes.

• for (100,21), $21\times 4 = 84$, and the remainder is $16$
• for (10,2.1), $2.1\times 4 = 8.4$, and the remainder should be $1.6$.

The problem is that "2.1" is not representable precisely as an IEEE double float. So ten times 2.1 it is not the same as 21 (in IEEE math).

But multiplying both sides by a power of ten (and then dividing the result by the same power) may be a good enough approximation. Or maybe not? I don't remember looking into implementations of fmod before in other places, so I don't know if there are techniques for getting more precise results.

[jpellegrini]

[ this comment is not that important; see later about rationalizing Y ]

And even Maxima does some weird things with large numbers:

(%i5) mod(1000000000000000000000000000000000000000000000000000000000000000000000000000000,21);
(%o5)                                  1
(%i6) mod(100000000000000000000000000000000000000000000000000000000000000000000000000000,2.1);
(%o6)                        1.2855504354071922e61

So for fmod I'm not sure it's feasible to extend the domain. See below!

[jpellegrini]

[ this comment is not that important; see later about rationalizing Y ]

Or, maybe:

1. Extract mantissa $m$ and exponent $e$ of Y

Y = m · 2^e ,   0.5 ≤ m < 1   (binary floating‑point)

Now shift X right by $e$ bits (this is integer division by $2^e$).
Result: an integer X' = floor(X / 2^e) and the remainder r = X rem 2^e -- both integers.

2. Do modulo with the mantissa:

Convert the mantissa $m$ to some rational $p/q$ where $p$ and $q$ (e.g., p = round(m·2^k), q = 2^k, with $k$ being enough bits to represent $m$).

Now compute t = (X' rem q) * p rem q. The fractional part is $t / q$.

3. Combine the two remainders

result = (r / 2^e) + (t / q)

If $result \geq Y$ then subtract Y once (this is possible due to rounding).

Not sure if this is correct.

[jpellegrini]

But if we could rationalize Y, then it would be perfect!

float_to_rational(Float) would return numerator and denominator.

fmod(X,Y) would be X*N rem D, where N and D are numerator and denominator of Y

Small proof-of-concept implementation of a new fmod with wider domain in Erlang.

-module(new_fmod).

-export([new_fmod/2,float_to_rational/1]).

-define(MIN_SMALL, - (1 bsl 27)).   % 32‑bit VM example
-define(MAX_SMALL,   (1 bsl 27) - 1).

float_to_rational(Float) ->
    Precision    = 10000,
    Fraction     = round(Float * Precision),
    Denominator  = Precision,
    Gcd          = gcd(Fraction, Denominator),
    Numerator    = Fraction div Gcd,
    Denominator1 = Denominator div Gcd,
    {Numerator, Denominator1}.

gcd(A, 0) -> A;
gcd(A, B) -> gcd(B, A rem B).

is_big_int(N) when is_integer(N),
                   (N < ?MIN_SMALL orelse N > ?MAX_SMALL) -> true;
is_big_int(_) -> false.

new_fmod(X,Y) ->
    IntX    = is_integer(X),
    BigintX = is_big_int(X), 
    FloatX  = is_float(X),
    IntY    = is_integer(Y),
    BigintY = is_big_int(Y),
    FloatY  = is_float(Y),
    case {IntX, BigintX, IntY, BigintY, FloatY} of
        {true, _, true, _, _}            -> float(X rem Y);
        {true, true, false, false, true} -> begin
                                                {N,D} = float_to_rational(Y),
                                                float(D*X rem N)
                                            end;
        {_,_,_,_,_}                      -> math:fmod(X,Y)
    end.

• The line Precision = 10000 would need to be adjusted
• MAX_SMALL and MIN_SMALL would need to be adjusted. Actually, the is_big_int function (which uses them) is a hack to get this example to work. I am sure the VM has this internally, with another name.

See that this is more precise than the current fmod for at least one case, even for integers! (See below)

20> math:fmod(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000,2.1).
0.532783491487582
21> math:fmod(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*10,21).
2.0
22> new_fmod:new_fmod(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000,2.1).
19.0

(The correct answer is 19)

Cool side-effect: a function to obtain the rational representation of a float could also be added to the OTP API. :)

[jpellegrini]

Forgot to mention: if there is interest, the rationalizing function could be re-written using continued fractions. Then a function to give the user a rational representation of a float would return a "simple" fraction.

For example, most Scheme and Lisp implementations will turn inexact (float) numbers into exact (rationals) using a standard algorithm:

> (exact 1.0001234)
2252077685782257/2251799813685248

One implementation (Gauche Scheme) uses continued fractions instead:

gosh$ (exact 1.0001234)
5000617/5000000

And guess? Those two rationals will evaluate to the SAME IEEE floating point number! :)

Just in case the Erlang implementors feel that having a float_to_rational function is a good idea.

Or, since "rational" is not an Erlang type, call it float_to_num_and_denom, or something...