install cephes on the ARM platform

[dimpase]

Currently, cephes spkg is only installed on Cygwin. We need it on ARM to take care of similar numerical noise.

Closing this as irrelevant. Using cephes is a wrong way, there is already GSL's gamma(), etc., that can be used.

[Snark]

On newton (x86_64 processor), I have :

jpuydt@newton:~/sage-4.8$ grep install.log -e Success |grep ceph
Successfully installed cephes-2.8

so could you please tell me more about this cephes issue?

I would like to work on those ARM numerical glitches now they're the last standing issues.

[Snark]

I'm trying to track down the numerical glitches related to the lgammal libm(in glibc) implementation, as discussed on the ubuntu glibc  ticket, as fixing the problem upstream looks like a better option in that case.

So far, I think I tracked down the implementation in either in sysdeps/ieee754/ldbl-128/e_lgammal_r.c or sysdeps/ieee754/ldbl-96/e_lgammal_r.c. The problem is that I don't know which is compiled in, as I don't even see those in any Makefile!

[dimpase]

Replying to Snark:

So far, I think I tracked down the implementation in either in sysdeps/ieee754/ldbl-128/e_lgammal_r.c or sysdeps/ieee754/ldbl-96/e_lgammal_r.c. The problem is that I don't know which is compiled in, as I don't even see those in any Makefile!

you can change the source .deb (to print stuff in these functions), install .deb from source, and see which one is invoked.

[dimpase]

Replying to dimpase:

Replying to Snark:

So far, I think I tracked down the implementation in either in sysdeps/ieee754/ldbl-128/e_lgammal_r.c or sysdeps/ieee754/ldbl-96/e_lgammal_r.c. The problem is that I don't know which is compiled in, as I don't even see those in any Makefile!

you can change the source .deb (to print stuff in these functions), install .deb from source, and see which one is invoked.

just in case: as the root, do the following: 1) add the line

deb-src http://ports.ubuntu.com/ubuntu-ports/ oneiric main restricted universe

to /etc/apt/sources.list

2) at the shell prompt

apt-get install devscripts
apt-get build-dep eglibc
apt-get source eglibc
cd eglibc-2.13/

now you can edit the source as you please. Then, in eglibc-2.13/

debuild -us -uc

builds the library creates .deb file(s) in ../. Probably there will be a library to link against in eglibc-2.13/, to test things. Else, one can just

cd ..
dpkg --install eglibc-2.13.deb

and then compile/link as usual.

[dimpase]

Replying to Snark:

I'm trying to track down the numerical glitches related to the lgammal libm(in glibc) implementation, as discussed on the ubuntu glibc  ticket, as fixing the problem upstream looks like a better option in that case.

So far, I think I tracked down the implementation in either in sysdeps/ieee754/ldbl-128/e_lgammal_r.c or sysdeps/ieee754/ldbl-96/e_lgammal_r.c. The problem is that I don't know which is compiled in, as I don't even see those in any Makefile!

no, the trouble is in sysdeps/ieee754/ldbl-64/e_lgammal_r.c It is actually the one that gets called. (long double is 8 bit on armel). You can check this as follows:

/* foo.c */
#include <stdio.h>
#include <math.h>

int
main (int argc, char* argv[])
{
  long double x = 6.0;
  int i;
  printf("lgammal (%.20Lf)=%.20Lf\n", x, lgammal(x));
  printf("__lgamma (%.20Lf)=%.20Lf\n", x, __ieee754_lgamma_r(x,&i));
  printf("tgammal (%.20Lf)=%.20Lf\n", x, tgammal(x));
  printf("__gamma (%.20Lf)=%.20Lf\n", x, __ieee754_gamma_r(x,&i));
  return 0;
}

Compile this with a static linking:

gcc foo.c /usr/lib/arm-linux/gnueabi/libm.a

Running it you see:

lgammal (6.00000000000000000000)=4.78749174278204581157
__lgamma (6.00000000000000000000)=4.78749174278204581157
tgammal (6.00000000000000000000)=119.99999999999997157829
__gamma (6.00000000000000000000)=119.99999999999997157829

[dimpase]

There is a loss of precision occurring in the algorithm used by eglibc to compute gamma(); it does exp(lgamma()), but this is dangerous; running the following on ARM

#include <stdio.h>
#include <math.h>
static long double xxx=4.78749174278204599545L;
int
main (int argc,
      char* argv[])
{
  long double x = 6.0;
  int i;
  printf("lgammal (%.20Lf)=%.20Lf\n", x, lgammal(x));
  printf("hex lgammal (%.20Lf)=%llx\n", x, lgammal(x));
  printf("hex xxx=%llx\n", xxx);
  return 0;
}

shows this:

lgammal (6.00000000000000000000)=4.78749174278204581157
hex lgammal (6.00000000000000000000)=401326643c4479c9
hex xxx=401326643c4479c9

so both "good" (kept in xxx) and "bad" (computed by lgamma) values are actually the same 64-bit fp number.

This gives:

sage: exp(4.78749174278204599545) # "good" lgamma value
120.00000000000000014
sage: exp(4.78749174278204581157) # "bad"  lgamma value
119.99999999999997808

explaining the "mystery" of 119.99999999999997.

[dimpase]

The last comment makes no real sense. In fact, the bug(s) are already in the Compare the outputs of the following on ARM and on x86:

#include <stdio.h>
#include <math.h>
static long double xxx=4.78749174278204599545L;
int main ()
{
  printf("xxx=%.20Lf\n", xxx);
  printf("exp(xxx)=%.20Lf\n", expl(xxx));
  return 0;
}

On ARM:

xxx=4.78749174278204581157
exp(xxx)=119.99999999999997157829

On x86_64:

xxx=4.78749174278204599545
exp(xxx)=120.00000000000000014572

[Snark]

I see two problems :

1. please check how big "long double" is on both boxes -- I'm pretty sure sizeof(long double) is 16 in a case and 8 in the other (so one is quadruple precision, while the test should be on double precision...) ;
2. 53 binary bits give less then 16 decimal digits, so your numbers should get cut on double-precision floats.

[dimpase]

Closing this as irrelevant. Using cephes is a wrong way, there is already GSL's gamma(), etc., that can be used.

[jdemeyer]

I don't understand why you say Cephes is "wrong". I don't see the problem. If it's wrong for ARM, then why isn't it wrong for Cygwin?

And Julien a.k.a. Snark should certainly give his opinion first.

[Snark]

Well, even if I didn't comment here, we discussed on irc ; and the conclusion was that using cephes was probably too complex : why add a new spkg on ARM if we can solve the problem by using an already installed package?

I'm still unsure about how to fix things properly, but I agree this bug report was premature : we should first look at the options and plan a solution before we open a ticket.

[jdemeyer]

Replying to Snark:

Well, even if I didn't comment here, we discussed on irc ; and the conclusion was that using cephes was probably too complex : why add a new spkg on ARM if we can solve the problem by using an already installed package?

Well, if you both agree to close this package, that's fine for me.

we should first look at the options and plan a solution before we open a ticket.

I totally disagree with this. Nothing wrong with opening a ticket and then realizing later the ticket makes no sense.