A math module better than C89

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

A math module better than C89

Michael Truog
There was a recent pull request at https://github.com/erlang/otp/pull/1298 which was unable to go further, due to it not being clear how bigint support should exist in the math module.  The absence of bigint support has been a long-term annoyance, and the pull request attempted to move custom code from elsewhere in Erlang/OTP to make a custom bigint pow/2 (and pow/3) exposed in the math module.  I like the idea of pursuing the principle of least astonishment, and allowing the existing pow function to handle bigints, similar how you don't have to think about the difference in python (and other programming languages that people enjoy math in).  However, after determining how to do the API, it would then require work to take place to support bigints, which ideally would just be calls to libgmp functions. So, anyone with strong opinions should respond with the hope that the next source code addition might not get refused.

Best Regards,
Michael

_______________________________________________
erlang-questions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlang-questions
Reply | Threaded
Open this post in threaded view
|

Re: A math module better than C89

Richard A. O'Keefe-2
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
[hidden email]
http://erlang.org/mailman/listinfo/erlang-questions
Reply | Threaded
Open this post in threaded view
|

Re: A math module better than C89

Michael Truog
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
[hidden email]
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!

_______________________________________________
erlang-questions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlang-questions
Reply | Threaded
Open this post in threaded view
|

Re: A math module better than C89

Richard A. O'Keefe-2


On 20/01/17 7:00 AM, Michael Truog wrote:

> 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.

That is good to hear.

As a minor documentation issue, because there will be people exposed to
C or Java who expect functions called round, trunc, floor, ceil to be
in the same place as sqrt and acosh, it would be helpful if the
manual page for the math module mentioned these functions, maybe at the
end, and said where to really find them.  I'm not talking about
repeating their descriptions, just pointing to them.




_______________________________________________
erlang-questions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlang-questions
Reply | Threaded
Open this post in threaded view
|

Re: A math module better than C89

Michael Truog
On 01/19/2017 03:30 PM, Richard A. O'Keefe wrote:

>
>
> On 20/01/17 7:00 AM, Michael Truog wrote:
>
>> 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.
>
> That is good to hear.
>
> As a minor documentation issue, because there will be people exposed to
> C or Java who expect functions called round, trunc, floor, ceil to be
> in the same place as sqrt and acosh, it would be helpful if the
> manual page for the math module mentioned these functions, maybe at the
> end, and said where to really find them.  I'm not talking about
> repeating their descriptions, just pointing to them.

I was wrong about this a bit.  There are ceil/floor functions added to the math function that use the same logic as the erlang module ceil/floor functions.  The main difference is that the math module return type is float while the erlang module return type is integer.  So, ceil/floor are present in the math module, though the math module doc probably wants a pointer to the erlang module ceil/floor by mentioning that it returns integer instead of float.

_______________________________________________
erlang-questions mailing list
[hidden email]
http://erlang.org/mailman/listinfo/erlang-questions