gamma function: add constants for 128-bit reals.

[redstar]

No description provided.

[tsbockman]

Is there a platform now where we can actually test D code with 128-bit reals?

[redstar]

I am working on an AArch64 port of LDC. This platform uses 128-bit reals.

[tsbockman]

Cool.

There has been some discussion by maintainers of std.math about doing, essentially, a re-write to cleanse it of unnecessary C-isms, over-use of inline assembler, etc. I had hoped this could be done before adding 128-bit real and cent support, to avoid duplication of effort.

But, on the other hand, already having a working 128-bit platform to test against would be really helpful...

[redstar]

It seems easy to add support for 128-bit reals. Only this PR, #4036 and an implementation of exp() are required.

[tsbockman]

At a minimum, there are many other places in std.math where 80-bit constants need to be upgraded to 128-bits (assuming you want the extra precision to actually be useful, anyway).

Also, a number of functions in std.math only really support real: using 128-bit floating-point operations for 32- and 64-bit values is likely to be unacceptably slow. Overloads will need to be created.

I also expect that once we're able to actually test the 128-bit code paths in std.math, problems will be revealed; I have personally fixed at least one such issue (which I only noticed because the code didn't quite work right for 80-bit, either).

Having said all that, a lot of my concern over 128-bit real actually has do with my previous investigation of what it would take to get imbExtended working. That's a far less sane format than the ieeeQuadruple one that you're working on, though, so I don't think you should be too worried.

[redstar]

My first goal is to compile all modules on AArch64. Without it I cannot run the test suite automatically (which will certainly reveal tons of other errors).

[tsbockman]

My first goal is to compile all modules on AArch64. Without it I cannot run the test suite automatically (which will certainly reveal tons of other errors).

Makes sense. Just remember that the std.math test suite, itself, will likely be somewhat broken until the constants have been expanded for 128-bit some places.

[DmitryOlshansky]

Can't really check the constants.. Anyhow LGTM
As this is related to math - @9il ?

[9il]

Can't really check the constants..

We should add sqrt2p constant to std.math, thought (hex representation from wolfram alpha is required).

[9il]

@redstar Could you please explain what is this constants and how they was computed?

[wilzbach]

According to wikipedia and the gnu's libquadmath the max value with quadruple precision is
1.18973149535723176508575932662800702e4932, thus:

11356.52340629414394949193107797076489126

afaict the current representation has less digits?

real.min is 3.36210314311209350626267781732175260e-4932 and thus

-11355.13711193302405887309661372784853821

[tsbockman]

@wilzbach

afaict the current representation has less digits?

That's because those are in hexadecimal.

[wilzbach]

@tsbockman yep I know, I tried to convert the real string, but obviously failed with low precision on my 64bit system.
Hence I wanted to point out that results with wolfram alpha differ, btw for max it gives me:

0x2c5c.85fdf473de6af278ece600fcbdab54a087542d48a06307294043f0285fee8215026a32a8c36b4f6906fde3d9c2c481034dfb0912176b3604e2b217c387b932e90332ca820b965b56bf4b36fd61860709a2b353606c41270a00931efe794d0e24013d67dd3fdb5125d25a5d48a125dcef5f0f23479ae4fdb6bb8893c201

[jpf91]

I just reverse engineered this, there are obviously too many things you can do wrong when calculating this value...

• Make sure log is the natural logarithm!
• real.max is really the max value, but you need the exact value (e.g. in binary representation, not some decimal approximation). For 64 bit: 2^1023·(1 + (1 − 2^−52)). Note: you can get these from the D compiler: writefln!"%a"(real.max or real.man_dig or real.epsilon); the returned value is a bit exact result. (But if there's a calculation involved, the calculation is of course affected by the FPUs range and rounding). Remember to convert for wolfram alpha: 0xf.fffffffffffffffp+16380 ==> (f.fffffffffffffff base 16) * 2^16380
• The hex format is X.XXXX * 2^y, so mantissa X needs to be base 16 and exponent y needs to be base 2! This is not the wolfram alpha base 16 format!
• You cannot simply calculate this with standard math functions in C/D/... because of limited precision!
• Calculate base 2 representation in wolfram alpha, then convert the base 2 mantissa to base 16
  64 bit: http://www.wolframalpha.com/input/?i=ln(2%5E1023%C2%B7(1+%2B+(1+%E2%88%92+2%5E%E2%88%9252)))+in+base+2
  80 bit: http://www.wolframalpha.com/input/?i=ln((f.fffffffffffffff+base+16)+*+2%5E16380)+in+base+2
  128 bit: http://www.wolframalpha.com/input/?i=ln((1.ffffffffffffffffffffffffffff+base+16)+*+2%5E16383)+in+base+2

void main()
{
    import std.algorithm, std.stdio, std.conv;
    string mantissa = "1.01100010111001000010111111101111101000111001111011110011010100111001001111000111011001110011000000000111111001011101110101011";
    write(mantissa[0 .. 2]); // always 1.
    mantissa = mantissa[2 .. $];
    for (size_t i = 0; i < mantissa.length/4; i++)
    {
        writef!"%x"(to!ubyte(mantissa[0..4], 2));
        mantissa = mantissa[4 .. $];
    }
    writeln();
}

To summarize: I hate the internal/math modules. It delayed the ARM port years ago, now it delayed the AArch64 port for 2 years. Nobody maintains this code, nobody knows what it's supposed to do, the referenced original C code doesn't match this code, so it isn't useful either.... I doubt this code is even correct. For example, in gammafunction.d we have MAXLOG for extended/double format. But in errorfunction.d we use only the extended MAXLOG for all float types. Is this correct? I guess not, I don't know and I don't know who could know. The whole math module situation is a mess. ping @andralex I think we need some long term solution for these modules....

For now, as this is the last thing stopping GDC (and probably LDC) from compiling phobos on AArch64 I'll post a revived version of this pull request tomorrow.

[wilzbach]

Ping @redstar

[wilzbach]

http://www.wolframalpha.com/input/?i=sqrt(2pi)

2.5066282746310005024157652848110452798109E0L
2.506628274631000502415765284811045253006986740609938316629

(latter is from Wolfram Alpha)

[9il]

wolfram allows to get hex representation, see constants in std.math

[wilzbach]

Oh nice :)
http://www.wolframalpha.com/input/?i=sqrt(2pi)+in+base+16
Cut wherever you want:

0x2.81b263fec4e0b2caf9483f5ce459dc5f19f3ea6416000b50dc2f412ddeeb948b5068337b65698726847d560b8818fb98307f0f37d42a17c1c2b86baec9ffa

[JackStouffer]

Closing due to author inactivity. @redstar please ping me if you want to reopen.

[jpf91]

New pull request: #5947