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


C言語による最新アルゴリズム事典
 奥村晴彦

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