samer@0: /*
samer@0:  *	Mathx.java	
samer@0:  *
samer@0:  *	Copyright (c) 2000, Samer Abdallah, King's College London.
samer@0:  *	All rights reserved.
samer@0:  *
samer@0:  *	This software is provided AS iS and WITHOUT ANY WARRANTY; 
samer@0:  *	without even the implied warranty of MERCHANTABILITY or 
samer@0:  *	FITNESS FOR A PARTICULAR PURPOSE.
samer@0:  */
samer@0: 
samer@0: package samer.maths;
samer@0: 
samer@0: /**
samer@0: 	A Class full of static functions, mainly for manipulating
samer@0: 	double arrays.
samer@0:   */
samer@0: 
samer@0: public class Mathx
samer@0: {
samer@0: 	private static double [] zeros = null;
samer@0: 	private static int zlength = 0;
samer@0: 
samer@0: 	public double [] doubleArray(int n) { return new double[n]; }
samer@0: 	
samer@0: 	/** set to zeros */
samer@0: 	public static void zero(double [] x)
samer@0: 	{
samer@0: 		if (zlength<x.length) {
samer@0: 			zlength = x.length;
samer@0: 			zeros = new double[zlength];
samer@0: 			for (int i=0; i<zlength; i++) zeros[i]=0;
samer@0: 		}
samer@0: 		System.arraycopy(zeros,0,x,0,x.length);
samer@0: 	}
samer@0: 
samer@0: 	/** set all to given value */
samer@0: 	public static void setAll(double [] x, double a) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=a;
samer@0: 	}
samer@0: 	
samer@0: 	/** log in place */
samer@0: 	public static void log(double [] x) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=Math.log(x[i]);
samer@0: 	}
samer@0: 
samer@0: 	/** negate in place */
samer@0: 	public static void negate(double [] x) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=-x[i];
samer@0: 	}
samer@0: 
samer@0: 	/** negate to y: x[i]=y[i] for all i */
samer@0: 	public static void negate(double [] x, double [] y) {
samer@0: 		for (int i=0; i<x.length; i++) y[i]=-x[i];
samer@0: 	}
samer@0: 
samer@0: 	/** subtract: x[i]=y[i]-z[i] for all i */
samer@0: 	public static void sub(double [] x, double [] y, double [] z) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=y[i]-z[i];
samer@0: 	}
samer@0: 
samer@0: 	public static void sub(int n, double [] x, double [] y, double [] z) {
samer@0: 		for (int i=0; i<n; i++) x[i]=y[i]-z[i];
samer@0: 	}
samer@0: 
samer@0: 	/** subtract in place y: x[i]-=y[i] for all i */
samer@0: 	public static void sub(double [] x, double [] y) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]-=y[i];
samer@0: 	}
samer@0: 	
samer@0: 	public static void sub(int n, double [] x, double [] y) {
samer@0: 		for (int i=0; i<n; i++) x[i]-=y[i];
samer@0: 	}
samer@0: 
samer@0: 	
samer@0: 	/** subtract scalar in place y: x[i]-=y for all i */
samer@0: 	public static void sub(double [] x, double y) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]-=y;
samer@0: 	}
samer@0: 
samer@0: 	/** add two arrays: x[i]=y[i]+z[i] for all i */
samer@0: 	public static void add(double [] x, double [] y, double [] z) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=y[i]+z[i];
samer@0: 	}
samer@0: 
samer@0: 	/** add in place y: x[i]+=y[i] for all i */
samer@0: 	public static void add(double [] x, double [] y) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]+=y[i];
samer@0: 	}
samer@0: 
samer@0: 	public static void add(int n, double [] x, double [] y) {
samer@0: 		for (int i=0; i<n; i++) x[i]+=y[i];
samer@0: 	}
samer@0: 
samer@0: 	/** add scalar in place y: x[i]+=y for all i */
samer@0: 	public static void add(double [] x, double y) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]+=y;
samer@0: 	}
samer@0: 
samer@0: 	/** return sum of elements in x */
samer@0: 	public static double sum(double [] x)
samer@0: 	{
samer@0: 		double total=0;
samer@0: 		for (int i=0; i<x.length; i++) total+=x[i];
samer@0: 		return total;
samer@0: 	}
samer@0: 
samer@0: 	public static void muladd(int n, double [] x, double [] y, double k) {
samer@0: 		for (int i=0; i<n; i++) x[i]+=k*y[i];
samer@0: 	}
samer@0: 
samer@0: 	/** return product of elements in x */
samer@0: 	public static double prod(double [] x)
samer@0: 	{
samer@0: 		double total=1;
samer@0: 		for (int i=0; i<x.length; i++) total*=x[i];
samer@0: 		return total;
samer@0: 	}
samer@0: 
samer@0: 	/** multiply by scalar in place: x[i] *=k */
samer@0: 	public static void mul(double [] x, double k) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]*=k;
samer@0: 	}
samer@0: 
samer@0: 	public static void mul(int n, double [] x, double k) {
samer@0: 		for (int i=0; i<n; i++) x[i]*=k;
samer@0: 	}
samer@0: 
samer@0: 	/** scalar multiply  x[i] = k*y[i] */
samer@0: 	public static void mul(double [] x, double [] y, double k) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=k*y[i];
samer@0: 	}
samer@0: 
samer@0: 	/** elementwise multiply in place: x[i] *=y[i] */
samer@0: 	public static void mul(double [] x, double [] y) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]*=y[i];
samer@0: 	}
samer@0: 
samer@0: 	/** elementwise multiply  x[i] =y[i]*z[i] */
samer@0: 	public static void mul(double [] x, double [] y, double [] z) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=y[i]*z[i];
samer@0: 	}
samer@0: 
samer@0: 	/** elementwise divide in place: x[i] /=y[i] */
samer@0: 	public static void div(double [] x, double [] y) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]/=y[i];
samer@0: 	}
samer@0: 
samer@0: 	/** elementwise divide: x[i]=y[i]/z[i] */
samer@0: 	public static void div(double [] x, double [] y, double [] z) {
samer@0: 		for (int i=0; i<x.length; i++) x[i]=y[i]/z[i];
samer@0: 	}
samer@0: 
samer@0: 	/** return maximum element of x */
samer@0: 	public static double max(double [] x)
samer@0: 	{
samer@0: 		double xmax=x[0];
samer@0: 		for (int i=1; i<x.length; i++) if (x[i]>xmax) xmax=x[i];
samer@0: 		return xmax;
samer@0: 	}
samer@0: 
samer@0: 	/** return minimum element of x */
samer@0: 	public static double min(double [] x)
samer@0: 	{
samer@0: 		double xmin=x[0];
samer@0: 		for (int i=1; i<x.length; i++) if (x[i]<xmin) xmin=x[i];
samer@0: 		return xmin;
samer@0: 	}
samer@0: 	
samer@0: 	/** return index of supremum of elements in x */
samer@0: 	public static int argmax(double [] x)
samer@0: 	{
samer@0: 		double xmax=x[0];
samer@0: 		int	 imax=0;
samer@0: 		for (int i=1; i<x.length; i++) {
samer@0: 			if (x[i]>xmax) { xmax=x[i]; imax=i; }
samer@0: 		}
samer@0: 		return imax;
samer@0: 	}
samer@0: 
samer@0: 	/** return dot product of a and b: sum_i a[i]*b[i] */
samer@0: 	public static double dot( double [] a, double [] b)
samer@0: 	{
samer@0: 		double t=0;
samer@0: 		for (int i=0; i<a.length; i++) t+=a[i]*b[i];
samer@0: 		return t;
samer@0: 	}
samer@0: 
samer@0: 	/** return 2-norm of a: sum_i a[i]*a[i] */
samer@0: 	public static double norm2( double [] a)
samer@0: 	{
samer@0: 		double t=0;
samer@0: 		for (int i=0; i<a.length; i++) t+=a[i]*a[i];
samer@0: 		return t;
samer@0: 	}
samer@0: 
samer@0: 	public static double norm2( int n, double [] a)
samer@0: 	{
samer@0: 		double t=0;
samer@0: 		for (int i=0; i<n; i++) t+=a[i]*a[i];
samer@0: 		return t;
samer@0: 	}
samer@0: 
samer@0: 	/** copy a to b: b[i]=a[i] */
samer@0: 	public static void copy( double [] a, double [] b) {
samer@0: 		System.arraycopy(a,0,b,0,a.length);
samer@0: 	}
samer@0: 
samer@0: 	/** return copy a in new array */
samer@0: 	public static double [] copy( double [] a) {
samer@0: 		double [] b = new double[a.length];
samer@0: 		System.arraycopy(a,0,b,0,a.length);
samer@0: 		return b;
samer@0: 	}
samer@0: 
samer@0: 	public static double tanh(double t) {
samer@0: //		if (t<-24) return -1;
samer@0: //		if (t>24) return 1;
samer@0: 		if (t<0) { double z=Math.exp(2*t); return (z-1)/(z+1); }
samer@0: 		else     { double z=Math.exp(-2*t); return (1-z)/(1+z); }
samer@0: 	}
samer@0: 	
samer@0: 	public static double atanh(double t) { return Math.log((1+t)/(1-t))/2; }
samer@0: 	public static double cosh(double t) { double z=Math.exp(-t); return (1/z+z)/2; }
samer@0: 	public static double sinh(double t) { double z=Math.exp(-t); return (1/z-z)/2; }
samer@0: 
samer@0: 	/** matrix multiplication: y[i] = sum_j A[i][j] x[j] */
samer@0: 	public static void mul( double [] y, double [][]A, double [] x)
samer@0: 	{
samer@0: 		int n=y.length;
samer@0: 		int m=x.length;
samer@0: 		double [] arow;
samer@0: 
samer@0: 		for (int i=0; i<n; i++) {
samer@0: 			double t=0;
samer@0: 			arow=A[i];
samer@0: 			for (int j=0; j<m; j++) t+=arow[j]*x[j];
samer@0: 			y[i]=t;
samer@0: 		}
samer@0: 	}
samer@0: 
samer@0: 
samer@0: 	/** set each elements of a by extracting successive values from g */
samer@0: 	public static void set( double[] a, Generator g) {
samer@0: 		g.next(a);
samer@0: 		// int n=a.length; for (int i=0; i<n; i++) a[i]=g.next();
samer@0: 	}
samer@0: 
samer@0: 	/** increment each element of a by a value extracted from g */
samer@0: 	public void add( double[] a, Generator g) {
samer@0: 		int n=a.length; for (int i=0; i<n; i++) a[i]+=g.next();
samer@0: 	}
samer@0: 
samer@0: 	/** set each element using a function of its index */
samer@0: 	public static void index( double[] a) {
samer@0: 		for (int i=0; i<a.length; i++) a[i]=(double)i;
samer@0: 	}
samer@0: 
samer@0: 	public static void indexFr( double[] a) {
samer@0: 		int N=a.length; for (int i=0; i<N; i++) a[i]=(double)i/N;
samer@0: 	}
samer@0: 
samer@0: 	/** set each element using a function of its index */
samer@0: 	public static void table( double[] a, Function f) {
samer@0: 		index(a); f.apply(a);
samer@0: 	}
samer@0: }
samer@0: 
samer@0: