SC10 Supplemental Material

4039 days ago by SC2010

Computational and Experimental Mathematics

Further Experimental Opportunities

Neil Calkin and Dan Warner

November 16, 2009

 

A Novel Approach to Computing $\pi$

In 1997 David Bailey, Peter Borwein, and Simon Plouffe used experimental mathematics, including the PSLQ algorithm, to determine a new series for computing $\pi$.  This series has a unique property.  It can be used to effectively compute the digits of $\pi$ starting at any specified digit without computing all the previous digits.  As we will see there is an initialization step that can be done very quickly - at least relative to the cost of computing all the initial digits, and it has extremely low memory requirements.

The series is given by

$\pi = \sum\limits_{k=0}^\infty \frac{1}{16^k}\left( \frac{4}{8*k+1} - \frac{2}{8*k+4} - \frac{1}{8*k+5} - \frac{1}{8*k+6}\right)$  

$\pi = \left( {4} - \frac{2}{4} - \frac{1}{5} - \frac{1}{6}\right) + \frac{1}{16}\left( \frac{4}{9} - \frac{2}{12} - \frac{1}{13} - \frac{1}{14}\right) + \frac{1}{16^2}\left( \frac{4}{17} - \frac{2}{20} - \frac{1}{21} - \frac{1}{22}\right) + \cdots$


# Some numerical values for the first few terms of the BBP Algorithm s = RDF(0) for k in range(5): t = (4/(8*k+1) - 2/(8*k+4) - 1/(8*k+5) - 1/(8*k+6))/(16^k) s += t print t print s 
       
# How about 10 terms? s = RDF(0) for k in range(11): t = (4/(8*k+1) - 2/(8*k+4) - 1/(8*k+5) - 1/(8*k+6))/(16^k) s += t print 'BBP(10) =', s.n(digits=15) print 'Pi = ', pi.n(digits=15) 
       
# How about 100 terms? Clearly we will need more precision. R350 = RealField(350) s = R350(0) for k in range(101): t = (4/(8*k+1) - 2/(8*k+4) - 1/(8*k+5) - 1/(8*k+6))/(16^k) s += t print 'BBP(100) =', s print 'Pi = ', pi.n(digits=105) 
       
def convert_hex(x,n): h = [floor(x)] x = x - floor(x) for k in range(n): y = 16*x z = floor(y) x = y - z h.append(z) return h u = 12.0/16^5 print u h = convert_hex(u,15) s = 0.0 for k in range(len(h)): s += h[k]/16^k print h print s u = RDF(pi) print u h = convert_hex(u,15) s = 0.0 for k in range(len(h)): s += h[k]/16^k print h print s print 'pi =',u.n(digits=15) 
       
R200 = RealField(200) p200 = R200(pi) print p200 h = convert_hex(p200,50) s = R200(0) for k in range(len(h)): s += h[k]/16^k print h print s h[47] 
       

Find the 47-th HEX digit without computing the earlier digits

In particular do the calculation in 64 bit arithmetic instead of 200 bit arithmetic

$16^{46} \pi = \sum\limits_{k=0}^{46} 16^{46-k}\left( \frac{4}{8*k+1} - \frac{2}{8*k+4} - \frac{1}{8*k+5} - \frac{1}{8*k+6}\right) + \sum\limits_{k=47}^\infty 16^{46-k}\left( \frac{4}{8*k+1} - \frac{2}{8*k+4} - \frac{1}{8*k+5} - \frac{1}{8*k+6}\right)$

Let

$S_j = \sum\limits_{k=0}^{\infty}  \frac{1}{8*k+j} $

Then

$\pi = 4S_1 - 2 S_4 - S_5 - S_6$

$\left\{ 16^d \pi \right\} = \left\{4\left\{ 16^d S_1\right\} - 2\left\{ 16^d S_4\right\}- \left\{16^d S_5\right\} -\left\{ 16^d S_6 \right\}\right\}$

def frac(x): return (x - floor(x)) def BBP0(n): nd = n-1 s = RDF(0) s1 = s4 = s5 = s6 = s for k in range(25): s1 += (1.0/(8*k+1))*16^(nd-k) s4 += (1.0/(8*k+4))*16^(nd-k) s5 += (1.0/(8*k+5))*16^(nd-k) s6 += (1.0/(8*k+6))*16^(nd-k) # print s1,s4,s5,s6 # print frac(s1),frac(s4),frac(s5),frac(s6) s = 4*frac(s1) - 2*frac(s4) - frac(s5) - frac(s6) s = frac(s) return s 
       
u = RDF(pi) print u hpi = convert_hex(u,15) print hpi[5:] mypi = BBP0(5) print mypi hd = convert_hex(mypi,10) print hd 
       
R200 = RealField(200) p200 = R200(pi) print p200 hpi = convert_hex(p200,50) print hpi[47:] mypi = BBP0(47) print mypi hd = convert_hex(mypi,10) print hd 
       

However, Bailey, Borwein, and Plouffe had figured out a neat trick around the problem of keeping track of all those leading digits - use modular arithmetic

$\left\{ 16^d \pi \right\} = \left\{4\left\{ 16^d S_1\right\} - 2\left\{ 16^d S_4\right\}- \left\{16^d S_5\right\} -\left\{ 16^d S_6 \right\}\right\}$

$\left\{ 16^d S_j \right\} = \left\{ \left\{ \sum\limits_{k=0}^{d} \frac{16^{d-k}}{8*k+j}\right\}  + \sum\limits_{k=d+1}^\infty  \frac{16^{d-k}}{8*k+j} \right\}$

$\left\{ 16^d S_j \right\} = \left\{ \left\{ \sum\limits_{k=0}^{d} \frac{16^{d-k}\mathrm{mod}(8k+j)}{8*k+j}\right\}  + \sum\limits_{k=d+1}^\infty  \frac{16^{d-k}}{8*k+j} \right\}$

# The Binary Exponentiation Algorithm # Consider the problem of calculating x^n # For example n = 45 # In binary 45 = 101101 = 1*32 + 0*16 + 1*8 + 1*4 + 0*2 + 1*1 # This means that # # x^45 = x^32 * x^8 * x^4 * x # # which leads to this simple algorithm s = 1 n = 45 xp = x while 0<n: if mod(n,2)==1: s = s*xp xp = xp*xp n = n//2 print s 
       
def int_mod(a,b): # mod(a,b) q = a // b r = a - q*b return r def power_mod(n,k): b = 16 r = 1 if n==0: return r else: t = 1 while t <= n: t = 2*t t = t/2 while 1 <= t: # print r,t,n if t <= n: r = int_mod(b*r,k) n = n-t t = t/2 if 1 <= t: r = int_mod(r^2,k) return r power_mod(4,3);power_mod(5,4);power_mod(4,7);power_mod(25,9);power_mod(10,8*15+9) 
       
mod(16^4,3);mod(16^5,4);mod(16^4,7);mod(16^25,9);mod(16^10,8*15+9) 
       
def BBP(n): nd = n-1 s0 = RDF(0) sv = [s0,s0,s0,s0] jv = [1,4,5,6] for j in range(4): sj = sv[j] for k in range(nd): k8j = k*8+jv[j] if 1 < k8j: sj += power_mod(nd-k,k8j)/k8j sj = frac(sj) for k in range(nd,nd+10): k8j = k*8+jv[j] sj += 1/(16^(k-nd)*k8j) sv[j] = frac(sj) s0 = frac(4*sv[0] - 2*sv[1] - sv[2] - sv[3]) return s0 
       
R200 = RealField(200) p200 = R200(pi) print p200 hpi = convert_hex(p200,50) print hpi[47:] nd = 47 s = BBP(nd) hd = convert_hex(s,15) print hd 
       
Re5 = RealField(100000) pe5 = Re5(pi) hpi = convert_hex(pe5,25000) print hpi[24990:] nd = 24990 s = BBP(nd) hd = convert_hex(s,15) print hd 
       

The PSLQ Algorithm

The PSLQ algorithm was developed in 1993 by Helaman Ferguson, who is also well known for his mathematical scultures.

 
       
 
       
r""" An implementation of the PSLQ algorithm developed by H. R. P. Ferguson, D. H. Bailey and S. Arno. AUTHORS: -- Agnes Jany (2009-05-10) EXAMPLES: sage: t=vector([11,27,31]) sage: p=pslq(t) Relation found: (-8, 9, -5) sage: p=pslq(t,iteration=true) Iteration no. 0 Iteration no. 1 Iteration no. 2 Iteration no. 3 Relation found: (-8, 9, -5) sage: p=pslq(t,digits=1000,iteration=true) Iteration no. 0 Iteration no. 1 Iteration no. 2 Iteration no. 3 Iteration no. 4 Relation found: (-6, -1, 3) sage: p=pslq(t,tau=1.01,iteration=true) Iteration no. 0 Iteration no. 1 Iteration no. 2 Iteration no. 3 Iteration no. 4 Iteration no. 5 Relation found: (-8, 9, -5) sage: p=pslq(t,tau=1.99,iteration=true) Iteration no. 0 Iteration no. 1 Relation found: (0, -31, 27) sage: p=pslq(t,tau=1,iteration=true) Choose tau such that 1<tau<2 ! sage: p=pslq(t,digits=-5,iteration=true) Set digits to a positive integer! sage: t=vector(RDF,[11,27,31]) sage: p=pslq(t) Relation found: (-6, -1, 3) sage: t=vector(RealField(213),[11,27,31]) sage: p=pslq(t) Relation found: (-1, 5, -4) sage: alp=vector(RealField(212),17) sage: for i in range(0,17,1): sage: alp[i]=((3^(1/4)-2^(1/4))^i) sage: p=pslq(alp,digits=62,tau=1.000006145) Relation found: (1, 0, 0, 0, -3860, 0, 0, 0, -666, 0, 0, 0, -20, 0, 0, 0, 1) sage:p=pslq(alp,tau=1.000006145) Relation found: (192, -1536, 162, 85, 287, 274, 814, -503, 1148, 786, -206, -464, -279, 83, 144, -316, -725) sage: p=pslq(alp,digits=100,tau=1.000006145) /usr/local/sage/local/bin/sage-sage: line 347: 3995 Killed python "$@" sage: t=vector([1,(pi),(sqrt(2)),(sqrt(3)),(ln(1+sqrt(2))),(ln(2+sqrt(3))), sage: 0.926390055174046729218163586547779014444960190107335046732521921271418504594036683829313473075349968212]) sage: p=pslq(t) Relation found: (4, -7, 17, -6, 21, 42, -75) sage: t=vector([gamma(1/2),gamma(7/2),gamma(13/2)]) sage: p=pslq(t) Relation found: (15, -8, 0) sage: t=vector([gamma(1/3),1,sqrt(5)]) sage: p=pslq(t,30) Relation found: (727989880380837, -562243113417127, -620731142335100) sage: p=pslq(t,100) /usr/local/sage/local/bin/sage-sage: line 347: 4092 Segmentation fault python "$@" """ #***************************************************************************** # Copyright (C) 2009 Agnes Jany <agnes.jany@yahoo.de> # # Distributed under the terms of the GNU General Public License (GPL) # http://www.gnu.org/licenses/ #***************************************************************************** class pslq(object): def __init__(self, vec ,digits=50,tau=(4/3)^0.5,iteration=false,printing=true): r""" Initialization step of the PSLQ algorithm. INPUT: vec -- input vector digits -- integer tau -- real number iteration -- boolean value NOTES: Define vectors x and s. Define matrices H,A and B. Perform first reduction. See [FBA] p.8. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ if tau<=1 or tau>=2: print "Choose tau such that 1<tau<2 !" else: if digits <=0: print "Set digits to a positive integer!" else: self.printing=printing sys.setrecursionlimit(100000) self.it=iteration self.d=digits self.gamma=sqrt((4*tau^2)/(4-tau^2)) self.count=-1 n=len(vec) RF=RealField(212) x=vector( RF,(1/norm(vec))*(vec)) s=self._set_s_(x,n) H=self._define_H_(s,x,n) A1=MatrixSpace(ZZ,n) A=A1.identity_matrix() B1=MatrixSpace(ZZ,n) B=B1.identity_matrix() H,A,B,x=self._Reduce_(H,A,B,x,n,false) self._iterate_(H,A,B,x,n) def _iterate_(self,H,A,B,x,n): r""" This function executes the iteration steps of PSLQ until it is halted by _terminate_(). INPUT: H -- matrix A,B -- matrices with integer values and A=B^-1 x -- vector n -- length of the input vector NOTES: This function uses the algorithm of [FBA] to find integer relations. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ if self._terminate_(H,A,B,x,n)==true: if self.printing: print self.reason[0] print self.reason[1] return self.reason else: self.count=self.count+1 if self.it: print "Iteration no.", self.count r=self._compare_gamma_(self.gamma,H,n) if r<n-2: self._predef_(H,n,r) self.R=self._define_R_(r,n) x=x*self.R H=self.R*H A=self.R*A B=B*self.R if r<n-2: Q=self._define_Q_(H,n,r) H=H*Q H,A,B,x = self._Reduce_(H,A,B,x,n,r) self._iterate_(H,A,B,x,n) def _set_s_(self,x,n): r""" Define a vector of partial sums. INPUT: x -- vector n -- length of the input vector OUTPUT: s -- vector of partial sums. NOTES: See [FBA] p.4. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ RF=RealField(212) j=0 s=vector(RF,copy(x)) while j<n: k=j sum=0 while k<n: sum=sum+x[k]^2 k=k+1 s[j]=sqrt(sum) j=j+1 return s def _nint_(self,x): r""" Nearest integer function. INPUT: x -- real number. OUTPUT: xnew -- nearest integer to x. For exact half-integer values: integer with greater absolute value. AUTHORS: - Agnes Jany (2009-05-10) """ if x/0.5 not in ZZ or x in ZZ: xnew=round(x) if x/0.5 in ZZ: if x>0: xnew=ceil(x) else: xnew=floor(x) return xnew def _define_H_(self,s,x,n): r""" Define a $n\times n-1$ lower trapezodial matrix H. INPUT: s -- vector of partial sums x -- vector n -- length of the input vector OUTPUT: H -- matrix NOTES: See [FBA] p.5. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ RF=RealField(212) H=matrix(RF,n,n-1) for i in range(0,n,1): for j in range(0,n-1,1): if i<j: H[i,j]=0 if i==j: H[i,j]=s[i+1]/s[i] if j<i: H[i,j]=(-x[i]*x[j])/(s[j]*s[j+1]) return H def _compare_gamma_(self,gamma,H,n): r""" Choose an integer j such that $\gamma^j |h_{jj}|\geq \gamma^i |h_{ii}|$ for all $1\leq i \leq n-1$. INPUT: gamma -- constant real number H -- matrix n -- length of the input vector OUTPUT: j -- integer NOTES: See [FBA] p.9. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ j=0 while j<=n-2: i=0 while i<=n-2: g1=(gamma^j)*abs(H[j,j]) g2=(gamma^i)*abs(H[i,i]) if g1<g2: break else: i=i+1 if i==n-1: return j break else: j=j+1 def _define_R_(self,j,n): r""" Define a permutation matrix R. INPUT: j -- integer n -- length of the input vector OUTPUT: R -- matrix NOTES: See [FBA] p.9. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ R1=MatrixSpace(ZZ,n) R=R1.identity_matrix() R[j,j]=0 R[j,j+1]=1 R[j+1,j+1]=0 R[j+1,j]=1 return R def _predef_(self,H,n,r): r""" Remember certain entries of H before defining Q. INPUT: H -- matrix r -- integer n -- length of the input vector NOTES: See [FBA] p.8-9. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ self.Qa=H[r,r] self.Qb=H[r+1,r] self.Qc=H[r+1,r+1] self.Qd=sqrt(self.Qc^2+self.Qb^2) def _define_Q_(self,H,n,r): r""" Define a $n-1\times n-1$ matrix Q. INPUT: H -- matrix r -- integer n -- length of the input vector OUTPUT: Q -- matrix NOTES: See [FBA] p.9. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ RF=RealField(500) b=self.Qb c=self.Qc d=self.Qd Q1=MatrixSpace(RF,n-1) Q=Q1.identity_matrix() if not (b/d).is_infinity(): Q[r,r]=b/d if not (-c/d).is_infinity(): Q[r,r+1]=-c/d if not (c/d).is_infinity(): Q[r+1,r]=c/d if not (b/d).is_infinity(): Q[r+1,r+1]=b/d return Q def _Reduce_(self,H,A,B,x,n,r): r""" Perform Hermite reduction. INPUT: H,A,B -- matrices x -- vector r -- integer n -- length of the input vector OUTPUT: H,A,B -- matrices x -- vector NOTES: See [FBA] p.8-9 and p.17-18. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ if self.count==-1: for i in range(1,n,1): for j in range(i-1,-1,-1): t=Integer(self._nint_((H[i,j]/H[j,j]))) x[j]=(copy(x[j])+t*x[i]) for k in range(0,j+1,1): H[i,k]=(copy(H[i,k])-t*H[j,k]) for k in range(0,n,1): A[i,k]=(copy(A[i,k])-t*A[j,k]) B[k,j]=(copy(B[k,j])+t*B[k,i]) return H,A,B,x else: for i in range(r+1,n,1): z=min(i-1,r+1) for j in range(z,-1,-1): t=Integer(self._nint_((H[i,j]/H[j,j]))) x[j]=(copy(x[j])+t*x[i]) for k in range(0,j+1,1): H[i,k]=(copy(H[i,k])-t*H[j,k]) for k in range(0,n,1): A[i,k]=(copy(A[i,k])-t*A[j,k]) B[k,j]=(copy(B[k,j])+t*B[k,i]) return H,A,B,x def _zero_H_(self,H,n): r""" Tests for zeros in the diagonal elements of H. INPUT: H -- matrices n -- length of the input vector OUTPUT: boolean value NOTES: See [FBA] p.9. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ i=0 while i<=n-2: if (abs(H[i,i]))<10^(-1000): break else: i=i+1 if i<n-1: return true else: return false def _ratio_x_(self,x,n): r""" Check if the ratio of the smallest and largest x vector entries is very small. INPUT: x -- vector n -- length of the input vector OUTPUT: boolean value NOTES: See [B] p.3. REFERENCES: [B] D. H. Bailey: "Finding New Mathematical Identities via Numerical Computations" ACM SIGNUM, vol. 33, no. 1 (Jan. 1998), pp. 17-22. AUTHORS: - Agnes Jany (2009-05-10) """ xmax=max(abs(x[i]) for i in range(0,n,1)) xmin=min((abs(x[i]) for i in range(0,n,1))) xratio=xmin/xmax if xratio<10^(-20): return true else: return false def _zero_x_(self,x,n,d): r""" Test for zero entries in x. INPUT: x -- vector n -- length of the input vector d -- integer OUTPUT: integer NOTES: Set precision to $10^{-d}$. If x has a zero-entry return its index else return -1. See [FBA] p.9. REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ i=0 while i<=n-2: if abs(x[i])<10^(-d): break else: i=i+1 if i==n-1: return -1 else: return i def _check_A_(self,A): r""" Test for too large entries in A. INPUT: A -- matrix OUTPUT: boolean value NOTES: See [FBA] p.18 . REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. AUTHORS: - Agnes Jany (2009-05-10) """ n=A.nrows() m=A.ncols() Amax=0 for i in range (0,n,1): z=max(abs(A[i,j]) for j in range(0,m,1)) if z>Amax: Amax=z if Amax>10^(10000): return true else: return false def _terminate_(self,H,A,B,x,n): r""" Check termination conditions. INPUT: H,A,B -- matrices x -- vector n -- length of the input vector OUTPUT: boolean value NOTES: See [FBA] p.9 and [B] p.3 . REFERENCES: [FBA] H. R. P. Ferguson, D. H. Bailey, S. Arno : "Analysis of PSLQ, An Integer Relation Finding Algorithm" Math. Comput. 68, 351-369, 1999. [B] D. H. Bailey: "Finding New Mathematical Identities via Numerical Computations" ACM SIGNUM, vol. 33, no. 1 (Jan. 1998), pp. 17-22. AUTHORS: - Agnes Jany (2009-05-10) """ if self._zero_H_(H,n)==true: #wieso zuerst das und dann erst zero x.. self.tex="\text{No relation found:}\quad H_{ii}=0." self.reason="No relation found:","H[i,i]=0" return true else: zerox=self._zero_x_(x,n,self.d) if zerox!=-1: self.tex=latex(B.column(zerox)) self.reason="Relation found:",B.column(zerox) return true else: if self._ratio_x_(x,n): zerox2=self._zero_x_(x,n,floor(self.d/4)) if zerox2!=-1: self.tex=latex(B.column(zerox2)) self.reason="Relation found:",B.column(zerox2) return true else: if self._check_A_(A): self.tex="\text{No relation found:}\quad A_{ij}\quad \text{too large.}" self.reason="No relation found:","A[i,j] too large" return true else: return false def _latex_(self): r""" Return \Latex representation of pslq output. EXAMPLES: sage: t=vector([11,27,31]) sage: a=pslq(t,digits=20) Relation found: (-1, 5, -4) sage: latex(a) \left(-1,5,-4\right) """ return self.tex 
       
t=vector([11,27,31]) p=pslq(t) 
       
 
       
R = RealField(350) s = R(5) phi = (1+ sqrt(s))/2 v = vector([R(1),phi,phi^2]) p = pslq(v) 
       
# In other words # # 1 + x - x^2 = 0 or 1 = x(x-1) or 1/(x-1) = x/1 # sol = solve(1+x-x^2 == 0,x) show(sol) 
       

The $\alpha(N)$ Function

$\alpha(N) = \frac{E^{\prime}(k_N)}{K(k_N)} -\frac{\pi}{4K^2(k_N)}

s = R(3) alpha3 = (sqrt(s)-1)/2 v = vector([1, alpha3, alpha3^2]) p = pslq(v) 
       
sol = solve(1 - 2*x - 2*x^2 == 0, x) show(sol) 
       
s = R(5) alpha5 = (sqrt(s)-sqrt(2*sqrt(s)-2))/2 v = vector([1, alpha5, alpha5^2]) p = pslq(v) 
       
v = vector([1, alpha5, alpha5^2, alpha5^3, alpha5^4]) p = pslq(v) 
       
sol = solve(29 - 80*x - 24*x^2 + 16*x^4 == 0, x) show(sol) 
       
mypi = R(pi) t = R(2) v = vector([arctan(t), arctan(1/t), mypi]) p = pslq(v) 
       
(arctan(t) + arctan(1/t)) 
       
mypi/2 
       
 
       

How about the inverse calculator?

http://glooscap.cs.dal.ca:8087/

print 'phi' print phi 
       
print 'alpha3' print alpha3 
       
print 'alpha5' print alpha5 
       
print 'arctan(2) + arctan(1/2)' print (arctan(t) + arctan(1/t)) 
       
mye = R(e) print '7*pi/11 + e/3' print 7*mypi/11 + mye/3 
       
mye = R(e) print '7*pi/11 + exp(1/3)' print 7*mypi/11 + mye^(1/3) 
       
 
       
# Some scraps for the curious 
       
# The Rabinowitz and Wagon spigot algorithm for pi nd = 20 nt = floor(10*nd/3) mypi = [] num = [k for k in range(nt)] num[0] = 1 den = [2*k+1 for k in range(nt)] dig = [2 for k in range(nt)] indx = range(nt) indx.reverse() for n in range(nd): carry = 0 for k in indx: dk = 10*dig[k] + carry q = dk // den[k] r = dk - q*den[k] dig[k] = r carry = q*num[k] q = carry // 10 dig[0] = carry - q*10 #if 10 <= q: # c = q // 10 # q = q - c*10 # mypi[-1] += c mypi.append(q) print mypi print pi.n(digits=(nd+1)) 
       
# Experimental project: how often do these glitches happen? # Is there a heuristic reason why it should be about 2% of the time? 
       
 
       
# Lambert's Continued Fraction p1 = 4/(1+(1^2/3)) print p1, p1.n(digits=10) p2 = 4/(1+(1^2/(3+2^2/5))) print p2, p2.n(digits=10) p3 = 4/(1+(1^2/(3+2^2/(5+3^2/7)))) print p3, p3.n(digits=10) p4 = 4/(1+(1^2/(3+2^2/(5+3^2/(7+4^2/9))))) print p4, p4.n(digits=10) 
       
z = var('z') L1 = Matrix(2,2,[0, 4, 1, 1]) L2 = Matrix(2,2,[0, 1^2, 1, 3]) print L1 print L2 zz = vector([z,1]) LL = L1*L2 p,q = LL*zz pq = (p/q)(z=0) print p/q, pq, pq.n(digits=10) L3 = Matrix(2,2,[0, 2^2, 1, 5]) LL = LL*L3 p,q = LL*zz pq = (p/q)(z=0) print p/q, pq, pq.n(digits=10) L4 = Matrix(2,2,[0, 3^2, 1, 7]) LL = LL*L4 p,q = LL*zz pq = (p/q)(z=0) print p/q, pq, pq.n(digits=10) L5 = Matrix(2,2,[0, 4^2, 1, 9]) LL = LL*L5 p,q = LL*zz pq = (p/q)(z=0) print p/q, pq, pq.n(digits=10) 
       
L3 s = 0.0 for k in range(1,10): s += 1/(k*2^k) print s