多倍長ライブラリ(更新中)

計算が遅かったり多少の誤差があるのも仕様です。

package jp.wshounen;

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;

public class BigMath {
	public static int precision = 20;
	public static MathContext context = new MathContext(24,
			RoundingMode.HALF_EVEN);
	private static BigDecimal ATAN0D1 = null;
	private static BigDecimal ATAN1 = null;
	private static BigDecimal LOG1D1 = null;
	private static BigDecimal LOG2 = null;
	public static BigDecimal PI = pi();
	public static BigDecimal E = e();

	public static BigDecimal round(final BigDecimal in, final int precision) {
		return in.round(new MathContext(precision, RoundingMode.HALF_EVEN));
	}

	public static BigDecimal round(final BigDecimal in) {
		return in.setScale(0, RoundingMode.HALF_EVEN);
	}

	public static BigDecimal degToRad(final BigDecimal arcsec) {
		return arcsec.multiply(ATAN1).divide(new BigDecimal(45), context);
	}

	public static BigDecimal radToDeg(final BigDecimal rad) {
		return rad.multiply(new BigDecimal(45)).divide(ATAN1, context);
	}

	public static BigDecimal floor(final BigDecimal in) {
		return in.setScale(0, RoundingMode.FLOOR);
	}

	public static BigDecimal ceil(final BigDecimal in) {
		return in.setScale(0, RoundingMode.CEILING);
	}

	private static BigDecimal pi() {
		ATAN0D1 = atan0d1();
		ATAN1 = atan(BigDecimal.ONE);
		return shift(ATAN1, 2);
	}

	private static BigDecimal e() {
		LOG1D1 = log1d1();
		LOG2 = log(new BigDecimal(2));
		System.out.println(BigMath.LOG1D1);
		System.out.println(BigMath.LOG2);
		return exp(BigDecimal.ONE);
	}

	public static void setPrecision(final int prec) {
		if (prec >= 20 && prec != precision) {
			precision = prec;
			setContext((int) Math.ceil(Math.log10(prec) * 3) + prec,
					RoundingMode.HALF_EVEN);
		}
	}

	public static void setContext(final int prec, final RoundingMode mode) {
		if (prec > precision && prec != context.getPrecision()) {
			context = new MathContext(prec, mode);
			if (ATAN0D1.precision() < prec) {
				LOG1D1 = null;
				LOG2 = null;
				ATAN0D1 = null;
				ATAN1 = null;
				E = e();
				PI = pi();
			}
		}
	}

	public static BigDecimal shift(final BigDecimal in, final int index) {
		if (index == 0) {
			return in;
		}
		BigInteger bi1 = in.unscaledValue();
		int num = in.scale();
		if (index > 0) {
			return new BigDecimal(bi1.shiftLeft(index), num);
		}
		BigInteger bi2 = bi1.shiftLeft(index);
		bi1 = bi2.shiftRight(index).negate().add(bi1);
		return (bi1.signum() == 0) ? new BigDecimal(bi2, num) : new BigDecimal(
				BigInteger.valueOf(5).pow(-index).multiply(bi1), num - index)
				.stripTrailingZeros().add(new BigDecimal(bi2, num));
	}

	public static BigDecimal shift(final BigDecimal in, final int index,
			final MathContext context) {
		return shift(in, index).plus(context);
	}

	public static BigDecimal asin(final BigDecimal in) {
		int num = in.abs().compareTo(BigDecimal.ONE);
		if (num == 0) {
			return (in.signum() > 0) ? shift(ATAN1, 1) : shift(ATAN1, 1)
					.negate();
		} else if (num == -1) {
			return (in.signum() == 0) ? BigDecimal.ZERO : shift(atan(in.divide(
					sqrt(in.pow(2).negate().add(BigDecimal.ONE)).add(
							BigDecimal.ONE), context)), 1);
		} else if (round(in.abs(), precision).compareTo(BigDecimal.ONE) == 0) {
			return (in.signum() > 0) ? shift(ATAN1, 1) : shift(ATAN1, 1)
					.negate();
		} else {
			return null;
		}
	}

	public static BigDecimal acos(final BigDecimal in) {
		int num = in.abs().compareTo(BigDecimal.ONE);
		if (num == 0) {
			return (in.signum() > 0) ? BigDecimal.ZERO : shift(ATAN1, 2);
		} else if (num == -1) {
			return (in.signum() == 0) ? shift(ATAN1, 1)
					: shift(atan(
							in.divide(
									sqrt(in.pow(2).negate().add(BigDecimal.ONE))
											.add(BigDecimal.ONE), context))
							.negate().add(ATAN1), 1);
		} else if (round(in.abs(), precision).compareTo(BigDecimal.ONE) == 0) {
			return (in.signum() > 0) ? BigDecimal.ZERO : shift(ATAN1, 2);
		} else {
			return null;
		}
	}

	private static BigDecimal atan0d1() {
		if (ATAN0D1 != null) {
			return ATAN0D1.plus();
		}
		BigDecimal bd = BigDecimal.ONE;
		BigDecimal[] temp = new BigDecimal[] { BigDecimal.ONE
				.scaleByPowerOfTen(-1) };
		BigDecimal coef = temp[0].plus();
		do {
			bd = new BigDecimal(2).add(bd);
			coef = coef.scaleByPowerOfTen(-2).negate();
			temp = new BigDecimal[] {
					temp[0].multiply(bd).add(coef).divide(bd, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal atan(final BigDecimal in) {
		BigDecimal perten = BigDecimal.ONE.scaleByPowerOfTen(-1);
		BigDecimal bd1, bd2, coef;
		BigDecimal[] temp;
		int num = in.signum();
		bd1 = in.abs(context);
		bd2 = BigDecimal.ZERO;
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (bd1.compareTo(perten) == 0) {
			return (num > 0) ? ATAN0D1.plus() : ATAN0D1.negate();
		} else if (bd1.compareTo(BigDecimal.ONE) == 0 && num == 1
				&& ATAN1 != null) {
			return ATAN1.plus();
		} else if (bd1.compareTo(BigDecimal.ONE) == 1) {
			bd1 = new BigDecimal(-1).add(bd1).divide(bd1.add(BigDecimal.ONE),
					context);
			bd2 = bd2.add(ATAN1);
		}
		perten = perten.negate();
		while (bd1.signum() == 1) {
			bd1 = perten.add(bd1).divide(
					bd1.scaleByPowerOfTen(-1).add(BigDecimal.ONE), context);
			bd2 = bd2.add(ATAN0D1);
		}
		if (num == -1) {
			bd1 = bd1.negate();
			bd2 = bd2.negate();
		}
		coef = bd1.plus();
		temp = new BigDecimal[] { coef.add(bd2, context) };
		bd2 = BigDecimal.ONE;
		bd1 = bd1.pow(2, context).negate();
		do {
			bd2 = new BigDecimal(2).add(bd2);
			coef = coef.multiply(bd1, context);
			temp = new BigDecimal[] {
					temp[0].multiply(bd2).add(coef).divide(bd2, context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal root(final BigDecimal in, final int index) {
		BigDecimal bd = in.round(context);
		int num = bd.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1 || index <= 0) {
			return null;
		} else if (bd.compareTo(BigDecimal.ONE) == 0 || index == 1) {
			return bd;
		}
		int i = index;
		while (i % 2 == 0) {
			i /= 2;
			bd = sqrt(bd);
		}
		if (i > 2) {
			bd = pow(bd, BigDecimal.ONE.divide(new BigDecimal(i), context));
		}
		return bd;
	}

	public static BigDecimal sqrt(final BigDecimal in) {
		BigDecimal bd = in.round(context);
		int num = bd.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			return null;
		} else if (bd.compareTo(BigDecimal.ONE) == 0) {
			return bd;
		}
		BigDecimal[] temp = new BigDecimal[] { bd };
		do {
			temp = new BigDecimal[] {
					temp[0].pow(2).add(bd).divide(shift(temp[0], 1), context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	private static BigDecimal log1d1() {
		if (LOG1D1 != null) {
			return LOG1D1.plus();
		}
		BigDecimal bd = BigDecimal.ONE;
		BigDecimal[] temp = new BigDecimal[] { BigDecimal.ONE
				.scaleByPowerOfTen(-1) };
		BigDecimal coef = temp[0].plus();
		do {
			bd = BigDecimal.ONE.add(bd);
			coef = coef.negate().scaleByPowerOfTen(-1);
			temp = new BigDecimal[] {
					temp[0].multiply(bd).add(coef).divide(bd, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal log(BigDecimal in) {
		if (in.signum() < 1) {
			return null;
		}
		BigDecimal nperten = new BigDecimal(11).scaleByPowerOfTen(-1);
		BigDecimal bd1 = in.round(context);
		int num = bd1.compareTo(BigDecimal.ONE);
		BigDecimal coef, bd2;
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			coef = BigDecimal.ONE;
		} else {
			bd1 = BigDecimal.ONE.divide(bd1, context);
			coef = new BigDecimal(-1);
		}
		bd2 = BigDecimal.ZERO;
		nperten = new BigDecimal(5).scaleByPowerOfTen(-1);
		while (bd1.compareTo(nperten) == -1) {
			bd1 = shift(bd1, 1);
			bd2 = bd2.add(LOG2);
		}
		if (bd1.plus(context).compareTo(nperten) == 0 && LOG2 != null) {
			return (num > 0) ? bd2.add(LOG2) : bd2.add(LOG2).negate();
		}
		nperten = new BigDecimal(11).scaleByPowerOfTen(-1);
		while (bd1.compareTo(BigDecimal.ONE) == -1) {
			bd1 = bd1.multiply(nperten);
			bd2 = bd2.add(LOG1D1);
		}
		if (num == -1) {
			bd2 = bd2.negate();
		}
		if (bd1.plus(context).compareTo(BigDecimal.ONE) == 0) {
			return bd2;
		}
		bd1 = sqrt(sqrt(bd1));
		coef = shift(coef, 2);
		bd1 = BigDecimal.ONE.add(new BigDecimal(-1).divide(bd1, context));
		coef = coef.multiply(bd1, context);
		BigDecimal[] temp = new BigDecimal[] { coef.add(bd2, context) };
		bd2 = BigDecimal.ONE;
		do {
			bd2 = BigDecimal.ONE.add(bd2);
			coef = coef.multiply(bd1, context);
			temp = new BigDecimal[] {
					temp[0].multiply(bd2).add(coef).divide(bd2, context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal sin(final BigDecimal in) {
		BigDecimal bd1, bd2;
		bd2 = shift(ATAN1, 2);
		bd1 = in.remainder(shift(bd2, 1));
		int num = bd1.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num < 0) {
			return sin(bd1.negate()).negate();
		}
		if (bd1.compareTo(bd2) > -1) {
			return sin(bd1.add(bd2)).negate();
		} else if (bd1.compareTo(ATAN1) > 0) {
			if (bd1.compareTo(shift(ATAN1, 1)) > 0) {
				return sin(bd1.negate().add(bd2));
			} else {
				return cos(bd1.negate().add(shift(ATAN1, 1)));
			}
		}
		bd1 = shift(bd1, 1).scaleByPowerOfTen(-1);
		BigDecimal[] temp = new BigDecimal[] { bd1.plus(context) };
		BigDecimal coef = temp[0].plus();
		bd1 = bd1.pow(2, context);
		bd2 = BigDecimal.ZERO;
		do {
			bd2 = new BigDecimal(2).add(bd2);
			coef = coef.multiply(bd1).divide(
					bd2.add(BigDecimal.ONE).multiply(bd2).negate(), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		bd1 = shift(temp[0], 1).pow(2, context);
		return bd1.add(new BigDecimal(-5)).multiply(bd1).add(new BigDecimal(5))
				.multiply(temp[0], context);
	}

	public static BigDecimal cos(final BigDecimal in) {
		BigDecimal bd1, bd2;
		bd2 = shift(ATAN1, 2);
		bd1 = in.remainder(shift(bd2, 1));
		int num = bd1.signum();
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num < 0) {
			return cos(bd1.negate());
		}
		if (bd1.compareTo(bd2) > -1) {
			return cos(bd1.add(bd2)).negate();
		} else if (bd1.compareTo(ATAN1) > -1) {
			if (bd1.compareTo(shift(ATAN1, 1)) > 0) {
				return cos(bd1.negate().add(bd2)).negate();
			} else {
				return sin(bd1.negate().add(shift(ATAN1, 1)));
			}
		}
		bd1 = shift(bd1, 1).pow(2, context).scaleByPowerOfTen(-2);
		BigDecimal[] temp = new BigDecimal[] { BigDecimal.ONE };
		BigDecimal coef = temp[0].plus();
		bd2 = BigDecimal.ZERO;
		do {
			bd2 = new BigDecimal(2).add(bd2);
			coef = coef.multiply(bd1).divide(
					bd2.negate().add(BigDecimal.ONE).multiply(bd2), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		bd1 = shift(temp[0], 1).pow(2, context);
		return bd1.add(new BigDecimal(-5)).multiply(bd1).add(new BigDecimal(5))
				.multiply(temp[0], context);
	}

	public static BigDecimal tan(final BigDecimal in) {
		BigDecimal bd = cos(in);
		return (bd.signum() == 0) ? null : sin(in).divide(bd, context);
	}

	public static BigDecimal exp(final BigDecimal in) {
		int num = in.signum();
		BigDecimal coef = BigDecimal.ONE;
		BigDecimal bi;
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num == 1) {
			bi = in.negate(context);
			if (bi.compareTo(LOG1D1.negate()) == 0) {
				return new BigDecimal(11).scaleByPowerOfTen(-1);
			} else if (bi.compareTo(LOG2.negate()) == 0) {
				return new BigDecimal(2);
			}
		} else {
			bi = in.plus(context);
			if (bi.compareTo(LOG1D1.negate()) == 0) {
				return BigDecimal.ONE.divide(new BigDecimal(11), context)
						.scaleByPowerOfTen(1);
			} else if (bi.compareTo(LOG2.negate()) == 0) {
				return new BigDecimal(5).scaleByPowerOfTen(-1);
			}
		}
		if (bi.compareTo(new BigDecimal(-1)) != 0) {
			while (bi.signum() == -1) {
				bi = bi.add(BigDecimal.ONE);
				coef = coef.multiply(E, context);
			}
		}

		if (num == 1) {
			bi = bi.negate();
		} else {
			coef = BigDecimal.ONE.divide(coef, context);
		}
		BigDecimal[] temp = new BigDecimal[] { coef.round(context) };
		num = 0;
		do {
			num++;
			coef = coef.multiply(bi).divide(new BigDecimal(num), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal pow(final BigDecimal in, final BigDecimal index) {
		BigDecimal bi = log(in);
		if (bi == null) {
			return null;
		} else {
			return exp(bi.multiply(index));
		}
	}

	/**
	 * Matrix of rotation Z to Y. R1(rad).
	 */
	public static BigDecimal[][] rotX(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] {
				{ BigDecimal.ONE, BigDecimal.ZERO, BigDecimal.ZERO },
				{ BigDecimal.ZERO, Cos, Sin.negate() },
				{ BigDecimal.ZERO, Sin, Cos } };
	}

	/**
	 * Matrix of rotation X to Z. R2(rad).
	 */
	public static BigDecimal[][] rotY(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, BigDecimal.ZERO, Sin },
				{ BigDecimal.ZERO, BigDecimal.ONE, BigDecimal.ZERO },
				{ Sin.negate(), BigDecimal.ZERO, Cos } };
	}

	/**
	 * Matrix of rotation Y to X. R3(rad).
	 */
	public static BigDecimal[][] rotZ(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin.negate(), BigDecimal.ZERO },
				{ Sin, Cos, BigDecimal.ZERO },
				{ BigDecimal.ZERO, BigDecimal.ZERO, BigDecimal.ONE } };
	}

	/**
	 * Matrix of rotation Y to X. R(rad).
	 */
	public static BigDecimal[][] rot2(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin.negate() }, { Sin, Cos } };
	}

	/**
	 * Rotation (1,0,0) to A[0] = (x,y,z). Rotation (0,1,0) to A[1] = (x,y,z).
	 * Rotation (0,0,1) to A[2] = (x,y,z). Ar = rotA(A,r).
	 */
	public static BigDecimal[] rotA(final BigDecimal[][] A,
			final BigDecimal[] star) {
		if (star.length != A.length) {
			return null;
		} else {
			for (int i = 1; i < A.length; i++) {
				if (A[i].length != A[0].length) {
					return null;
				}
			}
		}
		BigDecimal[] res = new BigDecimal[A[0].length];
		for (int i = 0; i < A[0].length; i++) {
			res[i] = BigDecimal.ZERO;
			for (int j = 0; j < star.length; j++) {
				res[i] = res[i].add(star[j].multiply(A[j][i]));
			}
		}
		return res;
	}

	/**
	 * AB = (rotA(A,B[0]),rotA(A,B[1]),rotA(A,B[2])).
	 */
	public static BigDecimal[][] rotAB(final BigDecimal[][] A,
			final BigDecimal[][] B) {
		BigDecimal[][] res = new BigDecimal[B.length][];
		for (int i = 0; i < B.length; i++) {
			res[i] = rotA(A, B[i]);
		}
		return res;
	}
}

JAVAを用いた多倍長ライブラリ(第2版)

多倍長ライブラリの一部を更新しました。以下の資料をご覧下さい。実際には5乗根の計算と0.2乗の計算は似たような速度でしたが計算桁数の関係で平方根以外は逆数乗で計算しています。

多倍長ライブラリの作り方と高次方程式の解の計算方法

通常の多倍長ライブラリです。回転行列を軸の回転に統一しました。

package jp.wshounen;

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;

public class BigMath {
	public static int precision = 20;
	public static MathContext context = new MathContext(24,
			RoundingMode.HALF_EVEN);
	public static BigDecimal PI = atan(BigDecimal.ONE).multiply(
			new BigDecimal(4));
	public static BigDecimal E = exp(BigDecimal.ONE);

	public static BigDecimal round(final BigDecimal in, final int precision) {
		return in.round(new MathContext(precision, RoundingMode.HALF_EVEN));
	}

	public static BigDecimal round(final BigDecimal in) {
		return in.setScale(0, RoundingMode.HALF_EVEN);
	}

	public static BigDecimal degToRad(final BigDecimal arcsec) {
		return arcsec.multiply(PI).divide(new BigDecimal(180), context);
	}

	public static BigDecimal radToDeg(final BigDecimal rad) {
		return rad.multiply(new BigDecimal(180)).divide(PI, context);
	}

	public static BigDecimal floor(final BigDecimal in) {
		return in.setScale(0, RoundingMode.FLOOR);
	}

	public static BigDecimal ceil(final BigDecimal in) {
		return in.setScale(0, RoundingMode.CEILING);
	}

	public static void setPrecision(final int prec) {
		if (prec >= 20 && prec != precision) {
			precision = prec;
			setContext((int) Math.ceil(Math.log10(prec) * 3) + prec,
					RoundingMode.HALF_EVEN);
		}
	}

	public static void setContext(final int prec, final RoundingMode mode) {
		if (prec >= precision && prec != context.getPrecision()) {
			context = new MathContext(prec, mode);
			if (PI.precision() < prec) {
				E = exp(BigDecimal.ONE);
				PI = atan(BigDecimal.ONE).multiply(new BigDecimal(4));
			}
		}
	}

	public static BigDecimal asin(final BigDecimal in) {
		int num = in.abs().compareTo(BigDecimal.ONE);
		if (num == 0) {
			return PI.divide(new BigDecimal(in.signum() * 2));
		} else if (num == -1) {
			BigDecimal bi = atan(in.divide(
					root(in.pow(2).negate().add(BigDecimal.ONE), 2).add(
							BigDecimal.ONE), context));
			return bi.add(bi, context);
		} else if (round(in.abs(), precision).compareTo(BigDecimal.ONE) == 0) {
			return PI.divide(new BigDecimal(in.signum() * 2));
		} else {
			return null;
		}
	}

	public static BigDecimal acos(final BigDecimal in) {
		if (in.round(context).compareTo(BigDecimal.ONE) == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bi = asin(in);
		if (bi == null) {
			return bi;
		} else {
			return PI.add(bi.add(bi).negate()).divide(new BigDecimal(2));
		}
	}

	public static BigDecimal atan(final BigDecimal in) {
		int num = in.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bd1 = in.abs(context);
		BigDecimal bd2 = BigDecimal.ZERO;
		BigDecimal perfive = new BigDecimal(2).scaleByPowerOfTen(-1);
		if (bd1.compareTo(perfive) == 1) {
			BigDecimal perfiveat = atan(perfive);
			while (bd1.signum() == 1) {
				bd1 = perfive
						.negate()
						.add(bd1)
						.divide(bd1.multiply(perfive).add(BigDecimal.ONE),
								context);
				bd2 = bd2.add(perfiveat);
			}
		}
		if (num == -1) {
			bd1 = bd1.negate();
			bd2 = bd2.negate();
		}
		BigDecimal coef = bd1.plus();
		BigDecimal[] temp = new BigDecimal[] { coef.add(bd2, context) };
		num = 1;
		bd1 = bd1.pow(2, context).negate();
		do {
			num += 2;
			bd2 = new BigDecimal(num);
			coef = coef.multiply(bd1, context);
			temp = new BigDecimal[] {
					temp[0].multiply(bd2).add(coef).divide(bd2, context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal root(final BigDecimal in, final int index) {
		BigDecimal bd = in.round(context);
		int num = bd.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1 || index <= 0) {
			return null;
		} else if (bd.compareTo(BigDecimal.ONE) == 0 || index == 1) {
			return bd;
		}
		int i = index;
		while (i % 2 == 0) {
			i /= 2;
			bd = sqrt(bd);
		}
		if (i > 2) {
			bd = pow(bd, BigDecimal.ONE.divide(new BigDecimal(i), context));
		}
		return bd;
	}

	public static BigDecimal sqrt(final BigDecimal in) {
		BigDecimal bd = in.round(context);
		int num = bd.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			return null;
		} else if (bd.compareTo(BigDecimal.ONE) == 0) {
			return bd;
		}
		BigDecimal[] temp = new BigDecimal[] { bd };
		do {
			temp = new BigDecimal[] {
					temp[0].multiply(temp[0]).add(bd)
							.divide(temp[0].add(temp[0]), context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal log(BigDecimal in) {
		if (in.signum() < 1) {
			return null;
		}
		BigDecimal bd1 = in.round(context);
		int num = bd1.compareTo(BigDecimal.ONE);
		BigDecimal coef, bd2;
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			coef = BigDecimal.ONE;
		} else {
			bd1 = BigDecimal.ONE.divide(bd1, context);
			coef = new BigDecimal(-1);
		}
		bd2 = BigDecimal.ZERO;
		do {
			bd1 = bd1.multiply(E, context);
			bd2 = bd2.add(BigDecimal.ONE);
		} while (bd1.compareTo(BigDecimal.ONE) == -1);
		if (num == -1) {
			bd2 = bd2.negate();
		}
		while (bd1.pow(16).compareTo(E) == 1) {
			bd1 = sqrt(bd1);
			coef = coef.add(coef);
		}
		bd1 = BigDecimal.ONE.add(new BigDecimal(-1).divide(bd1, context));
		coef = coef.multiply(bd1, context);
		BigDecimal[] temp = new BigDecimal[] { coef.add(bd2, context) };
		num = 1;
		do {
			num++;
			coef = coef.multiply(bd1, context);
			bd2 = new BigDecimal(num);
			temp = new BigDecimal[] {
					temp[0].multiply(bd2).add(coef).divide(bd2, context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal sin(final BigDecimal in) {
		BigDecimal bd1, bd2;
		bd2 = PI.add(PI);
		bd1 = in.remainder(bd2);
		int num = bd1.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num < 0) {
			return sin(bd1.negate()).negate();
		}
		if (bd1.compareTo(PI) > -1) {
			return sin(bd1.add(PI)).negate();
		} else if (bd1.compareTo(PI.divide(new BigDecimal(4))) > 0) {
			bd2 = PI.divide(new BigDecimal(2));
			if (bd1.compareTo(bd2) > 0) {
				return sin(bd1.negate().add(PI));
			} else {
				return cos(bd1.negate().add(bd2));
			}
		}
		bd1 = bd1.divide(new BigDecimal(3), context);
		BigDecimal[] temp = new BigDecimal[] { bd1.plus() };
		BigDecimal coef = temp[0].plus();
		bd1 = bd1.pow(2, context);
		num = 0;
		do {
			num++;
			coef = coef.multiply(bd1).divide(
					new BigDecimal(-2 * num * (2 * num + 1)), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0].multiply(temp[0]).multiply(new BigDecimal(-4))
				.add(new BigDecimal(3)).multiply(temp[0], context);
	}

	public static BigDecimal cos(final BigDecimal in) {
		BigDecimal bd1, bd2;
		bd2 = PI.add(PI);
		bd1 = in.remainder(bd2);
		int num = bd1.signum();
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num < 0) {
			return cos(bd1.negate());
		}
		if (bd1.compareTo(PI) > -1) {
			return cos(bd1.add(PI)).negate();
		} else if (bd1.compareTo(PI.divide(new BigDecimal(4))) > 0) {
			bd2 = PI.divide(new BigDecimal(2));
			if (bd1.compareTo(bd2) > 0) {
				return cos(bd1.negate().add(PI)).negate();
			} else {
				return sin(bd1.negate().add(bd2));
			}
		}
		bd1 = bd1.pow(2).divide(new BigDecimal(9), context);
		BigDecimal[] temp = new BigDecimal[] { BigDecimal.ONE };
		BigDecimal coef = temp[0].plus();
		num = 0;
		do {
			num++;
			coef = coef.multiply(bd1).divide(
					new BigDecimal(-2 * num * (2 * num - 1)), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0].multiply(temp[0]).multiply(new BigDecimal(4))
				.add(new BigDecimal(-3)).multiply(temp[0], context);
	}

	public static BigDecimal tan(final BigDecimal in) {
		BigDecimal bd = cos(in);
		if (bd.add(BigDecimal.ONE, context).compareTo(BigDecimal.ONE) == 0) {
			return null;
		} else {
			return sin(in).divide(bd, context);
		}
	}

	public static BigDecimal exp(final BigDecimal in) {
		int num = in.signum();
		BigDecimal coef = BigDecimal.ONE;
		BigDecimal bi;
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num == 1) {
			bi = in.negate(context);
		} else {
			bi = in.round(context);
		}
		if (bi.compareTo(new BigDecimal(-1)) != 0) {
			while (bi.signum() == -1) {
				bi = bi.add(BigDecimal.ONE);
				coef = coef.multiply(E, context);
			}
		}

		if (num == 1) {
			bi = bi.negate();
		} else {
			coef = BigDecimal.ONE.divide(coef, context);
		}
		BigDecimal[] temp = new BigDecimal[] { coef.round(context) };
		num = 0;
		do {
			num++;
			coef = coef.multiply(bi).divide(new BigDecimal(num), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal pow(final BigDecimal in, final BigDecimal index) {
		BigDecimal bi = log(in);
		if (bi == null) {
			return null;
		} else {
			return exp(bi.multiply(index));
		}
	}

	/**
	 * Matrix of rotation Z to Y. R1(rad).
	 */
	public static BigDecimal[][] rotX(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] {
				{ BigDecimal.ONE, BigDecimal.ZERO, BigDecimal.ZERO },
				{ BigDecimal.ZERO, Cos, Sin.negate() },
				{ BigDecimal.ZERO, Sin, Cos } };
	}

	/**
	 * Matrix of rotation X to Z. R2(rad).
	 */
	public static BigDecimal[][] rotY(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, BigDecimal.ZERO, Sin },
				{ BigDecimal.ZERO, BigDecimal.ONE, BigDecimal.ZERO },
				{ Sin.negate(), BigDecimal.ZERO, Cos } };
	}

	/**
	 * Matrix of rotation Y to X. R3(rad).
	 */
	public static BigDecimal[][] rotZ(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin.negate(), BigDecimal.ZERO },
				{ Sin, Cos, BigDecimal.ZERO },
				{ BigDecimal.ZERO, BigDecimal.ZERO, BigDecimal.ONE } };
	}

	/**
	 * Matrix of rotation Y to X. R(rad).
	 */
	public static BigDecimal[][] rot2(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin.negate() }, { Sin, Cos } };
	}

	/**
	 * Rotation (1,0,0) to A[0] = (x,y,z). Rotation (0,1,0) to A[1] = (x,y,z).
	 * Rotation (0,0,1) to A[2] = (x,y,z). Ar = rotA(A,r).
	 */
	public static BigDecimal[] rotA(final BigDecimal[][] A,
			final BigDecimal[] star) {
		if (star.length != A.length) {
			return null;
		} else {
			for (int i = 1; i < A.length; i++) {
				if (A[i].length != A[0].length) {
					return null;
				}
			}
		}
		BigDecimal[] res = new BigDecimal[A[0].length];
		for (int i = 0; i < A[0].length; i++) {
			res[i] = BigDecimal.ZERO;
			for (int j = 0; j < star.length; j++) {
				res[i] = res[i].add(star[j].multiply(A[j][i]));
			}
		}
		return res;
	}

	/**
	 * AB = (rotA(A,B[0]),rotA(A,B[1]),rotA(A,B[2])).
	 */
	public static BigDecimal[][] rotAB(final BigDecimal[][] A,
			final BigDecimal[][] B) {
		BigDecimal[][] res = new BigDecimal[B.length][];
		for (int i = 0; i < B.length; i++) {
			res[i] = rotA(A, B[i]);
		}
		return res;
	}
}

暦計算用の定数付き20桁精度のライブラリです。

package jp.wshounen;

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;

public class DEMath {
	public static final int precision = 20;
	public final static MathContext context = new MathContext(24,
			RoundingMode.HALF_EVEN);
	public static final BigDecimal au = new BigDecimal(149597870691l)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal emrat = new BigDecimal(8130056)
			.scaleByPowerOfTen(-5);
	public static final BigDecimal clight = new BigDecimal(2997922458l)
			.scaleByPowerOfTen(-4);
	public static final BigDecimal jultimeOfJ2000 = new BigDecimal(2451545);
	public static final BigDecimal equatorialRadius = new BigDecimal(6378137)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal polarRadius = new BigDecimal(6356752)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal solarRadius = new BigDecimal(696000);
	public static final BigDecimal lunarRadius = new BigDecimal(1738091)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal PI = atan(BigDecimal.ONE).multiply(
			new BigDecimal(4));
	public static final BigDecimal E = exp(BigDecimal.ONE);
	public static final BigDecimal[][] icrfRot = rotAB(
			rotAB(rotX(arcsecToRad(new BigDecimal(-68192).scaleByPowerOfTen(-7))),
					rotY(arcsecToRad(new BigDecimal(-166170)
							.scaleByPowerOfTen(-7)))),
			rotZ(arcsecToRad(new BigDecimal(146).scaleByPowerOfTen(-4))));
	public static final BigDecimal obliquityArcsecOfJ2000 = new BigDecimal(
			84381406).scaleByPowerOfTen(-3);

	public static BigDecimal round(final BigDecimal in, final int precision) {
		return in.round(new MathContext(precision, RoundingMode.HALF_EVEN));
	}

	public static long divideToFloor(final long value, final long divisor) {
		long remainder = value % divisor;
		return ((remainder < 0 && divisor > 0) || (remainder > 0 && divisor < 0)) ? (value - remainder)
				/ divisor - 1
				: (value - remainder) / divisor;
	}

	public static BigDecimal round(final BigDecimal in) {
		return in.setScale(0, RoundingMode.HALF_EVEN);
	}

	public static BigDecimal arcsecToRad(final BigDecimal arcsec) {
		return arcsec.multiply(PI).divide(new BigDecimal(648000), context);
	}

	public static BigDecimal radToArcsec(final BigDecimal rad) {
		return rad.multiply(new BigDecimal(648000)).divide(PI, context);
	}

	public static BigDecimal floor(final BigDecimal in) {
		return in.setScale(0, RoundingMode.FLOOR);
	}

	public static BigDecimal ceil(final BigDecimal in) {
		return in.setScale(0, RoundingMode.CEILING);
	}

	public static BigDecimal asin(final BigDecimal in) {
		int num = in.abs().compareTo(BigDecimal.ONE);
		if (num == 0) {
			return PI.divide(new BigDecimal(in.signum() * 2));
		} else if (num == -1) {
			BigDecimal bi = atan(in.divide(
					root(in.pow(2).negate().add(BigDecimal.ONE), 2).add(
							BigDecimal.ONE), context));
			return bi.add(bi, context);
		} else if (round(in.abs(), precision).compareTo(BigDecimal.ONE) == 0) {
			return PI.divide(new BigDecimal(in.signum() * 2));
		} else {
			return null;
		}
	}

	public static BigDecimal acos(final BigDecimal in) {
		if (in.round(context).compareTo(BigDecimal.ONE) == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bi = asin(in);
		if (bi == null) {
			return bi;
		} else {
			return PI.add(bi.add(bi).negate()).divide(new BigDecimal(2));
		}
	}

	public static BigDecimal atan(final BigDecimal in) {
		int num = in.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bd1 = in.abs(context);
		BigDecimal bd2 = BigDecimal.ZERO;
		BigDecimal perfive = new BigDecimal(2).scaleByPowerOfTen(-1);
		if (bd1.compareTo(perfive) == 1) {
			BigDecimal perfiveat = atan(perfive);
			while (bd1.signum() == 1) {
				bd1 = perfive
						.negate()
						.add(bd1)
						.divide(bd1.multiply(perfive).add(BigDecimal.ONE),
								context);
				bd2 = bd2.add(perfiveat);
			}
		}
		if (num == -1) {
			bd1 = bd1.negate();
			bd2 = bd2.negate();
		}
		BigDecimal coef = bd1.plus();
		BigDecimal[] temp = new BigDecimal[] { coef.add(bd2, context) };
		num = 1;
		bd1 = bd1.pow(2, context).negate();
		do {
			num += 2;
			bd2 = new BigDecimal(num);
			coef = coef.multiply(bd1, context);
			temp = new BigDecimal[] {
					temp[0].multiply(bd2).add(coef).divide(bd2, context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal root(final BigDecimal in, final int index) {
		BigDecimal bd = in.round(context);
		int num = bd.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1 || index <= 0) {
			return null;
		} else if (bd.compareTo(BigDecimal.ONE) == 0 || index == 1) {
			return bd;
		}
		int i = index;
		while (i % 2 == 0) {
			i /= 2;
			bd = sqrt(bd);
		}
		if (i > 2) {
			bd = pow(bd, BigDecimal.ONE.divide(new BigDecimal(i), context));
		}
		return bd;
	}

	public static BigDecimal sqrt(final BigDecimal in) {
		BigDecimal bd = in.round(context);
		int num = bd.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			return null;
		} else if (bd.compareTo(BigDecimal.ONE) == 0) {
			return bd;
		}
		BigDecimal[] temp = new BigDecimal[] { bd };
		do {
			temp = new BigDecimal[] {
					temp[0].multiply(temp[0]).add(bd)
							.divide(temp[0].add(temp[0]), context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal log(BigDecimal in) {
		if (in.signum() < 1) {
			return null;
		}
		BigDecimal bd1 = in.round(context);
		int num = bd1.compareTo(BigDecimal.ONE);
		BigDecimal coef, bd2;
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			coef = BigDecimal.ONE;
		} else {
			bd1 = BigDecimal.ONE.divide(bd1, context);
			coef = new BigDecimal(-1);
		}
		bd2 = BigDecimal.ZERO;
		do {
			bd1 = bd1.multiply(E, context);
			bd2 = bd2.add(BigDecimal.ONE);
		} while (bd1.compareTo(BigDecimal.ONE) == -1);
		if (num == -1) {
			bd2 = bd2.negate();
		}
		while (bd1.pow(16).compareTo(E) == 1) {
			bd1 = sqrt(bd1);
			coef = coef.add(coef);
		}
		bd1 = BigDecimal.ONE.add(new BigDecimal(-1).divide(bd1, context));
		coef = coef.multiply(bd1, context);
		BigDecimal[] temp = new BigDecimal[] { coef.add(bd2, context) };
		num = 1;
		do {
			num++;
			coef = coef.multiply(bd1, context);
			bd2 = new BigDecimal(num);
			temp = new BigDecimal[] {
					temp[0].multiply(bd2).add(coef).divide(bd2, context),
					temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal sin(final BigDecimal in) {
		BigDecimal bd1, bd2;
		bd2 = PI.add(PI);
		bd1 = in.remainder(bd2);
		int num = bd1.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num < 0) {
			return sin(bd1.negate()).negate();
		}
		if (bd1.compareTo(PI) > -1) {
			return sin(bd1.add(PI)).negate();
		} else if (bd1.compareTo(PI.divide(new BigDecimal(4))) > 0) {
			bd2 = PI.divide(new BigDecimal(2));
			if (bd1.compareTo(bd2) > 0) {
				return sin(bd1.negate().add(PI));
			} else {
				return cos(bd1.negate().add(bd2));
			}
		}
		bd1 = bd1.divide(new BigDecimal(3), context);
		BigDecimal[] temp = new BigDecimal[] { bd1.plus() };
		BigDecimal coef = temp[0].plus();
		bd1 = bd1.pow(2, context);
		num = 0;
		do {
			num++;
			coef = coef.multiply(bd1).divide(
					new BigDecimal(-2 * num * (2 * num + 1)), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0].multiply(temp[0]).multiply(new BigDecimal(-4))
				.add(new BigDecimal(3)).multiply(temp[0], context);
	}

	public static BigDecimal cos(final BigDecimal in) {
		BigDecimal bd1, bd2;
		bd2 = PI.add(PI);
		bd1 = in.remainder(bd2);
		int num = bd1.signum();
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num < 0) {
			return cos(bd1.negate());
		}
		if (bd1.compareTo(PI) > -1) {
			return cos(bd1.add(PI)).negate();
		} else if (bd1.compareTo(PI.divide(new BigDecimal(4))) > 0) {
			bd2 = PI.divide(new BigDecimal(2));
			if (bd1.compareTo(bd2) > 0) {
				return cos(bd1.negate().add(PI)).negate();
			} else {
				return sin(bd1.negate().add(bd2));
			}
		}
		bd1 = bd1.pow(2).divide(new BigDecimal(9), context);
		BigDecimal[] temp = new BigDecimal[] { BigDecimal.ONE };
		BigDecimal coef = temp[0].plus();
		num = 0;
		do {
			num++;
			coef = coef.multiply(bd1).divide(
					new BigDecimal(-2 * num * (2 * num - 1)), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0].multiply(temp[0]).multiply(new BigDecimal(4))
				.add(new BigDecimal(-3)).multiply(temp[0], context);
	}

	public static BigDecimal tan(final BigDecimal in) {
		BigDecimal bd = cos(in);
		if (bd.add(BigDecimal.ONE, context).compareTo(BigDecimal.ONE) == 0) {
			return null;
		} else {
			return sin(in).divide(bd, context);
		}
	}

	public static BigDecimal exp(final BigDecimal in) {
		int num = in.signum();
		BigDecimal coef = BigDecimal.ONE;
		BigDecimal bi;
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num == 1) {
			bi = in.negate(context);
		} else {
			bi = in.round(context);
		}
		if (bi.compareTo(new BigDecimal(-1)) != 0) {
			while (bi.signum() == -1) {
				bi = bi.add(BigDecimal.ONE);
				coef = coef.multiply(E, context);
			}
		}

		if (num == 1) {
			bi = bi.negate();
		} else {
			coef = BigDecimal.ONE.divide(coef, context);
		}
		BigDecimal[] temp = new BigDecimal[] { coef.round(context) };
		num = 0;
		do {
			num++;
			coef = coef.multiply(bi).divide(new BigDecimal(num), context);
			temp = new BigDecimal[] { temp[0].add(coef, context), temp[0] };
		} while (temp[0].compareTo(temp[1]) != 0);
		return temp[0];
	}

	public static BigDecimal pow(final BigDecimal in, final BigDecimal index) {
		BigDecimal bi = log(in);
		if (bi == null) {
			return null;
		} else {
			return exp(bi.multiply(index));
		}
	}

	/**
	 * Matrix of rotation Y to Z. R1(-rad).
	 */
	public static BigDecimal[][] rotX(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] {
				{ BigDecimal.ONE, BigDecimal.ZERO, BigDecimal.ZERO },
				{ BigDecimal.ZERO, Cos, Sin },
				{ BigDecimal.ZERO, Sin.negate(), Cos } };
	}

	/**
	 * Matrix of rotation X to Z. R2(rad).
	 */
	public static BigDecimal[][] rotY(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, BigDecimal.ZERO, Sin },
				{ BigDecimal.ZERO, BigDecimal.ONE, BigDecimal.ZERO },
				{ Sin.negate(), BigDecimal.ZERO, Cos } };
	}

	/**
	 * Matrix of rotation X to Y. R3(-rad).
	 */
	public static BigDecimal[][] rotZ(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin, BigDecimal.ZERO },
				{ Sin.negate(), Cos, BigDecimal.ZERO },
				{ BigDecimal.ZERO, BigDecimal.ZERO, BigDecimal.ONE } };
	}

	/**
	 * Matrix of rotation X to Y. R(-rad).
	 */
	public static BigDecimal[][] rot2(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin }, { Sin.negate(), Cos } };
	}

	/**
	 * Rotation (1,0,0) to A[0] = (x,y,z). Rotation (0,1,0) to A[1] = (x,y,z).
	 * Rotation (0,0,1) to A[2] = (x,y,z). Ar = rotA(A,r).
	 */
	public static BigDecimal[] rotA(final BigDecimal[][] A,
			final BigDecimal[] star) {
		if (star.length != A.length) {
			return null;
		} else {
			for (int i = 1; i < A.length; i++) {
				if (A[i].length != A[0].length) {
					return null;
				}
			}
		}
		BigDecimal[] res = new BigDecimal[A[0].length];
		for (int i = 0; i < A[0].length; i++) {
			res[i] = BigDecimal.ZERO;
			for (int j = 0; j < star.length; j++) {
				res[i] = res[i].add(star[j].multiply(A[j][i]));
			}
		}
		return res;
	}

	/**
	 * AB = (rotA(A,B[0]),rotA(A,B[1]),rotA(A,B[2])).
	 */
	public static BigDecimal[][] rotAB(final BigDecimal[][] A,
			final BigDecimal[][] B) {
		BigDecimal[][] res = new BigDecimal[B.length][];
		for (int i = 0; i < B.length; i++) {
			res[i] = rotA(A, B[i]);
		}
		return res;
	}

	public static BigDecimal[] getApproximation(final BigDecimal[] y,
			final BigDecimal unit) {
		if (y == null || y.length == 0) {
			return null;
		} else if (y.length == 1) {
			return new BigDecimal[] { y[0].plus() };
		} else if (unit == null || unit.signum() == 0) {
			return new BigDecimal[] { y[(y.length + 1) / 2].plus() };
		}
		BigDecimal[] res = new BigDecimal[y.length];
		for (int i = 0; i < y.length; i++) {
			res[i] = BigDecimal.ZERO;
		}
		for (int i = 0; i < y.length; i++) {
			ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
			temp.add(y[i].plus());
			for (int j = 0; j < y.length; j++) {
				if (i != j) {
					BigDecimal dif = new BigDecimal(j - (y.length + 1) / 2)
							.multiply(unit);
					temp.add(0, temp.get(0).multiply(dif));
					for (int k = 1; k < temp.size() - 1; k++) {
						temp.set(k,
								temp.get(k + 1).multiply(dif).add(temp.get(k)));
					}
					dif = new BigDecimal(j - i).multiply(unit);
					for (int k = 0; k < temp.size(); k++) {
						temp.set(k, temp.get(k).divide(dif, DEMath.context));
					}
				}
			}
			for (int j = 0; j < y.length; j++) {
				res[j] = res[j].add(temp.get(j));
			}
		}
		return res;
	}

	public static BigDecimal approximationR(final BigDecimal[] coef,
			final BigDecimal x) {
		if (coef == null || coef.length == 0) {
			return BigDecimal.ZERO;
		} else if (x == null || x.signum() == 0) {
			return coef[0].plus();
		}
		BigDecimal y = coef[coef.length - 1].plus();
		for (int i = coef.length - 2; i >= 0; i--) {
			y = y.multiply(x).add(coef[i]);
		}
		return y;
	}

	public static BigDecimal approximationV(final BigDecimal[] coef,
			final BigDecimal x) {
		if (coef == null || coef.length <= 1) {
			return BigDecimal.ZERO;
		} else if (x == null || x.signum() == 0) {
			return coef[1].plus();
		}
		BigDecimal y = new BigDecimal(coef.length - 1)
				.multiply(coef[coef.length - 1]);
		for (int i = coef.length - 2; i >= 1; i--) {
			y = y.multiply(x).add(new BigDecimal(i).multiply(coef[i]));
		}
		return y;
	}
}

JAVAを用いた多倍長ライブラリ

BigDecimalで多倍長計算するライブラリです。

package jp.wshounen;

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;

public class BigMath {
	private int precision;
	private MathContext context;
	private BigDecimal pi;
	private BigDecimal e;

	public BigMath(int precision) {
		this.precision = precision;
		context = new MathContext((int) Math.ceil(precision * 1.125),
				RoundingMode.HALF_EVEN);
		pi = null;
		e = null;
	}

	public BigDecimal round(BigDecimal in) {
		return in.round(new MathContext(precision, RoundingMode.HALF_EVEN));
	}

	public BigDecimal floor(BigDecimal in) {
		BigDecimal[] bi = in.divideAndRemainder(BigDecimal.ONE);
		if (bi[1].signum() == -1) {
			return bi[0].add(BigDecimal.ONE.negate());
		} else {
			return bi[0];
		}
	}

	public BigDecimal ceil(BigDecimal in) {
		BigDecimal[] bi = in.divideAndRemainder(BigDecimal.ONE);
		if (bi[1].signum() == 1) {
			return bi[0].add(BigDecimal.ONE);
		} else {
			return bi[0];
		}
	}

	public BigDecimal pi() {
		if (pi == null) {
			pi = atan(BigDecimal.ONE).multiply(new BigDecimal(4), context);
		}
		return pi;
	}

	public BigDecimal e() {
		if (e == null) {
			e = exp(BigDecimal.ONE);
		}
		return e;
	}

	public MathContext context() {
		return context;
	}

	public int precision() {
		return precision;
	}

	public BigDecimal asin(BigDecimal in) {
		int num = in.abs().compareTo(BigDecimal.ONE);
		if (num == 0) {
			return pi().divide(new BigDecimal(in.signum() * 2), context);
		} else if (num == -1) {
			BigDecimal bi = atan(in.divide(
					root(in.pow(2).negate().add(BigDecimal.ONE), 2).add(
							BigDecimal.ONE), context));
			return bi.add(bi, context);
		} else if (round(in.abs()).compareTo(BigDecimal.ONE) == 0) {
			return pi().divide(new BigDecimal(in.signum() * 2), context);
		} else {
			return null;
		}
	}

	public BigDecimal acos(BigDecimal in) {
		if (in.round(context).compareTo(BigDecimal.ONE) == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bi = asin(in);
		if (bi == null) {
			return bi;
		} else {
			return pi().divide(new BigDecimal(2), context).add(bi.negate(),
					context);
		}
	}

	public BigDecimal atan(BigDecimal in) {
		int num = in.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bi = in.abs().round(context);
		BigDecimal res = BigDecimal.ZERO;
		BigDecimal perfive = new BigDecimal("0.2");
		if (bi.compareTo(perfive) == 1) {
			BigDecimal perfiveat = atan(perfive);
			while (bi.signum() == 1) {
				bi = bi.add(perfive.negate()).divide(
						bi.multiply(perfive).add(BigDecimal.ONE), context);
				res = res.add(perfiveat);
			}
		}
		if (num == -1) {
			bi = bi.negate();
			res = res.negate();
		}
		BigDecimal coef = bi.round(context);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		temp.add(coef.add(res, context));
		bi = bi.pow(2).negate();
		do {
			num = temp.size();
			coef = coef.multiply(bi, context);
			temp.add(temp.get(num - 1).add(
					coef.divide(new BigDecimal(2 * num + 1), context), context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		return temp.get(num);
	}

	public BigDecimal root(BigDecimal in, int index) {
		int num = in.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			return null;
		} else if (index <= 0) {
			return null;
		} else if (in.compareTo(BigDecimal.ONE) == 0 || index == 1) {
			return in.round(context);
		}
		BigDecimal bi = in.round(context);
		int i = index * 1;
		int j;
		do {
			j = 2;
			while (i % j != 0 && i > j * j) {
				j++;
			}
			if (i == j * j) {
				j = i;
			} else if (i % j == 0) {
				i = i / j;
			} else {
				j = i;
			}
			if (j > 5) {
				bi = pow(bi, BigDecimal.ONE.divide(new BigDecimal(j), context));
			} else {
				ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
				temp.add(bi);
				do {
					num = temp.size();
					temp.add(temp
							.get(num - 1)
							.pow(j)
							.multiply(new BigDecimal(j - 1))
							.add(bi)
							.divide(temp.get(num - 1).pow(j - 1)
									.multiply(new BigDecimal(j)), context));
				} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
				bi = temp.get(num);
			}
		} while (i > j);
		return bi;
	}

	public BigDecimal log(BigDecimal in) {
		if (in.signum() < 1) {
			return null;
		}
		int num = in.compareTo(BigDecimal.ONE);
		BigDecimal coef, bi, res;
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			bi = in.round(context);
			coef = BigDecimal.ONE.negate();
		} else {
			bi = BigDecimal.ONE.divide(in, context);
			coef = BigDecimal.ONE;
		}
		res = BigDecimal.ZERO;
		do {
			bi = bi.multiply(e(), context);
			res = res.add(BigDecimal.ONE);
		} while (bi.compareTo(BigDecimal.ONE) == -1);
		if (bi.pow(5).compareTo(e()) == 1) {
			bi = root(bi, 5);
			coef = coef.multiply(new BigDecimal(5));
		}
		if (num == -1) {
			res = res.negate();
		}
		bi = BigDecimal.ONE.negate().add(BigDecimal.ONE.divide(bi, context));
		coef = coef.multiply(bi, context);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		temp.add(coef.add(res, context));
		bi = bi.negate();
		do {
			num = temp.size();
			coef = coef.multiply(bi, context);
			temp.add(temp.get(num - 1).add(
					coef.divide(new BigDecimal(num + 1), context), context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		return temp.get(num);
	}

	public BigDecimal sin(BigDecimal in) {
		BigDecimal[] bi = in.divideAndRemainder(pi(), context);
		bi[1] = bi[1].divide(new BigDecimal(3), context);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		if (bi[0].abs().remainder(new BigDecimal(2), context)
				.compareTo(BigDecimal.ONE) == 0) {
			temp.add(bi[1].negate(context));
		} else {
			temp.add(bi[1].round(context));
		}
		BigDecimal coef = temp.get(0);
		bi[1] = bi[1].pow(2);
		int num;
		do {
			num = temp.size();
			coef = coef.multiply(bi[1]).divide(
					new BigDecimal(-2 * num * (2 * num + 1)), context);
			temp.add(temp.get(num - 1).add(coef, context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		bi[1] = temp.get(num).pow(3).negate().multiply(new BigDecimal(4))
				.add(temp.get(num).multiply(new BigDecimal(3)), context);
		num = bi[1].signum();
		return bi[1].add(new BigDecimal(-num), context).add(
				new BigDecimal(num), context);
	}

	public BigDecimal cos(BigDecimal in) {
		BigDecimal[] bi = in.divideAndRemainder(pi(), context);
		bi[1] = bi[1].divide(new BigDecimal(3), context).pow(2);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		if (bi[0].abs().remainder(new BigDecimal(2), context)
				.compareTo(BigDecimal.ONE) == 0) {
			temp.add(BigDecimal.ONE.negate());
		} else {
			temp.add(BigDecimal.ONE);
		}
		BigDecimal coef = temp.get(0);
		int num;
		do {
			num = temp.size();
			coef = coef.multiply(bi[1]).divide(
					new BigDecimal(-2 * num * (2 * num - 1)), context);
			temp.add(temp.get(num - 1).add(coef, context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		bi[1] = temp
				.get(num)
				.pow(3)
				.multiply(new BigDecimal(4))
				.add(temp.get(num).negate().multiply(new BigDecimal(3)),
						context);
		num = bi[1].signum();
		return bi[1].add(new BigDecimal(-num), context).add(
				new BigDecimal(num), context);
	}

	public BigDecimal tan(BigDecimal in) {
		BigDecimal bi = cos(in);
		if (bi.add(BigDecimal.ONE, context).compareTo(BigDecimal.ONE) == 0) {
			return null;
		} else {
			return sin(in).divide(bi, context);
		}
	}

	public BigDecimal exp(BigDecimal in) {
		int num = in.signum();
		BigDecimal coef = BigDecimal.ONE;
		BigDecimal bi;
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num == 1) {
			bi = in.negate(context);
		} else {
			bi = in.round(context);
		}
		if (in.abs().round(context).compareTo(BigDecimal.ONE) != 0) {
			while (bi.signum() == -1) {
				bi = bi.add(BigDecimal.ONE);
				coef = coef.multiply(e(), context);
			}
		}
		if (num == 1) {
			bi = bi.negate();
		} else {
			coef = BigDecimal.ONE.divide(coef, context);
		}
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		temp.add(coef.round(context));
		do {
			num = temp.size();
			coef = coef.multiply(bi).divide(new BigDecimal(num), context);
			temp.add(temp.get(num - 1).add(coef, context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		return temp.get(num);
	}

	public BigDecimal pow(BigDecimal in, BigDecimal index) {
		BigDecimal bi = log(in);
		if (bi == null) {
			return null;
		} else {
			return exp(bi.multiply(index));
		}
	}
}

暦計算用の20桁精度の定数付きライブラリです。

package jp.wshounen;

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;
import java.util.ArrayList;

public class DEMath {
	public final static MathContext context = new MathContext(24,
			RoundingMode.HALF_EVEN);
	public static final int precision = 20;
	public static final BigDecimal au = new BigDecimal(149597870691l)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal emrat = new BigDecimal(8130056)
			.scaleByPowerOfTen(-5);
	public static final BigDecimal clight = new BigDecimal(2997922458l)
			.scaleByPowerOfTen(-4);
	public static final BigDecimal jultimeOfJ2000 = new BigDecimal(2451545);
	public static final BigDecimal equatorialRadius = new BigDecimal(6378137)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal polarRadius = new BigDecimal(6356752)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal solarRadius = new BigDecimal(695989245)
			.scaleByPowerOfTen(-3);
	public static final BigDecimal lunarRadius = new BigDecimal(17380908)
			.scaleByPowerOfTen(-4);
	public static final BigDecimal PI = atan(BigDecimal.ONE).multiply(
			new BigDecimal(4), context);
	public static final BigDecimal E = exp(BigDecimal.ONE);
	public static final BigDecimal[][] icrfRot = rotAB(
			rotAB(rotX(arcsecToRad(new BigDecimal(-68192).scaleByPowerOfTen(-7))),
					rotY(arcsecToRad(new BigDecimal(-166170)
							.scaleByPowerOfTen(-7)))),
			rotZ(arcsecToRad(new BigDecimal(146).scaleByPowerOfTen(-4))));
	public static final BigDecimal obliquityArcsecOfJ2000 = new BigDecimal(
			84381406).scaleByPowerOfTen(-3);

	public static BigDecimal round(final BigDecimal in, final int precision) {
		return in.round(new MathContext(precision, RoundingMode.HALF_EVEN));
	}

	public static long divideToFloor(final long value, final long divisor) {
		long remainder = value % divisor;
		return ((remainder < 0 && divisor > 0) || (remainder > 0 && divisor < 0)) ? (value - remainder)
				/ divisor - 1
				: (value - remainder) / divisor;
	}

	public static BigDecimal round(final BigDecimal in) {
		return in.setScale(0, RoundingMode.HALF_EVEN);
		// BigDecimal[] bi = in.divideAndRemainder(BigDecimal.ONE);
		// int num = bi[1].add(bi[1]).abs().compareTo(BigDecimal.ONE);
		// if (num == 1
		// || (num == 0 && bi[0].remainder(new BigDecimal(2)).signum() != 0)) {
		// return bi[0].add(new BigDecimal(bi[1].signum()));
		// } else {
		// return bi[0];
		// }
	}

	public static BigDecimal arcsecToRad(final BigDecimal arcsec) {
		return arcsec.multiply(PI).divide(new BigDecimal(648000), context);
	}

	public static BigDecimal radToArcsec(final BigDecimal rad) {
		return rad.multiply(new BigDecimal(648000)).divide(PI, context);
	}

	public static BigDecimal floor(final BigDecimal in) {
		return in.setScale(0, RoundingMode.FLOOR);
		// BigDecimal[] bi = in.divideAndRemainder(BigDecimal.ONE);
		// if (bi[1].signum() == -1) {
		// return bi[0].add(BigDecimal.ONE.negate());
		// } else {
		// return bi[0];
		// }
	}

	public static BigDecimal ceil(final BigDecimal in) {
		return in.setScale(0, RoundingMode.CEILING);
		// BigDecimal[] bi = in.divideAndRemainder(BigDecimal.ONE);
		// if (bi[1].signum() == 1) {
		// return bi[0].add(BigDecimal.ONE);
		// } else {
		// return bi[0];
		// }
	}

	public static BigDecimal asin(final BigDecimal in) {
		int num = in.abs().compareTo(BigDecimal.ONE);
		if (num == 0) {
			return PI.divide(new BigDecimal(in.signum() * 2), context);
		} else if (num == -1) {
			BigDecimal bi = atan(in.divide(
					root(in.pow(2).negate().add(BigDecimal.ONE), 2).add(
							BigDecimal.ONE), context));
			return bi.add(bi, context);
		} else if (round(in.abs(), 20).compareTo(BigDecimal.ONE) == 0) {
			return PI.divide(new BigDecimal(in.signum() * 2), context);
		} else {
			return null;
		}
	}

	public static BigDecimal acos(final BigDecimal in) {
		if (in.round(context).compareTo(BigDecimal.ONE) == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bi = asin(in);
		if (bi == null) {
			return bi;
		} else {
			return PI.divide(new BigDecimal(2), context).add(bi.negate(),
					context);
		}
	}

	public static BigDecimal atan(final BigDecimal in) {
		int num = in.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		}
		BigDecimal bi = in.abs().round(context);
		BigDecimal res = BigDecimal.ZERO;
		BigDecimal perfive = new BigDecimal(2).scaleByPowerOfTen(-1);
		if (bi.compareTo(perfive) == 1) {
			BigDecimal perfiveat = atan(perfive);
			while (bi.signum() == 1) {
				bi = bi.add(perfive.negate()).divide(
						bi.multiply(perfive).add(BigDecimal.ONE), context);
				res = res.add(perfiveat);
			}
		}
		if (num == -1) {
			bi = bi.negate();
			res = res.negate();
		}
		BigDecimal coef = bi.round(context);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		temp.add(coef.add(res, context));
		bi = bi.pow(2).negate();
		do {
			num = temp.size();
			coef = coef.multiply(bi, context);
			temp.add(temp.get(num - 1).add(
					coef.divide(new BigDecimal(2 * num + 1), context), context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		return temp.get(num);
	}

	public static BigDecimal root(final BigDecimal in, final int index) {
		int num = in.signum();
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			return null;
		} else if (index <= 0) {
			return null;
		} else if (in.compareTo(BigDecimal.ONE) == 0 || index == 1) {
			return in.round(context);
		}
		BigDecimal bi = in.round(context);
		int i = index;
		int j = 2;
		boolean b = true;
		while (i >= j && b) {
			while (i % j != 0 && i > j * j) {
				j++;
			}
			if (i % j == 0) {
				i /= j;
			} else {
				j = i;
				b = false;
			}
			if (j > 5) {
				bi = pow(bi, BigDecimal.ONE.divide(new BigDecimal(j), context));
			} else {
				ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
				temp.add(bi);
				do {
					num = temp.size();
					temp.add(temp
							.get(num - 1)
							.pow(j)
							.multiply(new BigDecimal(j - 1))
							.add(bi)
							.divide(temp.get(num - 1).pow(j - 1)
									.multiply(new BigDecimal(j)), context));
				} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
				bi = temp.get(num);
			}
		}
		return bi;
	}

	public static BigDecimal log(final BigDecimal in) {
		if (in.signum() < 1) {
			return null;
		}
		int num = in.compareTo(BigDecimal.ONE);
		BigDecimal coef, bi, res;
		if (num == 0) {
			return BigDecimal.ZERO;
		} else if (num == -1) {
			bi = in.round(context);
			coef = BigDecimal.ONE.negate();
		} else {
			bi = BigDecimal.ONE.divide(in, context);
			coef = BigDecimal.ONE;
		}
		res = BigDecimal.ZERO;
		do {
			bi = bi.multiply(E, context);
			res = res.add(BigDecimal.ONE);
		} while (bi.compareTo(BigDecimal.ONE) == -1);
		if (bi.pow(5).compareTo(E) == 1) {
			bi = root(bi, 5);
			coef = coef.multiply(new BigDecimal(5));
		}
		if (num == -1) {
			res = res.negate();
		}
		bi = BigDecimal.ONE.negate().add(BigDecimal.ONE.divide(bi, context));
		coef = coef.multiply(bi, context);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		temp.add(coef.add(res, context));
		bi = bi.negate();
		do {
			num = temp.size();
			coef = coef.multiply(bi, context);
			temp.add(temp.get(num - 1).add(
					coef.divide(new BigDecimal(num + 1), context), context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		return temp.get(num);
	}

	public static BigDecimal sin(final BigDecimal in) {
		BigDecimal[] bi = in.divideAndRemainder(PI, context);
		bi[1] = bi[1].divide(new BigDecimal(3), context);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		if (bi[0].abs().remainder(new BigDecimal(2), context)
				.compareTo(BigDecimal.ONE) == 0) {
			temp.add(bi[1].negate(context));
		} else {
			temp.add(bi[1].round(context));
		}
		BigDecimal coef = temp.get(0);
		bi[1] = bi[1].pow(2);
		int num;
		do {
			num = temp.size();
			coef = coef.multiply(bi[1]).divide(
					new BigDecimal(-2 * num * (2 * num + 1)), context);
			temp.add(temp.get(num - 1).add(coef, context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		bi[1] = temp.get(num).pow(3).negate().multiply(new BigDecimal(4))
				.add(temp.get(num).multiply(new BigDecimal(3)), context);
		num = bi[1].signum();
		return bi[1].add(new BigDecimal(-num), context).add(
				new BigDecimal(num), context);
	}

	public static BigDecimal cos(final BigDecimal in) {
		BigDecimal[] bi = in.divideAndRemainder(PI, context);
		bi[1] = bi[1].divide(new BigDecimal(3), context).pow(2);
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		if (bi[0].abs().remainder(new BigDecimal(2), context)
				.compareTo(BigDecimal.ONE) == 0) {
			temp.add(BigDecimal.ONE.negate());
		} else {
			temp.add(BigDecimal.ONE);
		}
		BigDecimal coef = temp.get(0);
		int num;
		do {
			num = temp.size();
			coef = coef.multiply(bi[1]).divide(
					new BigDecimal(-2 * num * (2 * num - 1)), context);
			temp.add(temp.get(num - 1).add(coef, context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		bi[1] = temp
				.get(num)
				.pow(3)
				.multiply(new BigDecimal(4))
				.add(temp.get(num).negate().multiply(new BigDecimal(3)),
						context);
		num = bi[1].signum();
		return bi[1].add(new BigDecimal(-num), context).add(
				new BigDecimal(num), context);
	}

	public static BigDecimal tan(final BigDecimal in) {
		BigDecimal bi = cos(in);
		if (bi.add(BigDecimal.ONE, context).compareTo(BigDecimal.ONE) == 0) {
			return null;
		} else {
			return sin(in).divide(bi, context);
		}
	}

	public static BigDecimal exp(final BigDecimal in) {
		int num = in.signum();
		BigDecimal coef = BigDecimal.ONE;
		BigDecimal bi;
		if (num == 0) {
			return BigDecimal.ONE;
		} else if (num == 1) {
			bi = in.negate(context);
		} else {
			bi = in.round(context);
		}
		if (in.abs().round(context).compareTo(BigDecimal.ONE) != 0) {
			while (bi.signum() == -1) {
				bi = bi.add(BigDecimal.ONE);
				coef = coef.multiply(E, context);
			}
		}
		if (num == 1) {
			bi = bi.negate();
		} else {
			coef = BigDecimal.ONE.divide(coef, context);
		}
		ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
		temp.add(coef.round(context));
		do {
			num = temp.size();
			coef = coef.multiply(bi).divide(new BigDecimal(num), context);
			temp.add(temp.get(num - 1).add(coef, context));
		} while (temp.get(num).compareTo(temp.get(num - 1)) != 0);
		return temp.get(num);
	}

	public static BigDecimal pow(final BigDecimal in, final BigDecimal index) {
		BigDecimal bi = log(in);
		if (bi == null) {
			return null;
		} else {
			return exp(bi.multiply(index));
		}
	}

	/**
	 * Matrix of rotation Y to Z. R1(-rad).
	 */
	public static BigDecimal[][] rotX(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] {
				{ BigDecimal.ONE, BigDecimal.ZERO, BigDecimal.ZERO },
				{ BigDecimal.ZERO, Cos, Sin },
				{ BigDecimal.ZERO, Sin.negate(), Cos } };
	}

	/**
	 * Matrix of rotation X to Z. R2(rad).
	 */
	public static BigDecimal[][] rotY(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, BigDecimal.ZERO, Sin },
				{ BigDecimal.ZERO, BigDecimal.ONE, BigDecimal.ZERO },
				{ Sin.negate(), BigDecimal.ZERO, Cos } };
	}

	/**
	 * Matrix of rotation X to Y. R3(-rad).
	 */
	public static BigDecimal[][] rotZ(final BigDecimal rad) {
		BigDecimal Cos = cos(rad);
		BigDecimal Sin = sin(rad);
		return new BigDecimal[][] { { Cos, Sin, BigDecimal.ZERO },
				{ Sin.negate(), Cos, BigDecimal.ZERO },
				{ BigDecimal.ZERO, BigDecimal.ZERO, BigDecimal.ONE } };
	}

	/**
	 * Rotation (1,0,0) to A[0] = (x,y,z). Rotation (0,1,0) to A[1] = (x,y,z).
	 * Rotation (0,0,1) to A[2] = (x,y,z). Ar = rotA(A,r).
	 */
	public static BigDecimal[] rotA(final BigDecimal[][] A,
			final BigDecimal[] star) {
		BigDecimal[] res = new BigDecimal[] { BigDecimal.ZERO, BigDecimal.ZERO,
				BigDecimal.ZERO };
		for (int i = 0; i < 3; i++)
			for (int j = 0; j < 3; j++)
				res[i] = res[i].add(star[j].multiply(A[j][i]));
		return res;
	}

	/**
	 * AB = (rotA(A,B[0]),rotA(A,B[1]),rotA(A,B[2])).
	 */
	public static BigDecimal[][] rotAB(final BigDecimal[][] A,
			final BigDecimal[][] B) {
		return new BigDecimal[][] { rotA(A, B[0]), rotA(A, B[1]), rotA(A, B[2]) };
	}

	public static BigDecimal[] getApproximation(final BigDecimal[] y,
			final BigDecimal unit) {
		if (y == null || y.length == 0) {
			return null;
		} else if (y.length == 1) {
			return new BigDecimal[] { y[0].plus() };
		} else if (unit == null || unit.signum() == 0) {
			return new BigDecimal[] { y[(y.length + 1) / 2].plus() };
		}
		BigDecimal[] res = new BigDecimal[y.length];
		for (int i = 0; i < y.length; i++) {
			res[i] = BigDecimal.ZERO;
		}
		for (int i = 0; i < y.length; i++) {
			ArrayList<BigDecimal> temp = new ArrayList<BigDecimal>(0);
			temp.add(y[i].plus());
			for (int j = 0; j < y.length; j++) {
				if (i != j) {
					BigDecimal dif = new BigDecimal(j - (y.length + 1) / 2)
							.multiply(unit);
					temp.add(0, temp.get(0).multiply(dif));
					for (int k = 1; k < temp.size() - 1; k++) {
						temp.set(k,
								temp.get(k + 1).multiply(dif).add(temp.get(k)));
					}
					dif = new BigDecimal(j - i).multiply(unit);
					for (int k = 0; k < temp.size(); k++) {
						temp.set(k, temp.get(k).divide(dif, DEMath.context));
					}
				}
			}
			for (int j = 0; j < y.length; j++) {
				res[j] = res[j].add(temp.get(j));
			}
		}
		return res;
	}

	public static BigDecimal approximationR(final BigDecimal[] coef,
			final BigDecimal x) {
		if (coef == null || coef.length == 0) {
			return BigDecimal.ZERO;
		} else if (x == null || x.signum() == 0) {
			return coef[0].plus();
		}
		BigDecimal y = coef[coef.length - 1].plus();
		for (int i = coef.length - 2; i >= 0; i--) {
			y = y.multiply(x).add(coef[i]);
		}
		return y;
	}

	public static BigDecimal approximationV(final BigDecimal[] coef,
			final BigDecimal x) {
		if (coef == null || coef.length <= 1) {
			return BigDecimal.ZERO;
		} else if (x == null || x.signum() == 0) {
			return coef[1].plus();
		}
		BigDecimal y = new BigDecimal(coef.length - 1)
				.multiply(coef[coef.length - 1]);
		for (int i = coef.length - 2; i >= 1; i--) {
			y = y.multiply(x).add(new BigDecimal(i).multiply(coef[i]));
		}
		return y;
	}
}

Macintoshで月影図1枚描くのに8分かかるので高速化の作業に入ります。だけどbyte配列の数を使わないといけないような計算はこの世に存在するのか？boolean配列の多倍長ライブラリもつくってみようかと思ったけどboolean配列だと無駄が多すぎて実装も大変なのでやめます。

このソースコードはMITライセンスですが万人に自由に使っていただくためJAVAを用いた多倍長計算の上記の機能の改良版について特許出願することを禁止いたします。

月影図を表示します

若くないと感じる時、それは計算の勘が鈍ったとき。最近数ヶ月使わなければすぐプログラミングのやり方を忘れて過去の作品が遺構になったりAndroidのプロジェクトが全部動かなくなったり。何事も日頃の行いがものを言いますね。

ブログネタ: 【賞品付き】もう若くないと感じた瞬間は？

来年の金環日蝕の月縁図のサイト

バスに大学時代の知人にそっくりな人が乗っていてびっくりしました。来年の金環日蝕の際の月縁図です。東京スカイツリーの直下の座標がセットされています。月縁付近を誇張して描いているため外向きに曲がる曲線は太陽の輪郭の一部です。

2012年5月21日の金環食の月縁図(東京スカイツリー直下)

このサイトは日本国内の座標で表示した場合だけ正確な月縁図を表示します。このURLは以下のようになっています。

http://wshounen.la.coocan.jp/20120521/contourImage.php?nd=緯度(度)&ed=経度(度)&ht=標高(m)

児童教育のために別ブログをやっていますが、重要度の高い記事が発生したので公開を一日遅らせます。向こうの記事も大事なのですが、緊急を要する状況ではないためです。自分の力で考えられないような大人のために記事を書くより全世界に貢献する仕事を遂行するための準備をするほうが重要です。

月縁図の進捗状況

今日ものどかな朝です。もしかしたらこのブログでは初公開かもしれませんが、来年の5月21日の金環食のときの月縁図を表示してみました。月縁を誇張して描くため月縁付近の距離を外側と内側に引き延ばしてみました。なので、月の輪郭が実際よりもギザギザが強調され、金環食の最初と最後の付近の時刻では月縁付近にかかる太陽の輪郭が外側に押し広げられています。太陽は日蝕のときに右から左へ動くので右側の曲線が金環食開始時刻付近、左側の曲線が金環食終了時刻付近です。

夏至の昼は節電チャレンジ

小鳥のさえずりとともに新しい朝が始まる。眠れない夜はずっと続くわけではない。今頃の早朝は明るくなるのが一年で最も早く外は小鳥の合唱ですがすがしい朝をむかえられます。今日は夏至の日。もう夏至は数時間過ぎていますが、昨日の夕方は晴れて蒸し暑く夏至に一番近いダイヤモンド富士はどうなんだろうと駅の空を眺めていました。今日は横浜では13時から15時まで節電チャレンジです。事業所、家庭などでどれくらい省エネできるか試される日です。

省エネと言えばLED電球の省エネ効果はかなり期待度BIGです。計画停電中もどうしても電灯が欲しくなったときに乾電池でLED電球を照らしたら太陽のようにすごく明るくてびっくりしました。それも、父の日曜大工で4個の電球をつくっていました。さすがに暑い日はちょっといやだけど涼しい日であればガスが止まるわけではないので乾電池LED電球で計画停電も乗り切れます。

ブログネタ: やってみようと思う節電対策は？

執念

休みの日はごろ寝・・・ではなく家で天文計算をしています。もしかしたら将来の論文のネタがここから生まれるかもしれませんよ。趣味の探求で大いに青春しています。最近はパソコンで時間のかかる数値計算をすることが多くなったので扇風機は欠かせません。網の上にパソコンを置いて真横から扇風機の風を連続的に当てるのです。完璧ではありませんが冷却効果が期待できパソコンのヒョーズが飛ぶのを防ぎます。科学計算は今でもフォートランが主流ですが私はJAVAで行います。多倍長ライブラリは自分で書きました。問題なく動いていますよ。JAVAで行うのは完成したときにアプリに加工できるからです。ある計算がどうしても時間の関数として表現できないのか試していたら朝になってしまいました。寝て起きたらうまく近似する方法が見つかり、時間軸の関数によるアニメーションの作成が不可能ではないことが分かりました。

ブログネタ: 【賞品付き】休みの日は何してる？

画像圧縮

去年面白いデータ圧縮方法を思いつきました。数値ばかりのデータをPNG画像として保存するのです。JAVAだとImageIOで数値データから作成されたPNG画像を可逆的に作成できるので科学計算の係数を画像として保持する事ができます。例えば、DEシリーズのデータはlong型の有効数字とbyte型の指数部分で表現できるので、一部long型の最大値を超えているデータもありますが1桁削って18桁にしても計算精度にはさほど影響しません。もしかしたらDE405は三倍精度実数である拡張倍精度型で計算されているのかもしれません。DE406は完全に倍精度実数の範囲内なので画像化したデータで完全に可逆的に展開できます。以前お見せした非公式ミラーサイトも実はそのしくみを使っているのです。有効数字と指数で9byteあるのでそれを3byteずつRGBカラーに押し込みます。ただし、画像なのでうっかり編集されないようにRGBカラーのチェックビットをつけます。もしチェックビットが検証値と異なる場合はNumberFormatExceptionを発するようにしてあります。Androidでも同様にできるか試してみたけどRGB展開の際に色空間が縮小されるためかうまくいきませんでした。

SELENEの月縁データ：地球は青かった

やる気を出す方法は月を拝む事です。見たまえ、この緑あふれる月を。

(c)JAXA,(c)WADA

地球は青かったというのはガガーリンの言葉の訳語ですが、ガガーリンが見た青い地球をかぐやの月面データで表現してみました。どう、ゲージツ的でしょう。この画像はJAXAから月面標高データをいただいて標高が基準面より低い場合は青く、高い場合は緑色に、ほぼ同じ場合は白く彩色しています。元データは緯度経度ともに1/16度間隔のデータリストでJAVAを使って元のデータから数値処理をして描画しているためこの絵は正確です。ただし、無限遠からの投影で標高は彩色以外では考慮していません。雲をまとった海と緑のある地球のように見えてきませんか。あれ、ウサギはどこでしょう、お休みでしょうか。いや、これは左が表側、右が裏側です。よく見ているとウサギが見えてきませんか。

やる気を出す方法はずばり今一番やりたい事をやること。忙しくてダメというときは歩くとやる気が出てくるので歩くのもいいでしょう。歩くと血行がよくなって適度に脳を刺激するので気分が前向きになります。会社員の生き方と就職活動の方法について書いている吉田典史様のtwitterでも歩く事がモチベーション向上に最適と書かれており、確かに歩くと多少疲れている時でもやる気が回復します。SEとかは平日は激務ですが週末は休養日で寝て過ごすのではなく外に出て心の栄養を充電しましょう。適度な運動は体の健康維持にも役立ちますので基礎体力の低下を防ぎ、それがモチベーションの向上にもつながります。病は気からと言いますが体の疲れは気の疲れ。時間が許せばパール富士などで自然の神秘を見つめるのもいいかもしれません。

ブログネタ: 【賞品付き】ヤル気を出す方法、教えて！