<html>
<head>
<meta content="text/html; charset=utf-8" http-equiv="Content-Type">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<div class="moz-cite-prefix">On 01/19/2017 12:28 AM, Richard A.
O'Keefe wrote:<br>
</div>
<blockquote
cite="mid:5d72f814-7dc6-bb56-96d8-fc68035ebbc8@cs.otago.ac.nz"
type="cite">I honestly don't think there's much difficulty
figuring
<br>
out what the math functions *should* do.
<br>
<br>
To start with, we cannot offer a math library better than
<br>
C89 until we offer a math library as good as C89, and we
<br>
just don't do that.
<br>
- there is no frexp and that turns out to matter a lot
<br>
- there is no ldexp
<br>
- there is no floor (round is provided elsewhere)
<br>
- there is no ceil (trunc is provided elsewhere)
<br>
- there is no fmod
<br>
- there is no modf
<br>
There's also a function that was in the Unix libraries but
<br>
didn't make it into the C standard until 1999:
<br>
- there is no hypot (hypot and atan2 are *partners*).
<br>
<br>
Since it's now 18 years since C99, shouldn't we be thinking
<br>
about a math module better than C99? (The real part anyway.)
<br>
<br>
Let's review what's in the math module, working from
<br>
<a class="moz-txt-link-freetext" href="http://erlang.org/doc/man/math.html">http://erlang.org/doc/man/math.html</a>
<br>
<br>
In working on my Smalltalk library, I wrote support
<br>
functions to convert a bignum to a normalised
<br>
fraction in float, double, or long double form and
<br>
an integer power of 2, rather like frexp. Putting
<br>
this in Erlang terms, if we have
<br>
frexp(Number) -> {Fraction, Scale}
<br>
then
<br>
log(Number) ->
<br>
{Fraction, Scale} = frexp(Number),
<br>
underlying:log(Fraction) + Scale * underlying:log(2.0).
<br>
<br>
This means not having to worry about floating point overflow
<br>
when dealing with logs of (positive) bignums.
<br>
<br>
As it happens, ANSI Smalltalk does not include a logarithm
<br>
operation on integers, so I *didn't* do this. Before going
<br>
to sleep tonight, I will!
<br>
<br>
acos(X)
<br>
X must be between -1 and 1 inclusive.
<br>
bignums must be reported as wrong.
<br>
acosh(X)
<br>
log(X + (X+1)sqrt((X-1)/(X+1))
<br>
(X-1)/(X+1) will be effectively 1,
<br>
sqrt((X-1)/(X+1)) will be effectively 1,
<br>
acosh(X) will be approximately log(2X+1).
<br>
See note about log above.
<br>
<br>
asin(X)
<br>
X must be between -1 and 1 inclusive.
<br>
bignums must be reported as wrong.
<br>
asinh(X)
<br>
log(X + sqrt(1 + X**2))
<br>
X**2 + 1 will be effectively 1,
<br>
sqrt(X**2 + 1) will be effectively |X|,
<br>
asinh(X) will be effectively log(2X).
<br>
See note about log above.
<br>
<br>
atan(X)
<br>
Bignums that would be converted to ±inf
<br>
may as well be allowed to,
<br>
atan(±inf) = ±pi
<br>
<br>
atan2(X, Y)
<br>
After dealing with X == 0 or Y == 0 specially,
<br>
remember the signs. We can work with
<br>
{FX, SX} = frexp(X),
<br>
{FY, SY} = frexp(Y),
<br>
FR = FY/FX, SR = SY-SX
<br>
R = if SR < lower threshold -> 0
<br>
; SR > upper threshold -> +infinity
<br>
; true -> ldexp(FR, SR)
<br>
end
<br>
Now compute the result from R and the saved signs.
<br>
<br>
atanh(X)
<br>
tanh(Y) lies between -1 and 1 inclusive, so
<br>
any integer, big or small, other than -1, 0, 1
<br>
is an error.
<br>
<br>
MISSING: cbrt/1.
<br>
See sqrt.
<br>
<br>
cos(X)
<br>
To be honest, I don't think this really makes sense.
<br>
I know that the Sun math library claimed to get as
<br>
good an answer as you could hope for using
<br>
"infinitely precise pi". I have no idea how that
<br>
works. The old UNIX libm interface used to report
<br>
doubles that were too big as basically not terribly
<br>
meaningful. (If an ulp is larger than pi, you
<br>
really cannot represent an angle in any meaningful
<br>
way. And that happens about 2.8e16. So frankly,
<br>
if a bignum would convert to infinity, I say let
<br>
it, and let the C library tell you that cos(inf) is NaN.
<br>
cosh(X)
<br>
(exp(X) + exp(-X))/2 will be near exp(|X|).
<br>
+infinity for a bignum. See exp.
<br>
<br>
exp(X)
<br>
The result is a float. For any float X greater than
<br>
about 709.78 the answer will be +infinity.
<br>
A negative bignum can give the answer 0.0.
<br>
A positive bignum should give the answer +infinity.
<br>
<br>
MISSING AND IMPORTANT: expm1/1
<br>
<br>
log(X)
<br>
Described above. We can do this!
<br>
<br>
log10(X)
<br>
is log(X)/log(10)
<br>
<br>
log2(X)
<br>
is log(X)/log(2)
<br>
<br>
MISSING AND IMPORTANT: log1p/1
<br>
<br>
pow(X, Y)
<br>
Smalltalk distinguishes between number raisedToInteger:
integer
<br>
(which typically uses the bit oriented square and multiply
<br>
approach) and base raisedTo: power which reverts to
<br>
raisedToInteger: if power is an integer otherwise uses
<br>
exp(log(base) * power).
<br>
This needs a careful case analysis.
<br>
pow(Bignum, Float) will of course only be allowed
<br>
when Bignum > 0. Since we can compute log(Bignum),
<br>
this reduces to exp(log(Bignum)*Float). Sometimes
<br>
that will give 0, sometimes a finite number, and
<br>
sometimes ±infinity.
<br>
It's all doable by careful case analysis.
<br>
<br>
sin(X)
<br>
See comment on cos(X).
<br>
sinh(X)
<br>
(exp(X) - exp(-X))/2.
<br>
If bignum X > 0, -> exp(X)/2 -> +infinity.
<br>
If bignum X < 0, -> -exp(-X)/2 -> -infinity.
<br>
<br>
sqrt(X)
<br>
Smalltalk actually defines sqrt on integers as integer
<br>
square root, which is a pain in the posterior. I have
<br>
wanted the integer square root function occasionally,
<br>
true. But I have much more often wanted
<br>
anInteger sqrt -> anInteger asFloat sqrt, which is what
<br>
Erlang does, hooray, hooray.
<br>
If bignum X < 0, this is an error.
<br>
Let {F, S} = frexp(X).
<br>
Let {G, T} = if S odd -> {sqrt(F*2.0), (S-1) div 2}
<br>
; S even -> (sqrt(F), S div 2)
<br>
end.
<br>
Answer ldexp(G, T).
<br>
This will give you good answers up to 3.2x10**616
<br>
and then infinity.
<br>
<br>
tan(X)
<br>
See comment on cos.
<br>
tanh(X)
<br>
(exp(X) - exp(-X))/(exp(X) + exp(-X))
<br>
For any X > about 709.78, answer +infinity.
<br>
For any X < about -709.78, answer -infinity.
<br>
<br>
erf(X)
<br>
From the graph of erf at <a class="moz-txt-link-freetext" href="http://mathworld.wolfram.com/Erf.html">http://mathworld.wolfram.com/Erf.html</a>
<br>
it is pretty obvious that
<br>
For bignum X > 0, erf(X) -> 1.0 and erfc(X) -> 0.0
<br>
For bignum X < 0, erf(X) -> -1.0 and erfc(X) -> 2.0
<br>
<br>
pi()
<br>
has no argument, so no bignum issue.
<br>
<br>
<br>
I presume pow(X, Y, M) computes (X mod M)**Y mod M?
<br>
This reduces to multiplication modulo M
<br>
(is that '*'(X, Y, M) ?) which is not problematic.
<br>
<br>
The core here is two operations that have been in the C library
<br>
(and in many languages under various guises) but are not in the
<br>
Erlang library: frexp and ldexp. I have a faint recollection
<br>
of libgmp offering frexp, if not I could donate my C code.
<br>
<br>
Of course many of these things are in other languages with
<br>
bignums. (Hands up everyone who realised I was staring at
<br>
page 308 of CLtL2? Which doesn't actually mention bignums.)
<br>
<br>
By the way, my HoD and HR are growling at me for not publishing
<br>
enough, and a friend this afternoon said "don't be ashamed of
<br>
publishing simple things". Could I get a paper out of this
<br>
if I finished the case analyses and wrapped a few more words
<br>
around it? (Seriously, I am quite stressed about the situation.)
<br>
<br>
_______________________________________________
<br>
erlang-questions mailing list
<br>
<a class="moz-txt-link-abbreviated" href="mailto:erlang-questions@erlang.org">erlang-questions@erlang.org</a>
<br>
<a class="moz-txt-link-freetext" href="http://erlang.org/mailman/listinfo/erlang-questions">http://erlang.org/mailman/listinfo/erlang-questions</a>
<br>
</blockquote>
<br>
<tt>We do have floor/ceil now (from
<a class="moz-txt-link-freetext" href="https://github.com/erlang/otp/pull/1145">https://github.com/erlang/otp/pull/1145</a>) when the source code is
released in Erlang/OTP 20.0 and these functions are guard
functions, though that means they both exist in the erlang module,
not the math module. </tt><code>crypto:mod_pow/3 exists, but the
math module would want its own implementation to avoid a
dependency on crypto. libgmp has pow functions for bigints though
the situation gets more complex when handling floating point usage
with bigints, as you have described.<br>
<br>
Thank you for identifying the missing pieces and the specific
problems!<br>
</code>
</body>
</html>