[erlang-questions] A math module better than C89

[Richard A. O'Keefe]

I honestly don't think there's much difficulty figuring
out what the math functions *should* do.

To start with, we cannot offer a math library better than
C89 until we offer a math library as good as C89, and we
just don't do that.
  - there is no frexp and that turns out to matter a lot
  - there is no ldexp
  - there is no floor (round is provided elsewhere)
  - there is no ceil (trunc is provided elsewhere)
  - there is no fmod
  - there is no modf
There's also a function that was in the Unix libraries but
didn't make it into the C standard until 1999:
  - there is no hypot (hypot and atan2 are *partners*).

Since it's now 18 years since C99, shouldn't we be thinking
about a math module better than C99?  (The real part anyway.)

Let's review what's in the math module, working from
http://erlang.org/doc/man/math.html

In working on my Smalltalk library, I wrote support
functions to convert a bignum to a normalised
fraction in float, double, or long double form and
an integer power of 2, rather like frexp.  Putting
this in Erlang terms, if we have
     frexp(Number) -> {Fraction, Scale}
then
     log(Number) ->
         {Fraction, Scale} = frexp(Number),
         underlying:log(Fraction) + Scale * underlying:log(2.0).

This means not having to worry about floating point overflow
when dealing with logs of (positive) bignums.

As it happens, ANSI Smalltalk does not include a logarithm
operation on integers, so I *didn't* do this.  Before going
to sleep tonight, I will!

acos(X)
     X must be between -1 and 1 inclusive.
     bignums must be reported as wrong.
acosh(X)
     log(X + (X+1)sqrt((X-1)/(X+1))
     (X-1)/(X+1) will be effectively 1,
     sqrt((X-1)/(X+1)) will be effectively 1,
     acosh(X) will be approximately log(2X+1).
     See note about log above.

asin(X)
     X must be between -1 and 1 inclusive.
     bignums must be reported as wrong.
asinh(X)
     log(X + sqrt(1 + X**2))
     X**2 + 1 will be effectively 1,
     sqrt(X**2 + 1) will be effectively |X|,
     asinh(X) will be effectively log(2X).
     See note about log above.

atan(X)
     Bignums that would be converted to ±inf
     may as well be allowed to,
     atan(±inf) = ±pi

atan2(X, Y)
     After dealing with X == 0 or Y == 0 specially,
     remember the signs.  We can work with
     {FX, SX} = frexp(X),
     {FY, SY} = frexp(Y),
     FR = FY/FX, SR = SY-SX
     R = if SR < lower threshold -> 0
          ; SR > upper threshold -> +infinity
          ; true                 -> ldexp(FR, SR)
         end
     Now compute the result from R and the saved signs.

atanh(X)
     tanh(Y) lies between -1 and 1 inclusive, so
     any integer, big or small, other than -1, 0, 1
     is an error.

MISSING: cbrt/1.
     See sqrt.

cos(X)
     To be honest, I don't think this really makes sense.
     I know that the Sun math library claimed to get as
     good an answer as you could hope for using
     "infinitely precise pi".  I have no idea how that
     works.  The old UNIX libm interface used to report
     doubles that were too big as basically not terribly
     meaningful.  (If an ulp is larger than pi, you
     really cannot represent an angle in any meaningful
     way.  And that happens about 2.8e16.  So frankly,
     if a bignum would convert to infinity, I say let
     it, and let the C library tell you that cos(inf) is NaN.
cosh(X)
     (exp(X) + exp(-X))/2 will be near exp(|X|).
     +infinity for a bignum.  See exp.

exp(X)
     The result is a float.  For any float X greater than
     about 709.78 the answer will be +infinity.
     A negative bignum can give the answer 0.0.
     A positive bignum should give the answer +infinity.

MISSING AND IMPORTANT: expm1/1

log(X)
     Described above.  We can do this!

log10(X)
     is log(X)/log(10)

log2(X)
     is log(X)/log(2)

MISSING AND IMPORTANT: log1p/1

pow(X, Y)
     Smalltalk distinguishes between number raisedToInteger: integer
     (which typically uses the bit oriented square and multiply
     approach) and base raisedTo: power which reverts to
     raisedToInteger: if power is an integer otherwise uses
     exp(log(base) * power).
     This needs a careful case analysis.
     pow(Bignum, Float) will of course only be allowed
     when Bignum > 0.  Since we can compute log(Bignum),
     this reduces to exp(log(Bignum)*Float).  Sometimes
     that will give 0, sometimes a finite number, and
     sometimes ±infinity.
     It's all doable by careful case analysis.

sin(X)
     See comment on cos(X).
sinh(X)
     (exp(X) - exp(-X))/2.
     If bignum X > 0, -> exp(X)/2 -> +infinity.
     If bignum X < 0, -> -exp(-X)/2 -> -infinity.

sqrt(X)
     Smalltalk actually defines sqrt on integers as integer
     square root, which is a pain in the posterior.  I have
     wanted the integer square root function occasionally,
     true.  But I have much more often wanted
     anInteger sqrt -> anInteger asFloat sqrt, which is what
     Erlang does, hooray, hooray.
     If bignum X < 0, this is an error.
     Let {F, S} = frexp(X).
     Let {G, T} = if S odd  -> {sqrt(F*2.0), (S-1) div 2}
                   ; S even -> (sqrt(F), S div 2)
                  end.
     Answer ldexp(G, T).
     This will give you good answers up to 3.2x10**616
     and then infinity.

tan(X)
     See comment on cos.
tanh(X)
     (exp(X) - exp(-X))/(exp(X) + exp(-X))
     For any X > about  709.78, answer +infinity.
     For any X < about -709.78, answer -infinity.

erf(X)
     From the graph of erf at http://mathworld.wolfram.com/Erf.html
     it is pretty obvious that
     For bignum X > 0, erf(X) -> 1.0   and erfc(X) -> 0.0
     For bignum X < 0, erf(X) -> -1.0  and erfc(X) -> 2.0

pi()
     has no argument, so no bignum issue.


I presume pow(X, Y, M) computes (X mod M)**Y mod M?
This reduces to multiplication modulo M
(is that '*'(X, Y, M) ?) which is not problematic.

The core here is two operations that have been in the C library
(and in many languages under various guises) but are not in the
Erlang library:  frexp and ldexp.  I have a faint recollection
of libgmp offering frexp, if not I could donate my C code.

Of course many of these things are in other languages with
bignums.  (Hands up everyone who realised I was staring at
page 308 of CLtL2?  Which doesn't actually mention bignums.)

By the way, my HoD and HR are growling at me for not publishing
enough, and a friend this afternoon said "don't be ashamed of
publishing simple things".  Could I get a paper out of this
if I finished the case analyses and wrapped a few more words
around it?  (Seriously, I am quite stressed about the situation.)