PSLQ Algorithm

PSLQ.java

package org.cab1729.algorithms;

import org.cab1729.matrixmp.*;
import jgmp.MPFloat;

import numbercruncher.matrix.MatrixException;

/**
 * @author jmenes
 * The PSLQ Algorithm
 * ref: 
 * http://www.cecm.sfu.ca/organics/papers/bailey/paper/html/node3.html#SECTION00030000000000000000
 * Borwein P., Computational Excursions in Analysis and Number Theory, Appendix B
 * 
 * numbercruncher package - 
 * http://www.amazon.com/Java-Number-Cruncher-Programmers-Numerical/dp/0130460419
 * 
 * fix 11/12/2011 - changed all computations to double precision
 * (package matrixdp is numbercruncher.matrix.* code with float types changed to double)
 * 
 * new version 12/14/11 - changed all computations to arbitrary multiple precision 
 * -- see MPFloat class
 * (package matrixmp is numbercruncher.matrix.* code with float types changed to MPFloat)
 *
 */
public class PSLQ {
	
	private SquareMatrix A, B = null;
	private Matrix H = null;
	private RowVector S, Y = null;
	
	private int n = 0;	// length of x (input vector)
	private static int e = -1; // index adjust factor
	// changed from static to instance since they need to be initialized in init()
	private MPFloat gamma;	// initialization moved to init() for precision
	private MPFloat mpfOne;	//
	
	private MPFloat prec; 	// detection threshold
	private int digits;     // precision
	
	protected int m = 0;
	
	protected Object result;	// can be MPFloat or MPFloat[] -- see termination test in iterate()

	public PSLQ () {
		
	}
	
	/**
	 * @param X		input vector
	 * @param T		digits of precision
	 */
	public void init (MPFloat[] X, int T) {
		
		n = X.length;
		digits = T;
		MPFloat.setDefaultPrec(digits);
		// initialization of static values moved here for precision
		// - changed from static to instance - see above
		gamma = new MPFloat(2/StrictMath.sqrt(3));
		mpfOne = new MPFloat((float)1);
		
		this.prec = new MPFloat(StrictMath.pow(10, (1-T)+5));
		
		// 1)
		A = new SquareMatrix(n);
		B = new SquareMatrix(n);
		IdentityMatrix.convert(A);
		IdentityMatrix.convert(B);
		
		S = new RowVector(n);
		Y = new RowVector(n);
		
		// 2)
		for (int k = 1+e; k <= n+e; k++) {
			MPFloat sum = new MPFloat(0.0);
			for (int j = k; j <= n+e; j++) {
				sum = sum.add(X[j].mult(X[j]));		// X[j]^2
			}
			S.set(k, sum.sqrt());
		}
		MPFloat t = S.at(1+e);
		
		for (int k = 1+e; k <= n+e; k++) {
			Y.set(k, X[k].div(t));
			S.set(k, S.at(k).div(t));
		}
		
		// 3)
		// H - lower trapezoidal matrix
		H = new Matrix(n, n - 1);
		for (int i = 1+e; i <= n+e; i++) {
			// Matrix class initializes every element to 0 value
//			for (int j = i+1; j <= (n+e)-1; j++) {
//				try {
//					H.set(i, j, new MPFloat());		// set to 0 by default
//				} catch (MatrixException e) {
//					// TODO Auto-generated catch block
//					e.printStackTrace();
//				}
//			}
			
			if (i <= (n+e)-1) {
				try {
					H.set(i, i, S.at(i+1).div(S.at(i)));
				} catch (MatrixException e) {
					// TODO Auto-generated catch block
					e.printStackTrace();
				}
			}
			for (int j = 1+e; j <= i - 1; j++) {
				try {
					H.set(i, j, ((Y.at(i).neg().mult(Y.at(j))).div((S.at(j).mult(S.at(j+1))))));
				} catch (MatrixException e) {
					// TODO Auto-generated catch block
					e.printStackTrace();
				}
			}
		}
		
		// 4)
		// full reduction
		//hermiteReduce(2);
		for (int i=(2+e)+1; i <= n+e; i++) {
			for (int j = i-1; j >= 1+e; j--) {
				try {
					t = (H.at(i, j).div(H.at(j, j))).trunc();
				} catch (MatrixException e) {
					// TODO Auto-generated catch block
					e.printStackTrace();
				}
				Y.set(j, Y.at(j).add(t.mult(Y.at(i))));
				for (int k = 1+e; k <= j; k++) {
					try {
						H.set(i, k, H.at(i, k).sub(t.mult(H.at(j, k))));
					} catch (MatrixException e) {
						// TODO Auto-generated catch block
						e.printStackTrace();
					}
				}
				for (int k = 1+e; k <= n+e; k++) {
					try {
						A.set(i, k, A.at(i, k).sub(t.mult(A.at(j, k))));
						B.set(k, j, B.at(k, j).add(t.mult(B.at(k, i))));
					} catch (MatrixException e) {
						// TODO Auto-generated catch block
						e.printStackTrace();
					}
				}
			}
		}
	}
	
	public void iterate() throws MatrixException {
		// Main Iteration
		
		MPFloat M = new MPFloat();
		
		// 1)
		m = 0;
		int hsize = H.columnCount();
		
		MPFloat g2i = new MPFloat();		// set to 0 by default
		MPFloat max[][] = new MPFloat[H.rowCount()][H.columnCount()];
		
		for (int i = 0; i < hsize; i++) {
			g2i = gamma.pow(i);
			try {
				max[i][i] = g2i.mult((H.at(i, i).abs()));
			} catch (MatrixException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}
		}
		
		MPFloat maxvalue = max[0][0];
		m = 0;
		for (int i = 0; i < hsize; i++) {
			if (max[i][i].compareTo(maxvalue) == 1) 	
			{
				// max[i][i] > maxvalue
				maxvalue = max[i][i];
				m = i;
			}
		}

		// TODO iteration starts here
		do {
		
			// 2)
			MPFloat valm = Y.at(m);
			MPFloat valmp1 = Y.at(m + 1);
			Y.set(m, valmp1);
			Y.set(m + 1, valm);

			RowVector rm = null;
			RowVector rmp1 = null;
			ColumnVector cm = null;
			ColumnVector cmp1 = null;

			try {
				rm = H.getRow(m);
				rmp1 = H.getRow(m + 1);
				H.setRow(rm, m + 1);
				H.setRow(rmp1, m);

				rm = A.getRow(m);
				rmp1 = A.getRow(m + 1);
				A.setRow(rm, m + 1);
				A.setRow(rmp1, m);

				cm = B.getColumn(m);
				cmp1 = B.getColumn(m + 1);
				B.setColumn(cm, m + 1);
				B.setColumn(cmp1, m);

			} catch (MatrixException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}

			// 3)
			MPFloat t0, t1, t2, t3, t4;  // set to 0 by default
			if (m <= (n+e) - 2) {
				try {
					t0 = (H.at(m, m).pow(2)).add((H.at(m, m+1).pow(2))).sqrt();
					t1 = H.at(m, m).div(t0);
					t2 = H.at(m, m+1).div(t0);

					for (int i = m; i < n; i++) {
						t3 = H.at(i, m);
						t4 = H.at(i, m+1);
						H.set(i, m, (t1.mult(t3)).add(t2.mult(t4)));
						H.set(i, m+1, (t2.neg().mult(t3).add(t1.mult(t4))));
					}

				} catch (MatrixException e) {
					// TODO Auto-generated catch block
					e.printStackTrace();
				}
			}	
			// 4)
			MPFloat t = new MPFloat();
			for (int i = m+1; i <= n+e; i++) {
				for (int j = StrictMath.min(i-1, m+1); j > 1; j--) {
					try {
						t = H.at(i, j).div(H.at(j, j)).trunc(); // round
						Y.set(j, Y.at(j).add(t.mult(Y.at(i))));
						for (int k = 0; k < j; k++) {
							H.set(i, k, H.at(i, k).sub(t.mult(H.at(j, k))));
						}
						for (int k = 0; k < n; k++) {
							A.set(i, k, A.at(i, k).sub(t.mult(A.at(j, k))));
							B.set(k, j, B.at(k, j).add(t.mult(B.at(k, i))));
						}
					} catch (MatrixException e) {
						// TODO Auto-generated catch block
						e.printStackTrace();
					}

				}
			}

			// 5)

			MPFloat maxnorm = new MPFloat((float)1);
			int hrows = H.rowCount();
			for (int j = 0; j < hrows; j++) {
				try {
					RowVector rj = H.getRow(j);
					MPFloat enorm = rj.norm();	// euclidian norm
					if (enorm.compareTo(maxnorm) == 1) {
						maxnorm = enorm;
					}

				} catch (MatrixException e) {
					// TODO Auto-generated catch block
					e.printStackTrace();
				}
			}
			M = mpfOne.div(maxnorm);
	
			// block reduction
			hermiteReduce(m+1);
			// get next m
			getMax();	
		
		//  iteration test
			if (maxmv(A).getPrecision() > digits || (minrv(Y).compareTo(prec) <1)) {
				break;
			}
		} while (true) ;
		// 5 - termination test
		// If the largest entry of A exceeds the level of numeric precision used, 
		// then precision is exhausted. 
		// If the smallest entry of the y vector is less than 
		// the detection threshold, a relation has been detected and is given in 
		// the corresponding column of B.

		if (maxmv(A).getPrecision() > digits) {
			result = M;
		} else {
			// TODO may not be getting the correct column index here
			result = B.getColumn(minrc(Y)).copyValues1D();
			// test only - print matrix values here to verify correct column
			System.out.println("A:");
			A.print(A.columnCount());
			System.out.println("B:");
			B.print(B.columnCount());
			System.out.println("H:");
			H.print(H.columnCount());
			System.out.println("S:");
			S.print(S.columnCount());
			System.out.println("Y:");
			Y.print(Y.columnCount());
		} 		
	}
	
	/**
	 * @param val determines if full or block reduction is performed
	 * -- update: this function is now called from iterate() only
	 */
	protected void hermiteReduce(int val) {
		
		MPFloat t = new MPFloat();
		
		for (int i=(val+e)+1; i <= n+e; i++) {
			for (int j = StrictMath.min(i-1, val+1); j >= 1+e; j--) {
				try {
					t = (H.at(i, j).div(H.at(j, j))).trunc();
				} catch (MatrixException e) {
					// TODO Auto-generated catch block
					e.printStackTrace();
				}
				Y.set(j, Y.at(j).add(t.mult(Y.at(i))));
				for (int k = 1+e; k <= j; k++) {
					try {
						H.set(i, k, H.at(i, k).sub(t.mult(H.at(j, k))));
					} catch (MatrixException e) {
						// TODO Auto-generated catch block
						e.printStackTrace();
					}
				}
				for (int k = 1+e; k <= n+e; k++) {
					try {
						A.set(i, k, A.at(i, k).sub(t.mult(A.at(j, k))));
						B.set(k, j, B.at(k, j).add(t.mult(B.at(k, i))));
					} catch (MatrixException e) {
						// TODO Auto-generated catch block
						e.printStackTrace();
					}
				}
			}
		}
	}
	
	protected void getMax() {
		m = 0;
		int hsize = H.columnCount();
		
		MPFloat g2i = new MPFloat();
		MPFloat max[][] = new MPFloat[H.rowCount()] [H.columnCount()];
		
		for (int i = 0; i < hsize; i++) {
			g2i = gamma.pow(i);
			try {
				max[i][i] = g2i.mult(H.at(i, i).floor());
			} catch (MatrixException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}
		}
		
		MPFloat maxvalue = max[0][0];
		m = 0;
		for (int i = 0; i < hsize; i++) {
			if (max[i][i].compareTo(maxvalue) == 1) {
				maxvalue = max[i][i];
				m = i + 1;
			}
		}

	}

	/**
	 * @param m
	 * @return diagonal of a matrix as RowVector
	 * @throws MatrixException
	 */
	protected RowVector diag(SquareMatrix m)
		throws MatrixException
	{
		
		int n = m.rowCount();
		RowVector rv = new RowVector(n);
		
		for (int i = 0; i<n; i++) {
			rv.set(i, m.at(i, i));
		}
		
		return rv;
	}
	
	/**
	 * @param m
	 * @return max value of a given matrix
	 */
	protected MPFloat maxmv(Matrix m) {
		
		MPFloat maxval = new MPFloat();
		
		int mrows = m.rowCount();
		int mcols = m.columnCount();
		MPFloat [][] mvalues = m.values();
		for (int r = 0; r < mrows; r++) {
			for (int c = 0; c < mcols; c++) {
				MPFloat mval = mvalues[r][c];
				if (!Double.isInfinite(mval.doubleValue()) && !Double.isNaN(mval.doubleValue())) {
					if (mval.compareTo(maxval) == 1) {
						maxval = mval;
					}
				}
			}
		}
		
		return maxval;
	}
	
	/**
	 * @param rv
	 * @return minimum value of a RowVector
	 */
	protected MPFloat minrv(RowVector rv) {
		
		MPFloat minval = new MPFloat();
		int rvcol = rv.columnCount();
		for (int c = 0; c < rvcol; c++) {
			MPFloat rvval = rv.at(c);
			if (!Double.isInfinite(rvval.doubleValue()) && !Double.isNaN(rvval.doubleValue())) {
				if (rvval.compareTo(minval) == -1) {
					minval = rvval;
				}
			}
		}
				
		return minval;
	}
	
	/**
	 * @param rv
	 * @return column index of minimum RowVector value
	 */
	protected int minrc(RowVector rv) {
		
		MPFloat minval = new MPFloat();
		int mincol = 0;
		int rvcol = rv.columnCount();
		for (int c = 0; c < rvcol; c++) {
			MPFloat rvval = rv.at(c);
			if (!Double.isInfinite(rvval.doubleValue()) && !Double.isNaN(rvval.doubleValue())) {
				if (rvval.compareTo(minval) == -1) {
					minval = rvval;
					mincol = c;
				}
			}
		}
		
		return mincol;
	}
	
	
	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		
		// test case
		// Experiment #1 from the paper Recognizing Numerical Constants 
		// by David H. Bailey and Simon Plouffe:
		// revised - see: 
		// Experimental Evaluation of Euler Sums
		// by David H. Bailey, Jonathan M. Borwein and Roland Gingersohn
		// RNR Technical report RNR-93-015 June 6, 1994
		
		int T = Integer.parseInt(args[0]);
		System.out.println("\nprecision selected: " + T);
		
		// values computed with GP/PARI
		MPFloat.setDefaultPrec(T);
		System.out.println("default prec set at: " + MPFloat.getDefaultPrec() + "\n");
		// TODO stub - need 135 prec value for Sa(2,3) 
		MPFloat C =
			new MPFloat("0.15616693338117691588103590968798819368577670984000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e0", 10, T);
		MPFloat zeta2 = 
			new MPFloat("1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890061975870530400431896233719067962872468700501e0", 10, T);
		MPFloat zeta3 = 
			new MPFloat("1.20205690315959428539973816151144999076498629234049888179227155534183820578631309018645587360933525814619915779526071941849199599867328e0", 10, T);
		MPFloat zeta4 =
			new MPFloat("1.08232323371113819151600369654116790277475095191872690768297621544412061618696884655690963594169991723299081390804274241458407157457005e0", 10, T);
		MPFloat zeta5 = 
			new MPFloat("1.03692775514336992633136548645703416805708091950191281197419267790380358978628148456004310655713333637962034146655660904280096177915597e0", 10, T);
		MPFloat log2 =
			new MPFloat("0.693147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996418687542001481020570685733685520235758131e0", 10, T);
		MPFloat plog4p5 = 
			new MPFloat("0.517479061673899386330758161898862945622377475141379258244319347977008281581864976936485777826560900647722864999993500396459287987951163e0", 10, T);
		MPFloat plog5p5 =
			new MPFloat("0.508400579242268707459108849258589941319541125664821648724497796352625394228780242619384210049344955062253148566177885373776251290109127e0", 10, T);
				
		MPFloat[] x = 
			new MPFloat[] {
				C, plog5p5, plog4p5.mult(log2), log2.pow(5), zeta5, zeta4.mult(log2), 
				zeta3.mult(log2.pow(2)), zeta2.mult(log2.pow(3)), zeta2.mult(zeta3)
				};
		
		System.out.println("Input vector:");
		for (MPFloat ve:x) {
			System.out.println("\nval: " + ve.toString() + "\ndval: " + ve.doubleValue() + "\nprec: " + ve.getPrecision());
		}
				
		PSLQ pslq = new PSLQ();
		pslq.init(x, T);
		try {
			pslq.iterate();
		} catch (MatrixException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
		// print result
		if (null != pslq.result) {
			if (MPFloat.class.isInstance(pslq.result))
			{
				System.out.println("\nprecision exausted");
				System.out.println("result: " + ((MPFloat)pslq.result).toString());
			}
			else
				if (pslq.result.getClass().isArray())
				{
					System.out.println("\nrelation found");
					for (MPFloat n:(MPFloat[])pslq.result)
						System.out.println("result: " + n.toString());
				}
		} else {
			System.out.println("error, no result available");
		}
		
	}

}

Java implementation of the PSLQ integer relation algorithm. This is still a work in progress. Should be used only as an example of using MPFloat (see: https://gist.github.com/1562150).