(* ::Package:: *)

BeginPackage["PSLQ`"];


Print["Basic PSLQ implementation by Armin Straub
	accompanying the paper \"A gentle introduction to PSLQ\"
	-- Tulane University -- Version 1.2 (2010/12/17)"];


PSLQ::usage=
	"PSLQ[x] tries to find integer relations between the real numbers x[1], x[2], ... using the PSLQ algorithm.
	If successful the method returns a (conjectural) integer relation; otherwise, a lower bound is given on the size of any potential relation.
	The following options are available:
		PrecisionGoal - controls the precision for deciding whether a relation holds (defaults to Automatic which currently is 5 decimal digits less than the precision of the input)
		MaxIterations - maximal number of iterations to be performed before giving up (defaults to $IterationLimit)
		StepMonitor - a function that is called after each iteration; see source code (defaults to None)
		Gamma - the value of the parameter \[Gamma] used in PSLQ (defaults to Sqrt[2])
		Verbose - if set to True extra output (number of iterations, relation bound, ...) is printed (defaults to False)
		FailSilently - if set then no warnings about precision will be printed if no relation is found (default: False)
		AlwaysReturnRelation - if set then the best found relation will be returned even though it is definitely incorrect (default: False)
	Usage:
		PSLQ[N[Join[{Pi},Table[Sum[1/(16^k(8k+j)),{k,0,Infinity}],{j,1,7}]],20]];
		PSLQ[N[EulerGamma^Range[0,10],1000], Verbose->True]";


Options[PSLQ]={
	PrecisionGoal->Automatic,
	MaxIterations:>$IterationLimit,
	StepMonitor->None,
	Gamma->Sqrt[2],
	Verbose->False,
	FailSilently->False,
	AlwaysReturnRelation->False};


Begin["Private`"];


PSLQ::lowprec="Precision too low to continue (`1` iterations performed).";
PSLQ::norel="No integer relation was found. The norm of any true integer relation is at least `1`.";
PSLQ::norelr="No integer relation was found (the returned one is valid only to `1` digits). The norm of any true integer relation is at least `2`.";
PSLQ[x_,OptionsPattern[]]:=Module[{prec,precg,monitor,\[Gamma]p,n=Length[x],i,j,t,it,A,B,s,y,H,r,a,b,d,rel={},relfound=False},
	(* Set up precision parameters. *)
	prec=Precision[x];
	precg=OptionValue[PrecisionGoal];
	If[precg===Automatic,precg=prec-5];
	(* \[Gamma] > 2/Sqrt[3], \[Tau] = 1/Sqrt[1/4+1/\[Gamma]^2] *)
	\[Gamma]p=N[Table[OptionValue[Gamma]^j,{j,1,n-1}]];
	monitor=OptionValue[StepMonitor];
	(* Initialize *)
	A=B=IdentityMatrix[n];
	s=Sqrt[Reverse[Accumulate[Reverse[x^2]]]];
	y=x/s[[1]];s=s/s[[1]];
	H=Table[If[i<j,0,If[i==j,s[[i+1]]/s[[i]],-y[[i]]y[[j]]/(s[[j]]s[[j+1]])]],{i,1,n},{j,1,n-1}];
	For[i=2,i<=n,i++,For[j=i-1,j>=1,j--,
		t=Round[H[[i,j]]/H[[j,j]]];
		y[[j]]=y[[j]]+t y[[i]];
		H[[i,1;;j]]=H[[i,1;;j]]-t H[[j,1;;j]];
		A[[i,1;;n]]=A[[i,1;;n]]-t A[[j,1;;n]];
		B[[1;;n,j]]=B[[1;;n,j]]+t B[[1;;n,i]];
	]];
	(* Iteration *)
	For[it=0,it<OptionValue[MaxIterations],it++,
		(* Monitor *)
		If[monitor=!=None,monitor[y,A,H]];
		(* Step 4: Termination *)
		If[Min[Abs[y]]<=10^-precg,
			(* Read off relation candidate from column of B. *)
			rel=Transpose[B][[First[Ordering[Abs[y],1]]]];
			relfound=Abs[Total[rel x]]<=10^-precg;
			If[relfound,Break[]];
		];
		If[Log[10,Max[Abs[A]]]>prec-5||Precision[Diagonal[H]]<=5,
			If[!OptionValue[FailSilently],Message[PSLQ::lowprec,it]];
			Break[]
		];
		(* Step 1: Exchange *)
		r=First[Ordering[\[Gamma]p Abs[Diagonal[H]],-1]];
		{y[[r]],y[[r+1]]}={y[[r+1]],y[[r]]};
		{A[[r]],A[[r+1]]}={A[[r+1]],A[[r]]};
		{H[[r]],H[[r+1]]}={H[[r+1]],H[[r]]};
		{B[[All,r]],B[[All,r+1]]}={B[[All,r+1]],B[[All,r]]};
		(* Step 2: Corner *)
		If[r<n-1,
			a=H[[r,r]];b=H[[r,r+1]];d=Sqrt[a^2+b^2];
			For[i=r,i<=n,i++,
				{H[[i,r]],H[[i,r+1]]}={a/d H[[i,r]]+b/d H[[i,r+1]],-b/d H[[i,r]]+a/d H[[i,r+1]]}];
		];
		(* Step 3: Reduction *)
		For[i=r+1,i<=n,i++,For[j=Min[i-1,r+1],j>=1,j--,
			t=Round[H[[i,j]]/H[[j,j]]];
			y[[j]]=y[[j]]+t y[[i]];
			H[[i,1;;j]]=H[[i,1;;j]]-t H[[j,1;;j]];
			A[[i,1;;n]]=A[[i,1;;n]]-t A[[j,1;;n]];
			B[[1;;n,j]]=B[[1;;n,j]]+t B[[1;;n,i]];
		]];
	];
	(* Read off best relation so far for diagnostic purposes. *)
	If[!relfound,
		rel=Transpose[B][[First[Ordering[Abs[y],1]]]]];
	(* Warn if integer relations is not valid. *)
	If[!relfound&&!OptionValue[FailSilently],
		If[OptionValue[AlwaysReturnRelation],
			Message[PSLQ::norelr,Floor[-Log[10,Abs[Total[rel x]]]],N[1/Max[Abs[Diagonal[H]]]]],
			Message[PSLQ::norel,N[1/Max[Abs[Diagonal[H]]]]]]];
	(* Print some statistics if requested. *)
	If[OptionValue[Verbose],
		Print[StringForm["`` iterations performed during which precision dropped to ``.",it,Round[Precision[Diagonal[H]]]]];
		Print["Relation norm bound: ",N[1/Max[Abs[Diagonal[H]]]]];
		Print["Relation residue: ",N[Total[rel x]]]];
	If[relfound||OptionValue[AlwaysReturnRelation],
		Return[rel],Return[{}]];
];


End[];


EndPackage[];
