多倍長整数の計算（5）

自乗の計算は、異なる積よりも速い。

f(x) = a₀ + a₁x + a₀x² + ...
f(x)² = a₀² + 2a₀a₁x + ( 2a₀a₂ + a₁²)x² + ...

ほとんどの場合で同じ乗算が2回行われる。
自乗にすると、8秒と1/2になった。

function square(x) {

    var x_length = x.length;

    var z_length = x_length * 2;

    var z = [ ];    // 返り値

    for(var i = 0; i < z_length; i++)

        z[i] = 0;

    

    for(var m = 0; m < z_length - 1; m++) {

        var z_m = 0;

        var x0, x1;     // z[m]を計算するときのxのindexの範囲

        if(m < x_length)

            x0 = 0;

        else

            x0 = m - x_length + 1;

        if(m & 1) {     // odd

            x1 = m >> 1;        // (m-1)/2

        }

        else {          // even

            x1 = (m >> 1) - 1;  // m/2-1

            

            var x_m2 = x[m>>1];

        }

        

        for(var k = x0; k <= x1; k++) {

            z_m += x[k] * x[m-k];

            if(z_m >= COUNT_MAX) {

                var mod = z_m % INT_MAX;

                var div = (z_m - mod) / INT_MAX;

                z_m = mod;

                var z_n;

                for(var n = m + 1; div != 0; n++) {

                    z_n = z[n] + div;

                    mod = z_n % INT_MAX;

                    div = (z_n - mod) / INT_MAX;

                    z[n] = mod;

                }

            }

        }

        

        z_m *= 2;

        if((m & 1) == 0)        // even

            z_m += x_m2 * x_m2;

        z_m += z[m];

        

        var n = m;

        while(true) {

            var mod = z_m % INT_MAX;

            var div = (z_m - mod) / INT_MAX;

            z[n] = mod;

            if(div == 0)

                break;

            z_m = z[n+1] + div;

            n++;

        }

    }

    

    if(z[z_length-1] == 0)

        z.length--;

    

    return z;

}