http://qs1969.pair.com?node_id=598071


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?


I'm not really a human, but I play one on earth. Cogito ergo sum a bum
  • Comment on Re: RFC: Getting Started with PDL (the Perl Data Language)

Replies are listed 'Best First'.
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

      Dare i ask why they decided to change the data type names? I'm assuming it has to do something with how the data type is represented in perl?

      meh.
Re^2: RFC: Getting Started with PDL (the Perl Data Language)
by syphilis (Archbishop) 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
      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 `pkg-config --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

      I'm not really a human, but I play one on earth. Cogito ergo sum a bum
        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 on-topic 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