[erlang-questions] A math module better than C89
Richard A. O'Keefe
ok@REDACTED
Thu Jan 19 09:28:14 CET 2017
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.)
More information about the erlang-questions
mailing list