Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl Monk, Perl Meditation
 
PerlMonks  

Portability of floor(log(N))

by toolic (Bishop)
on Mar 12, 2018 at 16:48 UTC ( [id://1210737]=perlquestion: print w/replies, xml ) Need Help??

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

I uploaded the latest version of Number::FormatEng (0.03) to CPAN four months ago, and all tests passed until last week. I am trying to understand why a test failed.

Only one test fails, and it happens on only one os/archname (OpenBSD.amd64-openbsd-thread-multi-ld) and on these two versions of Perl: 5.24.3 and 5.26.1. There are many passes (606) across all OS's.

The failing test is:

# Failed test at t/format_pref.t line 54. # got: '1000e-36' # expected: '1e-33'
This is when the test calls format_pref(1e-33). 1000e-36 is the same numeric value as 1e-33, but the point of this function is to display the number with a specific format.

I can't reproduce this because I don't have this OS or Perl version. But, from the tester report, I assume it fails on this line in FormatEng.pm:

my $e = floor( log($num) / log(1000) );
where floor is from POSIX.

When the test passes, $num=1e-33, and $e is -11. When the test fails, I assume $e is -12. Perhaps log(1e-33)/log(1000) evaluates to -11.0001, or -11.1, or -11.2, etc.

Is it possible for me to fix the code to make it more portable? Or, is there a known problem with this OS for these Perl versions?

Replies are listed 'Best First'.
Re: Portability of floor(log(N))
by pryrt (Abbot) on Mar 12, 2018 at 17:45 UTC

    Regarding how to look at the failing systems = sometimes, you can email the owner of the machine that gave you the failing test, and they are willing to manually run a debug version of your test suite and give you the results. When I haven't wanted to pester the owner, I've occasionally just burned a revision (or used an alpha revision) and submitted a debug revision to CPAN. (I've had problems with going to the next real rev after an alpha, which is why I've sometimes burned a normal revision). When I do that, I make sure everything will still pass, it will just give you extra debug information.

    For debugging, I would use some Test::More::diag() calls, especially looking at intermediate values, such as diag($num); diag(log($num)), diag(log($num)/log(1000)). You also might want to look at the hex representations, using the modern sprintf format %a1, which gives the hexadecimal-floating-point representation, or my Data::IEEE754::Tools's convertToHexString(), which will give a similar representation on systems that don't have that format.2

    My guess for the culprit was to look at the NVSIZE (use config; diag(nvsize => $Config{nvsize}); diag(nvtype => $Config{nvtype})) and some of the other types and sizes on the failing system(s)... I was wondering whether on that specific system, it was using a binary32 instead of binary64 for the NVTYPE (NVSIZE=4 instead of NVSIZE=8); or maybe it's gone the other way, and the failing systems are using binary128 (NVSIZE=16); either way, the change in precision can sometimes change which side of "exact" the binary floating-point representation lands on, which might change whether log(1e-33) would be slightly less or slightly greater than and the extra precision is causing the roundoff to be different.

    1: For binary32, %a would want %.6a, binary64 would want %.14a, and binary128 would want %.29a

    2: Currently, Data::IEEE754::Tools will coerce NVSIZE=16 (binary128) into NVSIZSE=8 (binary64), so it won't give the exact representation on a binary128 system. Given you're failing on 5.24 and 5.26, which are new enough to have %a, I'd probably recommend %a.

      I'd probably recommend %a

      Good advice - and the value I'd be looking at is $num itself.
      Perl is notorious for assigning values incorrectly - and I have a number of -Duselongdouble builds of perl (Windows and Linux) which assign the value 1e-33 incorrectly as 0xa.6274bbdd0fadd61p-113, instead of (correctly) 0xa.6274bbdd0fadd62p-113.
      This doesn't cause any test failures for me. It underestimates the value of '1e-33' by a measly 1 ULP, and I think we're looking for an overestimate for the test to fail.
      (I also don't discount the possibility that the error is in the calculation.)

      If the problem *is* in the value that perl assigns for '1e-33' then your solution is to make sure that a correct value is being assigned - for which you can use Math::MPFR or (I believe) Data::Float. There are perhaps other modules, too.
      Data::Float is pure perl and therefore makes for a more straightforward installation. I use Math::MPFR, which requires the MPFR C library.

      Afterthought: IIRC, some systems that implement the 'long double' type don't actually implement the full range of math functions for the long double.
      It's therefore worth determining whether the underlying libc for those 2 FAIL reports provides floorl() and logl(), or whether libc just uses the double precision floor() and log() functions.

      Cheers,
      Rob
      Thank you for the suggestions.

      I agree with you that, before I pester the tester, I should try to have a pretty good plan. Using diag for all the intermediate values is a great idea. I was unaware of %a. Now all I have to do is play with it until I have a grasp at what the hex fp format means. There's a ton I don't understand.

      I like the approach of releasing a new debug version of the module, too. However, in this case, it may not be practical since it could be months before I get any feedback (it did take 4 months to land on this one machine). I'd feel like an astronaut in one of those sci-fi movies transmitting a message back to earth and waiting 6 months for a response :)

      This is all more involved than I was expecting (and hoping). Creating portable software is hard! Based on your reply and that of syphilis as well, I'm struggling to decide whether my code has a fundamental flaw or whether this one machine is just badly broken in some way. Is my code bound to fail on other machines given enough variety in input values? Since the input is the set of all real numbers, it's impossible to exhaustively test it. I should try to test more thoroughly (but intelligently).

      Or, is there a way to change my code to avoid this error? Can I pre-process the value, say with something as simple as sprintf? Or, do I need to use one of syphilis's suggestions, like Data::Float? Or, can I post-process the value? And, if I do make this type of change, will it break backward-compatibility in a significant way? Or, will it be minor enough to justify the change? Lot's of questions for me to mull over.

        I guess I'd say the next steps depend on whether it bothers you (or your end users) to get 1000e-36 instead of 1e-33. If it does, post-processing is probably your best bet -- but I wouldn't use sprintf. I'd look at the number of digits before the "e" (or the decimal point), and if it's four digits, I'd divide by 1000 (and re-round, if necessary) and adjust the exponent. syphilis's suggestion of generating the 1e-33 using Data::Float's hex_float() would work to "fix" your test -- make the test pass; but if you or your end user would consider it wrong to ever see 1000e-36, then forcing a "correct" 1e-33 won't show the actual error.

        I was able to make it show 1000e-36 by running

        C:\WINDOWS\system32>perl -MNumber::FormatEng=format_pref -le "print fo +rmat_pref($_) for 0.9999999994e-33, 0.9999999995e-33, 0.9999999996e-3 +3" 999.999999e-36 999.999999e-36 1000e-36
        . Your code $num = 0 + sprintf '%.6f', ($num / $mult); only allows 6 digits after the decimal, and you don't do any more checks after that, so numbers that round from 999.999999x to 1000.000000 are going to show up wrong. In your library, you need to fix that. In your test suite, you need to make sure you test the edge cases of what your code does with that (for example, the %.6f introduces edge cases) -- you don't have to test every representable floating value in binary32, binary64, and binary128, just the edge cases that your specific code introduces. I'd probably do a simple test after that line that if(abs($num)>=1000) { ... $num /= 1000, $e++ } or similar. <code>
        Thank you for the suggestions.

        I think they were completely irrelevant to the problem you're facing.
        Silly me didn't think to actually check what the correct calculation of floor( log($num) / log(1000) ) should yield under 'long double' (64-bit) precision.

        When I do check, I find that the correct answer is -12.
        The error is in your test script ... either that, or in Math::MPFR:
        use strict; use warnings; use Math::MPFR qw(:mpfr); # Set precision for extended # precision 'long double' Rmpfr_set_default_prec(64); my $num = Math::MPFR->new('1e-33'); my $thou = Math::MPFR->new('1000'); # Convert $num & $thou to their # respective log() Rmpfr_log($num, $num, MPFR_RNDN); Rmpfr_log($thou, $thou, MPFR_RNDN); print $num / $thou; # prints -1.10000000000000000009e1
        This indicates that when the floating point precision is 64 bits, the correct floor() result is -12.
        And, IIUC, the test script should accept '1000e-36' as correct when NV precision is 64 bits (as is the case with those 2 FAIL reports).

        Interestingly, with full quad precision of 113 bits, the calculation reverts to returning -11.

        Cheers,
        Rob

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://1210737]
Approved by marto
Front-paged by davies
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others scrutinizing the Monastery: (5)
As of 2024-03-19 11:52 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found