syphilis has asked for the wisdom of the Perl Monks concerning the following question:

Hi,
The demo script (named 'double.pl'):
use warnings; use strict; use Inline C => Config => BUILD_NOISY => 1; use Inline C => <<'EOC'; SV * foo(SV * nv) { double x = (double)SvNV(nv); return newSVnv(x); } EOC print perl_foo(2.5), "\n"; print perl_foo(1.1), "\n"; # LINE 17 sub perl_foo { my $ret = foo($_[0]); if($_[0] - $ret) {warn "Loss of precision"} # LINE 21 return $ret; }
That script does pretty much what I want. On a perl built without long double support it simply prints out:
2.5 1.1
On a (linux, mdk-9.1) perl built with long double support it prints out (transcribed, not copy'n'pasted):
2.5 Loss of precision at double.pl line 21 1.10000000000000009
The key points to the exercise are:

1) That the warning is emitted for perl_foo(1.1) but not for perl_foo(2.5) - ie I want to see the warning on the -Duselongdouble build of perl only when there's an actual loss of precision;

2) The detection of the loss of precision is being handled by the perl code and not by the C code (over which, in the real life situation, I have no control);

3) The warning should never be seen on a perl that was not built with long double support.

The questions:

1) The warning is reported as coming from line 21. I would prefer to see line 17 reported as the line that generated the warning. How do I go about achieving that ? (I ask that question with a pained expression on my face, for I feel I should know the answer - or be able to find it for myself.)

2) Is simply subtracting the return value from the argument a reliable way of testing whether precision has been lost ? (It seems obvious to me that the answer to that question is "yes" ... but I'm always wary of the validity of answers that seem obvious :-)

3) Is there a better way that a -Duselongdouble build of perl could determine that foo(1.1) would lead to a loss of precision ? (It's only a minor issue, but I would prefer that a -Duselongdouble build of perl could tell in advance that foo(1.1) was going to lead to a loss of precision, rather than having to wait until foo(1.1) had executed.)

Cheers,
Rob

Replies are listed 'Best First'.
Re: XS, C doubles and -Duselongdouble
by BrowserUk (Patriarch) on Jun 10, 2007 at 12:52 UTC

    1. Carp provides for 'misreporting' the location of errors via $Carp::CarpLevel = 1;
    2. Seems right to me too. Take that with exactly as much authority as it deserves :)
    3. Seems to me there two possible ways:

      • It could yell about the possible loss of precision as a C compiler would when you assign a long double to a double.

        In that respect, when the XS code is compiled, the C compiler probably is yelling about it.

        I've looked for a proper way of pursuading XS to display warnings produced by the C compiler without success. I generally resort to adding an obviously fatal error at the end of the inline C--a # in column one of the last line works:). When the fatal error is reported, you also get to see all the warning errors that were previously kept hidden from you.

      • It could attempt to detect whether the number of significant bits set in the mantissa of the long double exceeds 53.

        *I'm not sure whether you want it to be your perl code, the XS-precompiler or the Inline generated wrappper stubs, or the C-code?

        For any given processor, or rather for a given real format, that would be a fairly simple operation. Just requiring a bit mask of the appropriate bits of the long double to check.

        For IEEE representation, you should be able to work out the appropriate bytes/bits to mask from IEEE 754 80-bit Extended double (long double) to 64-bit double unpack.

        Doing it such that it would be portable across platforms that don't use IEEE format reals could be a lot harder.


    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.
      Carp provides for 'misreporting' the location of errors via $Carp::CarpLevel = 1;

      Yep - that does it nicely. (After I posted the sopw I realised that particular question didn't have a real lot to do with "XS", "C doubles", or "-Duselongdouble" :-)

      In that respect, when the XS code is compiled, the C compiler probably is yelling about it

      Haven't noticed any complaints from gcc - even if the SvNV() is not explicitly cast to a double. On windows, with an MSVC-built perl, you aren't going to get a "proper" -Duselongdouble build of perl ... in that any Microsoft compiler that knows about long doubles will tell you that sizeof(long double) == sizeof(double) == 8.
      With MinGW, it should be theoretically possible to build a -Duse64bitint and -Duselongdouble build of perl .... but since MinGW uses the msvcrt.dll runtime library, I don't think it would be possible for such a perl to printf() the long double values correctly.

      I generally resort to adding an obviously fatal error at the end of the inline C

      The idea of the use Inline C => Config => BUILD_NOISY => 1; is that it provides verbosity (including any compiler warnings), thus eliminating the need to take such measures.

      It could attempt to detect whether the number of significant bits set in the mantissa of the long double exceeds 53

      Yes - that's what I was trying to think of. It's just a matter of finding the position of the least significant (set) bit in the argument that's given to perl_foo(). It's the perl code (the 'perl_foo' function to be precise) that needs to determine whether there has been a loss of precision.

      Thanks BrowserUk.

      Cheers,
      Rob
        It's just a matter of finding the position of the least significant (set) bit

        It's probably quicker to just test that there are no bits set in the lowest 8 bits (64-52) of the mantissa. Those 8 bits are stored in the last (10th) byte of the real*10 representation. So it should just be a case of (assuming you have the 10 bytes packed in scalar)

        if( substr $packedReal10, -1 ) { warn 'Loss of precision...'; }
        any Microsoft compiler that knows about long doubles will tell you that sizeof(long double) == sizeof(double) == 8.

        Gah. How stupid is that. You know, I have a funny feeling that I saw a set of 80bit real math functions kicking around inside one of the system DLLs. In theory, it would be possible to wrap a header file and an import library around that and get access to them from C/C++. I can't remember where I (think) I saw them and a quick grep didn't locate them. I'll have another look later.


        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
Re: XS, C doubles and -Duselongdouble
by ysth (Canon) on Jun 11, 2007 at 03:25 UTC
    I'm trying to figure out what your concern might actually be.

    In theory, the NV for the constant 1.1 will be the closest long double to the exact value for 1.1. As far as I can tell, the cast to double rounds to the nearest possible double value, which is almost always going to be the closest double value to the original exact value. So I don't see how you could be worried about those.

    The value returned by foo(), a long double laundered through a double, is almost always going to be different from the original value - given purely random long double input, something like 99.95 % of the time. Though you'll run across cases where no loss occurs more often than that, since many interesting values are exactly representable as a double or long double, e.g. 2.5, or 1.

    So, why does it matter to you that precision may be (rather, is probably going to be) lost? What do you want to do differently based on whether precision is lost?

    Update: maybe all you want is to use sprintf("%.15g", perl_foo($x)) instead of just perl_foo($x)?

      I'm trying to figure out what your concern might actually be

      The GMP C library can assign longs (signed and unsigned), doubles and strings to the mpz_t (integer) struct. But it has no inbuilt way of assigning long longs and long doubles to an mpz_t.

      I've merely been looking at how to deal with some of the gotchas that arise when you interface a -Duse64bitint build of perl (with/without -Duselongdouble) to that GMP library.

      Let's say I have a perl built with 64-bit int support (but no long double) and I've done:
      use warnings; use Math::GMPz qw(:mpz); $num = 144115188075868217;# this is an IV $obj = Math::GMPz->new($num);
      What should happen there ? Should it croak with the message "Hey ... the GMP library doesn't handle ints bigger than 32-bit!!" ? Should it silently assign the value as a 32-bit int ? I've chosen to have the new() XS function read $num as a string using SvPV_nolen(), and then have the GMP library's mpz_set_str() function assign from that string - and that seems to assign the correct 64-bit value.

      Unfortunately, it doesn't appear to be so simple when it comes to NV's and -Duselongdouble builds of perl:
      use warnings; use Math::GMPz qw(:mpz); $num = 2 ** 57 + 12345;# this is an NV $obj = Math::GMPz->new($num);
      If, under a long double build of perl, the XS function now reads $num using SvPV_nolen(), and the value gets assigned using mpz_set_str(), then it gets the value wrong. There is no simple way (that I can see) for $obj to be assigned the value represented by $num - using the existing assignment operations provided by the GMP library. Of course, if $num fits into a C double anyway, then there's no problem. So ... I started to think along the lines that, since the GMP library did not support long doubles, it was unreasonable to expect that Math::GMPz should support long doubles (even though the build of perl did support them) - and that a croak/warning should be emitted iff precision was being lost.

      But I'm now not so sure about that - your mention of sprintf() made me aware that I can probably make use of the C sprintf() function. That is, the new() XS function assigns the $num arg to a long double, sprintf() converts that long double to a string, and *that* string is then handed over to the GMP library's mpz_set_str() function ... and $obj ends up with the correct value and no loss of precision.

      I'll give that a go ... in the meantime, thoughts and comments welcome.

      Thanks ysth.

      Cheers,
      Rob
        It's really difficult to get XS code right that does different things for different types of input; if I were you, I'd document that Math::GMPz->new takes a string of digits, and let the caller worry about anything beyond that.

        Update: or provide a perl new_from_number that looks something like this:

        use Config; sub new_from_number { my $class = shift; my $num = shift; my $fmt = '%.0f'; if (exists($Config{nv_preserves_uv_bits}) && $Config{nv_preserves_uv_bits} < 8 * $Config{uvsize} ) { $fmt = $num < 0 ? '%d' : '%u'; } $class->new(sprintf($fmt, $num)); }
        I'm probably forgetting something important there, though.

        Update 2: yes, I was forgetting to revert to %.0f for numbers outside [IV_MIN, UV_MAX].