5mi11er has asked for the
wisdom of the Perl Monks concerning the following question:
Good Monks,
I've been puzzling over a question I've posed for myself involving floating point operations. My specific question was originally: given the number PI, if I plan to raise it to some power y, how many digits to the right of the decimal point would I need to ensure an answer that was correct to the least significant integer? Given small y, this is easy to play with in Excel, but I'm looking for the more general area of knowledge and how one would attack this problem for other numerical data and criteria. I am sure there is some branch of mathematics devoted to this type of analysis, but I don't even know what to search for. I've browsed a few of the large math web sites, but other than the extremely broad area of numerical analysis, I'm at a loss as to where else to head. Scott
Re: OT: How much float precision needed for operation? by Enlil (Parson) on Sep 02, 2005 at 19:11 UTC 
 [reply] 

Yes, indeed, significant digits is certainly one of the terms I need. I couldn't for the life of me remember that term, nor the closely related 'margin of error' when I first wrote my original post... Googling with "significant digits" hasn't yet led me to how to answer the question.I did find this link to a pdf discussing estimation of required significant digits, and I think it very likely applicable to my question, but the paper is far enough above my head that I'm unable to use that link without further assistance. Scott Update: Ack, I hadn't seen trammell's post below when I posted this, I'm studying it now...
 [reply] 
Re: OT: How much float precision needed for operation? by trammell (Priest) on Sep 02, 2005 at 20:11 UTC 
#!perl l
# Looking at function F = pi ^ y
#
# F = pi ^ y
# error = dF = y * pi ^ (y1) * dpi
# => dpi = error / ( y * x ^ (y1) )
#
# So if we want to limit the error in F to 0.1:
use constant ERROR => 0.1;
use constant PI => 4 * atan2(1,1);
warn "Check: pi is ", PI, "\n";
my @y = (1..9, map(10 * $_, 1..9), map(100*$_,1..10));
foreach my $y (@y) {
printf "y=%d dpi=%g\n", $y, dpi($y);
}
sub dpi {
my $y = $_[0];
my $denom = $y * PI ** ($y1);
return ERROR / $denom;
}
__END__
Check: pi is 3.14159265358979
y=1 dpi=0.1
y=2 dpi=0.0159155
y=3 dpi=0.00337737
y=4 dpi=0.000806288
y=5 dpi=0.00020532
y=6 dpi=5.44627e05
y=7 dpi=1.48594e05
y=8 dpi=4.13867e06
y=9 dpi=1.171e06
y=10 dpi=3.35468e07
y=20 dpi=1.79111e12
y=30 dpi=1.27507e17
y=40 dpi=1.02116e22
y=50 dpi=8.72341e28
y=60 dpi=7.76258e33
y=70 dpi=7.10495e38
y=80 dpi=6.6385e43
y=90 dpi=6.30114e48
y=100 dpi=6.05568e53
y=200 dpi=5.8364e103
y=300 dpi=7.5001e153
y=400 dpi=1.08428e202
y=500 dpi=1.67203e252
y=600 dpi=2.68581e302
y=700 dpi=0
y=800 dpi=0
y=900 dpi=0
y=1000 dpi=0
Update: yes, for example if we're calculating pi^y for @y=5..15, and we introduce an error to pi of 1e7:
#!perl l
use strict;
use warnings;
use constant PI => 4 * atan2(1,1);
use constant DPI => 1e7;
my @y = (5..15);
foreach my $y (@y) {
my $f1 = F(PI,$y);
my $f2 = F(PI+DPI,$y);
my $delta = $f2  $f1;
printf "y=%2d delta=%g\n", $y, $delta;
}
# F = pi ^ y
sub F {
my ($pi, $y) = @_;
return $pi ** $y,
}
__END__
y=5 delta=4.87045e05
y=6 delta=0.000183612
y=7 delta=0.000672972
y=8 delta=0.00241623
y=9 delta=0.00853968
y=10 delta=0.0298091
y=11 delta=0.103013
y=12 delta=0.353045
y=13 delta=1.20155
y=14 delta=4.06515
y=15 delta=13.6833
So an error of 1e7 in pi doesn't affect the output of F by more than 0.1 until we hit y=11.  [reply] [d/l] [select] 

So, I think what your data says is: To ensure an error of no more than plus/minus 0.1 we'd need ~7 digits to the right of the decimal point for y=10 and ~53 digits for y=100 Is this correct?Scott Update: trammel, thank you very much for your excellent examples above. So, for others who want to extend this to other functions: If we take the derivative of the function we're working with, multiplied by a delta (difference), and know what our tolerance for error is, we can solve for the maximum delta that will be within our tolerance for error.In math symbols, the process is:  given function F
 find the derative > dF
 maxdelta = errortolerance/dF
 [reply] 
Re: OT: How much float precision needed for operation? by spiritway (Vicar) on Sep 03, 2005 at 05:13 UTC 
You can figure it out by breaking your number into two parts, x (the true value of PI), and y (the error part). y is going to be some smallish figure, expressed as, say, 1e10, or whatever. Your problem then becomes calculating the difference between x**n, and (x+y)**n; you can estimate the error in this way. For n=2, you'd get x**2 + 2xy + y**2. Plug in the value for y to get your error. The error would be about 6.8e10 + 1e20.
You might consider using fractions instead of floats. If you are able to obtain a fraction that approximates PI, you can play around with exponents more safely. Using bigint, you could do this without concern for overflow, underflow, or other problems related to floats. Just a thought.
 [reply] [d/l] [select] 
Re: OT: How much float precision needed for operation? by shemp (Deacon) on Sep 02, 2005 at 22:59 UTC 
One approach that could work would be to do error analysis as if you were approximating an infinite series. This is often covered in a secondsemester calculus course (in the US), or at least touched on. Its been many years since i've actually done anything like this, but it was the first approach that came to mind.
In general, the idea is that you determine an reasonable upper bound that including another term (or digit in your case, which can be considered a term) will modify your answer by.
I will introduce a nonstandard notation here:
Let PI(n) denote an approximation of the constant PI, truncated after n decimal places.
For example if you're dealing with PI, and have some approximation that has 'n' decimal places (n digits after the decimal), you know that including the next digit will affect the value of PI as:
Error <= PI(n+1)  PI(n) < (1/10)^n
So if you consider PI^2, you can do set something up like this:
Error <= [ PI(n+1)  PI(n) ]^2 < ((1/10)^n)^2
I think that doing that analysis directly is not exactly how you want to proceed, but i can't remember exactly what to do.
When i get home from work i'll have my math books available, instead of only these computer books. If i find anything more concrete, i'll post it from home.
I guess i can include this. If you find a series error analysis technique, you could use this series for PI:
PI = sum( i = 0 .. INF ) (1 ^ i) / (2 * i + 1)
although that converges extremely slowly, and there are other better series that i cant think of right now.
I use the most powerful debugger available: print!
 [reply] [d/l] [select] 

