peter revised this gist 3 hours ago. Go to revision
1 file changed, 157 insertions
pi_spigot.cpp(file created)
| @@ -0,0 +1,157 @@ | |||
| 1 | + | // pi_spigot.cpp — compute N decimal digits of pi using the | |
| 2 | + | // Rabinowitz–Wagon spigot algorithm. | |
| 3 | + | // | |
| 4 | + | // "A Spigot Algorithm for the Digits of Pi" | |
| 5 | + | // Stanley Rabinowitz and Stan Wagon, 1995. | |
| 6 | + | // | |
| 7 | + | // The algorithm exploits the series | |
| 8 | + | // | |
| 9 | + | // pi = sum_{i=0..inf} (i! * 2^(i+1)) / (2i+1)!! | |
| 10 | + | // | |
| 11 | + | // rewritten in a mixed-radix form so that decimal digits can be | |
| 12 | + | // extracted one at a time using only integer arithmetic on a fixed | |
| 13 | + | // array of length L = floor(10N/3) + 1. No big-integer library is | |
| 14 | + | // required and every intermediate value fits comfortably in a 64-bit | |
| 15 | + | // integer for the problem sizes a single machine cares about. | |
| 16 | + | // | |
| 17 | + | // Output is "3." followed by N digits after the decimal point, and a | |
| 18 | + | // timing line on stderr. | |
| 19 | + | // | |
| 20 | + | // Build: c++ -std=c++20 -O2 -o pi_spigot pi_spigot.cpp | |
| 21 | + | // Run: ./pi_spigot 1000 | |
| 22 | + | ||
| 23 | + | #include <charconv> | |
| 24 | + | #include <chrono> | |
| 25 | + | #include <cstdint> | |
| 26 | + | #include <cstdlib> | |
| 27 | + | #include <iostream> | |
| 28 | + | #include <string> | |
| 29 | + | #include <string_view> | |
| 30 | + | #include <vector> | |
| 31 | + | ||
| 32 | + | namespace { | |
| 33 | + | ||
| 34 | + | // Parse a non-negative integer from argv. Returns nullopt-style sentinel | |
| 35 | + | // (-1) on failure so main() can print a usage message. | |
| 36 | + | long parse_count(std::string_view s) { | |
| 37 | + | long value = 0; | |
| 38 | + | auto [ptr, ec] = std::from_chars(s.data(), s.data() + s.size(), value); | |
| 39 | + | if (ec != std::errc{} || ptr != s.data() + s.size() || value < 0) { | |
| 40 | + | return -1; | |
| 41 | + | } | |
| 42 | + | return value; | |
| 43 | + | } | |
| 44 | + | ||
| 45 | + | // Compute `digits` decimal digits of pi using the Rabinowitz–Wagon | |
| 46 | + | // spigot. The leading "3" is the first digit emitted; the remaining | |
| 47 | + | // digits-1 are the fractional part. The result is returned as a plain | |
| 48 | + | // decimal string of length `digits`. | |
| 49 | + | // | |
| 50 | + | // Carry handling: each generated digit `q` may be 0..10. A value of 10 | |
| 51 | + | // means the previously buffered digit must be incremented (and any run | |
| 52 | + | // of intervening 9s flipped to 0s) before being flushed — the classic | |
| 53 | + | // "carry across a streak of nines" trick that makes the spigot work. | |
| 54 | + | std::string compute_pi(long digits) { | |
| 55 | + | if (digits <= 0) return {}; | |
| 56 | + | ||
| 57 | + | // Array length comes from the convergence bound in the paper: we | |
| 58 | + | // need at least ceil(10N/3) terms to be sure of the first N digits. | |
| 59 | + | const long len = (10 * digits) / 3 + 1; | |
| 60 | + | std::vector<std::int64_t> a(len, 2); | |
| 61 | + | ||
| 62 | + | std::string out; | |
| 63 | + | out.reserve(static_cast<std::size_t>(digits)); | |
| 64 | + | ||
| 65 | + | long predigit = 0; // digit waiting to be flushed | |
| 66 | + | long nines = 0; // count of pending 9s between predigit and current | |
| 67 | + | long produced = 0; // digits actually written to `out` | |
| 68 | + | ||
| 69 | + | auto emit = [&](char c) { | |
| 70 | + | if (produced < digits) { | |
| 71 | + | out.push_back(c); | |
| 72 | + | ++produced; | |
| 73 | + | } | |
| 74 | + | }; | |
| 75 | + | ||
| 76 | + | for (long j = 0; j < digits && produced < digits; ++j) { | |
| 77 | + | std::int64_t q = 0; | |
| 78 | + | ||
| 79 | + | // Sweep from the high end of the array down to index 0, | |
| 80 | + | // performing the mixed-radix reduction. Each cell i uses the | |
| 81 | + | // odd modulus (2i+1); the multiplier on the way down is i. | |
| 82 | + | for (long i = len; i > 0; --i) { | |
| 83 | + | std::int64_t x = 10 * a[i - 1] + q * i; | |
| 84 | + | a[i - 1] = x % (2 * i - 1); | |
| 85 | + | q = x / (2 * i - 1); | |
| 86 | + | } | |
| 87 | + | ||
| 88 | + | // Lowest cell holds the next decimal digit (mod 10); the rest | |
| 89 | + | // of q is the carry into the digit stream. | |
| 90 | + | a[0] = q % 10; | |
| 91 | + | q /= 10; | |
| 92 | + | ||
| 93 | + | if (q == 9) { | |
| 94 | + | // Buffer this 9 — it might still flip to 0 from a future | |
| 95 | + | // carry. | |
| 96 | + | ++nines; | |
| 97 | + | } else if (q == 10) { | |
| 98 | + | // Carry propagates: bump predigit, all buffered 9s become 0s. | |
| 99 | + | emit(static_cast<char>('0' + predigit + 1)); | |
| 100 | + | for (long k = 0; k < nines && produced < digits; ++k) emit('0'); | |
| 101 | + | predigit = 0; | |
| 102 | + | nines = 0; | |
| 103 | + | } else { | |
| 104 | + | // No carry — flush predigit and the run of 9s, start a new | |
| 105 | + | // pending digit. | |
| 106 | + | if (j != 0) { | |
| 107 | + | emit(static_cast<char>('0' + predigit)); | |
| 108 | + | for (long k = 0; k < nines && produced < digits; ++k) emit('9'); | |
| 109 | + | } | |
| 110 | + | predigit = q; | |
| 111 | + | nines = 0; | |
| 112 | + | } | |
| 113 | + | } | |
| 114 | + | ||
| 115 | + | // Flush the final pending digit and any trailing 9s, truncated to | |
| 116 | + | // the requested length. | |
| 117 | + | if (produced < digits) emit(static_cast<char>('0' + predigit)); | |
| 118 | + | for (long k = 0; k < nines && produced < digits; ++k) emit('9'); | |
| 119 | + | ||
| 120 | + | return out; | |
| 121 | + | } | |
| 122 | + | ||
| 123 | + | } // namespace | |
| 124 | + | ||
| 125 | + | int main(int argc, char** argv) { | |
| 126 | + | if (argc != 2) { | |
| 127 | + | std::cerr << "usage: " << (argc ? argv[0] : "pi_spigot") << " <digits>\n"; | |
| 128 | + | return EXIT_FAILURE; | |
| 129 | + | } | |
| 130 | + | ||
| 131 | + | const long n = parse_count(argv[1]); | |
| 132 | + | if (n < 0) { | |
| 133 | + | std::cerr << "error: <digits> must be a non-negative integer\n"; | |
| 134 | + | return EXIT_FAILURE; | |
| 135 | + | } | |
| 136 | + | ||
| 137 | + | const auto t0 = std::chrono::steady_clock::now(); | |
| 138 | + | const auto digits = compute_pi(n); | |
| 139 | + | const auto t1 = std::chrono::steady_clock::now(); | |
| 140 | + | ||
| 141 | + | // Format as "3.14159..." — first digit is the integer part, the | |
| 142 | + | // remaining n-1 sit after the decimal point. For n == 0 print | |
| 143 | + | // nothing; for n == 1 just print "3". | |
| 144 | + | if (n >= 1) { | |
| 145 | + | std::cout << digits.front(); | |
| 146 | + | if (n >= 2) { | |
| 147 | + | std::cout << '.' << std::string_view(digits).substr(1); | |
| 148 | + | } | |
| 149 | + | std::cout << '\n'; | |
| 150 | + | } | |
| 151 | + | ||
| 152 | + | const auto ms = std::chrono::duration<double, std::milli>(t1 - t0).count(); | |
| 153 | + | std::cerr << "computed " << n << " digit" << (n == 1 ? "" : "s") | |
| 154 | + | << " in " << ms << " ms\n"; | |
| 155 | + | ||
| 156 | + | return EXIT_SUCCESS; | |
| 157 | + | } | |
Newer
Older