PHPのMath関数を再発明してみた

言うまでもなくdankogai氏の1ヶ月位前のブログに触発されて。
普段から使うMath系の関数もあるけど、中で何をやってるのか知らなかったから再発明してみた。
乱数とX進数からX進数へ変換以外の関数は全て。
四角い車輪。
　
極力よく見かける公式とかに合わせた。
速度は気にしない。
Math::atan(1) が一向に収束しないのはどうにかしたいけど、公式をそのまま使うとそうだからとりあえずそれでいい。
　
これ作るのに、2週間位かかった。
数学やったのなんて中学生以来なんじゃないかと。
今までろくに勉強してこなかったのが悔やまれる。
これを作れて、どんな事をやってるのかが知れたからとりあえず満足。
一般教養的な意味でこれを知れたのは有意義だった。
　
Mathクラスと平行して「プログラマーでなくても名前ぐらい覚えておきたいアルゴリズムx10」もPHPで作ってみてる。
Mathクラスはまぁどうにかなっても、これをPHPで作るのはまずかった。
でも途中まではやってあるから、そのまま完成させる予定。
　
参考にしたページはブックマークしといたつもりだったけど、結構な数が行方不明になってしまった。
　
Math関数の再発明をしてて思った事は、数学は何千年も前から存在するオープンソースプロジェクトみたいだと言う事。
オープンソースというよりはRFCとかそっちのような気もするけど。
昔の人はすごいなぁ、なんて思いつつ、自分もいつか提供する側になりたいなぁ、なんて思いつつ、寝る。
　

再発明した関数一覧

abs
acos
acosh
asin
asinh
atan
atan2
atanh
ceil
cos
cosh
exp
expm1
deg2rad
floor
fmod
hypot
log
log10
log1p
max
min
pi
pow
round
sin
sinh
sqrt
tan
tanh
rad2deg
　
PHP: Math 関数 - Manual
http://php.net/manual/ja/ref.math.php
　

ソース

<?php
class Math
{
	// 絶対値
	public static function abs($a)
	{
		return ($a<0) ? $a * -1 : $a;
	}

	// 逆余弦(アークコサイン)
	public static function acos($x)
	{
		return self::pi() / 2 - self::asin($x);
	}

	// 逆双曲線余弦(アークハイパボリックコサイン)
	public static function acosh($x)
	{
		return self::log($x + self::sqrt($x * $x - 1));
	}

	// 逆正弦(アークサイン)
	public static function asin($x)
	{
		return self::atan($x / self::sqrt(1 - $x * $x));
	}

	// 逆双曲線正弦(アークハイパボリックサイン)
	public static function asinh($x)
	{
		return self::log($x + self::sqrt($x * $x + 1));
	}

	// 逆正接(アークタンジェント)
	public static function atan($a)
	{
		$pi = self::pi();
		$phase = (self::abs($a) / $pi * 2) % 4;
		$x = self::fmod($a, $pi);
		//printf("phrase: %d, a:%f, x: %f\n", $phase, $a, $x);
		if (1<self::abs($x)) {
			//printf("!!! atan abs %f\n", $a);
			return $pi / 2 - self::atan(1 / $a);
		}
		$result = $x;
		$tmp = 1;
		$minus = 1;
		for ($i=1; $result!=$result+$tmp; $i++) {
			$n = $i * 2 + 1;
			$minus *= -1;
			$tmp = self::pow($x, $n) / $n * $minus;
			$result += $tmp;
			//printf("i:%d, a:%d, m:%d, n:%d, t:%.30f, x:%.30f\n", $i, $i*2+1, $minus, $n, $tmp, $result);
		}
		if (2<=$phase) {
			$result *= -1;
		}
		return $result;
	}

	// 2変数のアークタンジェント
	public static function atan2($y, $x)
	{
		$minus_y = 1;
		$minus_x = 1;
		if ($y<0) {
			$y = self::abs($y);
			$minus_y = -1;
		}
		if ($x<0) {
			$x = self::abs($x);
			$minus_x = -1;
		}
		if ($y==0) {
			return 0;
		}
		if ($x==0) {
			return self::pi() / 2 * $minus_y;
		}
		return self::atan($y / $x * $minus_x) * $minus_y;
	}

	// 逆双曲線正接(アークハイパボリックタンジェント)
	public static function atanh($x)
	{
		if (1-$x==0) {
			return INF;
		}
		return self::log((1+$x) / (1-$x)) / 2;
	}

	// 切り上げ
	public static function ceil($a)
	{
		return (int)(((int)$a!=$a) ? ($a<0 ? $a : $a + 1) : $a);
	}

	// 余弦(コサイン)
	public static function cos($x)
	{
		return self::sin($x + self::pi() / 2);
	}

	// 双曲線余弦(ハイパボリックコサイン)
	public static function cosh($x)
	{
		return (Math::exp($x) + Math::exp($x * -1)) / 2;
	}

	// 累乗
	public static function exp($x)
	{
		if (!is_finite($x)) {
			return NAN;
		}
		$tmp = 1;
		$result = 1;
		for ($i=1; $result!=$result+$tmp; $i++) {
			$tmp = $tmp * $x / $i;
			$result += $tmp;
			//printf("x: %f, tmp: %f, result: %f\n", $x, $tmp, $result);
		}
		return $result;
	}

	// exp(number) - 1
	public static function expm1($x)
	{
		return self::exp($x) - 1;
	}

	// 度単位の数値をラジアン単位に変換
	public static function deg2rad($deg)
	{
		return $deg / 180 * self::pi();
	}

	// 切り捨て
	public static function floor($a)
	{
		return (int)(((int)$a!=$a) ? ($a<0 ? $a - 1 : $a) : $a);
	}

	// 剰余
	public static function fmod($x, $y)
	{
		return $x - $y * (int)($x / $y);
	}

	// 直角三角形の斜辺の長さ
	public static function hypot($x, $y)
	{
		return sqrt($x*$x + $y*$y);
	}

	// 自然対数を返す
	public static function log($a, $base=null)
	{
		if (isset($base) && $base<=0) {
			throw new Exception("log(): base must be greater than 0");
		}
		if ($a<=0) {
			return -INF;
		}

		$x = ($a - 1) / ($a + 1);
		$result = 0;
		$tmp = 1;
		for ($i=0; $result!==$result+$tmp; $i++) {
			$tmp = 2 * (1 / (2 * $i + 1)) * self::pow($x, 2 * $i + 1);
			$result += $tmp;
			//printf("%d: %.30f\n", $i, $result);
		}
		//printf("loop: %d\n", $i);

		if (isset($base)) {
			$base = self::log($base);
			if ($base==0) {
				return NAN;
			}
			return $result / $base;
		}
		return $result;
	}

	// 常用対数を返す
	public static function log10($a)
	{
		return self::log($a, 10);
	}

	// log(1 + number)
	public static function log1p($a)
	{
		return self::log(1 + $a);
	}

	// 最大値
	public static function max($a, $b)
	{
		return $a<$b ? $b : $a;
	}

	// 最小値
	public static function min($a, $b)
	{
		return $a<$b ? $a : $b;
	}

	// 円周率
	public static function pi()
	{
		return self::pi_by_gauss();
	}

	// 円周率 ヴィエートの公式
	public static function pi_by_viete()
	{
		$n = self::sqrt(1/2);
		$result = $n;
		//printf("%f\n", $n);
		while ($result!=$result*$n) {
			$n = self::sqrt(1/2 + 1/2 * $n);
			$result *= $n;
			//printf("%f, %.35f\n", $n, $result);
		}
		return 2 / $result;
	}

	// 円周率 モンテカルロ法
	public static function pi_by_montecalro($n=1000)
	{
		$p = 0;
		for ($i=$n; 0<$i; $i--) {
			$x = mt_rand(0, 10000) / 10000;
			$y = mt_rand(0, 10000) / 10000;
			if ($x*$x+$y*$y<=1) {
				$p++;
			}
			//printf("x:%f y:%f\n", $x, $y);
		}
		return $p / $n * 4;
	}

	// 円周率 ガウス＝ルジャンドルのアルゴリズム
	public static function pi_by_gauss()
	{
		$a = 1;
		$b = 1 / self::sqrt(2);
		$s = 1;
		$t = 4;
		for ($i=0; $i<3; $i++) {
			$last = $a;
			$a = ($a + $b) / 2;
			$b = self::sqrt($last * $b);
			$s -= $t * ($a - $last) * ($a - $last);
			$t *= 2;
			//printf("%.20f\n", ($a + $b) * ($a + $b) / $s);
		}
		return ($a + $b) * ($a + $b) / $s;
	}

	// 指数表現
	public static function pow($base, $exp)
	{
		if ($base==0) {
			return 0.0;
		}

		// integer
		if (is_int($exp)) {
			$result = 1;
			if ($exp<0) {
				$exp = self::abs($exp);
				$base = 1 / $base;
			}
			while ($exp--) {
				$result *= $base;
			}
			return $result;

		// float
		} else {
			return self::exp(self::log($base) * $exp);
		}
	}

	// 四捨五入
	public static function round($a)
	{
		return self::floor($a + 0.5);
	}

	// 正弦(サイン)
	public static function sin($x)
	{
		$pi = self::pi();
		$phase = (self::abs($x) / $pi * 2) % 4;
		$x = self::fmod($x, $pi);
		//printf("phrase: %d, x: %f\n", $phase, $x);
		$result = $x;
		$tmp = 1;
		$n = 1;
		$minus = 1;
		for ($i=1; $result!=$result+$tmp; $i++) {
			$n *= ($i * 2) * ($i * 2 + 1);
			$minus *= -1;
			$tmp = self::pow($x, $i*2+1) / $n * $minus;
			$result += $tmp;
			//printf("i:%d, a:%d, m:%d, n:%d, x:%f\n", $i, $i*2+1, $minus, $n, $result);
		}
		if (2<=$phase) {
			$result *= -1;
		}
		return $result;
	}

	// 双曲線正弦(ハイパボリックサイン)
	public static function sinh($x)
	{
		return (Math::exp($x) - Math::exp($x * -1)) / 2;
	}

	// 平方根
	public static function sqrt($a)
	{
		return self::sqrt_by_newton($a);
	}

	// 平方根 二分探索
	public static function sqrt_by_binarysearch($a)
	{
		$minus = 1;
		if ($a<0) {
			$minus = -1;
			$a *= -1;
		}
		if ($a==1) {
			return $a * $minus;
		}
		$low = 1;
		$high = $a;
		if ($a<1) {
			$low = $a;
			$high = 1;
		}
		$middle = 0;
		while ($middle*0.000000000000001<$high-$low) {
			$middle = 0.5 * ($low + $high);
			if ($middle * $middle < $a) {
				$low = $middle;
			} else {
				$high = $middle;
			}
			//printf("middle: %f\n", $middle);
		}
		return $middle * $minus;
	}

	// 平方根 ニュートン法
	public static function sqrt_by_newton($a)
	{
		if ($a==0) {
			return 0.0;
		}
		$a = self::abs($a);
		$x = $a;
		//var_dump($x);
		while (true) {
			$last = $x;
			$x = ($x + $a / $x) / 2;
			//var_dump($x);
			if ($last==$x) {
				break;
			}
		}
		return $x;
	}

	// 正接(タンジェント)
	public static function tan($x)
	{
		return self::sin($x) / self::cos($x);
	}

	// 双曲線正接(ハイパボリックタンジェント)
	public static function tanh($x)
	{
		return self::sinh($x) / self::cosh($x);
	}

	// ラジアン単位の数値を度単位に変換
	public static function rad2deg($rad)
	{
		return $rad * 180 / self::pi();
	}
}

参考URL

404 Blog Not Found:javascript - Mathを再発明してみた
http://blog.livedoor.jp/dankogai/archives/51518565.html
　
正弦関数 sin(x)をプログラムで実装するにあ？（３） - udongeinの日記
http://d.hatena.ne.jp/udongein/20080725/1216967329
　
AWK Users JP :: 自然対数の計算
http://gauc.no-ip.org/awk-users-jp/blis.cgi/DoukakuAWK_227
　
C言語の二分法のプログラムについて | OKWave
http://okwave.jp/qa/q3745693.html
　
三角関数 - Wikipedia
http://ja.wikipedia.org/wiki/%E4%B8%89%E8%A7%92%E9%96%A2%E6%95%B0
　
双曲線関数 - Wikipedia
http://ja.wikipedia.org/wiki/%E5%8F%8C%E6%9B%B2%E7%B7%9A%E9%96%A2%E6%95%B0
　
モンテカルロ法による円周率の計算
http://hp.vector.co.jp/authors/VA014765/pi/montecalro.html
　
指数関数のマクローリン展開
http://www.geocities.jp/yousukehonda_shuzan/calc/japanese/a0402exp.html
　
ケータイ(Vodafone)＠らるるのわがまま
http://www.rarul.com/k-tai/java/math/