in reply to
RFC: Getting Started with PDL (the Perl Data Language)
I think this is great. I do have a question/suggestion though. I have just been looking at the computational accuracy of the c language, specifically how accurate an angle can you precisely work with. It turns out that there is a 15 decimal point limit to PI, even if you declare it as a long double. So I was wondering, wouldn't it be good to mention the computational accuracy of PDL. For instance, can it use PI to 50 decimal places and still be accurate? Can a PDL data type hold huge numbers, or can piddles be a Math::BigInt, etc. I've been reading that Fortran (on which PDL is based) is superior to c in this regard, but is PDL?
Re^2: RFC: Getting Started with PDL (the Perl Data Language) by lin0 (Curate) on Feb 04, 2007 at 07:47 UTC 
Hi zentara
The core of PDL is written in the C programming language (see PDL::Internals). Fortran is only used in some libraries such as PGPLOT (used for plotting) and Slatec (used for manipulating matrices, calculate Fast Fourier Transforms, fit data using polynomials, interpolate data, etc.). About the piddles' data types, here is a table that will show you the data types already defined (in the file Basic/Core/Types.pm.PL):
PDL Type Real C Type
 
byte unsigned char
short short
ushort unsigned short
long long int
longlong long long int
float float
double double
Note that you can add new types as explained in the PDL documentation. You might also want to read the documentation on PDL::PP to learn how to add your own routines.
In short, I believe the computational accuracy of PDL to be closer to that of C than to the computational accuracy of Fortran.
Cheers,
lin0  [reply] [d/l] 

 [reply] 
Re^2: RFC: Getting Started with PDL (the Perl Data Language) by syphilis (Canon) on Feb 04, 2007 at 14:47 UTC 
It turns out that there is a 15 decimal point limit to PI, even if you declare it as a long double
How does that come about ? 8 bytes are allocated for a double, and 12 bytes for a long double  I would therefore have envisaged that extra precision could be attained using the "long double".
Cheers, Rob  [reply] 

I'm not a c compiler expert, but have a look at double precision pi It may be the difference between computational precision and "%Lf" display precision. I started looking into it, when I questioned the way GLib defined pi as a constant. They defined it out to 50 decimal places, but any use of it in it's long double is limited to 15 decimal places. One of the c gurus in that thread said that c's precision is only gauranteed to 10 decimal places, but with IEEE standards, it's common to see 15. I see 15. Eventually, I found mpfr which lets you set how many digits of precision you want to use. If you can show me a simple c program that computes and displays values out to 50 decimal places, with normal c, I would be greatful. Everything I see truncates it (pi) to 3.(15 decimal places).
For example, in this code, the value of pi is correctly printed as a string on the first line of output, but the subsequent lines all have garbage after the 15th decimal place.
#include <gtk/gtk.h>
/* gcc o test test.c `pkgconfig cflags libs gtk+2.0` */
int main (){
/* how the headers define G_PI */
/* #define G_PI 3.1415926535897932384626433832795028841971693993751
+
*/
long double PI = G_PI;
g_print("3.1415926535897932384626433832795028841971693993751\n");
g_print("%0.50e\n",G_PI);
g_print("%0.50Lf\n",G_PI);
g_print("%0.50Lf\n",PI);
g_print("%0.50Lg\n",PI);
g_print("%0.50Le\n",PI);
return 0;
}
Output:
^^^^^^^^^^^ > unstable after 15 dec
+imal
3.1415926535897932384626433832795028841971693993751
3.14159265358979311599796346854418516159057617187500e+00
0.00000000000000000000000000000000000000000000000000
3.14159265358979311599796346854418516159057617187500
3.141592653589793115997963468544185161590576171875
3.14159265358979311599796346854418516159057617187500e+00
^^^^^^^^^^^^ > unstable after 15 dec
+imal
 [reply] [d/l] [select] 

Eventually, I found mpfr which lets you set how many digits of precision you want to use.
Then there's really no need to look further  it's an excellent library. Incidentally, I placed a module (Math::MPFR) on CPAN that wraps the documented mpfr functions. I think I'm the only person that uses it. If you ever take a look at it, I'd be interested on any feedback. I believe it functions well, but there are aspects which bother me  eg, the documentation, the naming of the functions, and the fact that I wrote all of the XS functions "longhand" (as opposed to using a typemap).
Anyway ... back to the consideration of long doubles. Recent versions of mpfr allow you to convert directly from an mpfr_t to a long double, so one could do something like the following:
mpfr_t mpfr_pi;
long double ld_pi;
mpfr_set_default_prec(sizeof(long double) * 8);
mpfr_init(mpfr_pi);
mpfr_const_pi(mpfr_pi, GMP_RNDN);
ld_pi = mpfr_get_ld(mpfr_pi, GMP_RNDN);
I would expect that ld_pi and mpfr_pi would contain the same value, though there's no guarantee that I'm correct. I was actually hoping to have a bit of a play with this at work tonight, but a last minute breakdown made a complete mess of that plan. (I may yet take a closer look and send you an /msg, as it's not really ontopic for either this forum or this particular thread).
Perhaps long doubles don't provide the precision that I expect  I've never checked. Or perhaps it's just an issue with printf or a general problem with C (as you suggested). Or perhaps it's just an issue with Glib. It would be interesting to know ....
Cheers, Rob  [reply] [d/l] 



