#56　任意精度で√2の近似値を二分法で計算してみた

概要

この記事は前回の続きです

前回、二分法で√2の近似値を計算してみましたが、倍精度浮動小数の精度の限界で16桁程度までしか計算できませんでしたので、今回は任意精度での計算をします

任意精度の実現

JavaScript/TypeScriptにはBigIntという任意精度の整数を表現する方法があり、BigIntを利用して実数を表現します

具体的には、仮数部と指数部をBigIntで表現し、仮数部 * 2 ^ 指数部 = 表現したい値とすることで表現します

※この方法の場合、循環小数や無理数などは厳密に表現できませんが今回のテーマでは登場しないので考慮しないものとします

実装

bigNumber(小数を表現するクラス)

任意精度の実数を表現するためbigNumberクラスを作成します

• constructor(f: bigint, e: bigint = 0n): コンストラクタ。仮数部は必須で、指数部はオプション(デフォルト値は0)です
• normalize(): 正規化を行います。仮数部の最下位ビットが1になるよう仮数部と指数部を調整します
  • 仮数部の最下位ビットが0の場合は偶数なので2で除算(右シフト)し、指数部の値を1小さくします。この操作を繰り返します
• toString(): 文字列としてを出力します

class bigNumber {
    e = BigInt(0);
    f = BigInt(0);

    /**
     * f * (2 ^ e)
     * @param bigint 仮数部 (require)
     * @param bigint 指数部 (option)
     */
    constructor(f: bigint, e: bigint = 0n) {
        this.f = f;
        this.e = e;
    }

    /**
     * 正規化を行う
     * @returns bigNumber
     */
    normalize() {
        while (!(this.f & 0x01n)) {
            this.f >>= 1n;
            this.e++;
        }

        return this;
    }

    /**
     * @returns string
     */
    toString(): string {
        this.normalize();
        
        return `${this.f} * (2 ^ ${this.e})`;
    }
}

BigNumber(演算を行うメソッド群)

• BigNumber.add(operand1: bigNumber, operand2: bigNumber): 引数として渡された数値を加算します
  • 指数部をそろえて仮数部同士の加算を行い、最後に正規化を行います
• BigNumber.mul(operand1: bigNumber, operand2: bigNumber): 引数として渡された数値を乗算します
  • 仮数部同士を乗算し指数部同士を加算します。最後に正規化を行います
• BigNumber.sub(operand1: bigNumber, operand2: bigNumber): 1つ目の引数から2つ目の引数を減算します

namespace BigNumber {
    //足し算
    export function add(operand1: bigNumber, operand2: bigNumber): bigNumber {
        let o1 = new bigNumber(operand1.f, operand1.e);
        let o2 = new bigNumber(operand2.f, operand2.e);

        //指数部を小さいほうに合わせる
        if (o1.e > o2.e) {
            //(f1 * (2 ^ (e1 - e2)) + f2) * (2 ^ e2)
            o1.f <<= (o1.e - o2.e);
            o1.e = o2.e;
        }
        else {
            //(f1 + f2 * (2 ^ (e2 - e1))) * (2 ^ e1)
            o2.f <<= (o2.e - o1.e);
            o2.e = o1.e;
        }

        //仮数部の足し算
        o1.f += o2.f;

        //正規化
        o1.normalize();

        return o1;
    }

    //引き算
    export function sub(operand1: bigNumber, operand2: bigNumber): bigNumber {
        return BigNumber.add(BigNumber.mul(operand1, new bigNumber(-1n, 0n)), operand2);
    }

    //掛け算
    export function mul(operand1: bigNumber, operand2: bigNumber): bigNumber {
        let o1 = new bigNumber(operand1.f, operand1.e);
        let o2 = new bigNumber(operand2.f, operand2.e);

        //(f1 * f2) * (2 ^ (e1 + e2))
        o1.f *= o2.f;
        o1.e += o2.e;

        //正規化
        o1.normalize();

        return o1;
    }
}

前回作成したコードに組み込む

今回は1000回ループを回して√2の近似値を計算してみます

let targetNumber = new bigNumber(2n); // = 2
let min = new bigNumber(1n); // = 1
let max = new bigNumber(2n); // = 2

let half = new bigNumber(1n, -1n); // = 0.5

for (let i = 0; i < 1000; i++) {
    //最小値と最大値の中間を求める
    let sqrt = BigNumber.mul(BigNumber.add(min, max), half); //(min + max) * 0.5

    //√2の候補を2乗する
    let t2 = BigNumber.mul(sqrt, sqrt); //sqrt * sqrt

    //もし中間値が目的値より小さいなら範囲の最小値を更新、そうでないなら最大値を更新
    if (BigNumber.sub(targetNumber, t2).f < 0) {
        min = sqrt;
    }
    else {
        max = sqrt;
    }

    count++;
}
console.log(`min: ${min.toString()}\nmax: ${max.toString()}`);
console.log(`count: ${count}`);

実行・結果

npx ts-node ./src/main.tsを実行します

$ npx ts-node ./src/main.ts
min: 473544376400726394228831039451913480123688799751094630616578363333661533455686216226143010758906440540743599383075282477937325116004184859362870439857750473316274124854227871716121348224479943162567953976022762080355598781683844459583873884188875089875732987543720916348829425487472888348035796043929 * (2 ^ -995)
max: 15153420044823244615322593262461231363958041592035028179730507626677169070581958919236576344285006097303795180258409039293994403712133915499611854075448015146120771995335291894915883143183358181202174527232728386571379161013883022706683964294044002876023455601399069323162541615599132427137145473405729 * (2 ^ -1000)
count: 1000

計算結果は以下の通りになり、およそ300桁程度の精度まで計算することができました

min: 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140798964…
max: 1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851741864088919860955232923048430871432145083976260362799525140798973…

まとめ

BitIntを用いることで倍精度浮動小数よりも大きな桁数の実数が扱えるようになり、√2の近似値もより精度良く計算することができました