- 注册时间
- 2007-12-27
- 最后登录
- 1970-1-1
- 威望
- 星
- 金币
- 枚
- 贡献
- 分
- 经验
- 点
- 鲜花
- 朵
- 魅力
- 点
- 上传
- 次
- 下载
- 次
- 积分
- 50402
- 在线时间
- 小时
|
使用我上面的思路可以有下面的代码(我的输出是第一项是1，和OEIS保持一致），主要计算原则就是对于上界，总是向上舍入，对于下界，总是向下舍入。
- #include <gmpxx.h>
- #include <cstdlib>
- #include <iostream>
- #include <sstream>
- #include <stdexcept>
- #include <string>
- namespace {
- struct Options {
- unsigned long keep_digits = 1000;
- unsigned long steps = 20;
- bool verbose = false;
- };
- struct DecimalState {
- mpz_class len10;
- unsigned long kept = 1;
- mpz_class prefix_lo;
- mpz_class prefix_hi_exclusive;
- };
- // Directed integer bound mant * 2^shift.
- // For lower bounds, truncation always rounds mantissa downward.
- // For upper bounds / exclusive upper bounds, truncation always rounds upward.
- struct BinBound {
- mpz_class mant;
- mpz_class shift;
- };
- [[noreturn]] void Fail(const std::string& message) {
- throw std::runtime_error(message);
- }
- unsigned long BitLen(const mpz_class& x) {
- if (x <= 0) {
- return 0;
- }
- return mpz_sizeinbase(x.get_mpz_t(), 2);
- }
- mpz_class Pow10(unsigned long n) {
- mpz_class out = 1;
- mpz_ui_pow_ui(out.get_mpz_t(), 10, n);
- return out;
- }
- std::string BinaryStringFixed(const mpz_class& x, unsigned long width) {
- std::string s = x.get_str(2);
- if (s.size() > width) {
- Fail("binary width overflow while formatting prefix");
- }
- if (s.size() < width) {
- s.insert(s.begin(), width - s.size(), '0');
- }
- return s;
- }
- mpz_class DecimalFromBitString(const std::string& bits) {
- return mpz_class(bits, 10);
- }
- mpz_class CeilShiftRight(const mpz_class& x, unsigned long bits) {
- if (bits == 0) {
- return x;
- }
- mpz_class add = (mpz_class(1) << bits) - 1;
- return (x + add) >> bits;
- }
- void TruncateDown(BinBound& v, unsigned long keep_bits) {
- if (v.mant <= 0) {
- Fail("invalid lower binary bound");
- }
- unsigned long bits = BitLen(v.mant);
- if (bits <= keep_bits) {
- return;
- }
- unsigned long drop = bits - keep_bits;
- v.mant >>= drop;
- v.shift += drop;
- }
- void TruncateUp(BinBound& v, unsigned long keep_bits) {
- if (v.mant <= 0) {
- Fail("invalid upper binary bound");
- }
- unsigned long bits = BitLen(v.mant);
- if (bits <= keep_bits) {
- return;
- }
- unsigned long drop = bits - keep_bits;
- v.mant = CeilShiftRight(v.mant, drop);
- v.shift += drop;
- }
- BinBound ExactIntegerDown(const mpz_class& n, unsigned long keep_bits) {
- if (n <= 0) {
- Fail("exact lower bound requires a positive integer");
- }
- BinBound out{n, 0};
- TruncateDown(out, keep_bits);
- return out;
- }
- BinBound ExactIntegerUp(const mpz_class& n, unsigned long keep_bits) {
- if (n <= 0) {
- Fail("exact upper bound requires a positive integer");
- }
- BinBound out{n, 0};
- TruncateUp(out, keep_bits);
- return out;
- }
- BinBound MultiplyDown(const BinBound& a, const BinBound& b,
- unsigned long keep_bits) {
- BinBound out;
- out.mant = a.mant * b.mant;
- out.shift = a.shift + b.shift;
- TruncateDown(out, keep_bits);
- return out;
- }
- BinBound MultiplyUp(const BinBound& a, const BinBound& b,
- unsigned long keep_bits) {
- BinBound out;
- out.mant = a.mant * b.mant;
- out.shift = a.shift + b.shift;
- TruncateUp(out, keep_bits);
- return out;
- }
- BinBound PowSmallBaseDown(unsigned long base, const mpz_class& exp,
- unsigned long keep_bits) {
- BinBound result = ExactIntegerDown(1, keep_bits);
- BinBound cur = ExactIntegerDown(base, keep_bits);
- mpz_class e = exp;
- while (e > 0) {
- if (mpz_odd_p(e.get_mpz_t())) {
- result = MultiplyDown(result, cur, keep_bits);
- }
- e >>= 1;
- if (e > 0) {
- cur = MultiplyDown(cur, cur, keep_bits);
- }
- }
- return result;
- }
- BinBound PowSmallBaseUp(unsigned long base, const mpz_class& exp,
- unsigned long keep_bits) {
- BinBound result = ExactIntegerUp(1, keep_bits);
- BinBound cur = ExactIntegerUp(base, keep_bits);
- mpz_class e = exp;
- while (e > 0) {
- if (mpz_odd_p(e.get_mpz_t())) {
- result = MultiplyUp(result, cur, keep_bits);
- }
- e >>= 1;
- if (e > 0) {
- cur = MultiplyUp(cur, cur, keep_bits);
- }
- }
- return result;
- }
- mpz_class LowerBitLength(const BinBound& v) {
- if (v.mant <= 0) {
- Fail("invalid lower bound for bit length");
- }
- return v.shift + BitLen(v.mant);
- }
- mpz_class UpperLastBitLength(const BinBound& upper_exclusive) {
- if (upper_exclusive.mant <= 0) {
- Fail("invalid upper bound for bit length");
- }
- if (upper_exclusive.shift == 0) {
- return BitLen(upper_exclusive.mant - 1);
- }
- if (upper_exclusive.mant == 1) {
- return upper_exclusive.shift;
- }
- return upper_exclusive.shift + BitLen(upper_exclusive.mant);
- }
- mpz_class ShiftedLowerQuotient(const BinBound& lower, const mpz_class& drop) {
- if (drop < 0) {
- Fail("negative drop for lower quotient");
- }
- if (lower.shift >= drop) {
- mpz_class delta = lower.shift - drop;
- if (!delta.fits_ulong_p()) {
- Fail("left shift too large while extracting lower prefix bits");
- }
- return lower.mant << delta.get_ui();
- }
- mpz_class delta = drop - lower.shift;
- if (!delta.fits_ulong_p()) {
- return 0;
- }
- return lower.mant >> delta.get_ui();
- }
- mpz_class ShiftedUpperExclusiveQuotient(const BinBound& upper_exclusive,
- const mpz_class& drop) {
- if (drop < 0) {
- Fail("negative drop for upper quotient");
- }
- if (upper_exclusive.shift >= drop) {
- mpz_class delta = upper_exclusive.shift - drop;
- if (!delta.fits_ulong_p()) {
- Fail("left shift too large while extracting upper prefix bits");
- }
- return upper_exclusive.mant << delta.get_ui();
- }
- mpz_class delta = drop - upper_exclusive.shift;
- if (!delta.fits_ulong_p()) {
- return 1;
- }
- return CeilShiftRight(upper_exclusive.mant, delta.get_ui());
- }
- mpz_class ClampKeptDigits(const mpz_class& len10, unsigned long keep_digits) {
- if (len10 <= keep_digits) {
- return len10;
- }
- return mpz_class(keep_digits);
- }
- unsigned long AsUnsignedLongChecked(const mpz_class& x, const std::string& what) {
- if (!x.fits_ulong_p()) {
- Fail(what + " does not fit into unsigned long");
- }
- return x.get_ui();
- }
- void PrintState(unsigned long index, const DecimalState& s) {
- std::cout << "n=" << index
- << " len=" << s.len10.get_str()
- << " kept=" << s.kept
- << " prefix_lo=" << s.prefix_lo.get_str()
- << " prefix_hi_exclusive=" << s.prefix_hi_exclusive.get_str()
- << "\n";
- }
- DecimalState InitialState() {
- return DecimalState{1, 1, 2, 3};
- }
- DecimalState NextState(const DecimalState& cur, unsigned long keep_digits,
- bool verbose) {
- if (cur.prefix_hi_exclusive <= cur.prefix_lo) {
- Fail("invalid decimal interval");
- }
- mpz_class kept_mpz = ClampKeptDigits(cur.len10, keep_digits);
- unsigned long kept = AsUnsignedLongChecked(kept_mpz, "kept digits");
- mpz_class e10 = cur.len10 - kept_mpz;
- const unsigned long keep_bits = keep_digits + 32;
- BinBound pow5_lo = PowSmallBaseDown(5, e10, keep_bits);
- BinBound pow5_hi = PowSmallBaseUp(5, e10, keep_bits);
- BinBound prefix_lo = ExactIntegerDown(cur.prefix_lo, keep_bits);
- BinBound prefix_hi = ExactIntegerUp(cur.prefix_hi_exclusive, keep_bits);
- BinBound x_lo = MultiplyDown(prefix_lo, pow5_lo, keep_bits);
- x_lo.shift += e10;
- BinBound x_hi_exclusive = MultiplyUp(prefix_hi, pow5_hi, keep_bits);
- x_hi_exclusive.shift += e10;
- mpz_class next_len_lo = LowerBitLength(x_lo);
- mpz_class next_len_hi = UpperLastBitLength(x_hi_exclusive);
- if (next_len_lo != next_len_hi) {
- std::ostringstream oss;
- oss << "precision exhausted: next length is only known to be in ["
- << next_len_lo.get_str() << ", " << next_len_hi.get_str() << "]";
- Fail(oss.str());
- }
- mpz_class next_len = next_len_lo;
- mpz_class next_kept_mpz = ClampKeptDigits(next_len, keep_digits);
- unsigned long next_kept =
- AsUnsignedLongChecked(next_kept_mpz, "next kept digits");
- mpz_class drop = next_len - next_kept_mpz;
- mpz_class next_prefix_lo_bits = ShiftedLowerQuotient(x_lo, drop);
- mpz_class next_prefix_hi_bits_exclusive =
- ShiftedUpperExclusiveQuotient(x_hi_exclusive, drop);
- mpz_class upper_limit = mpz_class(1) << next_kept;
- if (next_prefix_hi_bits_exclusive > upper_limit) {
- next_prefix_hi_bits_exclusive = upper_limit;
- }
- std::string lo_bits = BinaryStringFixed(next_prefix_lo_bits, next_kept);
- mpz_class next_prefix_lo = DecimalFromBitString(lo_bits);
- mpz_class next_prefix_hi_exclusive;
- if (next_prefix_hi_bits_exclusive == upper_limit) {
- next_prefix_hi_exclusive = Pow10(next_kept);
- } else {
- std::string hi_bits =
- BinaryStringFixed(next_prefix_hi_bits_exclusive, next_kept);
- next_prefix_hi_exclusive = DecimalFromBitString(hi_bits);
- }
- if (verbose) {
- std::cout << " e10=" << e10.get_str()
- << " next_len_lo=" << next_len_lo.get_str()
- << " next_len_hi=" << next_len_hi.get_str()
- << " next_prefix_lo_bits=" << next_prefix_lo_bits.get_str(2)
- << " next_prefix_hi_bits_exclusive="
- << next_prefix_hi_bits_exclusive.get_str(2) << "\n";
- }
- if (next_prefix_hi_exclusive <= next_prefix_lo) {
- Fail("next decimal prefix interval collapsed");
- }
- return DecimalState{next_len, next_kept, next_prefix_lo,
- next_prefix_hi_exclusive};
- }
- Options ParseOptions(int argc, char** argv) {
- Options opt;
- for (int i = 1; i < argc; ++i) {
- std::string arg = argv[i];
- auto need_value = [&](const std::string& name) -> std::string {
- if (i + 1 >= argc) {
- Fail("missing value for " + name);
- }
- return argv[++i];
- };
- if (arg == "--keep-digits") {
- opt.keep_digits = std::stoul(need_value(arg));
- } else if (arg == "--steps") {
- opt.steps = std::stoul(need_value(arg));
- } else if (arg == "--verbose") {
- opt.verbose = true;
- } else if (arg == "--help" || arg == "-h") {
- std::cout
- << "Usage: bin2dec [--keep-digits M] [--steps N] [--verbose]\n";
- std::exit(0);
- } else {
- Fail("unknown argument: " + arg);
- }
- }
- if (opt.keep_digits == 0) {
- Fail("--keep-digits must be positive");
- }
- if (opt.steps == 0) {
- Fail("--steps must be positive");
- }
- return opt;
- }
- } // namespace
- int main(int argc, char** argv) {
- try {
- const Options opt = ParseOptions(argc, argv);
- DecimalState state = InitialState();
- PrintState(1, state);
- for (unsigned long i = 2; i <= opt.steps; ++i) {
- state = NextState(state, opt.keep_digits, opt.verbose);
- PrintState(i, state);
- }
- return 0;
- } catch (const std::exception& ex) {
- std::cerr << "error: " << ex.what() << "\n";
- return 1;
- }
- }
比如输入精度为1000时，可以计算到第602项。 |
|
楼主: KeyTo9_Fans
IP卡
狗仔卡
发表于 
显身卡
发表于 
发表于