Beefy Boxes and Bandwidth Generously Provided by pair Networks
There's more than one way to do things
 
PerlMonks  

Fix floats like you do in your head

by tachyon (Chancellor)
on Dec 24, 2002 at 02:24 UTC ( [id://222013]=CUFP: print w/replies, xml ) Need Help??

I have just completed work on Math::Trig::Units with Math::Trig::Degree, Math::Trig::Radian and Math::Trig::Gradian subclasses. Something that has always anoyed me is how float operations return 0.499999999998 or 5.00000000001 when the actual value is 0.5. The approx sub was written for this module.

Because of the limit on the accuracy of the vaule of Pi that is easily supported via a float you will get values like dsin(30) = 0.49999999999999945 when using degrees (or gradians). This can be fixed using the approx() function.

By default the approx sub will modify numbers so if we have a number like 0.499999945 with 6 9s or 0.50000012 with 6 0s the number will be rounded to 0.5. It also works on numbers like 5.250000001. This is useful when using degrees or gradians. In degrees these functions will return 0.5 as expected

approx(dsin(30)) approx(dcos(30))
The approx sub takes a second optional argument that specifies how many 0s or 9s in a row will trigger rounding. The default is 6.
approx($num, 7); # will return 0.5 for 0.500000001 but 0.50000001 + if # that is passed as it only has 6 zeros.

Numbers that do not fulfill the requisite criteria are returned unchanged. For example 0.5000001 will not be rounded to 0.5 as it only has 5 0s.

sub approx { my ( $num, $dp ) = @_; $dp ||= 6; if ( $num =~ m/\d*\.(\d*?)(9{$dp,})\d*/ ) { my $exp = 10** (length $1 + length $2); return int(($num * $exp) +1 )/$exp; } elsif ( $num =~ m/\d*\.(\d*?)(0{$dp,})\d*/ ) { my $exp = 10** (length $1 + length $2); return int($num * $exp)/$exp; } else { return $num; } }

Replies are listed 'Best First'.
Re: Fix floats like you do in your head
by mojotoad (Monsignor) on Dec 24, 2002 at 21:32 UTC
    This is a clever way to deal with the problem, but not very consistent (in that some numbers are treated differently than others) -- this inconsistency will make most mathematicians and scientists quake in their boots.

    The proper solution is to define an epsilon that represents the precision used to compensate for accumulated rounding errors.

    Indeed, this is the driving force behind the IEEE 754 standard.

    Also of interest: What Every Computer Scientist Should Know About Floating-Point Arithmetic (1991) by David Golberg, is a great discussion of the IEEE 754 standard and its practical application.

    Matt

      Great links. As you say this is not a ceil() or a round() function which behaves in strict mathematical terms, nor does it attempt to solve any of the problems associated with doing decimal fp arithmetic in binary.

      The beauty/purpose of it is that is does what a human tends to do. I see 0.9999999 and think 1 or 0.00999999 and think 0.01 which is what you get with the approx() function. I did call it approx() rather than exact(), actually() or really() ;-). I had it in mind for simple display functions not anything serious.

      It also only works on the decimal fraction but could be extended to work on the integer component so 999999.9 or 1000000.1 would become 1,000,000.

      cheers

      tachyon

      s&&rsenoyhcatreve&&&s&n.+t&"$'$`$\"$\&"&ee&&y&srve&&d&&print

        No worries. After I went and perused your Math::Trig::Units work, it became evident that this wasn't an oversight on your part. Ah well, the links are there for posterity, in case anyone else wants to learn about it.

        As for names...approx() works. How about fudge()?

        At any rate, at least you know that we're dealing with a psychological tendency, here. With 6 sig figs, for example, 1.99999 is no more special than 1.46793 on a mathematical basis. On a psychological basis it is somewhat irritating because of our preference for nice, neat numbers.

        Matt

Re: Fix floats like you do in your head
by toma (Vicar) on Dec 26, 2002 at 17:52 UTC
    Numerical methods are deceptively difficult. They require exquisite attention to detail to make them work properly. For example, if you are getting dsin(30)=0.49999999999999945 then you have a serious bug.

    For most hardware out there, you cannot properly initialize a variable such as pi by entering lots of digits in your program. This is because the internal precision of the machine exceeds the precision that you are normally allowed to print or initialize. So if you catch yourself typing in lots of digits of pi, you are doing something wrong. Here is a better way to do it:

    use Math::Trig; my $pi= 4*atan2(1,1); my $deg_to_radian= $pi/180; print sin(30*$deg_to_radian),"\n";
    This prints 0.5, as you would hope that it would.

    As far as the snippet goes, it is seriously borked. Among other problems, it has problems with the implicit string/numeric conversion that perl uses. This causes it to give different results for

    approx(".00000000001")
    and
    approx(".00000000001"+0)

    The recent CPAN module Math::SnapTo version 0.02 also has serious errors. For good rounding, use Math::Round. It doesn't do exactly what you want (that is, round like a human), but neither does this snippet or Math::SnapTo.

    It should work perfectly the first time! - toma

      Yeah, I agree it is seriously borked. As you note Math::Round does a better job than Math::SnapTo so I have scheduled it for deletion :-(

      That said the pi thing is interesting. As your snippet shows using the value of pi from atan2 provides an accurate result. You are not actually using Math::Trig BTW as it has no sin.

      Theoretically speaking as it is impossible to express pi exactly it is also impossible to express 30 degrees in radians exactly so by definition you should not be able to get 0.5 unless there is internal rounding going on.

      As a demonstration of the vagaries of fp decimal math on binary based systems dsin(30) from Math::Trig::Units-0.02 with the hard coded pi returns 0.5 on the system I am on now, whereas it returns the 0.4999..... value on a similar win2k (obvoiusly not identical). Interestingly that system has and Athlon CPU whereas this one has PIII Xeons.

      It appears as though Math::Trig::Units will become Math::Trig in the near future. I have reverted to using atan2(1,1)*4 to get a value for pi (as Math::Trig does) and deleted the kludge subs from it.

      tachyon

      s&&rsenoyhcatreve&&&s&n.+t&"$'$`$\"$\&"&ee&&y&srve&&d&&print

        I agree with the Math::SnapTo decision, it didn't seem consistent with your standard of coding excellence!

        IIRC some ancient Intel hardware used to keep something like 80 bits of internal precision, and this was typically rounded to 64 bit IEEE by most compilers in some but not all operations. In this case, internal precision was a huge win!

        When testing multi-platform floating point test code, it is often convenient to account for rounding differences on different platforms. A form of diff called ndiff has been created to "compare putatively similar files, ignoring small numeric differences." I think ndiff would be a useful perl module!

        It should work perfectly the first time! - toma

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: CUFP [id://222013]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others surveying the Monastery: (3)
As of 2024-06-22 06:28 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found

    Notices?
    erzuuli‥ 🛈The London Perl and Raku Workshop takes place on 26th Oct 2024. If your company depends on Perl, please consider sponsoring and/or attending.