[erlang-questions] A math module better than C89

Michael Truog mjtruog@REDACTED
Thu Jan 19 19:00:35 CET 2017


On 01/19/2017 12:28 AM, Richard A. O'Keefe wrote:
> 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.)
>
> _______________________________________________
> erlang-questions mailing list
> erlang-questions@REDACTED
> http://erlang.org/mailman/listinfo/erlang-questions

We do have floor/ceil now (from https://github.com/erlang/otp/pull/1145) 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. |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.

Thank you for identifying the missing pieces and the specific problems!
|
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://erlang.org/pipermail/erlang-questions/attachments/20170119/7e95f5f7/attachment.htm>


More information about the erlang-questions mailing list