/*
 * Copyright 2009 Yuichiro Moriguchi
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package net.morilib.lisp;

import java.math.BigInteger;

public final class LispMath {
	
	private static final LispNumber TWO_I =
		LispComplex.newComplex(0.0, 2.0);

	private static final LispNumber N_INF_PI_HALF_I =
		LispComplex.newComplex(Double.NEGATIVE_INFINITY, -Math.PI / 2);
	
	private static final Double P_ZERO = new Double( 0.0);
	
	private static final Double M_ZERO = new Double(-0.0);
	
	private LispMath() {
		// do nothing
	}

	private static LispNumber exptC(LispNumber nm1, LispNumber nm2) {
		double a1 = nm1.getRealDouble();
		double b1 = nm1.getImagDouble();
		double r1 = Math.hypot(a1, b1);
		double t1 = Math.atan2(b1, a1);
		double a2 = nm2.getRealDouble();
		double b2 = nm2.getImagDouble();
		double rr = a2 * Math.log(r1) - t1 * b2;
		double tr = b2 * Math.log(r1) + t1 * a2;
		
		return LispComplex.newComplex(
				Math.exp(rr) * Math.cos(tr),
				Math.exp(rr) * Math.sin(tr));
	}

	public static boolean isPlusZero(double r) {
		return (r == 0.0) && P_ZERO.equals(new Double(r));
	}

	public static boolean isMinusZero(double r) {
		return (r == 0.0) && M_ZERO.equals(new Double(r));
	}
	
	public static LispNumber expt(LispNumber nm1, LispNumber nm2) {
		if(nm1.isNaN() || nm2.isNaN()) {
			return LispDouble.NaN;
		} else if(nm2.isZero()) {
			if(nm1.isExact() && nm2.isExact()) {
				return LispInteger.ONE;
			} else {
				return LispDouble.ONE;
			}
		} else if(nm1.isZero()) {
			if(nm1.isExact() && nm2.isExact()) {
				return LispInteger.ZERO;
			} else {
				return LispDouble.ZERO;
			}
		} else if(nm1.isOne()) {
			if(nm1.isExact() && nm2.isExact()) {
				return LispInteger.ONE;
			} else {
				return LispDouble.ONE;
			}
		} else if(nm2.isOne()) {
			if(nm1.isExact() && nm2.isExact()) {
				return nm1;
			} else {
				return nm1.toInexact();
			}
		} else if(!nm2.isReal()) {
			return exptC(nm1, nm2);
		} else if(!nm1.isReal()) {
			return LispMath.expt(nm1, nm2.getRealDouble());
		} else if(!nm2.isInteger()) {
			if(nm1.getRealDouble() > 0) {
				return new LispDouble(Math.pow(
						nm1.getRealDouble(), nm2.getRealDouble()));
			} else {
				return LispMath.expt(nm1, nm2.getRealDouble());
			}
		} else if(!nm1.isExact() || !nm2.isExact()) {
			return new LispDouble(Math.pow(
					nm1.getRealDouble(), nm2.getRealDouble()));
		} else if(nm1.isInteger()) {
			LispReal rm2 = nm2.getReal();
			BigInteger b = rm2.getBigInteger();
			Integer b2 = IntLispUtils.toIntExact(b);
			BigInteger a = nm1.getReal().getBigInteger();
			
			if(b2 == null) {
				return new LispDouble(Math.pow(
						nm1.getRealDouble(), nm2.getRealDouble()));
			}
			
			if(b.signum() > 0) {
				return LispInteger.valueOf(a.pow(b2.intValue()));
			} else {
				return LispRational.newRational(
						BigInteger.ONE,
						a.pow(-b2.intValue()));
			}
		} else if(nm1.isRational()) {
			LispReal rm2 = nm2.getReal();
			BigInteger b = rm2.getBigInteger();
			Integer b2 = IntLispUtils.toIntExact(b);
			BigInteger an = nm1.getReal().getNumerator();
			BigInteger ad = nm1.getReal().getDenominator();
			
			if(b2 == null) {
				return new LispDouble(Math.pow(
						nm1.getRealDouble(), nm2.getRealDouble()));
			}
			
			if(b.signum() > 0) {
				return LispRational.newRational(
						an.pow(b2.intValue()),
						ad.pow(b2.intValue()));
			} else {
				return LispRational.newRational(
						ad.pow(-b2.intValue()),
						an.pow(-b2.intValue()));
			}
		}
		return exptC(nm1, nm2);
	}

	public static LispNumber expt(LispNumber nm1, double r) {
		double a1 = nm1.getRealDouble();
		double b1 = nm1.getImagDouble();
		
		if(nm1.isNaN() || Double.isNaN(r)) {
			return LispDouble.NaN;
		} else if(b1 == 0.0 && a1 == 1.0) {
			return LispDouble.ONE;
		} else if(r == Double.POSITIVE_INFINITY) {
			double r1 = Math.hypot(a1, b1);
			
			if(r1 < 1.0 && r1 > -1.0) {
				return LispDouble.ZERO;
			} else if(b1 == 0.0 && a1 > 1.0) {
				return LispDouble.POSITIVE_INFINITY;
			} else {
				return LispDouble.NaN;
			}
		} else if(r == Double.NEGATIVE_INFINITY) {
			double r1 = Math.hypot(a1, b1);
			
			if(r1 > 1.0 || r1 < -1.0) {
				return LispDouble.ZERO;
			} else if(b1 == 0.0 && a1 < 1.0 && a1 > 0.0) {
				return LispDouble.POSITIVE_INFINITY;
			} else {
				return LispDouble.NaN;
			}
		} else if(b1 == 0.0 && a1 > 0.0) {
			return new LispDouble(Math.pow(a1, r));
		} else {
			double r1 = Math.hypot(a1, b1);
			double t1 = Math.atan2(b1, a1);
			double a2 = r;
			double rr = a2 * Math.log(r1);
			double tr = t1 * a2;
			
			return LispComplex.newComplex(
					Math.exp(rr) * Math.cos(tr),
					Math.exp(rr) * Math.sin(tr));
		}
	}

	public static LispNumber sqrt(LispNumber n) {
		if(n.isReal()) {
			if(n.getRealDouble() == 0.0) {
				return new LispDouble(0.0);
			} else if(n.getRealDouble() > 0.0) {
				return new LispDouble(Math.sqrt(n.getRealDouble()));
			} else {
				return LispComplex.newComplex(
						0.0, Math.sqrt(-n.getRealDouble()));
			}
		} else {
			return LispMath.expt(n, 0.5);
		}
	}

	public static LispNumber exp(LispNumber n) {
		if(n.isReal()) {
			return new LispDouble(Math.exp(n.getRealDouble()));
		} else {
			double a = n.getRealDouble();
			double b = n.getImagDouble();
			
			return LispComplex.newComplex(
					Math.exp(a) * Math.cos(b),
					Math.exp(a) * Math.sin(b));
		}
	}

	public static LispNumber log(LispNumber n) {
		if(n.getReal().isZero() && isMinusZero(n.getImagDouble())) {
			return N_INF_PI_HALF_I;
		} else if(n.isReal() && n.getReal().signum() > 0) {
			return new LispDouble(Math.log(n.getRealDouble()));
		} else {
			double a = n.getRealDouble();
			double b = n.getImagDouble();
			double r = Math.hypot(a, b);
			double t = Math.atan2(b, a);
			
			return LispComplex.newComplex(Math.log(r), t);
		}
	}

	public static LispNumber sin(LispNumber n) {
		if(n.isReal()) {
			return new LispDouble(Math.sin(n.getRealDouble()));
		} else if(n.getReal().isZero()) {
			double y = n.getImagDouble();
			double b = Math.sinh(y);
			
			return LispComplex.newComplex(0.0, b);
		} else {
			double x = n.getRealDouble();
			double y = n.getImagDouble();
			double a = Math.sin(x) * Math.cosh(y);
			double b = Math.cos(x) * Math.sinh(y);
			
			return LispComplex.newComplex(a, b);
		}
	}

	public static LispNumber cos(LispNumber n) {
		if(n.isReal()) {
			return new LispDouble(Math.cos(n.getRealDouble()));
		} else if(n.getReal().isZero()) {
			double y = n.getImagDouble();
			double a = Math.cosh(y);
			
			return LispComplex.newComplex(a, 0.0);
		} else {
			double x = n.getRealDouble();
			double y = n.getImagDouble();
			double a = Math.cos(x) * Math.cosh(y);
			double b = -(Math.sin(x) * Math.sinh(y));
			
			return LispComplex.newComplex(a, b);
		}
	}

	public static LispNumber tan(LispNumber n) {
		if(n.isReal()) {
			return new LispDouble(Math.tan(n.getRealDouble()));
		} else if(n.getReal().isZero()) {
			double y = n.getImagDouble();
			
			return LispComplex.newComplex(0.0, Math.tanh(y));
		} else if(n.getImagDouble() == Double.POSITIVE_INFINITY) {
			double x = n.getRealDouble();
			LispNumber nm, dn;
			
			nm = LispComplex.newComplex(Math.tan(x), 1);
			dn = LispComplex.newComplex(1, Math.tan(x));
			return nm.div(dn);
		} else if(n.getImagDouble() == Double.NEGATIVE_INFINITY) {
			double x = n.getRealDouble();
			LispNumber nm, dn;
			
			nm = LispComplex.newComplex(Math.tan(x), -1);
			dn = LispComplex.newComplex(1, -Math.tan(x));
			return nm.div(dn);
		} else {
			return sin(n).div(cos(n));
		}
	}

	public static LispNumber asin(LispNumber n) {
		double x = n.getRealDouble();
		double y = n.getImagDouble();
		
		if(n.isReal() && x >= -1.0 && x <= 1.0) {
			return new LispDouble(Math.asin(n.getRealDouble()));
		} else if(Double.isInfinite(x) && Double.isInfinite(y)) {
			return LispDouble.NaN;
		} else if(x == Double.POSITIVE_INFINITY) {
			if(y <= 0) {
				x = -x;
			}
			return LispComplex.newComplex(Math.PI / 2, x);
		} else if(x == Double.NEGATIVE_INFINITY) {
			if(y >= 0) {
				x = -x;
			}
			return LispComplex.newComplex(-Math.PI / 2, x);
		} else if(y == Double.POSITIVE_INFINITY) {
			return LispComplex.newComplex(0, y);
		} else if(y == Double.NEGATIVE_INFINITY) {
			return LispComplex.newComplex(0, y);
		} else {
			LispNumber z1 = LispComplex.I.mul(n);
			LispNumber z2 = sqrt(LispDouble.ONE.sub(n.mul(n)));
			LispNumber z3 = log(z1.add(z2));
			
			//return LispComplex.MINUS_I.mul(z3);
			return LispComplex.newComplex(
					z3.getImagDouble(), -z3.getRealDouble());
		}
	}

	public static LispNumber acos(LispNumber n) {
		double x = n.getRealDouble();
		
		if(n.isReal() && x >= -1.0 && x <= 1.0) {
			return new LispDouble(Math.acos(n.getRealDouble()));
		} else {
			LispNumber pi2 = new LispDouble(Math.PI / 2);
			
			return pi2.sub(asin(n));
		}
	}

	public static LispNumber atan(LispNumber n) {
		if(n.isReal()) {
			return new LispDouble(Math.atan(n.getRealDouble()));
		} else if(Double.isInfinite(n.getImagDouble())) {
			double im = n.getImagDouble();
			
			if(Double.isInfinite(n.getRealDouble())) {
				return LispDouble.NaN;
			} else if(im > 0) {
				return new LispDouble(Math.PI / 2);
			} else {
				return new LispDouble(-Math.PI / 2);
			}
		} else {
			LispNumber z1 = LispDouble.ONE.add(LispComplex.I.mul(n));
			LispNumber z2 = LispDouble.ONE.sub(LispComplex.I.mul(n));
			LispNumber z3 = log(z1).sub(log(z2));
			
			return z3.div(TWO_I);
		}
	}
	
}
