diff --git a/src/math/Maths.java b/src/math/Maths.java
new file mode 100644
index 0000000..6af8a60
--- /dev/null
+++ b/src/math/Maths.java
@@ -0,0 +1,1457 @@
+package math;
+
+import java.math.BigInteger;
+import java.util.Arrays;
+
+import math.matrix.Matrix;
+import math.statistics.Probability;
+import stdlib.StdArrayIO;
+import stdlib.StdOut;
+
+/**
+ * Class to operate real numbers, x∈ℝ.
+ * @author Daniel Weschke
+ * @version October 2013
+ */
+final public class Maths{
+	
+	/**
+	 * The double value that is closer than any other to π,
+	 * the ratio of the circumference of a circle to its diameter.
+	 */
+	public static final double π = java.lang.Math.PI;
+	
+	/**
+	 * The double value that is closer than any other to e,
+	 * the base of the natural logarithms.
+	 */
+	public static final double e = java.lang.Math.E;
+	
+	/**
+	 * ε small number as relative error tolerance
+	 */
+	public static final double ε = 1.0e-15;
+	
+	/**
+	 * small number replacing zero, e. g. in LU decomposition
+	 */
+	public static final double tiny = 1.0e-100;
+	
+	/**
+	 * @param a number
+	 * @return the absolute value
+	 */
+	public static double abs(double a){ return a >= 0 ? a : -a; }
+	public static float abs(float a){ return a >= 0 ? a : -a; }
+	public static int abs(int a){ return a >= 0 ? a : -a; }
+	public static long abs(long a){ return a >= 0 ? a : -a; }
+
+	/**
+	 * Pythagoras
+	 * (a<sup>2</sup> + b<sup>2</sup>)<sup>1/2</sup> without under- or overflow.
+	 * @param a leg
+	 * @param b leg
+	 * @return hypotenuse
+	 */
+	public static double hypot(double a, double b) {
+		double r;
+		double aa = Maths.abs(a);
+		double ab = Maths.abs(b);
+		if(aa > ab){
+			r = b/a;
+			r = aa*Math.sqrt(1+r*r);
+		} else if(b != 0){
+			r = a/b;
+			r = ab*Math.sqrt(1+r*r);
+		} else
+			r = 0.0;
+		return r;
+	}
+	
+	/**
+	 * a<sup>2</sup>
+	 * @param a number
+	 * @return square of the number
+	 */
+	public static double sqare(double a){ return a*a; }
+	
+	/**
+	 * Computes the square root of a nonnegative number c using
+	 * Newton's method:<br />
+	 * First initialize t = a,<br />
+	 * than replace t with the average of c/t and t<br />
+	 * and last repeat until desired accuracy reached (ε=10<sup>-15</sup>).
+	 * @param a number
+	 * @return the estimate of the square root of a
+	 * @see #sqrt(double, double)
+	 */
+	public static double sqrt(double a){
+        return sqrt(a, ε);
+	}
+
+	/**
+	 * @param a number
+	 * @param ε relative error tolerance
+	 * @return the estimate of the square root of a
+	 * @see #sqrt(double)
+	 */
+    public static double sqrt(double a, double ε) {
+		if(a<0) return Double.NaN;
+        double t = a; // estimate of the square root of c
+        // repeatedly apply Newton update step until desired precision is achieved
+        while(Math.abs(t - a/t) > ε*t)
+            t = (a/t + t) / 2.0;
+        return t;
+    }
+	
+	public static double cube(double a){ return sqare(a)*a; }
+
+	/**
+	 * n-th power of 2.<br />
+	 * Only works if 0 <= n < 31 since 2^31 overflows an int.
+	 * @param n exponent
+	 * @return 2<sup>n</sup>
+	 */
+	public static int powOfTwo(double n){
+        int i = 0;
+        int powerOfTwo = 1;
+        while(i <= n){
+            powerOfTwo *= 2;
+            i += 1;
+        }
+        return powerOfTwo;
+	}
+	
+	/**
+	 * k-th power of n.
+	 * @param base n
+	 * @param exponent k
+	 * @return n<sup>k</sup>
+	 */
+	public static long pow(long base, long exponent){
+	    if(exponent == 0) return 1;
+	    else return base * pow(base, exponent -1);
+	}
+
+	/**
+	 * k-th power of n.
+	 * @param base n
+	 * @param exponent k
+	 * @return n<sup>k</sup>
+	 */
+	public static double pow(double base, double exponent){
+	    if(exponent == 0) return 1;
+	    else return base * pow(base, exponent -1);
+	}
+
+	/**
+	 * Determine the sign; i.e. even numbers have sign + and odds -.
+	 * used in matrix, chessboard rule.
+	 * @param i number
+	 * @return +1 for even number, -1 for odd number
+	 */
+	public static int changeSign(int i){ return i%2==0?1:-1; }
+	
+	public static double cos(double a){ return java.lang.Math.cos(a); }
+	public static double cosd(double a){ return cos(toRad(a)); }
+	
+	public static double cosh(double a){ return java.lang.Math.cosh(a); }
+    public static double cosh2(double x) {
+        return (Math.exp(x) + Math.exp(-x)) / 2.0;
+    }
+	public static double coshd(double a){ return cosh(toRad(a)); }
+	
+	public static double acos(double a){ return java.lang.Math.acos(a); }
+	public static double acosd(double a){ return toDeg(acos(a)); }
+	
+	public static double sin(double a){ return java.lang.Math.sin(a); }
+	public static double sind(double a){ return sin(toRad(a)); }
+	
+	public static double sinh(double a){ return java.lang.Math.sinh(a); }
+    public static double sinh2(double x) {
+        return (Math.exp(x) - Math.exp(-x)) / 2.0;
+    }
+	public static double sinhd(double a){ return sinh(toRad(a)); }
+	
+	public static double asin(double a){ return java.lang.Math.asin(a); }
+	public static double asind(double a){ return toDeg(asin(a)); }
+	
+	public static double tan(double a){ return java.lang.Math.tan(a); }
+	public static double tand(double a){ return tan(toRad(a)); }
+	
+	public static double tanh(double a){ return java.lang.Math.tanh(a); }
+    public static double tanh2(double x) {
+        return sinh(x) / cosh(x);
+    }
+	public static double tanhd(double a){ return tanh(toRad(a)); }
+	
+	public static double atan(double a){ return java.lang.Math.atan(a); }
+	public static double atand(double a){ return toDeg(atan(a)); }
+	
+	public static double atan(double y, double x){ return java.lang.Math.atan2(y, x); }
+	public static double atand(double y, double x){ return toDeg(atan(y, x)); }
+	
+	/**
+	 * Number of decimal places, the number of digits after the decimal point.
+	 * 0.0123 = 4
+	 * @param number
+	 * @return number of decimal places
+	 */
+	public static int decimalPlaces(double number){
+		String s = Double.toString(number);
+//		System.out.println(number+" "+s+" ");
+		s = s.substring(s.indexOf(".") + 1);
+		if(s.length()==1 && s.equals("0")) s = "";
+//		System.out.println(s.length()+" ");
+		return s.length();
+	}
+	
+	/**
+	 * Number of decimal places, the number of digits after the decimal point.
+	 * 0.0123 = 4
+	 * @param numbers
+	 * @return number of decimal places
+	 */
+	public static int decimalPlaces(Matrix numbers){
+		int i, j;
+		int decimalPlaces = 0, currentDecimalPlaces;
+		for(i=0; i<numbers.getM(); i++)
+			for(j=0; j<numbers.getN(); j++){
+				currentDecimalPlaces = decimalPlaces(numbers.get(i, j));
+				if(currentDecimalPlaces > decimalPlaces) decimalPlaces = currentDecimalPlaces;
+			}
+		return decimalPlaces;
+	}
+	
+	/**
+	 * Reverses the digits of a positive integer a using arithmetic.
+	 * @param a positive integer
+	 * @return reverse digits
+	 */
+	public static int digitReverse(int a){
+		if(a<0) throw new IllegalArgumentException();
+		int b = 0;
+		while(a != 0) {
+			b = (10 * b) + (a % 10);
+			a /= 10;
+		}
+		return b;
+	}
+	
+	/**
+	 * sum of decimal digits. 701=7+0+1=8
+	 * @param number
+	 * @return the digit sum
+	 */
+	public static int digitSum(int number){
+		int sum = 0;
+		while(number > 0) {
+			sum = sum + (number % 10);
+			number /= 10;
+		}
+		return sum;
+	}
+	
+    /**
+     * @param n
+     * @return n!
+     */
+    public static BigInteger factorial(int n) {
+        BigInteger ansBI = new BigInteger("1");
+        if(n == 0) return ansBI;
+        return new BigInteger(String.valueOf(n)).multiply(factorial(n-1));
+    }
+	
+    /**
+     * assuming n ≥ 0
+     * @param n
+     * @return n!
+     */
+    public static long factorial(long n) {
+    	if(n>20) throw new RuntimeException("cause overflow, only factorial till 20 can be calculated.");
+        if(n == 0) return 1;
+        return n * factorial(n-1);
+    }
+
+	/**
+	 * @param x double value
+	 * @return the digits after the decimal point
+	 */
+	public static double frac(double x){
+		if(x > 0.0) return x - Math.floor(x);
+		else        return x - Math.ceil(x);
+	}
+	
+	/**
+	 * Computes the greatest common divisor of p and q using Euclid's algorithm.
+	 * @param p
+	 * @param q
+	 * @return greatest common divisor
+	 */
+    public static int gcd(int p, int q) {
+        if(q == 0) return p;
+        else return gcd(q, p % q);
+    }
+    
+    /**
+     * φ = (1+√5)/2 = 1.6180339887498948482...
+     * @return golden ratio
+     */
+    public static double golden() {
+        return (1.0 + Math.sqrt(5)) / 2.;
+    }
+    
+    /**
+     * Computes an approximation to the golden ratio using the recursive
+     * formula f<sub>0</sub> = 1, f<sub>n</sub> = 1 + 1 / f<sub>n-1</sub> if n > 0.
+     * <p>
+     * n = 5: 1.625<br />
+     * n = 10: 1.6179775280898876<br />
+     * n = 20: 1.618033985017358<br />
+     * n = 30: 1.6180339887496482<br />
+     * @param n
+     * @return golden ratio
+     */
+    public static double golden(int n) {
+        if(n == 0) return 1;
+        return 1.0 + 1.0 / golden(n-1);
+    }
+	
+	/**
+	 * Computes the number of primes less than or equal to a
+	 * using the Sieve of Eratosthenes.
+	 * @param a number
+	 * @return the number of primes <= a
+	 */
+	public static int primeSieve(int a){
+		int i, j;
+		// initially assume all integers are prime
+        boolean[] isPrime = new boolean[a + 1];
+        for(i=2; i<=a; i++)
+            isPrime[i] = true;
+        // mark non-primes <= N using Sieve of Eratosthenes
+        for(i=2; i*i<=a; i++) {
+            // if i is prime, then mark multiples of i as nonprime
+            // suffices to consider mutiples i, i+1, ..., N/i
+            if (isPrime[i])
+                for(j=i; i*j<=1; j++)
+                    isPrime[i*j] = false;
+        }
+        // count primes
+        int primes = 0;
+        for(i=2; i<=a; i++)
+            if (isPrime[i]) primes++;
+        return primes;
+	}
+	
+	public static boolean kroneckerDelta(int x, int y){
+		return (x == y);
+	}
+
+	/**
+	 * @param a array
+	 * @return maximum value in array, -∞ if no such value
+	 */
+	public static int max(int... a){
+		int i;
+		int max = Integer.MIN_VALUE;
+		for(i=0; i<a.length; i++)
+		     if(a[i] > max)
+		    	 max = a[i];
+		return max;
+	}
+
+	/**
+	 * @param a double array
+	 * @return maximum value in double array, -∞ if no such value
+	 */
+	public static int max(int[][] a){
+		int i, j;
+		int max = Integer.MIN_VALUE;
+		for(i=0; i<a.length; i++)
+			for(j=0; j<a[0].length; j++)
+				if(a[i][j] > max)
+					max = a[i][j];
+		return max;
+	}
+	
+	/**
+	 * @param a array
+	 * @return maximum value in array, -∞ if no such value.
+	 */
+	public static double max(double... a){
+		int i;
+		double max = Double.NEGATIVE_INFINITY;
+		for(i=0; i<a.length; i++)
+		     if(a[i] > max)
+		    	 max = a[i];
+		return max;
+	}
+	
+	/**
+	 * @param is first index
+	 * @param ie last index
+	 * @param a array
+	 * @return maximum value in subarray a[i<sub>s</sub>...i<sub>e</sub>], -∞ if no such value.
+	 */
+	public static double max(int is, int ie, double... a) {
+		if(!inBounds(a, is, ie))
+			throw new RuntimeException("Subarray indices out of bounds");
+		int i;
+		double max = Double.NEGATIVE_INFINITY;
+		for(i=is; i<=ie; i++)
+			if(a[i] > max)
+				max = a[i];
+		return max;
+	}
+	
+	/**
+	 * @param a double array
+	 * @return maximum value in double array, -∞ if no such value.
+	 */
+	public static double max(double[][] a){
+		int i, j;
+		double max = Double.NEGATIVE_INFINITY;
+		for(i=0; i<a.length; i++)
+			for(j=0; j<a[0].length; j++)
+				if(a[i][j] > max)
+					max = a[i][j];
+		return max;
+	}
+	
+	public static int maxOrderOfMagnitude(double[] a){
+		int i, log;
+		int maxlog = Integer.MIN_VALUE;
+		for(i=0; i<a.length; i++){
+			log = (int) Math.log10(Math.abs(a[i]));
+			if(log > maxlog)
+				maxlog = log;
+			}
+		return maxlog;
+	}
+
+	/**
+	 * Sum of values of a data set divided by number of values.
+	 * @param a data field
+	 * @return average value in array, NaN if no such value
+	 */
+	public static double mean(int... a) {
+		if(a.length == 0) return Double.NaN;
+        double sum = sum(a);
+        return sum / a.length;
+	}
+	
+	/**
+	 * Sum of values of a data set divided by number of values.
+	 * @param a data field
+	 * @return average value in array, NaN if no such value
+	 */
+	public static double mean(double... a){
+        if(a.length == 0) return Double.NaN;
+        double sum = sum(a);
+        return sum / a.length;
+	}
+
+	/**
+	 * Sum of values of a data set divided by number of values.
+	 * @param is first index
+	 * @param ie last index
+	 * @param a data field
+	 * @return average value in subarray a[i<sub>s</sub>...i<sub>e</sub>], NaN if no such value.
+	 */
+	public static double mean(int is, int ie, double... a) {
+		if(!inBounds(a, is, ie))
+			throw new RuntimeException("Subarray indices out of bounds");
+		int length = ie - is + 1;
+		if(length == 0) return Double.NaN;
+		double sum = sum(a, is, ie);
+		return sum / length;
+	}
+	
+	/**
+	 * Middle value separating the greater and lesser halves of a data set.
+	 * @param a data filed
+	 * @return the median value of a data field
+	 */
+	public static double median(double... a){
+		Arrays.sort(a);
+		if((a.length+1)%2==0) return a[(a.length+1)/2-1];
+		else{
+			int pos = (int) (a.length+1)/2;
+			return (a[pos-1]+a[pos])/2;
+		}
+	}
+
+	/**
+	 * @param a array
+	 * @return minimum value in array, +∞ if no such value
+	 */
+	public static int min(int... a){
+		int i;
+		int min = Integer.MAX_VALUE;
+		for(i=0; i<a.length; i++)
+		     if(a[i] < min)
+		    	 min = a[i];
+		return min;
+	}
+
+	/**
+	 * @param a double array
+	 * @return minimum value in double array, +∞ if no such value
+	 */
+	public static int min(int[][] a){
+		int i, j;
+		int min = Integer.MAX_VALUE;
+		for(i=0; i<a.length; i++)
+			for(j=0; j<a[0].length; j++)
+				if(a[i][j] < min)
+					min = a[i][j];
+		return min;
+	}
+
+	/**
+	 * @param a array
+	 * @return minimum value in array, +∞ if no such value
+	 */
+	public static double min(double... a){
+		int i;
+		double min = Double.POSITIVE_INFINITY;
+		for(i=0; i<a.length; i++)
+		     if(a[i] < min)
+		    	 min = a[i];
+		return min;
+	}
+	
+	/**
+	 * @param is first index
+	 * @param ie last index
+	 * @param a array
+	 * @return minimum value in subarray a[i<sub>s</sub>...i<sub>e</sub>], +∞ if no such value
+	 */
+    public static double min(int is, int ie, double... a) {
+		if(!inBounds(a, is, ie))
+            throw new RuntimeException("Subarray indices out of bounds");
+        int i;
+        double min = Double.POSITIVE_INFINITY;
+        for(i=is; i<=ie; i++)
+            if (a[i] < min)
+            	min = a[i];
+        return min;
+    }
+
+	/**
+	 * @param a double array
+	 * @return minimum value in double array, +∞ if no such value
+	 */
+	public static double min(double[][] a){
+		int i, j;
+		double min = Double.POSITIVE_INFINITY;
+		for(i=0; i<a.length; i++)
+			for(j=0; j<a[0].length; j++)
+				if(a[i][j] < min)
+					min = a[i][j];
+		return min;
+	}
+	
+	/**
+	 * @param a array
+	 * @return minimum and maximum value in array, +∞ for minimum or -∞ for maximum if no such value
+	 */
+	public static double[] minmax(double... a){
+		int i;
+		double min = Double.POSITIVE_INFINITY;
+		double max = Double.NEGATIVE_INFINITY;
+		for(i=0; i<a.length; i++){
+		     if(a[i] < min)
+		    	 min = a[i];
+		     if(a[i] > max)
+		    	 max = a[i];
+		}
+		return new double[]{min,max};
+	}
+	
+	/**
+	 * @param a array
+	 * @return minimum and maximum value in array, +∞ for minimum or -∞ for maximum if no such value
+	 */
+	public static double[] minmax(double[][] a){
+		int i, j;
+		double min = Double.POSITIVE_INFINITY;
+		double max = Double.NEGATIVE_INFINITY;
+		for(i=0; i<a.length; i++)
+			for(j=0; j<a[0].length; j++){
+				if(a[i][j] < min)
+					min = a[i][j];
+				if(a[i][j] > max)
+					max = a[i][j];
+			}
+		return new double[]{min,max};
+	}
+	
+	/**
+	 * TODO:
+	 * Most frequent value in a data set.
+	 * @param a array
+	 * @return the mode of a data field
+	 */
+	public static double mode(double... a){
+//		MATLAB:
+//		X = sort(x);
+//		indices   =  find(diff([X; realmax]) > 0); % indices where repeated values change
+//		[modeL,i] =  max (diff([0; indices]));     % longest persistence length of repeated values
+//		mode      =  X(indices(i));
+		return 0;
+	}
+	
+    /**
+     * Modular exponentiation a^e mod m.
+     * Computes a^e mod m using repeated squaring.
+     * 17^5 mod 10 = 7.
+     * 43417^53535 mod 34310 = 12053.
+     * 345343417^542323535 mod 324244310 = 23504373.
+     * @param a base
+     * @param e exponent
+     * @param m modulus
+     * @return Modular exponentiation a^e mod m
+     */
+    public static int modExp(int a, int e, int m) {
+        if(e == 0) return 1;
+        long t = modExp(a, e/2, m);  // use long for intermediate computations to eliminate overflow
+        long c = (t * t) % m;
+        if(e % 2 == 1)
+           c = (c * a) % m;
+        return (int) c;
+    }
+	
+	public static double random(){ return java.lang.Math.random(); }
+	
+	public static double[] randomSeq(int n){
+		int i;
+		double[] result = new double[n];
+        for(i=0; i<n; i++)
+            result[i] = random();
+        return result;
+	}
+	
+	/**
+	 * @param a
+	 * @return pseudo-random number between 0 and a.
+	 */
+	public static double random(double a){ return random()*a; }
+	
+	public static float random(float a){ return (float) random()*a; }
+	
+	/**
+	 * @param a
+	 * @return pseudo-random number between 0 and a-1.
+	 */
+	public static int random(int a){ return (int) (random()*a); }
+	
+	public static long random(long a){ return (long) (random()*a); }
+	
+	/**
+	 * Simulate a fair coin flip.
+	 * @return true or false
+	 */
+	public static boolean flip(){ return random() < 0.5; }
+	
+	/**
+	 * Binomial coefficient. (n choose k) = n!/((n-k)!k!)
+	 * Combinations refer to the combination of n things taken k at a time without repetition.
+	 * @param n total number of things
+	 * @param k taken thinks
+	 * @return combinations without repetition
+	 * @see math.statistics.Probability#combinations(int, int)
+	 */
+	public static long binomialCoefficient(int n, int k){
+		return Probability.combinations(n,k);
+	}
+	
+	/**
+	 * @param a number of elements
+	 * @param b range from 0, 1, ..., a-1
+	 * @return a random sample of a of the integers from 0 to b-1
+	 */
+	public static int[] sample(int a, int b){
+		int i;
+        // create permutation 0, 1, ..., N-1
+        int[] perm = new int[b];
+        for(i=0; i<b; i++)
+            perm[i] = i;
+        // create random sample in perm[0], perm[1], ..., perm[M-1]
+        for(i=0; i<a; i++){
+            // random integer between i and N-1
+            int r = i + (int) random(b-1);
+            // swap elements at indices i and r
+            int t = perm[r];
+            perm[r] = perm[i];
+            perm[i] = t;
+        }
+        int[] result = new int [a];
+        for(i=0; i<a; i++)
+            result[i] = perm[i];
+        return result;
+	}
+	
+	/**
+	 * @param n number of distinct card types
+	 * @return the number of how many random cards do you need to collect before you have at least one of each type
+	 */
+	public static int couponCollector(int n){
+		boolean[] found = new boolean[n]; // found[i] = true if card i has been collected
+        int cardcnt = 0;                  // total number of cards collected
+        int valcnt = 0;                   // number of distinct cards
+        // repeatedly choose a random card and check whether it's a new one
+        while(valcnt < n) {
+            int val = (int) random(n);    // random card between 0 and N-1
+            cardcnt++;                    // we collected one more card
+            if (!found[val]) valcnt++;    // it's a new card type
+            found[val] = true;            // update found[]
+        }
+        return cardcnt; // return the total number of cards collected
+	}
+	
+	/**
+	 * @param a array
+	 * @return sample standard deviation of array, NaN if no such value.
+	 */
+    public static double stddev(int[] a) {
+        return Math.sqrt(var(a));
+    }
+    
+    /**
+     * @param a array
+     * @return sample standard deviation of array, NaN if no such value.
+     */
+    public static double stddev(double[] a) {
+        return Math.sqrt(var(a));
+    }
+    
+    /**
+     * @param a array
+     * @param is first index
+     * @param ie last index
+     * @return sample standard deviation of subarray a[i<sub>s</sub>...i<sub>e</sub>], NaN if no such value.
+     */
+    public static double stddev(double[] a, int is, int ie) {
+        return Math.sqrt(var(a, is, ie));
+    }
+    
+    /**
+     * @param a array
+     * @return population standard deviation of array, NaN if no such value.
+     */
+    public static double stddevp(double[] a) {
+        return Math.sqrt(varp(a));
+    }
+    
+    /**
+     * @param a array
+     * @param is first index
+     * @param ie last index
+     * @return population standard deviation of subarray a[i<sub>s</sub>...i<sub>e</sub>], NaN if no such value.
+     */
+    public static double stddevp(double[] a, int is, int ie) {
+        return Math.sqrt(varp(a, is, ie));
+    }
+
+	/**
+	 * @param a
+	 * @return sum of all values in array.
+	 */
+	public static int sum(int[] a) {
+		int i, sum = 0;
+		for(i=0; i<a.length; i++)
+			sum += a[i];
+		return sum;
+	}
+	
+	/**
+	 * @param a
+	 * @return sum of all values in array.
+	 */
+	public static double sum(double[] a) {
+		int i;
+		double sum = 0.0;
+		for(i=0; i<a.length; i++)
+			sum += a[i];
+		return sum;
+	}
+
+	/**
+	 * @param a array
+     * @param is first index
+     * @param ie last index
+	 * @return sum of all values in subarray a[i<sub>s</sub>...i<sub>e</sub>]
+	 */
+	public static double sum(double[] a, int is, int ie) {
+		if(!inBounds(a, is, ie))
+			throw new RuntimeException("Subarray indices out of bounds");
+		int i;
+		double sum = 0.0;
+		for(i=is; i<=ie; i++)
+			sum += a[i];
+		return sum;
+	}
+	
+	/**
+	 * @param a
+	 * @return sum of all values in array.
+	 */
+	public static double sumAbs(double[] a) {
+		int i;
+		double sum = 0.0;
+		for(i=0; i<a.length; i++)
+			sum += abs(a[i]);
+		return sum;
+	}
+	
+	/**
+	 * @param a array
+	 * @return sample variance of array, NaN if no such value
+	 */
+	public static double var(int[] a) {
+    	int n = a.length;
+		if(n == 0) return Double.NaN;
+		int i;
+		double avg = mean(a);
+		double sum = 0.0;
+		for(i=0; i<n; i++)
+			sum += (a[i] - avg) * (a[i] - avg);
+		return sum / (n - 1);
+	}
+	
+	/**
+	 * @param a array
+	 * @return sample variance of array, NaN if no such value
+	 */
+    public static double var(double[] a) {
+    	int n = a.length;
+        if(n == 0) return Double.NaN;
+		int i;
+        double avg = mean(a);
+        double sum = 0.0;
+        for(i=0; i<n; i++)
+            sum += (a[i] - avg) * (a[i] - avg);
+        return sum / (n - 1);
+    }
+    
+    /**
+     * @param a array
+     * @param is first index
+     * @param ie last index
+     * @return sample variance of subarray a[i<sub>s</sub>...i<sub>e</sub>], NaN if no such value
+     */
+    public static double var(double[] a, int is, int ie) {
+        int length = ie - is + 1;
+        if(!inBounds(a, is, ie))
+            throw new RuntimeException("Subarray indices out of bounds");
+        if(length == 0) return Double.NaN;
+        int i;
+        double avg = mean(is, ie, a);
+        double sum = 0.0;
+        for(i=is; i<=ie; i++)
+            sum += (a[i] - avg) * (a[i] - avg);
+        return sum / (length - 1);
+    }
+    
+    /**
+     * @param a array
+     * @return population variance of array, NaN if no such value
+     */
+    public static double varp(double[] a) {
+    	int n = a.length;
+        if(n == 0) return Double.NaN;
+        int i;
+        double avg = mean(a);
+        double sum = 0.0;
+        for(i = 0; i < n; i++)
+            sum += (a[i] - avg) * (a[i] - avg);
+        return sum / n;
+    }
+    
+    /**
+     * @param a array
+     * @param is first index
+     * @param ie last index
+     * @return population variance of subarray a[i<sub>s</sub>...i<sub>e</sub>],  NaN if no such value.
+     */
+    public static double varp(double[] a, int is, int ie) {
+        int length = ie - is + 1;
+        if(!inBounds(a, is, ie))
+            throw new RuntimeException("Subarray indices out of bounds");
+        if(length == 0) return Double.NaN;
+        double avg = mean(is, ie, a);
+        double sum = 0.0;
+        for (int i = is; i <= ie; i++) {
+            sum += (a[i] - avg) * (a[i] - avg);
+        }
+        return sum / length;
+    }
+    
+    /**
+     * Bubble sort
+     * @param a array
+     */
+    public static void bubbleSort(double[] a){ 
+        int i;
+        double k; 
+        for(i=0; i<a.length-1; i++){ 
+            if(a[i]<a[i+1])
+                continue;
+            k = a[i];
+            a[i] = a[i+1]; 
+            a[i+1] = k; 
+            bubbleSort(a); 
+        }
+    }
+    
+    /**
+     * Insert sort
+     * @param a array
+     */
+    public static void insertSort(double[] a){
+    	int i, j;
+        double k; 
+        for(i=0; i<a.length; i++){ 
+            for(j=a.length-1; j>0; j--){ 
+                if(a[j-1] > a[j]){ 
+                    k = a[j]; 
+                    a[j] = a[j-1]; 
+                    a[j-1] = k; 
+                } 
+            } 
+        } 
+    }
+    
+    /**
+     * Merge sort
+     * @param a array
+     * @param l left index
+     * @param r right index
+     */
+    public static void mergeSort(double[] a, int l, int r) { 
+        if (l < r){
+            int q = (l + r) / 2;
+            mergeSort(a, l, q);
+            mergeSort(a, q + 1, r);
+            merge(a, l, q, r);
+        }
+    }
+    
+    // @see #mergeSort
+    private static void merge(double[] a, int l, int q, int r) { 
+        double[] arr = new double[a.length]; 
+        int i, j, k; 
+        for(i=l; i<=q; i++)
+            arr[i] = a[i];
+        for(j=q+1; j<=r; j++)
+            arr[r+q+1-j] = a[j];
+        i = l;
+        j = r;
+        for(k=l; k<=r; k++){
+            if(arr[i] <= arr[j]){
+                a[k] = arr[i];
+                i++;
+            } else {
+                a[k] = arr[j];
+                j--;
+            }
+        }
+    }
+    
+    /**
+     * OETsort. Odd-Even-Transposition-Sort
+     * @param a array
+     */
+    public static void oetSort(double[] a){
+    	int i, j;
+        double k;
+        for(i=0; i<a.length; i++){ 
+            for(j=0+(i%2); j<a.length-1; j+=2){ 
+                if(a[j]>a[j+1]){ 
+                    k = a[j]; 
+                    a[j] = a[j+1]; 
+                    a[j+1] = k; 
+                } 
+            } 
+        }
+    }
+    
+    /**
+     * quick sort.
+     * @param a array
+     */
+    public static void quickSort(double[] a){
+    	quickSort(a, 0, a.length-1);
+    }
+    
+    /**
+     * quick sort
+     * @param a array
+     * @param l left index
+     * @param r right index
+     */
+    public static void quickSort(double[] a, int l, int r){ 
+        int q; 
+        if(l < r){ 
+            q = quickSortPartition(a, l, r); 
+            quickSort(a, l, q); 
+            quickSort(a, q + 1, r); 
+        }
+    }
+
+    /**
+     * define partition index of quick sort algorithm.
+     * @param a array
+     * @param l left index
+     * @param r right index
+     * @return partitioned index
+     */
+    static int quickSortPartition(double[] a, int l, int r) { 
+        int i, j;
+        double x = a[(l+r)/2]; 
+        i = l - 1;
+        j = r + 1;
+        do{ i++; } while(a[i] < x);
+        do{ j--; } while(a[j] > x);
+        if(i<j){
+            double k = a[i]; 
+            a[i] = a[j];
+            a[j] = k;
+        } else
+            return j;
+        return -1;
+    }
+    
+    /**
+     * Ripple sort
+     * @param a array
+     */
+    public static void rippleSort(double[] a){ 
+        boolean switched;
+        int i;
+        double k;
+        do { 
+            switched = false; 
+            for(i=0; i<a.length-1; i++){ 
+                if(a[i]>a[i+1]){ 
+                    k = a[i]; 
+                    a[i] = a[i+1]; 
+                    a[i+1] = k; 
+                    switched = true; 
+                } 
+            } 
+        } while (switched == true);
+    }
+    
+    /**
+     * Select sort
+     * @param a array
+     */
+    public static void selectSort(double[] a){
+    	int i, j, q;
+        double k;
+        for(i=a.length-1; i>=1; i--){
+            q = 0;
+            for(j=1; j<=i; j++)
+                if(a[j] > a[q])
+                    q = j;
+            k = a[q];
+            a[q] = a[i];
+            a[i] = k;
+        }
+    }
+    
+    /**
+     * Shaker sort
+     * @param a array
+     */
+    public static void shakerSort(double[] a){
+        int i = 0, l = a.length; 
+        while(i < l) { 
+            shaker1(a, i, l); 
+            l--; 
+            shaker2(a, i, l); 
+            i++; 
+        }
+    } 
+
+    private static void shaker1(double[] a, int i, int l){
+    	int j;
+    	double k;
+        for(j=i; j<l-1; j++){  // increment from lower index
+            if(a[j] > a[j+1]){
+                k = a[j]; 
+                a[j] = a[j+1]; 
+                a[j+1] = k; 
+            } 
+        } 
+    } 
+
+    private static void shaker2(double[] a, int i, int l) { 
+    	int j;
+    	double k;
+        for(j=l-1; j>=i; j--){ // decrement from upper index
+            if(a[j] > a[j+1]){ 
+                k = a[j]; 
+                a[j] = a[j+1]; 
+                a[j+1] = k; 
+            } 
+        } 
+    } 
+    
+    /**
+     * complexity ~ O(N<sup>3/2</sup>)
+     * @param a array
+     */
+    public static void shellSort(double[] a){
+    	int i, j, k, N = a.length;
+    	double v;
+    	for(k=1; k<=N/9; k=3*k+1) ;
+    	for( ; k>0; k/=3)
+    		for(i=k; i<=N; i++){
+    			v = a[i-1];
+    			j = i-1;
+    			while(j>k-1 && a[j-k]>v){
+    				a[j] = a[j-k];
+    				j -= k;
+    			}
+    			a[j] = v;
+    		}
+    }
+    
+    /**
+     * Simple sort
+     * @param a array
+     */
+    public static void simpleSort(double[] a){
+    	int i, j;
+        double k;
+        for(i=a.length-1; i>1; i--){
+            for(j=0; j<=i-1; j++){
+                if(a[j] >= a[i]){
+                    k = a[i];
+                    a[i] = a[j];
+                    a[j] = k;
+                }
+            }
+        }
+    }
+
+    public static void print(double[] a){
+    	System.out.print("[");
+    	int i;
+    	for(i=0; i<a.length; i++){
+    		System.out.print(a[i]);
+    		if(i!=(a.length-1)) System.out.print(", ");
+    		else System.out.println("]");
+    	}
+    }
+    public static void print(double[][] a){
+    	System.out.println(a.length+" "+a[0].length);
+    	System.out.print("[");
+    	int i, j;
+    	for(i=0; i<a.length; i++){
+    		for(j=0; j<a[i].length; j++){
+    			if(i!=0 && j==0) System.out.print(" ");
+    			System.out.print(a[i][j]);
+    			if(j!=(a[i].length-1)) System.out.print(", ");
+    			else if(i!=(a.length)-1) System.out.println(";");
+    			else System.out.println("]");
+    		}
+    	}
+    }
+    public static void printF(double[] a){
+    	System.out.println();
+    	int i;
+    	for(i=0; i<a.length; i++)
+    		System.out.println("a["+i+"]= "+a[i]);
+    }
+    public static void printF(double[][] a){
+    	System.out.println();
+    	int i, j;
+    	for(i=0; i<a.length; i++)
+        	for(j=0; j<a[i].length; j++)
+        		System.out.println("a["+i+"]["+j+"]= "+a[i][j]);
+    }
+    
+    /**
+     * Determines the check digit of an ISBN number given the first 9 digits.
+     * @param isbn
+     * @return check digit
+     * @see #isbn(String)
+     */
+    public static int isbnCheckDigit(String isbn){
+    	int n = Integer.parseInt(isbn);
+        // compute the weighted sum of the digits, from right to left
+        int i, sum = 0, digit;
+        for(i=2; i<=10; i++){
+           digit = n % 10; // rightmost digit
+           sum = sum + i * digit;
+           n = n / 10;
+        }
+        // return check digit, 10 (X)
+        if     (sum % 11 == 1) return 10;
+        else if(sum % 11 == 0) return 0;
+        else                   return (11 - (sum % 11));
+    }
+    
+    /**
+     * An ISBN number is legal if it consists of 10 digits and
+     * d<sub>1</sub> +2*d<sub>2</sub> + 3*d<sub>3</sub> + ... + 10*d<sub>10</sub>
+     * is a multiple of 11.
+     * For example, 0-201-31452-5 is legal since
+     * 1*5 + 2*2 + 3*5 + 4*4 + 5*1 + 6*3 + 7*1 + 8*0 + 9*2 + 10*0 = 88
+     * and 88 is a multiple of 11.
+     * @param isbn number given the first 9 digits
+     * @return full isbn number, with check digit
+     * @see #isbnCheckDigit(String)
+     */
+    public static String isbn(String isbn){
+    	System.out.println(isbn);
+        int checkDigit = isbnCheckDigit(isbn);
+        // return check digit, use X for 10
+        String result = "" + isbn;
+        if(checkDigit == 10) result += "X";
+        else                 result += checkDigit;
+        return result;
+    }
+    
+    /**
+     * Determines if a VIN number is valid by computing its check digit.
+     * Do aggressive error checking.
+     * <p>
+     * 1B4YEM9P4KP186543 : Invalid<br />
+     * 1FA-CP45E-X-LF192944 : Valid<br />
+     * 1FA-CP45E-6-LF192944 : Invalid<br />
+     * QFA-CP45E-X-LF192944 :
+     * Exception in thread "main" java.lang.RuntimeException: Illegal character: Q<br />
+     * 1FA-CP45E-G-LF192944 :
+     * Exception in thread "main" java.lang.RuntimeException: Illegal check digit: G<br />
+     * 1FA-CP45E-X-LF19294 :
+     * Exception in thread "main" java.lang.RuntimeException: VIN number must be 17 characters<br />
+     * @param s VIN number
+     * @return boolean true: valid, false: invalid
+     */
+    public static boolean vin(String s){
+        int[] values = { 1, 2, 3, 4, 5, 6, 7, 8, 0, 1,
+                         2, 3, 4, 5, 0, 7, 0, 9, 2, 3,
+                         4, 5, 6, 7, 8, 9 };
+        int[] weights = { 8, 7, 6, 5, 4, 3, 2, 10, 0, 9,
+                          8, 7, 6, 5, 4, 3, 2 };
+        s = s.replaceAll("-", "");
+        s = s.toUpperCase();
+        if(s.length() != 17)
+            throw new RuntimeException("VIN number must be 17 characters");
+        int sum = 0;
+        int i;
+        for(i=0; i<17; i++){
+            char c = s.charAt(i);
+            int value;
+            int weight = weights[i];
+            // letter
+            if(c >= 'A' && c <= 'Z') {
+                value = values[c - 'A'];
+                if(value == 0)
+                    throw new RuntimeException("Illegal character: " + c);
+            }
+            // number
+            else if(c >= '0' && c <= '9') value = c - '0';   
+            // illegal character
+            else throw new RuntimeException("Illegal character: " + c);
+            sum = sum + weight * value;
+        }
+        // check digit
+        sum = sum % 11;
+        char check = s.charAt(8);
+        if(check != 'X' && (check < '0' || check > '9'))
+            throw new RuntimeException("Illegal check digit: " + check);
+        if     (sum == 10 && check == 'X') return true; // System.out.println("Valid");
+        else if(sum == check - '0')        return true; // System.out.println("Valid");
+        else                               return false; // System.out.println("Invalid");
+    }
+    
+    /**
+     * @param n oder of matrix
+     * @return Hadamard matrix of order N. Assumes N is a power of 2.
+     */
+    public static boolean[][] hadarmad(int n){
+    	int i, j, k;
+        boolean[][] H = new boolean[n][n];
+        // initialize Hadamard matrix of order N
+        H[0][0] = true;
+        for(k=1; k<n; k+=k){
+           for(i=0; i<k; i++){
+              for(j=0; j<k; j++){
+                 H[i+k][j]   =  H[i][j];
+                 H[i][j+k]   =  H[i][j];
+                 H[i+k][j+k] = !H[i][j];
+              }
+           }
+        }
+        return H;
+    }
+    
+    public static boolean isTrue(int i){ return (i!=0); }
+    public static int isTrue(boolean b){ return b?1:0; }
+    
+    public static boolean isOdd(int n) {
+        return n % 2 != 0;
+    }
+
+	/**
+	 * @param a array
+	 * @param is first index of array
+	 * @param ie last index of array
+	 * @return boolean true if i<sub>s</sub> is greater than i<sub>e</sub> and 0 and
+	 * if i<sub>e</sub> is lower than array.length
+	 */
+	public static boolean inBounds(int[] a, int is, int ie){
+        return(is>0 && ie<a.length && is<ie);
+	}
+
+	/**
+	 * @param a array
+	 * @param is first index of array
+	 * @param ie last index of array
+	 * @return boolean true if i<sub>s</sub> is greater than i<sub>e</sub> and 0 and
+	 * if i<sub>e</sub> is lower than array.length
+	 */
+	public static boolean inBounds(double[] a, int is, int ie){
+        return(is>0 && ie<a.length && is<ie);
+	}
+	
+	/**
+	 * @param year
+	 * @return true if year corresponds to a leap year, and false otherwise.
+	 * Assumes year >= 1582, corresponding to a year in the Gregorian calendar.
+	 */
+	public static boolean isLeapYear(double year){
+        // divisible by 4 and not 100 unless divisible by 400
+        return (((year % 4 == 0) && (year % 100 != 0)) || (year % 400 == 0));
+	}
+	
+	/**
+	 * @param i value
+	 * @param j value
+	 * @return true or false if i divides j or j divides i
+	 */
+	public static boolean itDivides(double i, double j){
+		return (i % j == 0 || j % i == 0);
+	}
+
+	public static boolean isPrime(int number){
+        boolean isPrime = true;
+        if (number < 2) isPrime = false;
+
+        // try all possible factors i of N
+        // if if N has a factor, then it has one less than or equal to sqrt(N),
+        // so for efficiency we only need to check i <= sqrt(N) or equivalently i*i <= N
+        for(long i = 2; i*i <= number; i++) {
+
+            // if i divides evenly into N, N is not prime, so break out of loop
+            if(number % i == 0) {
+                isPrime = false;
+                break;
+            }
+        }
+        return isPrime;
+	}
+	
+	/**
+	 * Angle from radians to degrees 
+	 * @param angrad an angle, in radians 
+	 * @return degrees
+	 */
+	public static double toDeg(double angrad){ return angrad*180/π; }
+	
+	/**
+	 * Angle from degrees to radians 
+	 * @param angdeg - an angle, in degrees 
+	 * @return radians
+	 */
+	public static double toRad(double angdeg){ return angdeg*π/180; }
+
+	/**
+	 * @param a number
+	 * @return a in binary
+	 * @see Integer#toBinaryString(int)
+	 */
+	public static int toBinary(int a){
+		boolean isNegative = a<0;
+		if(isNegative) a *= -1;
+        // set v to the largest power of two that is <= a
+        int v = 1;
+        String result = "";
+        while(v <= a/2)
+            v = v * 2;
+        // check for presence of powers of 2 in a, from largest to smallest
+        while(v > 0){
+            if(a < v) // v is not present in n 
+            	result += "0";
+            else { // v is present in n, so remove v from n
+            	result += "1";
+                a = a - v;
+            }
+            v = v / 2; // next smallest power of 2
+        }
+		return Integer.valueOf(result);
+	}
+	
+	/**
+	 * @param args
+	 */
+	public static void main(String[] args) {
+		double[] numbers = new double[]{26.1, 25.6, 25.7, 25.2, 25.0};
+        System.out.println(Arrays.toString(numbers));
+        Arrays.sort(numbers);
+        System.out.println(Arrays.toString(numbers));
+		double[] numbers2 = new double[]{26.1, 25.6, 25.7, 25.2, 25.0, 24.7};
+
+        System.out.println(toBinary(16));
+        
+        System.out.println(Arrays.toString(sample(6, 49)));
+        
+        System.out.println("sqrt(5) = "+sqrt(5));
+        System.out.println("sqrt(5) = "+Math.sqrt(5));
+        
+        System.out.println("hypo/sqrt(3^2+4^2) = "+hypot(3, 4)+" : "+(hypot(3, 4)==5.));
+        
+        System.out.println("The full ISBN number of 020131452 is " + isbn("020131452"));
+        
+		double mean;
+        mean = Maths.mean(numbers);
+        System.out.println(mean);
+        
+        double median;
+        median = Maths.median(numbers);
+        System.out.println(median);
+        median = Maths.median(numbers2);
+        System.out.println(median);
+
+        System.out.println(decimalPlaces(10));
+        System.out.println(decimalPlaces(0.01));
+        
+        System.out.println(factorial(20));
+        System.out.println(factorial(40));
+        
+    	double[] a = new double[]{0.1, 0.4, 0.3};
+    	// Convert command-line arguments to array of doubles and call various methods.
+    	if(args.length!=0) a = StdArrayIO.readDouble1D(); // Maths.java < 0.1 0.4 0.3
+        StdOut.printf("       min %7.3f\n", min(a));
+        StdOut.printf("      mean %7.3f\n", mean(a));
+        StdOut.printf("       max %7.3f\n", max(a));
+        StdOut.printf("   std dev %7.3f\n", stddev(a));
+
+        System.out.println("\nShellSort demo");
+        java.util.Random rg = new java.util.Random();
+        a = new double[30];
+        int i;
+        for(i=0; i<a.length; i++)
+        	a[i] = rg.nextDouble();
+        printF(a);
+        shellSort(a);
+        printF(a);
+	}
+}