pi_spigot.cpp
· 5.3 KiB · C++
Raw
// pi_spigot.cpp — compute N decimal digits of pi using the
// Rabinowitz–Wagon spigot algorithm.
//
// "A Spigot Algorithm for the Digits of Pi"
// Stanley Rabinowitz and Stan Wagon, 1995.
//
// The algorithm exploits the series
//
// pi = sum_{i=0..inf} (i! * 2^(i+1)) / (2i+1)!!
//
// rewritten in a mixed-radix form so that decimal digits can be
// extracted one at a time using only integer arithmetic on a fixed
// array of length L = floor(10N/3) + 1. No big-integer library is
// required and every intermediate value fits comfortably in a 64-bit
// integer for the problem sizes a single machine cares about.
//
// Output is "3." followed by N digits after the decimal point, and a
// timing line on stderr.
//
// Build: c++ -std=c++20 -O2 -o pi_spigot pi_spigot.cpp
// Run: ./pi_spigot 1000
#include <charconv>
#include <chrono>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <string>
#include <string_view>
#include <vector>
namespace {
// Parse a non-negative integer from argv. Returns nullopt-style sentinel
// (-1) on failure so main() can print a usage message.
long parse_count(std::string_view s) {
long value = 0;
auto [ptr, ec] = std::from_chars(s.data(), s.data() + s.size(), value);
if (ec != std::errc{} || ptr != s.data() + s.size() || value < 0) {
return -1;
}
return value;
}
// Compute `digits` decimal digits of pi using the Rabinowitz–Wagon
// spigot. The leading "3" is the first digit emitted; the remaining
// digits-1 are the fractional part. The result is returned as a plain
// decimal string of length `digits`.
//
// Carry handling: each generated digit `q` may be 0..10. A value of 10
// means the previously buffered digit must be incremented (and any run
// of intervening 9s flipped to 0s) before being flushed — the classic
// "carry across a streak of nines" trick that makes the spigot work.
std::string compute_pi(long digits) {
if (digits <= 0) return {};
// Array length comes from the convergence bound in the paper: we
// need at least ceil(10N/3) terms to be sure of the first N digits.
const long len = (10 * digits) / 3 + 1;
std::vector<std::int64_t> a(len, 2);
std::string out;
out.reserve(static_cast<std::size_t>(digits));
long predigit = 0; // digit waiting to be flushed
long nines = 0; // count of pending 9s between predigit and current
long produced = 0; // digits actually written to `out`
auto emit = [&](char c) {
if (produced < digits) {
out.push_back(c);
++produced;
}
};
for (long j = 0; j < digits && produced < digits; ++j) {
std::int64_t q = 0;
// Sweep from the high end of the array down to index 0,
// performing the mixed-radix reduction. Each cell i uses the
// odd modulus (2i+1); the multiplier on the way down is i.
for (long i = len; i > 0; --i) {
std::int64_t x = 10 * a[i - 1] + q * i;
a[i - 1] = x % (2 * i - 1);
q = x / (2 * i - 1);
}
// Lowest cell holds the next decimal digit (mod 10); the rest
// of q is the carry into the digit stream.
a[0] = q % 10;
q /= 10;
if (q == 9) {
// Buffer this 9 — it might still flip to 0 from a future
// carry.
++nines;
} else if (q == 10) {
// Carry propagates: bump predigit, all buffered 9s become 0s.
emit(static_cast<char>('0' + predigit + 1));
for (long k = 0; k < nines && produced < digits; ++k) emit('0');
predigit = 0;
nines = 0;
} else {
// No carry — flush predigit and the run of 9s, start a new
// pending digit.
if (j != 0) {
emit(static_cast<char>('0' + predigit));
for (long k = 0; k < nines && produced < digits; ++k) emit('9');
}
predigit = q;
nines = 0;
}
}
// Flush the final pending digit and any trailing 9s, truncated to
// the requested length.
if (produced < digits) emit(static_cast<char>('0' + predigit));
for (long k = 0; k < nines && produced < digits; ++k) emit('9');
return out;
}
} // namespace
int main(int argc, char** argv) {
if (argc != 2) {
std::cerr << "usage: " << (argc ? argv[0] : "pi_spigot") << " <digits>\n";
return EXIT_FAILURE;
}
const long n = parse_count(argv[1]);
if (n < 0) {
std::cerr << "error: <digits> must be a non-negative integer\n";
return EXIT_FAILURE;
}
const auto t0 = std::chrono::steady_clock::now();
const auto digits = compute_pi(n);
const auto t1 = std::chrono::steady_clock::now();
// Format as "3.14159..." — first digit is the integer part, the
// remaining n-1 sit after the decimal point. For n == 0 print
// nothing; for n == 1 just print "3".
if (n >= 1) {
std::cout << digits.front();
if (n >= 2) {
std::cout << '.' << std::string_view(digits).substr(1);
}
std::cout << '\n';
}
const auto ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
std::cerr << "computed " << n << " digit" << (n == 1 ? "" : "s")
<< " in " << ms << " ms\n";
return EXIT_SUCCESS;
}
| 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 | } |
| 158 |