W少年のお散歩日記: JAVA

2012年4月
日	月	火	水	木	金	土
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30

携帯URL

お気に入りサイト

ブログパーツ

Please Comment!

Your Mail-Address

Your Mail-Address again.

Comment
天気予報と郵便番号検索

上のブログパーツの天気予報検索に必要な7桁の郵便番号を調べます。郵便を送るときに郵便番号を調べるのにも便利です。
うまく翻訳できるでしょうか
ジャンケンゲーム
地球上の距離の計算

JAVA

JAVAを用いた多倍長ライブラリがようやく完成しました。でもこれは始まりでしかないのです。ライブラリには完成はなくそのときの時代に応じて最適化されていくものです。なので、逐次更新できるような仕組みを用意しました。まずはAbstractMathからです。多倍長ライブラリの主計算メソッドを集めてユーザーの好みに合わせてチューニングできるようにしてあります。このクラスを拡張して多倍長ライブラリにします。

多倍長ライブラリ本体であるBigMathです。BigDecimalを使って計算するようになっています。計算可能桁数は現実的に全く制約を感じないくらい多いです。

2011年11月12日 (土) 11時20分 JAVA, プログラマー日記 | 固定リンク
Tweet

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

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

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;
	}
}

2011年10月 4日 (火) 23時48分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

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;
	}
}

2011年10月 2日 (日) 09時00分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

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を用いた多倍長計算の上記の機能の改良版について特許出願することを禁止いたします。

2011年9月30日 (金) 22時46分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

DKA法で日蝕要素の計算をやってみました

もしロボットが家にいたら私に代わってコーディングしてほしいかな。面倒な計算のアルゴリズムを考えて変数などの処理もしてアプリをつくってくれるロボットがいたらいいな。

今日は日蝕要素の計算だけです。緯度経度を指定できるオンデマンドな計算は今回は行っていません。あらかじめ用意された日蝕の日付で1分おきに日出入の限界線、南北限界線、中心線の瞬時の座標を計算するだけです。

ブログネタ: もしロボットが家にいたら何をしてもらう？参加数

2011年6月 5日 (日) 09時46分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

DE406の疑似アスキーファイル提供サイト

DE406のアスキーファイルは圧縮ファイルで提供されていて本家からダウンロードするといったん圧縮ファイルを解凍しないといけないので非公式ミラーサイトをつくってしまいました。ここで提供するアスキーファイルの数値は本家の数値と値は同じですが本家で期間が隣り合うファイルで重複している期間の片方を削除しているのでデータセットの冒頭の番号が変わる可能性があります。ここで提供しているアスキーファイルを生成するデータファイルは私のホームページのDE406のデータファイルそのものなので、生成されたアスキーファイルの使用条件はDE406本家の使用条件と同じですがデータファイルそのものを探したりコピーして使用したりする事はしないでください。

非公式ミラーサイトを起動する

2011年6月 2日 (木) 18時47分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

DE406の簡易計算(Sample of DE406 prediction)

DE406をJAVAで計算します。入力可能な日付はDE406で用意された範囲で西暦-3000年頃から西暦3000年頃までです。ただし、日付はユリウス日で入力するようになっています。初期値はJ2000、つまり力学時の西暦2000年1月1日正午になっております。計算では24桁使っていますが答えの有効数字は12桁です。結果の座標は位置、速度の順に並んでいます。単位はkm、km/dayです。

Execute DE406 with JAVA. Allowed date is about between B.C. 3000 and A.D. 3000 what DE406 equipped with. Please input date as Julian day. Default value is J2000 that is Jan. 1 A.D. 2000 at noon TT. Prediction shall execute 24 digits, but precision of result is 12 digits. Coordinate of result is listed Location and Velocity. Units are km and km/day.

2011年5月 9日 (月) 00時25分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

DE405の簡易計算(Sample of DE405 prediction)

DE405をJAVAで計算します。入力可能な日付はDE405で用意された範囲で西暦1600年頃から西暦2200年頃までです。ただし、日付はユリウス日で入力するようになっています。初期値はJ2000、つまり力学時の西暦2000年1月1日正午になっております。計算では24桁使っていますが答えの有効数字は12桁です。結果の座標は位置、速度の順に並んでいます。章動(Nutation)と秤動(Libration)の単位はrad、rad/dayでその他はkm、km/dayです。

Execute DE405 with JAVA. Allowed date is about between A.D. 1600 and A.D. 2200 what DE405 equipped with. Please input date as Julian day. Default value is J2000 that is Jan. 1 A.D. 2000 at noon TT. Prediction shall execute 24 digits, but precision of result is 12 digits. Coordinate of result is listed Location and Velocity. Units of Nutation and Libration are rad and rad/day. Others are km and km/day.

2011年5月 5日 (木) 08時48分 JAVA, プログラマー日記, 天文 | 固定リンク
Tweet

使えない電卓

Cocolog_oekaki_2011_04_28_22_37 つくった本人も完成品なのか分からない。入力アシスト付き16進多倍長電卓。2進数、4進数、8進数、10進数、16進数、桁数無制限の最強の電卓。算木の作成もけっこうやっかいだったけど入力アシストも難儀。何か数字を入力するまでどんな計算ができるのか分からない、ミステリー電卓。世の中割り切れないことだらけだけど、この電卓の中だけはきれいに割り切れる。嘘です、余分に計算して割り切れない部分は四捨五入しました。そのため割り切れない数でも割り切れる数のようにきれいに計算できます。入力アシストにバグがあったとしてもつくった本人は気付いていません。今日は疲れたのでこの電卓を投稿して寝ます。

ブログネタ: あなたの「コネタ」を募集中です。参加数

2011年4月28日 (木) 22時50分 JAVA, プログラマー日記 | 固定リンク
Tweet

足し算と引き算だけの16進数電卓

足し算と引き算しかしません。それは算木が不要な完璧な直列計算で誤差の出る余地がないからです。結果は完璧に誤差のない正確な計算を行います。というのは2進数も4進数も8進数も16進数も10進数で正確に計算できるからです。計算は10進数のBigDecimalにお任せして結果と入力の部分だけごにょごにょっとやりました。四則計算の電卓の公表予定は明日ではなく未定です。算木と表示精度の扱いにけっこう苦労しています。

2011年4月27日 (水) 00時00分 JAVA, プログラマー日記 | 固定リンク
Tweet

より以前の記事一覧

その他のカテゴリー

W少年のお散歩日記

最近の記事