16ビット演算を使用して級数を評価することによりπを計算するときにオーバーフローを回避しますか?

15
比尔盖子 2019-05-08 02:25.

πから1000桁以上の10進数を計算するプログラムを書こうとしています。

楽しみのために低レベルのプログラミングを練習するために、最終的なプログラムは、乗算または除算がなく、16ビットの加算のみを実行する8ビットCPU上でアセンブリで記述されます。実装を容易にするために、16ビットの符号なし整数演算のみを使用でき、反復アルゴリズムを使用できることが望ましいです。速度は大きな問題ではありません。また、高速の乗算と除算はこの質問の範囲を超えているため、これらの問題も考慮しないでください。

アセンブリに実装する前に、デスクトップコンピューターのCで使用可能なアルゴリズムを見つけようとしています。これまでのところ、次のシリーズはかなり効率的で、比較的簡単に実装できることがわかりました。

この式は、収束加速法を使用してライプニッツ級数から導出されます。導出するには、Carl D. Offnerによる「πでの桁の計算」を参照してください(https://cs.umb.edu/~offner/files/pi.pdf)、19-26ページ。最終的な数式を26ページに示します。私が書いた最初の数式にはいくつかのタイプミスがありました。ページを更新して、固定された数式を確認してください。2最大項の定数項については、54ページで説明しています。この論文では、高度な反復アルゴリズムについても説明していますが、ここでは使用しませんでした。

多くの(たとえば5000)項を使用して級数を評価すると、数千桁のπを簡単に取得できます。このアルゴリズムを使用しても、この級数を繰り返し評価するのは簡単です。

アルゴリズム

  1. まず、数式を並べ替えて、配列から定数項を取得します。

  1. 配列に2を入力して最初の反復を開始します。したがって、新しい数式は元の数式に似ています。

  2. しましょうcarry = 0

  3. 最大の用語から始めます。配列から1つの項(2)を取得し、項にPRECISIONを掛けて、に対して固定小数点除算を実行し2 * i + 1、リマインダーを新しい項として配列に保存します。次に、次の用語を追加します。ここでデクリメントiし、次の項に進み、まで繰り返しi == 1ます。最後に最後の項を追加しx_0ます。

  4. 16ビット整数が使用されているため、PRECISION10、したがって10進数の2桁が取得されますが、有効なのは最初の桁のみです。2桁目をキャリーとして保存します。最初の桁とキャリーを表示します。

  5. x_0 は整数2であり、連続する反復では追加しないでください。クリアしてください。

  6. 手順4に進み、必要な桁がすべて揃うまで、次の10進数を計算します。

実装1

このアルゴリズムをCに変換する:

#include <stdio.h>
#include <stdint.h>

#define N 2160
#define PRECISION 10

uint16_t terms[N + 1] = {0};

int main(void)
{
    /* initialize the initial terms */
    for (size_t i = 0; i < N + 1; i++) {
        terms[i] = 2;
    }

    uint16_t carry = 0;
    for (size_t j = 0; j < N / 4; j++) {
        uint16_t numerator = 0;
        uint16_t denominator;
        uint16_t digit;

        for (size_t i = N; i > 0; i--) {
            numerator += terms[i] * PRECISION;
            denominator = 2 * i + 1;

            terms[i] = numerator % denominator;
            numerator /= denominator;
            numerator *= i;
        }
        numerator += terms[0] * PRECISION;
        digit = numerator / PRECISION + carry;
        carry = numerator % PRECISION;

        printf("%01u", digit);

        /* constant term 2, only needed for the first iteration. */
        terms[0] = 0;
    }
    putchar('\n');
}

コードは、エラーが発生するまで、πから31桁の10進数を計算できます。

31415926535897932384626433832794
10 <-- wrong

時には、digit + carry9より大きいので、余分なキャリーが必要です。運が悪ければ、ダブルキャリーやトリプルキャリーなどが発生する可能性があります。リングバッファを使用して最後の4桁を格納します。余分なキャリーが検出された場合は、バックスペースを出力して前の桁を消去し、キャリーを実行して再印刷します。これは、概念実証の醜い解決策であり、オーバーフローに関する私の質問とは無関係ですが、完全を期すために、ここにあります。将来的にはもっと良いものが実装されるでしょう。

キャリーを繰り返す実装2

#include <stdio.h>
#include <stdint.h>

#define N 2160
#define PRECISION 10

#define BUF_SIZE 4

uint16_t terms[N + 1] = {0};

int main(void)
{
    /* initialize the initial terms */
    for (size_t i = 0; i < N + 1; i++) {
        terms[i] = 2;
    }

    uint16_t carry = 0;
    uint16_t digit[BUF_SIZE];
    int8_t idx = 0;

    for (size_t j = 0; j < N / 4; j++) {
        uint16_t numerator = 0;
        uint16_t denominator;

        for (size_t i = N; i > 0; i--) {
            numerator += terms[i] * PRECISION;
            denominator = 2 * i + 1;

            terms[i] = numerator % denominator;
            numerator /= denominator;
            numerator *= i;
        }
        numerator += terms[0] * PRECISION;
        digit[idx] = numerator / PRECISION + carry;

        /* over 9, needs at least one carry op. */
        if (digit[idx] > 9) {
            for (int i = 1; i <= 4; i++) {
                if (i > 3) {
                    /* allow up to 3 consecutive carry ops */
                    fprintf(stderr, "ERROR: too many carry ops!\n");
                    return 1;
                }
                /* erase a digit */
                putchar('\b');

                /* carry */
                digit[idx] -= 10;
                idx--;
                if (idx < 0) {
                    idx = BUF_SIZE - 1;
                }
                digit[idx]++;            
                if (digit[idx] < 10) {
                    /* done! reprint the digits */
                    for (int j = 0; j <= i; j++) {
                        printf("%01u", digit[idx]);
                        idx++;
                        if (idx > BUF_SIZE - 1) {
                            idx = 0;
                        }
                    }
                    break;
                }
            }
        }
        else {
            printf("%01u", digit[idx]);
        }

        carry = numerator % PRECISION;
        terms[0] = 0;

        /* put an element to the ring buffer */
        idx++;
        if (idx > BUF_SIZE - 1) {
            idx = 0;
        }
    }
    putchar('\n');
}

これで、エラーが発生するまで、プログラムは534桁のπを正しく計算できるようになりました。

3141592653589793238462643383279502884
1971693993751058209749445923078164062
8620899862803482534211706798214808651
3282306647093844609550582231725359408
1284811174502841027019385211055596446
2294895493038196442881097566593344612
8475648233786783165271201909145648566
9234603486104543266482133936072602491
4127372458700660631558817488152092096
2829254091715364367892590360011330530
5488204665213841469519415116094330572
7036575959195309218611738193261179310
5118548074462379962749567351885752724
8912279381830119491298336733624406566
43086021394946395
22421 <-- wrong

16ビット整数オーバーフロー

最初の最大の項の計算中に、最初の除数が約4000の範囲にあるため、誤差項が非常に大きくなることがわかります。系列を評価するとき、numerator実際にはすぐに乗算でオーバーフローし始めます。

最初の500桁を計算する場合、整数オーバーフローは重要ではありませんが、誤った結果が得られるまで、ますます悪化し始めます。

に変更uint16_t numerator = 0するuint32_t numerator = 0と、この問題を解決し、πから1000桁以上を計算できます。

ただし、前述したように、ターゲットプラットフォームは8ビットCPUであり、16ビット操作しかありません。1つ以上のuint16_tのみを使用して、ここで見られる16ビット整数オーバーフローの問題を解決するためのトリックはありますか?多倍長演算を回避できない場合、ここでそれを実装する最も簡単な方法は何ですか?どういうわけか、16ビットの「拡張ワード」を追加する必要があることはわかっていますが、どうすれば実装できるかわかりません。

そして、ここでの長い文脈を理解するためにあなたの忍耐に前もって感謝します。

3 answers

2
Spektre 2019-05-08 21:07.

関連するQAをご覧ください。

  • ベーキングパイチャレンジ-理解と改善

Wikiを使用しています:Bailey–Borwein–Plouffe_formulaは、整数演算により適しています。

ただし、実際の課題は次のとおりです。

  • 非常に長い2進数を10進数に変換するにはどうすればよいですか?。

おそらく12月ベースで数値を印刷したいので...

また、asmよりも高水準言語で運ぶ必要がある場合は、これを見てください。

  • キャリーを通じて価値を広めることができない

必要な数のキャリービットを処理するように変更できます(データ型のビット幅よりも小さい場合)。

[編集1] C ++ / VCLでのBBPの例

私はこの式を使用しました(上記のリンク先のWikiページから取得):

固定小数点に変換...

//---------------------------------------------------------------------------
AnsiString str_hex2dec(const AnsiString &hex)
    {
    char c;
    AnsiString dec="",s;
    int i,j,l,ll,cy,val;
    int  i0,i1,i2,i3,sig;
    sig=+1; l=hex.Length();
    if (l) { c=hex[l]; if (c=='h') l--; if (c=='H') l--; }
    i0=0; i1=l; i2=0; i3=l;
    for (i=1;i<=l;i++)      // scan for parts of number
        {
        char c=hex[i];
        if (c=='-') sig=-sig;
        if ((c=='.')||(c==',')) i1=i-1;
        if ((c>='0')&&(c<='9')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        if ((c>='A')&&(c<='F')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        if ((c>='a')&&(c<='f')) { if (!i0) i0=i; if ((!i2)&&(i>i1)) i2=i; }
        }

    l=0; s=""; if (i0) for (i=i0;i<=i1;i++)
        {
        c=hex[i];
             if ((c>='0')&&(c<='9')) c-='0';
        else if ((c>='A')&&(c<='F')) c-='A'-10;
        else if ((c>='a')&&(c<='f')) c-='A'-10;
        for (cy=c,j=1;j<=l;j++)
            {
            val=(s[j]<<4)+cy;
            s[j]=val%10;
            cy  =val/10;
            }
        while (cy>0)
            {
            l++;
            s+=char(cy%10);
            cy/=10;
            }
        }
    if (s!="")
        {
        for (j=1;j<=l;j++) { c=s[j]; if (c<10) c+='0'; else c+='A'-10; s[j]=c; }
        for (i=l,j=1;j<i;j++,i--) { c=s[i]; s[i]=s[j]; s[j]=c; }
        dec+=s;
        }
    if (dec=="") dec="0";
    if (sig<0) dec="-"+dec;

    if (i2)
        {
        dec+='.';
        s=hex.SubString(i2,i3-i2+1);
        l=s.Length();
        for (i=1;i<=l;i++)
            {
            c=s[i];
                 if ((c>='0')&&(c<='9')) c-='0';
            else if ((c>='A')&&(c<='F')) c-='A'-10;
            else if ((c>='a')&&(c<='f')) c-='A'-10;
            s[i]=c;
            }
        ll=((l*1234)>>10);  // num of decimals to compute
        for (cy=0,i=1;i<=ll;i++)
            {
            for (cy=0,j=l;j>=1;j--)
                {
                val=s[j];
                val*=10;
                val+=cy;
                s[j]=val&15;
                cy=val>>4;
                }
            dec+=char(cy+'0');
            for (;;)
                {
                if (!l) break;;
                if (s[l]) break;
                l--;
                }
            if (!l) break;;
            }
        }

    return dec;
    }
//---------------------------------------------------------------------------
AnsiString pi_BBP() // https://en.wikipedia.org/wiki/Bailey–Borwein–Plouffe_formula
    {
    const int N=100;        // 32*N bit uint arithmetics
    int sh;
    AnsiString s;
    uint<N> pi,a,b,k,k2,k3,k4;

    for (pi=0,sh=(N<<5)-8,k=0;sh>=0;k++,sh-=4)
        {
        k2=k*k;
        k3=k2*k;
        k4=k3*k;
        a =k2* 120;
        a+=k * 151;
        a+=     47;
        b =k4* 512;
        b+=k3*1024;
        b+=k2* 712;
        b+=k * 194;
        b+=     15;
        a<<=sh;
        pi+=a/b;
        }
    pi<<=4;
    s=pi.strhex();
    s=s.Insert(".",2);
    return str_hex2dec(s);
    }
//---------------------------------------------------------------------------

コードは、自己割り当て文字列であるVCL と、鉱山ALU32に基づくビット幅の符号なし整数演算であるAnsiString鉱山uint<N>テンプレートを使用しています。ご覧のとおり、これには大きな整数の除算の加算と乗算のみが必要です(他のすべての処理は通常の整数で実行できます)。32*N

ここでは、10年ごとの結果と1000桁の円周率の参照を比較しています。

ref: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989
BPP: 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778048187

計算されたbigint値は16進文字列にエクスポートされstr_hex2dec、上記のリンクを使用して10進ベースに変換されます。反復回数は、ターゲットのビット幅によって異なります。

コードはまだ最適化されていません...

1
Yves Daoust 2019-05-08 03:39.

32ビット演算の実装はどうですか?

さらに、上位2ワード(16ビット)を追加し、次に下位2ワードを追加し、オーバーフロービットをテストし、必要に応じて上位結果にキャリーします。

オーバーフローがいつ発生するかを予測できる場合は、必要に応じて16ビットから32ビットの演算に切り替えることができます。


オーバーフロービットのテストは純粋なCでは実行できません。インラインアセンブリまたは組み込み関数が必要になります。

そうでなければ、あなたはこの答えに触発されることができます: https://codereview.stackexchange.com/a/37178/39646

0
alx 2019-05-08 03:03.

トリックがあります:

分子には配列を使用し、分母には別の配列を使用することを検討してください。各位置は、実際の数を取得するためにその数が乗算される回数を表します。

例:

(1 * 2 * 3 * 7 * 7) / (3 * 6 * 8)

次のように表されます:

num[] = {1, 1, 1, 0, 0, 0, 2};
denom[] = {0, 0, 1, 0, 0, 1, 0, 1};

次に、数値を格納する前に、数値ごとに素数に因数分解することを検討してください。そうすれば、数値が小さくなります。次に、すべての素数を格納するための別の配列が必要になります。

primes[] = {2, 3, 5, 7};
num[] = {1, 1, 0, 2};
denom[] = {4, 2, 0, 0};

これにより、想像を絶するほど大きな数を格納できますが、遅かれ早かれそれらを数に戻したいので、最初にこれを単純化する必要があります。これを行う方法factors[i] += num[i] - denom[i]は、配列のすべてのフィールド、シリーズのすべての分数を減算することです。オーバーフローのリスクを最小限に抑えるために、各反復後に単純化する必要があります。

factors[] = {-3, -1, 0, 2};

数値が必要なnum *= pow(primes[i], factors[i]);場合はnum /= pow(primes, -factors[i]);、配列のすべてのフィールドに対して、係数が正の場合、または負の場合に実行します。(0の場合は何もしません。


numdenomは分数を格納するために使用される一時的な配列であり、結果が格納される配列はfactorsです。memset使用する前に、一時配列を覚えておいてください。


この説明は、大きな部分に役立ちます。特定の問題に適応させるには、整数のべき関数を使用し、10 ^何かを掛けて小数部分を整数部分に変換する必要がある場合があります。それはあなたの使命です、あなたがそれを受け入れるべきですか:)

Related questions

MORE COOL STUFF

Reba McEntire は、彼女が息子の Shelby Blackstock と共有する「楽しい」クリスマスの伝統を明らかにしました:「私たちはたくさん笑います」

Reba McEntire は、彼女が息子の Shelby Blackstock と共有する「楽しい」クリスマスの伝統を明らかにしました:「私たちはたくさん笑います」

Reba McEntire が息子の Shelby Blackstock と共有しているクリスマスの伝統について学びましょう。

メーガン・マークルは、自然な髪のスタイリングをめぐってマライア・キャリーと結ばれました

メーガン・マークルは、自然な髪のスタイリングをめぐってマライア・キャリーと結ばれました

メーガン・マークルとマライア・キャリーが自然な髪の上でどのように結合したかについて、メーガンの「アーキタイプ」ポッドキャストのエピソードで学びましょう.

ハリー王子は家族との関係を修復できるという「希望を持っている」:「彼は父親と兄弟を愛している」

ハリー王子は家族との関係を修復できるという「希望を持っている」:「彼は父親と兄弟を愛している」

ハリー王子が家族、特にチャールズ王とウィリアム王子との関係について望んでいると主張したある情報源を発見してください。

ワイノナ・ジャッドは、パニックに陥った休暇の瞬間に、彼女がジャッド家の家長であることを認識しました

ワイノナ・ジャッドは、パニックに陥った休暇の瞬間に、彼女がジャッド家の家長であることを認識しました

ワイノナ・ジャッドが、母親のナオミ・ジャッドが亡くなってから初めての感謝祭のお祝いを主催しているときに、彼女が今では家長であることをどのように認識したかを学びましょう.

セントヘレナのジェイコブのはしごを登るのは、気弱な人向けではありません

セントヘレナのジェイコブのはしごを登るのは、気弱な人向けではありません

セント ヘレナ島のジェイコブズ ラダーは 699 段の真っ直ぐ上る階段で、頂上に到達すると証明書が発行されるほどの難易度です。

The Secrets of Airline Travel Quiz

The Secrets of Airline Travel Quiz

Air travel is far more than getting from point A to point B safely. How much do you know about the million little details that go into flying on airplanes?

Where in the World Are You? Take our GeoGuesser Quiz

Where in the World Are You? Take our GeoGuesser Quiz

The world is a huge place, yet some GeoGuessr players know locations in mere seconds. Are you one of GeoGuessr's gifted elite? Take our quiz to find out!

バイオニック読書はあなたをより速く読むことができますか?

バイオニック読書はあなたをより速く読むことができますか?

BionicReadingアプリの人気が爆発的に高まっています。しかし、それは本当にあなたを速読術にすることができますか?

作家のアンバー・ラフィンとジェニー・ヘーゲルが上司のセス・マイヤーズを引き継ぐのを見る

作家のアンバー・ラフィンとジェニー・ヘーゲルが上司のセス・マイヤーズを引き継ぐのを見る

深夜のアンバー・ラフィンとジェニー・ヘーゲルが繰り返し「ジョーク・セス・カント・テル」で戻ってきました。これらのジョークの多くは観客をがっかりさせますが、最初から最後まで素晴らしいです。ルフィンとヘーゲルは黒人女性として自己紹介します。とゲイの女性、それぞれ、したがって、セスマイヤーズが10フィートのポールで触れることができない主題について賢明にクラックすることができます。

ジョンウィック:第3章は2019年5月に劇場への道を容赦なく殺します

ジョンウィック:第3章は2019年5月に劇場への道を容赦なく殺します

(写真:ライオンズゲート)この「キアヌ・リーブスはダッパースーツを着て人々を殺害する」というモチーフ全体が手元にあることをはっきりと知っているライオンズゲートは、スタイリッシュで復讐に燃えるジョン・ウィックのフランチャイズで3回目のリリース日を設定しました。犬をベースにした復讐のためのババ・ヤガの果てしない十字軍を支えるバットシット神話をより深く掘り下げることを約束する3番目のジョン・ウィック映画は、2019年5月17日に設定されました。これまでのところ、それはその日に上陸した唯一の映画です。

このパイロットは、This IsUsの残りの部分に高い基準を設定します

このパイロットは、This IsUsの残りの部分に高い基準を設定します

写真:NBCパイロットは良すぎるのでしょうか?ありそうもないようですが、This IsUsのファンの場合はそうかもしれません。クレイジー、バカ、ラブライターのダン・フォーゲルマンからの待望の新シリーズは、ツイストエンディングを中心に展開しています。シリーズを適切に設定しますが、非常に巧妙に行われているため、改善の余地はあまりありません。

ああ、GIFがついにFacebookで機能する

ああ、GIFがついにFacebookで機能する

ここにいくつかのニュースがあります:あなたは今FacebookにGIFを埋め込むことができます。まあ、技術的には、GIFへのリンクを投稿することができ、Facebookは、他のほとんどすべてのソーシャルネットワークが何年も行ってきたようにアニメーションを作成します。

米国のフィギュア スケートは、チーム イベントでの最終決定の欠如に「苛立ち」、公正な裁定を求める

米国のフィギュア スケートは、チーム イベントでの最終決定の欠如に「苛立ち」、公正な裁定を求める

ロシアのフィギュアスケーター、カミラ・バリエバが関与したドーピング事件が整理されているため、チームは2022年北京冬季オリンピックで獲得したメダルを待っています。

Amazonの買い物客は、わずか10ドルのシルクの枕カバーのおかげで、「甘やかされた赤ちゃんのように」眠れると言っています

Amazonの買い物客は、わずか10ドルのシルクの枕カバーのおかげで、「甘やかされた赤ちゃんのように」眠れると言っています

何千人ものAmazonの買い物客がMulberry Silk Pillowcaseを推奨しており、現在販売中. シルクの枕カバーにはいくつかの色があり、髪を柔らかく肌を透明に保ちます。Amazonで最大46%オフになっている間にシルクの枕カバーを購入してください

パデュー大学の教授が覚醒剤を扱った疑いで逮捕され、女性に性的好意を抱かせる

パデュー大学の教授が覚醒剤を扱った疑いで逮捕され、女性に性的好意を抱かせる

ラファイエット警察署は、「不審な男性が女性に近づいた」という複数の苦情を受けて、12 月にパデュー大学の教授の捜査を開始しました。

コンセプト ドリフト: AI にとって世界の変化は速すぎる

コンセプト ドリフト: AI にとって世界の変化は速すぎる

私たちの周りの世界と同じように、言語は常に変化しています。以前の時代では、言語の変化は数年または数十年にわたって発生していましたが、現在では数日または数時間で変化する可能性があります。

SF攻撃で91歳のアジア人女性が殴られ、コンクリートに叩きつけられた

犯罪擁護派のオークランドが暴力犯罪者のロミオ・ロレンゾ・パーハムを釈放

SF攻撃で91歳のアジア人女性が殴られ、コンクリートに叩きつけられた

認知症を患っている 91 歳のアジア人女性が最近、47 番街のアウター サンセット地区でロメオ ロレンゾ パーハムに襲われました。伝えられるところによると、被害者はサンフランシスコの通りを歩いていたところ、容疑者に近づき、攻撃を受け、暴行を受けました。

Precios accesibles, nuestro aprendizaje desde la perspectiva iOS

Precios accesibles, nuestro aprendizaje desde la perspectiva iOS

Cómo mejoramos la accesibilidad de nuestro componente de precio, y cómo nos marcó el camino hacia nuevos saberes para nuestro sistema de diseño. Por Ana Calderon y Laura Sarmiento Leer esta historia en inglés.

ℝ

“And a river went out of Eden to water the garden, and from thence it was parted and became into four heads” Genesis 2:10. ? The heart is located in the middle of the thoracic cavity, pointing eastward.

Language