The Theory of Rational Approximations to Irrational Numbers

(*    Any irrational can be approximated by fractions. The function, Phoenix, will deliver fractional approxomations with virtually unlimited accuracy for any real number.  E.g., it zeros in on two fractions (one on its left and one to its right) that are accurate for π to over 11830 decimal places in 8 minutes by using the command "Phoenix[π,8000,8000]," defined within
Mathematica7, below. The earth moving function plows through 8000 Partial Sequences, cultivating 78819 fractions to uncover the last two π impersonators. See Below. 
   *)

(*   We will define the term Partial Sequence to refer to a finite sequence of 2-D vectors, where the first element corresponds to the numerator of a given fraction and the 2'nd element corresponds to its denominator. When each 2-D vector of a Partial sequence is converted to a fraction, the sequence of fractions approaches a target number, usually irrational. To qualify as an entry to the P.S.,each vector (as a fraction) must be closer to the irrational target than the previous one. Also, all first elements of the 2-D vectors (the yin sequence) are in ascending order likewise for all second elements (the yang sequence). All three sequences are limited for a rational target and unlimited for an irrational target.
 
The first vector of each new Partial Sequence arises from operations on the last vector in each of the two previous Partial Sequences. In the function below Phoenix these two particular vectors are called CrimsonAsh and GoldenAsh, exploiting the saga of the God-like Bird, Phoenix, from Greek mythology, who sprang up again and again, from the ashes of previous incarnations.
 Once the first vector of a new P.S. is found, CrimsonAsh is recursively added until the 
newborn vector (converted to a fraction) is further from the target than the previous converted 
vector. At this point, the Partial Sequence dies. The lengths of Partial Sequences for π vary widely, 
from 1 (in many P.S.'s) to 10,389 (in P.S. 292).

In our use of 2-D vectors, we mean that they have the same operations of addition and scalar multiplication. 
We say that the vectors encountered here correspond to fractions. Thus, for all integers a, b, c, d and k, noting that for our purposes,  k = -1  or  k > 1, we have:
																		   Phoenix - Mathematica7
				Vector Addition:  [a, b] + [c, d] = [a + c, b + d]        GoldenAsh + CrimsonAsh
 
				Scalar Multiplication:    k[a, b] = [ka, kb].             n*GoldenAsh, sgn*CrimsonAsh
-------------------------------------------------------------------------------------------------------
					   U/D fcns.
						   Conversion:  [a, b]' = a/b                         (VP[v])

		  Comparison with target, t:  [a, b]'(t) = a/b - t                   (VecPrm[v,t])   

		   If negative, fraction is to left of target, otherwise to right.
         
		    Distance to the target, t:   [a, b]"(t) = |a/b - t|.               (VecPrm2[v,t])
	--------------------------------------------

					  The Fundamental Identity in Rational Approximation Theory:\\n\\n											   [ka, kb]' = [a, b]'				
	
	The FIRAT above involves conversion with scalsr multiplication of any vector.

    Most of the auxillary functions below are very simple ones. For example, the function VP[v] merely returns the fraction v[[1]]/v[[2]]. Each Partial Sequence, when their elements are converted to fractions represent rational numbers that get closer and closer to the target of the function Phoenix. VecPrm computes the difference of v with the target, while VecPrm2 forms the absolute value of this difference for comparison purposes. FirstVec returns the first vector of a P.S. LastVec returns the last one and Victor returns the n'th vector 
     of the k'th P.S. TestPolRev tests (T/F) to see if there has been a polarity reversal at the end of a given P.S. IndexSeqWithLength searches for the first P.S. that has the length specified. Phoenix[π, k, dp] returns the k'th sequence that converges to π, while the GoldenFireBird  only flies after the PhoenixMatx function has been used. This function keeps everything, including the matrix of two dimensional vectors. If the target is rational,
    the number of Partial Sequences is finite. The function stores the vectors into a (2 x n) matrix, where n is the total number of fractions found by the function, approaching the target. It's good  practice to execute the Exit[] command when finishing a session using PhoenixMatx in order to clear the memory.
                                                                                                                    *)

(*  Living in India and migrating to Egypt every 500 years, the Crimson and Golden Bird flies from Sunrise to Sunset.
  His name is legend:
                                       Φοινιξ                            *)

Phoenix[RealNo_,k_,dp_]:=  (*  Key to discover the (k^th) P.S. for any real number, with dp decimal places.   *)(Module[{CrimsonAsh,GoldenAsh,Joaquin,NewLeaf,JoaquinPhoenix,LSK,NumSeq,target,best,test,n,sgn,LenSeq},

CrimsonAsh={1,0};(*  Let the Force (1/0) be with you:  Floor[RealNo] < RealNo < ∞.  *)    GoldenAsh={Floor[RealNo],1};NumSeq=1; If[k==1,Return[{GoldenAsh,1}]];target=N[RealNo,Max[dp,4*k]];best=VecPrm2[GoldenAsh,target];LenSeq={1};

Label[Ra];LSK=1;n=1;sgn=1;If[CrimsonAsh≠{1,0},sgn=-Sign[VecPrm[CrimsonAsh,target]*VecPrm[GoldenAsh,target]]];

Label[Heliopolis];n=n+1;Joaquin=n*GoldenAsh+sgn*CrimsonAsh;test=VecPrm2[Joaquin,target];(*  Discovery Loop  *)   If[test>best,(*  Sojourn back to Sun City. In order to move forward, Joaquin must be closer to his target than he was before. *)  Goto[Heliopolis]];
(*  Otherwise, Joaquin turns over a New Leaf.  *)  NewLeaf=Joaquin;CrimsonAsh=GoldenAsh;GoldenAsh=NewLeaf;best=test;  Label[FrankincenseAndMyrrh];Joaquin=GoldenAsh+CrimsonAsh;test=VecPrm2[Joaquin,target];(*  Recursive Loop  *)If[test<best,best=test;LSK=LSK+1;GoldenAsh=Joaquin;Goto[FrankincenseAndMyrrh]];

(* The Spices have been prepared for the Nest of the Golden Fire Bird. *)
LenSeq=Append[LenSeq,LSK];NumSeq=NumSeq+1;If[NumSeq<k,(*  Return to Egyptian Sun God  *) Goto[Ra]];

(*  Eureka, the Great Bird flies again as the (k^th) P.S.:  (Subscript[S, k]) = {(Subscript[v, 2]) + (N - 1)(Subscript[v, 1])}  *)  JoaquinPhoenix={NewLeaf};
Do[NewLeaf=NewLeaf+CrimsonAsh;JoaquinPhoenix=Append[JoaquinPhoenix,NewLeaf],{var,LSK-1}];Print[" Φοινιξ - The Great Crimson and Golden Bird Soars Above"];Return[Append[JoaquinPhoenix,LenSeq]]]);

VP[v_]:=(v[[1]]/v[[2]]);                                                     (*     [a, b]' = a/b           *)
VecPrm[v_,target_]:=(N[VP[v]-target]);                  (*  [a, b'](t) = [a, b]' - t   *)
VecPrm2[v_,target_]:=(Abs[VecPrm[v,target]])      (*  [a, b]"(t) = |[a, b]'(t)|  *)

PhoenixMatx[RealNumber_,LimNoSeq_,dp_]:=  (*  Provide the target, no. of P.S.'s & decimal places desired.  *)(st=AbsoluteTime[];target=N[RealNumber,Max[dp,4*LimNoSeq]];RevPol=0;LenSeq={1};LSK=1;NumSeq=1;RPSeqs={};it=Floor[target];V1={it,1};V2={it+1,1};bestL=N[target-it,100];bestR=N[it+1-target,100];best=Min[bestL,bestR];If[bestL<bestR,Matx={{V1[[1]]},{V1[[2]]}},Matx={{V2[[1]]},{V2[[2]]}}];Label[initialize];V3=V1+V2;test=Abs[V3[[1]]/V3[[2]]-target];If[test<best,  (*  then  *)  If[bestL<bestR,Matx={{Matx[[1,1]],V3[[1]]},{Matx[[2,1]],V3[[2]]}},Matx={{V3[[1]],Matx[[1,1]]},{V3[[2]],Matx[[2,1]]}}],If[bestL<bestR,V2=V3,(* else *) V1=V3];Goto[initialize]];If[bestL<bestR,V2=V3,(* else *) V1=V3];numFrac=2;
Label[continue];LSK=1;n=1;sgn=1;Label[RecursiveLoop];V3=V1+V2;test=VecPrm2[V3,target];If[test<best,  (*  then  *)  best=test;LSK=LSK+1;Matx=Join[{Join[Matx[[1]],{V3[[1]]}]},{Join[Matx[[2]],{V3[[2]]}]}];numFrac=numFrac+1;V2=V3;Goto[RecursiveLoop]];(*  Otherwise, the partial sequence ends.  *)  LenSeq=Append[LenSeq,LSK];NumSeq=NumSeq+1;If[best<2*(10^(-dp)),Goto[fin]];If[NumSeq==LimNoSeq,Goto[fin]];If[VecPrm[V1,target]*VecPrm[V2,target]>0,sgn=-1;RevPol=RevPol+1];Label[DiscoveryLoop];n=n+1;V3=n*V2+sgn*V1;test=VecPrm2[V3,target];If[test>best,  (*  then go back until true  *)  Goto[DiscoveryLoop],  
(*  else the first element, V3, that is, of the next partial sequence has been found  *)  best=test;numFrac=numFrac+1;Matx=Join[{Join[Matx[[1]],{V3[[1]]}]},{Join[Matx[[2]],{V3[[2]]}]}];V1=V2;V2=V3;Goto[continue]];Label[fin];tt=AbsoluteTime[]-st;Print[" Φοινιξ - The Great Crimson and Golden Bird Soars Above"];If[NumSeq<LimNoSeq,Print["No. of P.S.'s = ",NumSeq]]);

   (*  This function delivers the k'th Partial Sequence from any Matx found in the PhoenixMatx function.  *)
GoldenFireBird[Matx_,k_,LenSeq_]:=   (*  Arguments are the matrix, the index of the P.S. and the LenSeq array.  *)(Module[{Vec1,Vec2,Vec3,LSK,GFB,NewBird},If[k==1,Return[{Matx[[1,1]],Matx[[2,1]]}]];Vec1=LastVec[Matx,k-1,LenSeq];Vec2=FirstVec[Matx,k,LenSeq];Vec3=LastVec[Matx,k,LenSeq];LSK=(Vec3[[2]]-Vec2[[2]])/Vec1[[2]]+1;GFB={Vec2};NewBird=Vec2;Print[" Φοινιξ - The Great Crimson and Golden Bird Soars Above"];
Do[NewBird=NewBird+Vec1;GFB=Append[GFB,NewBird],{n,LenSeq[[k]]-1}];Return[GFB]]);

GFBirdShort[Matx_,k_,LenSeq_]:=   (*  The output shortened form of the GoldenFireBird  *)(Module[{Vec1,Vec2,Vec3,LSK},If[k==1,Return[{Matx[[1,1]],Matx[[2,1]]}]];Vec1=LastVec[Matx,k-1,LenSeq];Vec2=FirstVec[Matx,k,LenSeq];Vec3=LastVec[Matx,k,LenSeq];LSK=(Vec3[[2]]-Vec2[[2]])/Vec1[[2]]+1;Print[" Φοινιξ - The Great Crimson and Golden Bird Soars Above"];If[LSK==1,Print[" Single Vector = ",Vec2];Return[]];
Print[" First Vector = ",Vec2];Print[" Incrementing by: ",Vec1];Print[" Last Vector =  ",Vec3];Print[" Number of Vectors in Partial Sequence = ",LSK];Return[]]);

FirstVec[Matx_,k_,LenSeq_]:=  (*  Returns the first vector in Sequence No. k.  *)(Module[{Col},Col=Sum[LenSeq[[i]],{i,k-1}];Return[{Matx[[1,Col+1]],Matx[[2,Col+1]]}]]);
LastVec[Matx_,k_,LenSeq_]:=  (*  Returns the last vector in Sequence No. k.  *)(Module[{Col},Col=Sum[LenSeq[[i]],{i,k}];Return[{Matx[[1,Col]],Matx[[2,Col]]}]]);
Victor[Matx_,n_,k_,LenSeqMtx_]:=  (*  What's your vector, Victor? In Matx, the n-th vector in P.S. k.  *) (If[n>LenSeqMtx[[k]],Return["Undefined"]];Module[{V1,V2},V1=LastVec[Matx,k-1,LenSeqMtx];V2=FirstVec[Matx,k,LenSeqMtx];Return[V2+(n-1)*V1]]);
IndexSeqWithLength[LenSeq_,Len_]:=(Module[{k},Do[If[LenSeq[[i]]==Len,k=i;Goto[fin]],{i,Length[LenSeq]}];Return["Not Found"];Label[fin];Return[k]]);

TestPolRev[Matx_,IdxLS_,target_]:=   (*  Tests for Polarity Reversal. IdxLS is the index of the last sequence,
thus, IdxLS - 1 is the index of the sequence before the last one.  *)(Module[{Colm1,Colm2,test},Colm2=Sum[LenSeq[[i]],{i,IdxLS}];Colm1=Sum[LenSeq[[i]],{i,IdxLS-1}];test=Sign[(Matx[[1,Colm2]]/Matx[[2,Colm2]]-target)*(Matx[[1,Colm1]]/Matx[[2,Colm1]]-target)];If[test<0,Return[False],Return[True]]]);

RatMat[Matx_,col_,dp_]:=(N[Abs[Matx[[1,col]]/Matx[[2,col]]-target],dp]);
RatMat2[Matx_,col_,dp_]:=(*  This returns the ratio of the accuracy of Matx in column, col, to its previous column, col - 1. The accuracy is computed to dp decimal places.  *)
(N[Abs[(Matx[[1,col-1]]/Matx[[2,col-1]]-target)/(Matx[[1,col]]/Matx[[2,col]]-target)],dp]);

Test[u_,v_,target_]:=(Module[{w},w=u+v;VecPrm2[w,target]]);
PhoenixNoBird[RealNo_,k_,dp_]:=  (*  Key to discover the (k^th) P.S. for any real number, with dp decimal places.   *)(Module[{CrimsonAsh,GoldenAsh,Joaquin,NewLeaf,JoaquinPhoenix,LSK,NumSeq,target,best,test,n,sgn,LenSeq},

CrimsonAsh={1,0};(*  Let the Force (1/0) be with you:  Floor[RealNo] < RealNo < ∞.  *)    GoldenAsh={Floor[RealNo],1};NumSeq=1; If[k==1,Return[{GoldenAsh,1}]];target=N[RealNo,Max[dp,4*k]];best=VecPrm2[GoldenAsh,target];LenSeq={1};

Label[Ra];LSK=1;n=1;sgn=1;If[CrimsonAsh≠{1,0},sgn=-Sign[VecPrm[CrimsonAsh,target]*VecPrm[GoldenAsh,target]]];

Label[Heliopolis];n=n+1;Joaquin=n*GoldenAsh+sgn*CrimsonAsh;test=VecPrm2[Joaquin,target];(*  Discovery Loop  *)   If[test>best,(*  Sojourn back to Sun City. In order to move forward, Joaquin must be closer to his target than he was before. *)  Goto[Heliopolis]];
(*  Otherwise, Joaquin turns over a New Leaf.  *)  NewLeaf=Joaquin;CrimsonAsh=GoldenAsh;GoldenAsh=NewLeaf;best=test;  Label[FrankincenseAndMyrrh];Joaquin=GoldenAsh+CrimsonAsh;test=VecPrm2[Joaquin,target];(*  Recursive Loop  *)If[test<best,best=test;LSK=LSK+1;GoldenAsh=Joaquin;Goto[FrankincenseAndMyrrh]];

(* The Spices have been prepared for the Nest of the Golden Fire Bird. *)
LenSeq=Append[LenSeq,LSK];NumSeq=NumSeq+1;If[NumSeq<k,(*  Return to Egyptian Sun God  *) Goto[Ra]];

(*  Eureka, the Great Bird flies again as the (k^th) P.S.:  (Subscript[S, k]) = {(Subscript[v, 2]) + (N - 1)(Subscript[v, 1])}  *)  JoaquinPhoenix={NewLeaf};
Do[NewLeaf=NewLeaf+CrimsonAsh;JoaquinPhoenix=Append[JoaquinPhoenix,NewLeaf],{var,LSK-1}];Return[Append[JoaquinPhoenix,LenSeq]]]);

VishnusVecs[Seq1_,Indx1_,Seq2_,Indx2_,target_]:=(Module[{u,v,m,z,LenPS1,LenPS2,PartSeq1,PartSeq2,OL},PartSeq1=PhoenixNoBird[target,Seq1,1000];LenPS1=Length[PartSeq1]-1;OL="Angler, you'll be fined for exceeding your fishing limit.";If[Indx1>LenPS1,Print["First selected Partial Sequence # ",Seq1,", only has ",LenPS1," vectors."];Return[OL]];PartSeq2=PhoenixNoBird[target,Seq2,1000];LenPS2=Length[PartSeq2]-1;If[Indx2>LenPS2,Print["Second selected Partial Sequence # ",Seq2,", only has ",LenPS2," vectors."];Return[OL]];u=PartSeq1[[Indx1]];v=PartSeq2[[Indx2]];m=Transpose[{u,v}];

uu={{u[[1]]},{u[[2]]}};vv={{v[[1]]},{v[[2]]}};Print["We have the following information about the selected vectors"];

Print["Vishnu's Vectors are:  ",MatrixForm[uu]," and  ",MatrixForm[vv]];Print["The matrix formed by this pair of vectors is ",MatrixForm[m]];Print["The determinant of the matrix above has the value:  ",Det[m]];z=Eigenvalues[m];
Print["The product of the eigenvalues of m is:  ",N[z[[1]]*z[[2]],55]];Print["The approximate eigenvalues of the matrix above are:"];
Print[N[z,33]];
Print["The sum of the eigenvalues of m is:   ",N[z[[1]]+z[[2]],33]];Print["The sum of its main diagonal pair is:  ",m[[1,1]]+m[[2,2]]]]);

VishnusDet[Seq1_,Indx1_,Seq2_,Indx2_,target_]:=(Module[{u,v,m},u=PhoenixNoBird[target,Seq1,1000][[Indx1]];v=PhoenixNoBird[target,Seq2,1000][[Indx2]];m=Transpose[{u,v}];Return[Det[m]]]);

(*  Two irrational constants are set, first Euler' little gamma  *)    γ=EulerGamma;
(*  "Riemann's Zeta Function," Edwards, p.67:   The infinite sum, where ρ represents an imaginary zero
 of the Zeta Fcn,  Σ 1/ρ =  *)    RiemannsGarden=1+ (1/2)γ + Log[ (1/2)(1/(Gamma[(1/2)])) ];    (*  The dueces are wild in this constant. RH is half full of 'em, by Gamma.     *)

st=AbsoluteTime[];v=Phoenix[π,8000,8000];AbsoluteTime[]-st

%/60

Length[v]

VecPrm[v[[1]],N[Pi,15000]]     (* fraction to the right of π, accurate to 11830 decimal places.  *)

StringLength[ToString[VP[v[[1]]]]]    (*  This is a boat-load of digits, top and bottom.  *)

VecPrm[v[[2]],N[Pi,15000]]      (* fraction to the left of π, accurate to 11831 decimal places.  *)

N[{Mean[v[[3]]],StandardDeviation[v[[3]]],Max[v[[3]]],Median[v[[3]]]}]

v=Phoenix[π,292,200];

v[[10390]]      (*  The last vector contains the lengths of the first 292 P.S.'s  *)

v2=Phoenix[Pi,2,20]

v3=Phoenix[Pi,3,20]

v4=Phoenix[Pi,4,20]

v5=Phoenix[Pi,5,20]

v6=Phoenix[Pi,6,20]

v7=Phoenix[Pi,7,20]       (*  Calculator accuracy  *)

v8=Phoenix[Pi,8,20]

v9=Phoenix[Pi,9,20]

v10=Phoenix[Pi,10,20]

VecPrm2[v10[[2]],N[Pi,25]]     (*  Accurate to 16 decimal places  *)

v20=Phoenix[Pi,20,20]

VecPrm2[v20[[2]],N[Pi,35]]     (*  Accurate to 29 decimal places  *)

v50=Phoenix[Pi,50,100]

VecPrm2[v50[[2]],N[Pi,99]]     (*  Accurate to 70 decimal places  *)

v150=Phoenix[Pi,150,200]

VecPrm2[v150[[1]],N[Pi,300]]     (*  Accurate to 225 decimal places  *)

v292=Phoenix[π,292,200];

v292[[10390]]      (*  The last vector contains the lengths of the first 292 P.S.'s  *)

VecPrm2[v292[[10389]],N[Pi,500]]     (*  Accurate to 440 decimal places  *)

st=AbsoluteTime[];v=Phoenix[π,8000,8000];AbsoluteTime[]-st

%/60       (*  8 minutes  *)

Length[v]

VecPrm[v[[1]],N[Pi,15000]]     (* fraction to the right of π, accurate to 11830 decimal places.  *)

StringLength[ToString[VP[v[[1]]]]]    (*  This is a boat-load of digits, top and bottom.  *)

VecPrm[v[[2]],N[Pi,15000]]      (* fraction to the left of π, accurate to 11831 decimal places.  *)

N[{Mean[v[[3]]],StandardDeviation[v[[3]]],Max[v[[3]]],Median[v[[3]]]}]

vE50=Phoenix[E,50,200]     (*  Unlike other trancendentals LenSeq for E exhibits a pattern  *)

(*  Place any real number in the first argument of Phoenix, the index of the Partial Sequence desired
into the second argument. If you get funny results, like all 1's after a certain point in LenSeq, just
enlarge the decimal place accuracy with the third argument of Phoenix.  *)