2010年09月14日 06:30 [Edit]

javascript - Mathを再発明してみた

「基本というからには四則演算で三角関数実装しないとねー」と思いつつ書いていたら…

Math.random()を除いてMathを全部再発明しおえたので。

多倍長演算バージョンを作る時の下ごしらえにもなるかも。


下ごしらえ

仕様は

アンチョコはもはや最新というにはあまりに古い、しかし代わりなき「C言語による最新アルゴリズム事典」。低レベルな車輪を再発明する人必携!

初期化と定数

定数の精度はおおげさに。

MyMath = {};

MyMath.E       = 2.718281828459045235360287471352662497757;
MyMath.LN2     = 0.6931471805599453094172321214581765680755;
MyMath.LN10    = 2.302585092994045684017991454684364207601;
MyMath.LOG2E   = 1.442695040888963407359924681001892137427;
MyMath.LOG10E  = 0.4342944819032518276511289189166050822944;
MyMath.PI      = 3.141592653589793238462643383279502884197;
MyMath.SQRT1_2 = 0.7071067811865475244008443621048490392848;
MyMath.SQRT2   = 1.41421356237309504880168872420969807857;

最大、最小

引数なしだと ±Infinity ∓Infinityになるってご存知でした?



MyMath.min = function() {
    if (!arguments.length) return 1 / 0; /* Infinity */
    var r = +arguments[0], i, l, n;
    for (i = 1, l = arguments.length; i < l; i++) {
        n = +arguments[i];
        if (isNaN(n)) return n;
        if (n < r) r = n;
    }
    return r;
};

MyMath.max = function() {
  if (!arguments.length) return -1 / 0; /* Infinity */
    var r = +arguments[0], i, l, n;
    for (i = 1, l = arguments.length; i < l; i++) {
        n = +arguments[i];
        if (isNaN(n)) return n;
        if (n > r) r = n;
    }
    return r;
};

整数化

%の右側って整数でなくてもいいんですね。



MyMath.float = function(n) { return n % 1 };
MyMath.int   = function(n) { return (n - (n % 1)) };
MyMath.abs   = function(n) { return n < 0 ? n * -1 : n };
MyMath.floor = function(n) { return n % 1 ? MyMath.int(n < 0 ? n - 1 : n) : n };
MyMath.ceil  = function(n) { return n % 1 ? MyMath.int(n < 0 ? n : n + 1) : n };
MyMath.round = function(n) {
    return n < 0 ? MyMath.int(n % 1 < -0.5 ? n - 1 : n )
                 : MyMath.int(n % 1 <  0.5 ? n : n + 1 );
};

平方根

けれんみのないニュートン法



MyMath.sqrt = function(n){
    if (n < 0)   return 0/0; /* NaN */;
    if (n === 0) return 0;
    var s, t;
    s = n < 1 ? 1 : n;
    do {
        t = s; s = ( n / s + s ) / 2;
    } while (s < t);
    return s;
};

指数

そのままテイラー展開で。コメントアウトされている部分をアンコメントすると気持ち速くなるが精度も落ちる。



MyMath.exp = function(n) {
    if (isNaN(n)) return 0/0;
    var i, s = 0, e, a, ep;
    /*
    k = MyMath.int(n / MyMath.LN2 + (n < 0 ? -0.5 : 0.5));
    n -= k * MyMath.LN2;
    */
    if (n < 0) s = 1, n *= -1;
    for( e = 1 + n, a = n, i = 2, ep = -1; e !== ep; i++ ) {
         ep = e; a *= n / i; e += a;
    };
    /* return MyMath.ldexp(s ? 1 / e : e, k); */
    return s ? 1 / e : e
};

対数

級数展開をそのまま使うと遅いので「C言語によるアルゴリズム事典」の方法で。



MyMath.frexp = function(n){
    var s = 1, x = 0;
    if (n < 0) s = -1, n *= -1;
    if (n === 0) {
        s = 0;
    } else if (n < 0.5) {
        while (n < 0.5) n *= 2, x--;
    } else if (n >= 1) {
        while (n >= 1) n /= 2, x++;
    }
    return { s:s, m:n, x:x };
};

MyMath.ldexp = function(x, n) {
    var b = n >= 0 ? 2 : 1/2;
    if (n < 0) n *= -1;
    for (; n > 0 ; b *= b, n >>>= 1) if (n & 1) x *= b;
    return x;
};

MyMath.log = function(n){
    if (n <= 0) return 0/0;
    var i, k, n2, s, sp;
    k = MyMath.frexp(n / MyMath.SQRT2).x;
    n /= MyMath.ldexp(1, k);
    n = (n - 1) / (n + 1); n2 = n * n; i = 1; s = n;
    do {
        n *= n2; i += 2; sp = s; s += n / i;
    } while (sp != s);
    return MyMath.LN2 * k + 2 * s;
};

ベキ乗

nが整数ならそのままn回かけ、そうでなければ en log bを。ただしn回乗算する代わりに、404 Blog Not Found:アルゴリズム - 同じ文字列のn回繰り返しをlog n回で作る方法と同じ手口で乗算回数はlog nに抑えてある。



MyMath.pow = function(b, n){
    if (n === MyMath.int(n)){
        var r = 1, s = 0;
        if (n < 0) n *= -1, s = 1;
        for (n *= 1; n > 0; n >>>= 1, b *= b) if (n & 1) r *= b;
        return s ? 1 / r : r;
    }else{
        return MyMath.exp(MyMath.log(b) * n);
    }
};

三角関数

cosのみきちんとかいてあとはそこから導出。%の使い方に注目。見ての通りπ/2の倍数は特別扱いしているので、tan(π/2)がちゃんとInfinityになる。あと、誤差を抑えるため、計算する範囲を(0, π/2)に留めてある。



MyMath.cos = function(n) {
    if (!isFinite(n)) return 0/0; /* to be safe */
    n %= 2 * Math.PI;
    if (n < 0) n *= -1;
    var r = MyMath.PI / 2;
    if (n % r === 0) return [1, 0, -1, 0][n / r];
    if (n > r) {
        return (n < 2 * r) ?  arguments.callee(2 * r - n)
            :  (n < 3 * r) ? -arguments.callee(n - 2 * r)
                           : -arguments.callee(4 * r - n);
    }
    var i, c, cp, a;
    for (i = 2, c = a = 1, cp = 2; cp !== c; i += 2) {
        cp = c; a *= - n * n / (i - 1) / i; c += a;
    }
    return c;
};

MyMath.sin = function(n) { return MyMath.cos(n - Math.PI / 2) };
MyMath.tan = function(n) { return MyMath.sin(n) / MyMath.cos(n) };

逆三角関数

「C言語によるアルゴリズム事典」の連分数展開で。atan2()の仕様には+0と-0が分けて書いてあるのだけど、JavaScriptでは-0.0は現れないようだ。+0.0 === -0.0===でも区別できず、1/+0.0 !== 1/-0.0で区別しなければならない。



MyMath.__ATAN_ITER__ = 24; /* good enough for long double */
MyMath.atan = function(n) {
    if (!n) return 0;
    if (isNaN(n)) return n;
    if (!isFinite(n)) return MyMath.PI / 2 * (n < 0 ? -1 : 1);
    var i, s = 0, a = 0;
    if (n > 1) { s = 1; n = 1 / n }
    else if (n < -1) { s = -1; n = 1 / n }
    for (i = MyMath.__ATAN_ITER__; i >= 1; i--){
        a = (i * i * n * n) / (2 * i + 1 + a);
    }
    return s ? (s > 0) ?   MyMath.PI / 2 - n / (1 + a)
                       : - MyMath.PI / 2 - n / (1 + a)
             : n / (1 + a);
};
MyMath.atan2 = function(y, x){ 
    if (isFinite(y)) {
        if (isFinite(x)){
            if (x !== 0){
                if (y !== 0){
                    return y > 0 ? x > 0 ? MyMath.atan(y/x)
                                         : +MyMath.PI + MyMath.atan(y/x)
                                 : x < 0 ? -MyMath.PI + MyMath.atan(y/x)
                                         : MyMath.atan(y/x);
                } else {
                    return x < 0 ? MyMath.PI : 0;
                }
            } else {
                return y !== 0 ? y < 0 ? -MyMath.PI/2 : MyMath.PI/2
                               : 1/x === 1/0 ? 0 : MyMath.PI;
            }
        } else {
            return x < 0 ? MyMath.PI : 0;
        }
    } else {
        return (
            isFinite(x) ? y < 0 ? -1/2 : 1/2
                        : y < 0 ? x < 0 ? -3/4 : -1/4
                                : x < 0 ?  3/4 :  1/4
        ) * MyMath.PI;             
    }          
};                                 
MyMath.asin = function(n){ return MyMath.atan(n / MyMath.sqrt(1 - n * n)) };
MyMath.acos = function(n){ return MyMath.PI/2 - MyMath.asin(n) };

乱数

これのみ実装していない。代わりにここでは

の実装を借用する形にしてある。



/*
 * <script src="http://homepage2.nifty.com/magicant/sjavascript/mt.js"></script>
 */
MyMath.__MT__ = new MersenneTwister;
MyMath.random = function() { return MyMath.__MT__.next() };

まとめ

これで

delete Math;

されても安心!?

Dan the Wheel Reinventor


この記事へのトラックバックURL

この記事へのトラックバック
言うまでもなくdankogai氏の1ヶ月位前のブログに触発されて。 普段から使うMath系の関数もあるけど、中で何をやってるのか知らなかったから再発明してみた。 乱数とX進数からX進数へ変換以外の関数は全て。 四角い車輪。   極力よく見かける公式とかに合わせた。 速度は気
[アルゴリズム][PHP]PHPのMath関数を再発明してみた【eth0jpの日記】at 2010年10月04日 05:18