| Title: | Mathematics at DEC |
| Moderator: | RUSURE::EDP |
| Created: | Mon Feb 03 1986 |
| Last Modified: | Fri Jun 06 1997 |
| Last Successful Update: | Fri Jun 06 1997 |
| Number of topics: | 2083 |
| Total number of notes: | 14613 |
For your recreation,
I tried to program this problem
" In 250 B.C.E the Greek mathematician Archimedes estimated the number
Pie (3.14) as follows, He looked at circle with diameter 1, and hence the
circumference was Pie. Inside the Circle He inscribed a square, the perimeter
of the square is smaller than the circumference of the circle and so it is a
lower bound for Pie, Archimedes then considered and inscribed Polygon, and
produced and even better approximation for Pie, using 96-sided inscribed and
circumscribed polygons, he was able to show that 223/71 < Pie < 22/7 .
there is recursive formula for these estimates.
n
Let P(n) be the PERIMETER of inscribed polygon with 2 sides.
i.e. P(2) is the perimeter of inscribed SQUARE, i.e. P(2)= 2 * SQRT(2) .
in general
---------------------------------
/ ---------------------
n / / n 2
P(n+1) = 2 \/ 2 (1- \/ 1- ( P(n) / 2 ) )
compute p(n) for n=3,4..,60
My problem is that at n=12 the term inside the second SQRT goes to 1,
12
since p(12)/ 2 --> 0 , and so the whole expression goes to 0 !
the following is my 'C' program, I was wondering if someone can show me
a way out of this, ( and show us their program !). may be it is MAPLE time?
I guess the trick is to simplify first and try to produce a formula where
this problem will not occur, I did not get much luck.
Thank You,
/naser
---------------------------------------------------------------------------
#include <stdio.h>
#include <math.h>
main()
{
int i=2;
float PN;
float j;
float p;
char t[1];
PN = (float)2.0 * (sqrt( (double)2.0) );
do
{
TWO(i,&p);
printf("number of sides= %f , permmeter= %f\n",p,PN);
j= PN/p;
j *=j;
j= (float)1.0 -j; /* this becomes 0 and causes the problem*/
j= sqrt(j);
j= (float)1.0 -j;
j= (float)2.0 *j;
j= sqrt(j);
PN= p*j;
i++;
}while (i<60);
}
/******/
TWO(int i,float *p)
{
int j;
int k=1;
for(j=0;j<i;j++)
k= 2*k;
*p= (float)k;
}
$run test
number of sides= 4.000000 , permmeter= 2.828427
number of sides= 8.000000 , permmeter= 3.061467
number of sides= 16.000000 , permmeter= 3.121444
number of sides= 32.000000 , permmeter= 3.136546
number of sides= 64.000000 , permmeter= 3.140333
number of sides= 128.000000 , permmeter= 3.141286
number of sides= 256.000000 , permmeter= 3.141519
number of sides= 512.000000 , permmeter= 3.141208
number of sides= 1024.000000 , permmeter= 3.142451
number of sides= 2048.000000 , permmeter= 3.142451
number of sides= 4096.000000 , permmeter= 3.162278
number of sides= 8192.000000 , permmeter= 3.162278
number of sides= 16384.000000 , permmeter= 2.828427 <-- WRONG
number of sides= 32768.000000 , permmeter= 0.000000
number of sides= 65536.000000 , permmeter= 0.000000
| T.R | Title | User | Personal Name | Date | Lines |
|---|---|---|---|---|---|
| 1299.1 | try this... | ALLVAX::JROTH | It's a bush recording... | Sun Sep 23 1990 02:19 | 92 |
The problem is loss of significance when subtracting nearly equal
quantities.
A solution is to replace the root with a power series expansion
which lets you cancel out the troublesome "1" - see below.
It would probably be more accurate to iterate on a variable
q(n) = (p(n)/2^n)^2, which would eliminate the square rooting and
squaring.
- Jim
#include <stdio.h>
#include <math.h>
#define PI 3.14159265358979323846264338
main()
{
int i=2;
double PN;
double j;
double p;
double s, t;
int k;
PN = 2.0 * sqrt(2.0);
do
{
TWO(i,&p);
printf("sides = %12.0f, perimeter = %22.17f\n", p, PN);
j= PN/p;
j *=j;
#if 0
j= 1.0 -j; /* this becomes 0 and causes the problem*/
j= sqrt(j);
j= 1.0 -j;
#endif
/* 1-sqrt(1-j) = 0.5*j*(1-((0.5-1)/2)*j*(1-((0.5-2)/3)*j*(1-... */
for (k = 2, t = 0.5*j, s = 0.0;
s != s+t;
t *= j*(k-1.5)/k, k++) s += t;
j= 2.0 * s;
j= sqrt(j);
PN= p*j;
i++;
}while (i<32);
}
TWO(int i, double *p)
{
int j;
for (*p = 1.0, j = 0; j < i; j++) *p += *p;
}
$ run pie
sides = 4, perimeter = 2.82842712474619012
sides = 8, perimeter = 3.06146745892071831
sides = 16, perimeter = 3.12144515225805247
sides = 32, perimeter = 3.13654849054593943
sides = 64, perimeter = 3.14033115695475312
sides = 128, perimeter = 3.14127725093277310
sides = 256, perimeter = 3.14151380114430129
sides = 512, perimeter = 3.14157294036709162
sides = 1024, perimeter = 3.14158772527715996
sides = 2048, perimeter = 3.14159142151120030
sides = 4096, perimeter = 3.14159234557011807
sides = 8192, perimeter = 3.14159257658487301
sides = 16384, perimeter = 3.14159263433856328
sides = 32768, perimeter = 3.14159264877698602
sides = 65536, perimeter = 3.14159265238659169
sides = 131072, perimeter = 3.14159265328899318
sides = 262144, perimeter = 3.14159265351459355
sides = 524288, perimeter = 3.14159265357099365
sides = 1048576, perimeter = 3.14159265358509371
sides = 2097152, perimeter = 3.14159265358861872
sides = 4194304, perimeter = 3.14159265358950002
sides = 8388608, perimeter = 3.14159265358972034
sides = 16777216, perimeter = 3.14159265358977546
sides = 33554432, perimeter = 3.14159265358978929
sides = 67108864, perimeter = 3.14159265358979273
sides = 134217728, perimeter = 3.14159265358979356
sides = 268435456, perimeter = 3.14159265358979378
sides = 536870912, perimeter = 3.14159265358979384
sides = 1073741824, perimeter = 3.14159265358979384
sides = 2147483648, perimeter = 3.14159265358979384
| |||||
| 1299.2 | thanks | SMAUG::ABBASI | Sun Sep 23 1990 11:07 | 7 | |
Thanks Jim !
I guess this is a lesson in how to keep the computation stable.
for polygon with 536870912 the perimeter came out to be
3.141592653589793837642929474895936436951160430908203125000000
using printf("%2.50f")
| |||||
| 1299.3 | precision and accuracy | CSSE::NEILSEN | I used to be PULSAR::WALLY | Fri Sep 28 1990 12:57 | 6 |
I think your value for pi becomes inaccurate right here
for polygon with 536870912 the perimeter came out to be
3.141592653589793837642929474895936436951160430908203125000000
^
238646
| |||||
| 1299.4 | precision and nit-picking :-) | ALLVAX::JROTH | It's a bush recording... | Sat Sep 29 1990 01:56 | 11 |
<<< Note 1299.3 by CSSE::NEILSEN "I used to be PULSAR::WALLY" >>>
-< precision and accuracy >-
I think your value for pi becomes inaccurate right here
for polygon with 536870912 the perimeter came out to be
3.141592653589793837642929474895936436951160430908203125000000
^
238646
^
46264338
| |||||