円周率の計算 (9)
8 月 20 日の記事で検討した spigot アルゴリズムの変更点を、8 月 15 日の記事の Arduino 用のスケッチに適用し、乗除算部分をインライン・アセンブラによる 8 ビット演算に置き換えたプログラムを下に示します。 (Spigot 関数の部分のみ)
全体のスケッチのコンパイル後のバイナリ・スケッチ・サイズは 5574 バイトです。
void Spigot( mp8_t num ) // MP work { const uint8_t base = 100; // n-ary number const uint8_t nterms = 128; // number of terms const uint8_t n_step = 4; // calculation interval register uint8_t i; // loop index register uint8_t n; // loop index register uint16_t temp; // temporary work / carry register uint8_t out; // output digits / carry adjust register uint8_t denom; // denominator of series #if defined(ARDUINO) register uint16_t temp2; #endif out = 2; for( i = 0 ; i < nterms ; i++ ) { num[i] = out; // initalize array } PRINT_STR ("Spigot (nterms = "); PRINT_DEC (nterms); PRINTLN_STR(")"); PRINT_STR ("pi = "); for( n = nterms ; n >= (nterms-2*Ndigits_Spigot); n -= n_step ) { temp = 0; for( i = n - 1 ; i >= 2 ; i-- ) { denom = 2 * i - 1; #if defined(ARDUINO) // product, multiplicand, lier, carry mul_add_u8(temp2, (uint8_t)temp, i, 0); mul_add_u8(temp, num[i], base, 0); temp2 += temp; temp = div_u16_u8(temp2, denom); // temp = rem8:quo8 num[i] = (temp >> 8); // remainder temp &= 0xff; // quotient #else temp = temp * i + num[i] * base; num[i] = temp % denom; temp /= denom; #endif } // for (i ... #if defined(ARDUINO) temp = div_u16_u8(temp, base); // temp = rem8:quo8 i = out + (0xff & temp); // quitient out = (temp >> 8); // remainder temp = div_u16_u8(i, 10); // temp = rem8:quo8 if (nterms != n) { PRINT_DEC(0xff & temp); // i / 10 } // if (nterms != n) ... PRINT_DEC(temp >> 8); // i % 10 if (0 == (div_u16_u8((nterms-n),(n_step * 5)))>>8) { PRINT_STR(((nterms) == n) ? ".\r\xa" : " "); } // if (0 == ... #else i = (out + temp / base); if (nterms != n) { PRINT_DEC(i / 10); } // if (nterms != n) ... PRINT_DEC(i % 10); if (0 == ((nterms-n) % (n_step * 5))) { PRINT_STR(((nterms) == n) ? ".\r\xa" : " "); } // if (0 == ... out = temp % base; #endif } // for (n ... PRINTLN_STR(""); PRINTLN_STR(""); } // void Spigot()
Spigot アルゴリズムの分の出力結果を下に示します。
整数部が 100 進法の「31」から「3」になったため、小数点以下の有効桁数が 38 から 37 に減りました。
Spigot (nterms = 128) pi = 3. 1415926535 8979323846 2643383279 5028841624