The Archimedes alogorithm for pi
James Van Buskirk
Archimedes algorithm is attractive. Here is an attempt to
get a little better estimate without using any calculus.
We need to prove a few trigonometric inequalities.
For any first-quadrant angle theta,
sec(theta) = sec(theta)**2*cos(theta)
= (1+tan(theta)**2)*(1-2*sin(theta/2)**2)
= (1+4*cos(theta/2)**2*sin(theta/2)**2/cos(theta)**2)* &
(1-2*sin(theta/2)**2)
= 1+2*sin(theta/2)**2+4*sec(theta)*sin(theta/2)**4
>= 1+2*sin(theta/2)**2
sec(theta/2) = (1+tan(theta/2)**2)*(1-2*sin(theta/4)**2)
= (1+sin(theta/2)**2/cos(theta/2)**2)*(1-sin(theta/2)**2/ &
(2*cos(theta/4)**2))
= 1+sin(theta/2)**2/2+(sec(theta/2)**2- &
sec(theta/4)**2/2*sec(theta/4)**2/4- &
sec(theta/2)**2*sec(theta/4)**2/2)*sin(theta/2)**4
= 1+sin(theta/2)**2/2+(4*sec(theta/2)+2)*sin(theta/4)**4
>= 1+sin(theta/2)**2/2
sin(theta/2) = sin(theta)*sec(theta/2)/2
>= sin(theta)/2*(1+sin(theta/2)**2/2)
>= sin(theta)/2*(1+sin(theta)**2/8)
And now we can prove Lemma 1:
For any first-quadrant angle theta and natural number n,
sin(theta/2**n) >= sin(theta)/2**n*(1+(1-1/4**n)/6*sin(theta)**2)
Proof: We have already established this for n = 1, and if true
for any n >= 1, then
sin(theta/2**(n+1)) = sin((theta/2)/2**n)
>= sin(theta/2)/2**n*(1+(1-1/4**n)/6*sin(theta/2)**2)
>= sin(theta)/2*(1+sin(theta)**2/8)/2**n* &
(1+(1-1/4**n)/6*(sin(theta)/2)**2)
= sin(theta)/2**(n+1)*(1+(1-1/4**(n+1))/6*sin(theta)**2+ &
(1-1/4**n)/192*sin(theta)**4)
>= sin(theta)/2**(n+1)*(1+(1-1/4**(n+1))/6*sin(theta)**2)
Now we can refine the Archimedes estimates for pi a little.
In figure 1 of:
http://home.comcast.net/~kmbtib/StupidJava/figsa.html
the first quadrant, NOP, of the unit circle is depicted.
Ray OQ is drawn from the origin making angle theta with
the x-axis OP and line segment QR is draw at right angles
to x-axis OP intersecting the unit circle at Q. Comparing
triangle ORQ and sector POQ, we see that the area of
half-segment RPQ is theta/2-sin(theta)*cos(theta)/2.
In figure 2, chord QP is drawn and bisector of angle
POQ intersects chord QP at T and the unit circle at S.
Considering half segments TPS and TSQ and triangle RPQ,
we find that
theta/2-sin(theta)*cos(theta)/2 = &
2*(theta/4-sin(theta/2)*cos(theta/2)/2)+ &
sin(theta)*(1-cos(theta))/2
Substituting theta/2 for theta and cleaning up, we get
theta-sin(theta) = 2*(theta/2-sin(theta/2))+ &
2*sin(theta/2)*(1-cos(theta/2))
Now,
2*sin(theta/2)*(1-cos(theta/2)) = &
2*sin(theta/2)*2*sin(theta/4)**2
= sin(theta/2)**3*sec(theta/4)**2
= sin(theta/2)**3*(1+sin(theta/4)**2*sec(theta/4)**2)
= sin(theta/2)**3*(1+sin(theta/2)**2*sec(theta/4)**4/4)
= sin(theta/2)**3*(1+sin(theta/2)**2/4* &
(1+sin(theta/2)**2*sec(theta/4)**4/2+tan(theta/4)**4))
>= sin(theta/2)**3*(1+sin(theta/2)**2/4+sin(theta/2)**4/8)
So
theta-sin(theta) >= 2*(theta/2-sin(theta/2))+ &
sin(theta/2)**3*(1+sin(theta/2)**2/4+sin(theta/2)**4/8)
Repeating this process a total of n times we get
theta-sin(theta) >= 2**n*(theta/2**n-sin(theta/2**n))+ &
sum([(2**k*(sin(theta/2**k)**3/2+sin(theta/2**k)**5/8+ &
sin(theta/2**k)**7/16),k=1,n)])
From Lemma 1 we have
sin(theta/2**n)**3 >= sin(theta)**3/8**n*(1+ &
(1-1/4**n)/2*sin(theta)**2)
sin(theta/2**n)**5 >= sin(theta)**5/32**n
sin(theta/2**n)**7 >= sin(theta)**7/128**n
theta-sin(theta) >= sum([((1/4**k)/2*sin(theta)**3+ &
(1/4**k-1/16**k)/4*sin(theta)**5+(1/16**k)/8*sin(theta)**5+ &
(1/64**k)/16*sin(theta)**7,k=1,n)])
= (1/6)*sin(theta)**3+(1/12-1/120)*sin(theta)**5+ &
(1/1008)*sin(theta)**7-(1/4**n)/6*sin(theta)**3- &
((1/4**n)/12-(1/16**n)/120)*sin(theta)**5- &
((1/64**n)/1008)*sin(theta)**7
We are assured by the Archimedean property of the real numbers
that for a given first-quadrant angle theta there exists some
n such that
(1/1008)*sin(theta)**7 >= (1/4**n)/6*sin(theta)**3+ &
((1/4**n)/12-(1/16**n)/120)*sin(theta)**5+ &
((1/64**n)/1008)*sin(theta)**7
So we have estalished
theta-sin(theta) >= (1/6)*sin(theta)**3+(3/40)*sin(theta)**5
for any first-quadrant angle theta.
To get an upper bound for theta, consider figure 3 where again
the first quadrant of the unit circle, NOP, is depecited.
Line PU is drawn tangent to the unit circle at the x-axis at
P and ray OU is drawn from the origing at angle theta to the
x-axis OP, intersecting the unit circle at Q and line PU at U.
Comparing triangle OPU and sector OPQ we see that the area of
the shape (what is the name of this kind of shape?) PUQ is
tan(theta)/2-theta/2.
In figure 4, bisector of angle POQ is drawn intersecting the
unit circle at W and line PU at V. Also line segment QV
joining points Q and V is drawn. Considering shapes PVW
and WVQ and triangle VUQ, we find that:
tan(theta)/2-theta/2 = 2*(tan(theta/2)/2-theta/4)+ &
tan(theta/2)*(sec(theta)-1)/2
Or,
tan(theta)-theta = 2*(tan(theta/2)-theta/2)+ &
sin(theta/2)*sec(theta/2)*2*sin(theta/2)**2*sec(theta)
>= 2*(tan(theta/2)-theta/2)+2*sin(theta/2)**3* &
(1+sin(theta/2)**2/2)*(1+2*sin(theta/2)**2)
= 2*(tan(theta/2)-theta/2)+2*sin(theta/2)**3+ &
5*sin(theta/2)**5+2*sin(theta/2)**7
From the secant inequalities given earlier. Repeating the
process a total of n times we have:
tan(theta)-theta >= 2**n*(tan(theta/2**n)-theta/2**n)+ &
sum([(2**k*(sin(theta/2**k)**3+(5/2)*sin(theta/2**k)**5+ &
sin(theta/2**k)**7),k=1,n)])
From Lemma 1 we have
sin(theta/2**n)**3 >= sin(theta)**3/8**n*(1+ &
(1-1/4**n)/2*sin(theta)**2+(1-2/4**n+1/16**n)/12*sin(theta)**4)
sin(theta/2**n)**5 >= sin(theta)**5/32**n* &
(1+(5-5/4**n)/6*sin(theta)**2)
sin(theta/2**n)**7 >= sin(theta)**7/128**n
tan(theta)-theta >= sum([((1/4**k)*sin(theta)**3+ &
(1/4**k-1/16**k)/2*sin(theta)**5+ &
(1/4**k-2/16**k+1/64**k)/12*sin(theta)**7+ &
(5/16**k)/2*sin(theta)**5+(25/16**k-25/64**k)/12*sin(theta)**7+ &
(1/64**k)*sin(theta)**7),k=1,n)])
= (1/3)*sin(theta)**3+(1/6+2/15)*sin(theta)**5+ &
(1/36+23/180-1/63)*sin(theta)**7-(1/4**n)/3*sin(theta)**3- &
((1/4**n)/6+(2/16**n)/15)*sin(theta)**5- &
((1/4**n)/36+(23/16**n)/180-(1/64**n)/63)*sin(theta)**7
Again, we know that for a given first-quadrant angle theta
there is an n large enough that
37*cos(theta)**2*sin(theta)**7/(2520*cos(theta/2)**4)
>= (1/4**n)/3*sin(theta)**3+((1/4**n)/6+ &
(2/16**n)/15)*sin(theta)**5+ &
((1/4**n)/36+(23/16**n)/180-(1/64**n)/63)*sin(theta)**7
so we have established
tan(theta)-theta >= (1/3)*sin(theta)**3+(3/10)* &
sin(theta)**5+(44/315)*sin(theta)**7- &
37*cos(theta)**2*sin(theta)**7/(2520*cos(theta/2)**4)
To render this inequality in a more palatable form, note that
sin(x)**3*(1/3+3*sin(x)**2/10+44*sin(x)**4/315)
= tan(x)**3*cos(x)*(1/3-sin(x)**2/30+ &
(44*cos(x)**2/315-3/10)*sin(x)**4)
= tan(x)**3*(1/3-sin(x)**2/5+(-105+352*cos(x/2)**4*cos(x)**2- &
756*cos(x/2)**4-554*cos(x)**2*cos(x/2)**2+420*cos(x/2)**2+ &
176*cos(x)**4*cos(x/2)**2)*sin(x)**4/(2520*cos(x/2)**4))
= tan(x)**3*(1/3-tan(x)**2/5+(504*cos(x/2)**4-105*cos(x)**2+ &
352*cos(x/2)**4*cos(x)**4-756*cos(x/2)**4*cos(x)**2- &
554*cos(x)**4*cos(x/2)**2+420*cos(x/2)**2*cos(x)**2+ &
176*cos(x)**6*cos(x/2)**2)*sin(x)**4/(2520*cos(x)**2*cos(x/2)**4))
= tan(x)**3*(1/3-tan(x)**2/5+((126+504*cos(x)+924*cos(x)**2+ &
1176*cos(x)**3+1050*cos(x)**4+823*cos(x)**5)*(1-cos(x))**2+ &
772*cos(x)**6*(1-cos(x))+37*cos(x)**7)*sin(x)**4/ &
(2520*cos(x)**2*cos(x/2)**4))
So we have found the expression we wanted, viz.
tan(theta)**3/3-tan(theta)**5/5
and the quantity we borrowed to dominate the remainder term
for the geometric series and even a little more and so have
proved that for any first quadrant angle theta:
tan(theta)-theta >= tan(theta)**3/3-tan(theta)**5/5
Combining lower and upper bounds:
sin(theta)+sin(theta)**3/6+3*sin(theta)**5/40
<= theta <= tan(theta)-tan(theta)**3/3+tan(theta)**5/5
OK, so it's hardly news that any partial sum of the
Maclaurin series for ASIN(x) is an underestimate nor that
that the sum of the first odd number of terms of the
Maclaurin series for ATAN(x) is an overestimate, but I have
tried to use concepts that were known in the time of
Archimedes, even if the technical execution of the algebra
was more strenuous than seen then.
Going back to the Archimedes method, he knew that
csc(pi/6) = 2 and cot(pi/6) = sqrt(3)
and used the identities
csc(theta) = sqrt(1+cot(theta)**2)
and
cot(theta/2) = cot(theta)+csc(theta)
to get to bounds on cot(pi/96) and csc(pi/96) and so to find
bounds on pi.
In keeping with the spirit of Archimedes, we will do the same,
using only integer arithmetic, but retaining more accuracy
throughout and applying the enhancements noted above:
C:\gfortran\james\archpi>type piday4.f90
module pifuncs
implicit none
integer,parameter :: ik8 = selected_int_kind(18)
integer,parameter :: pref = selected_int_kind(38)
integer,parameter :: ik16 = merge(pref,ik8,pref>=0)
integer,parameter :: MaxPiDen = 96
contains
subroutine surd(x,num,den,minden,hi)
integer(ik8),value :: x
integer(ik8) num,den,minden
logical,value :: hi
integer(ik8) isq,isqsq,t1,t2,delta
integer n
integer(ik8) n0,n1,n2,d0,d1,d2
integer(ik8) a,r,s
if(x <= 0) then
num = x
den = 1
return
end if
n = (bit_size(x)-leadz(x)-1)/2
isq = ishft(1_ik8,n)
isqsq = ishft(1_ik8,2*n)
delta = ishft(isq,-1)
t1 = isqsq
t2 = ishft(t1,-2)
do n = n-1,0,-1
if(isqsq+t1+t2 <= x) then
isqsq = isqsq+t1+t2
isq = isq+delta
t1 = ishft(t1,-1)+t2
else
t1 = ishft(t1,-1)
end if
t2 = ishft(t2,-2)
delta = ishft(delta,-1)
end do
if(isqsq == x) then
num = isq
den = 1
return
end if
n0 = 0
n1 = 1
d0 = 1
d1 = 0
r = 0
s = 1
do
hi = .NOT.hi
a = (isq+r)/s
n2 = n0+a*n1
n0 = n1
n1 = n2
d2 = d0+a*d1
d0 = d1
d1 = d2
if(hi .AND. d1>=minden) exit
r = a*s-r
s = (x-r**2)/s
end do
num = n1
den = d1
end subroutine surd
subroutine rat(x,y,num,den,minden,hi)
integer(ik8),value :: x,y
integer(ik8) num,den,minden
logical,value :: hi
integer(ik8) n0,n1,n2,d0,d1,d2
integer(ik8) a,y1
if(y<=0 .OR. x<=0) then
num = x
den = y
return
end if
n0 = 0
n1 = 1
d0 = 1
d1 = 0
do
hi = .NOT.hi
a = x/y
n2 = n0+a*n1
n0 = n1
n1 = n2
d2 = d0+a*d1
d0 = d1
d1 = d2
if(hi .AND. d1>=minden) exit
y1 = x-a*y
if(y1==0) exit
x = y
y = y1
end do
num = n1
den = d1
end subroutine rat
end module pifuncs
program piday4
use pifuncs
implicit none
integer(ik8) x,y
integer(ik8) CotNumLo0,CotDenLo0,CotNumHi0,CotDenHi0
integer(ik8) CscNumLo,CscDenLo,CscNumHi,CscDenHi
integer(ik8) CotNumLo1,CotDenLo1,CotNumHi1,CotDenHi1
integer(ik8) MinDen,PiDen
integer(ik8) SerNum1,SerDen1,SerNum0,SerDen0
minden = 200000
PiDen = 6
x = 3
call surd(x,CotNumLo0,CotDenLo0,MinDen,hi=.FALSE.)
call surd(x,CotNumHi0,CotDenHi0,MinDen,hi=.TRUE.)
write(*,'(5(i0,a))') CotNumLo0,'/',CotDenLo0,' < cot(pi/', &
PiDen,') < ',CotNumHi0,'/',CotDenHi0
write(*,'(3(a,i0))') 'Check: ',CotNumLo0**2,' < ',3*CotDenLo0**2
write(*,'(3(a,i0))') 'Check: ',3*CotDenHi0**2,' < ',CotNumHi0**2
x = 2
y = 1
call rat(x,y,CscNumLo,CscDenLo,MinDen,hi=.FALSE.)
call rat(x,y,CscNumHi,CscDenHi,MinDen,hi=.TRUE.)
write(*,'(5(i0,a))') CscNumLo,'/',CscDenLo,' <= csc(pi/', &
PiDen,') <= ',CscNumHi,'/',CscDenHi
do
PiDen = 2*PiDen
x = CotNumLo0*CscDenLo+CotDenLo0*CscNumLo
y = CscDenLo*CotDenLo0
call rat(x,y,CotNumLo1,CotDenLo1,MinDen,hi=.FALSE.)
x = CotNumHi0*CscDenHi+CotDenHi0*CscNumHi
y = CscDenHi*CotDenHi0
call rat(x,y,CotNumHi1,CotDenHi1,MinDen,hi=.True.)
write(*,'(7(i0,a)/6(a,i0))') CotNumLo1,'/',CotDenLo1,' <= ', &
CotNumLo0,'/',CotDenLo0,'+',CscNumLo,'/',CscDenLo, &
' < cot(pi/',PiDen,')',' < ',CotNumHi0,'/',CotDenHi0, &
'+',CscNumHi,'/',CscDenHi,' <= ',CotNumHi1,'/',CotDenHi1
write(*,'(3(a,i0))') 'Check: ',1_ik16*CotNumLo1*CotDenLo0*CscDenLo, &
' <= ',1_ik16*CotDenLo1*(CotNumLo0*CscDenLo+CotDenLo0*CscNumLo)
write(*,'(3(a,i0))') 'Check: ',1_ik16*CotDenHi1*(CotNumHi0*CscDenHi+ &
CotDenHi0*CscNumHi),' <= ',1_ik16*CotNumHi1*CotDenHi0*CscDenHi
CotNumLo0 = CotNumLo1
CotDenLo0 = CotDenLo1
CotNumHi0 = CotNumHi1
CotDenHi0 = CotDenHi1
x = CotNumLo0**2+CotDenLo0**2
call surd(x,x,y,MinDen,hi=.FALSE.)
y = y*CotDenLo0
call rat(x,y,CscNumLo,CscDenLo,MinDen,hi=.FALSE.)
x = CotNumHi0**2+CotDenHi0**2
call surd(x,x,y,MinDen,hi=.TRUE.)
y = y*CotDenHi0
call rat(x,y,CscNumHi,CscDenHi,MinDen,hi=.TRUE.)
write(*,'(5(i0,a)/4(a,i0))') CscNumLo,'/',CscDenLo,' <= sqrt((', &
CotNumLo0,'/',CotDenLo0,')**2+1) < csc(pi/',PiDen, &
')',' < sqrt((',CotNumHi0,'/',CotDenHi0,')**2+1) <= ', &
CscNumHi,'/',CscDenHi
write(*,'(3(a,i0))') 'Check: ',1_ik16*CscNumLo**2*CotDenLo0**2, &
' <= ',1_ik16*CscDenLo**2*(CotNumLo0**2+CotDenLo0**2)
write(*,'(3(a,i0))') 'Check: ',1_ik16*CscDenHi**2*(CotNumHi0**2+ &
CotDenHi0**2),' <= ',1_ik16*CscNumHi**2*CotDenHi0**2
if(PiDen >= MaxPiDen) exit
end do
x = 20*CscNumHi**2+9*CscDenHi**2
y = 120*CscNumHi**2
call rat(x,y,SerNum0,SerDen0,MinDen,hi=.FALSE.)
write(*,'(5(i0,a))') SerNum0,'/',SerDen0,' <= 1/6+3*(', &
CscDenHi,'/',CscNumHi,')**2/40 < 1/6+3*sin(pi/', &
PiDen,')**2/40'
write(*,'(3(a,i0))') 'Check: ',SerNum0*120*CscNumHi**2,' <= ', &
SerDen0*(20*CscNumHi**2+9*CscDenHi**2)
call rat(CscDenHi**2,CscNumHi**2,x,y,MinDen,hi=.FALSE.)
x = y*SerDen0+x*SerNum0
y = SerDen0*y
call rat(x,y,SerNum1,SerDen1,MinDen,hi=.FALSE.)
write(*,'(6(i0,a)/3(a,i0))') SerNum1,'/',SerDen1,' < 1+(', &
SerNum0,'/',SerDen0,')*(',CscDenHi,'/',CscNumHi, &
')**2',' < 1+sin(pi/',PiDen,')**2/6+3*sin(pi/',PiDen, &
')**4/40'
write(*,'(3(a,i0))') 'Check: ',1_ik16*SerNum1*SerDen0*CscNumHi**2,' < ',
&
SerDen1*(1_ik16*SerDen0*CscNumHi**2+1_ik16*SerNum0*CscDenHi**2)
x = PiDen*SerNum1*CscDenHi
y = SerDen1*CscNumHi
call rat(x,y,SerNum0,SerDen0,MinDen,hi=.FALSE.)
write(*,'(7(i0,a)/5(a,i0))') SerNum0,'/',SerDen0,' < ',PiDen, &
'*(',SerNum1,'/',SerDen1,')*(',CscDenHi,'/', &
CscNumHi,')',' < ',PiDen,'*(sin(pi/',PiDen,')+sin(pi/', &
PiDen,')**3/6+3*sin(pi/',PiDen,')**5/40) < pi'
write(*,'(3(a,i0))') 'Check: ',1_ik16*SerNum0*SerDen1*CscNumHi, &
' < ',1_ik16*SerDen0*PiDen*SerNum1*CscDenHi
x = 5*CotNumLo0**2-3*CotDenLo0**2
y = 15*CotNumLo0**2
call rat(x,y,SerNum0,SerDen0,MinDen,hi=.FALSE.)
write(*,'(5(a,i0))') '-1/3+tan(pi/',PiDen, &
')**2/5 < -1/3+(',CotDenLo0,'/',CotNumLo0, &
')**2/5 < ',-SerNum0,'/',SerDen0
write(*,'(3(a,i0))') 'Check: ',-1_ik16*SerDen0*(5*CotNumLo0**2- &
3*CotDenLo0**2),' < ',-15*1_ik16*CotNumLo0**2*SerNum0
call rat(CotDenHi0**2,CotNumHi0**2,x,y,MinDen,hi=.FALSE.)
x = SerNum0*x
y = SerDen0*y
call rat(x,y,SerNum1,SerDen1,MinDen,hi=.FALSE.)
write(*,'(6(a,i0),a/3(a,i0))') '-tan(pi/',PiDen,')**2/3+tan(pi/', &
PiDen,')**4/5 < -(',SerNum0, &
'/',SerDen0,')*(',CotDenHi0,'/',CotNumHi0,')**2',' < -', &
SerNum1,'/',SerDen1
write(*,'(3(a,i0))') 'Check: ',-1_ik16*SerDen1*SerNum0*CotDenHi0**2, &
' < ',-1_ik16*SerNum1*SerDen0*CotNumHi0**2
call rat(PiDen*CotDenLo0*(SerDen1-SerNum1),CotNumLo0*SerDen1, &
SerNum0,SerDen0,minden,hi=.TRUE.)
write(*,'(4(a,i0),a/7(a,i0))') 'pi < ',PiDen,'*(tan(pi/',PiDen, &
')-tan(pi/',PiDen,')**3/3+tan(pi/',PiDen, &
')**5/5)',' < ',PiDen,'*(',CotDenLo0,'/',CotNumLo0, &
')*(1-(',SerNum1,'/',SerDen1,') < ',SerNum0,'/',SerDen0
write(*,'(3(a,i0))') 'Check: ',PiDen*SerDen0*(1_ik16*CotDenLo0)* &
(SerDen1-SerNum1),' < ',1_ik16*CotNumLo0*SerDen1*SerNum0
end program piday4
C:\gfortran\james\archpi>gfortran piday4.f90 -opiday4
C:\gfortran\james\archpi>piday4
716035/413403 < cot(pi/6) < 978122/564719
Check: 512706121225 < 512706121227
Check: 956722646883 < 956722646884
2/1 <= csc(pi/6) <= 2/1
1542841/413403 <= 716035/413403+2/1 < cot(pi/12)
< 978122/564719+2/1 <= 2107560/564719
Check: 637815097923 <= 637815097923
Check: 1190179175640 <= 1190179175640
5874077/1520323 <= sqrt((1542841/413403)**2+1) < csc(pi/12)
< sqrt((2107560/564719)**2+1) <= 847399/219323
Check: 5896937408734549301348961 <= 5896937408734829649426010
Check: 229002748087119097400569 <= 229002748088750784806161
2021891/266187 <= 1542841/413403+5874077/1520323 < cot(pi/24)
< 2107560/564719+847399/219323 <= 2484518/327093
Check: 1270770805135998579 <= 1270770805137367038
Check: 307722126586738773 <= 307722126586900766
3015571/393611 <= sqrt((2021891/266187)**2+1) < csc(pi/24)
< sqrt((2484518/327093)**2+1) <= 14498523/1892437
Check: 644336597784810018141729 <= 644336597784899016036850
Check: 22490029431936219891214837 <= 22490029431937479648882321
8740414/572877 <= 2021891/266187+3015571/393611 < cot(pi/48)
< 2484518/327093+14498523/1892437 <= 11724022/768433
Check: 915769283676520398 <= 915769283676644106
Check: 7257203566558184165 <= 7257203566558788102
6019819/393715 <= sqrt((8740414/572877)**2+1) < csc(pi/48)
< sqrt((11724022/768433)**2+1) <= 14611518/955639
Check: 11892951275787962737643169 <= 11892951275793988599443125
Check: 126067368959957200040117333 <= 126067368959968205216354436
8560813/280252 <= 8740414/572877+6019819/393715 < cot(pi/96)
< 11724022/768433+14611518/955639 <= 8288213/271328
Check: 1930893666918728715 <= 1930893666919152796
Check: 6086403993465395456 <= 6086403993465661331
15937580/521463 <= sqrt((8560813/280252)**2+1) < csc(pi/96)
< sqrt((8288213/271328)**2+1) <= 16762542/548455
Check: 19949967692034661274425600 <= 19949967692037454820430537
Check: 20685641095184936883734825 <= 20685641095186254161842176
111801/670483 <= 1/6+3*(548455/16762542)**2/40 < 1/6+3*sin(pi/96)**2/40
Check: 6563363573433506016 <= 6563363591954758251
5989562/5988493 < 1+(111801/670483)*(548455/16762542)**2
< 1+sin(pi/96)**2/6+3*sin(pi/96)**4/40
Check: 1128398743026399603293058744 < 1128398743026590856777176241
729016/232053 < 96*(5989562/5988493)*(548455/16762542)
< 96*(sin(pi/96)+sin(pi/96)**3/6+3*sin(pi/96)**5/40) < pi
Check: 73180350515738041296 < 73180350515878620480
-1/3+tan(pi/96)**2/5 < -1/3+(280252/8560813)**2/5 < -170442/511655
Check: -187369070267287251115 < -187369070265905974470
-tan(pi/96)**2/3+tan(pi/96)**4/5 < -(170442/511655)*(271328/8288213)**2
< -224/627453
Check: -7873123227541116585984 < -7873123209213229115680
pi < 96*(tan(pi/96)-tan(pi/96)**3/3+tan(pi/96)**5/5)
< 96*(280252/8560813)*(1-(224/627453) < 4908616/1562461
Check: 26366669127711685248 < 26366669127714774024
So by these methods we have been able to show that
729016/232053 < pi < 4908616/1562461.
If Archimedes had known that his results would be analyzed for
millenia, he may have come up with a somewhat better estimate!
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
Transfer Principle
> Even though it's the oldest algorithm, the simplicity of the
> Archimedes algorithm is attractive. Here is an attempt to
> get a little better estimate without using any calculus.
> 729016/232053 < pi < 4908616/1562461.
Nice! The decimal equivalents of those are:
3.14159265340... < pi < 3.14159265415...
So that gives us almost nine decimal digits of accuracy.
> If Archimedes had known that his results would be analyzed for
> millenia, he may have come up with a somewhat better estimate!
Happy Pi Day, James!