洛谷 P1045 [NOIP 2003 普及组] 麦森数

// P1045 麦森数 —— 正解：位数用对数算，末 500 位用**高精度快速幂**
//
// ★★★ 这道题是「快速幂只要结合律」的第二个现场（第一个是 [P3390](/sol/p3390/)）：
//   这里把「取模」换成了「**只留末 500 位**」—— 而那本来就是 mod 10⁵⁰⁰。
//   它成立的理由和取模一模一样：**乘积的末 500 位，只由两个乘数的末 500 位决定。**
//   ⇒ 那五行一个字不改。
//
// 两件事各走各的路：
//   ① **位数**：2^P − 1 和 2^P 位数相同（2^P 不是 10 的幂）⇒ 位数 = ⌊P·log₁₀2⌋ + 1。
//      ⚠ 这一半**根本不用高精度**。
//   ② **末 500 位**：高精度快速幂算 2^P，再减 1。
//      ★ 减 1 不会借位：P ≥ 1 时 2^P 的末位是 2/4/6/8 ⇒ 只改最后一位。
//
// ⚠ 压位：一个 int 存 4 位十进制 ⇒ 500 位只要 125 个 limb，
//   一次乘法 125² = 15 625 次 —— 而快速幂一共只做四十几次乘法。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

static const int L = 125;          // 125 个 limb × 4 位 = 500 位
static const int BASE = 10000;

struct Big { int a[L]; Big() { memset(a, 0, sizeof a); } };

/** 只保留末 500 位的乘法（也就是 mod 10⁵⁰⁰） */
static Big mul(const Big& x, const Big& y) {
    static ll tmp[L];
    memset(tmp, 0, sizeof tmp);
    for (int i = 0; i < L; i++) {
        if (!x.a[i]) continue;
        for (int j = 0; i + j < L; j++)           // ★ i + j ≥ L 的那些位，正是被「取模」丢掉的
            tmp[i + j] += (ll)x.a[i] * y.a[j];
    }
    Big z;
    ll carry = 0;
    for (int i = 0; i < L; i++) { ll v = tmp[i] + carry; z.a[i] = (int)(v % BASE); carry = v / BASE; }
    return z;                                      // carry 溢出去的也丢掉 —— 同样是取模
}

int main() {
    int p;
    if (scanf("%d", &p) != 1) return 0;

    printf("%d\n", (int)(p * log10(2.0)) + 1);     // ① 位数

    Big res, base;                                  // ② 末 500 位：快速幂
    res.a[0] = 1;
    base.a[0] = 2;
    for (int b = p; b; b >>= 1) {
        if (b & 1) res = mul(res, base);
        base = mul(base, base);
    }
    res.a[0] -= 1;                                  // ★ 2^P 末位是偶数 ⇒ 不借位

    char s[501];
    for (int i = 0; i < L; i++) {                   // 铺成 500 个字符，高位在前
        int v = res.a[i];
        for (int d = 0; d < 4; d++) { s[500 - (i * 4 + d) - 1] = char('0' + v % 10); v /= 10; }
    }
    for (int r = 0; r < 10; r++) { fwrite(s + r * 50, 1, 50, stdout); putchar('\n'); }
    return 0;
}

// ⚠ 第 ① 版：不用快速幂 —— 老老实实乘 P 次 2（每次只留末 500 位）
//
// ★★ 它比你以为的强得多，而理由很具体：**乘以 2 是 O(len) 的，不是 O(len²)。**
//   快速幂省下来的是「四十几次 125×125 的乘法」，
//   而这一版做的是「三百多万次 125 长度的扫描」——
//   ⇒ 基本操作多了几百倍，可每一次都便宜得多。
// ⇒ 这一版到底过不过得去，见页面第 1 步那张表（结论和直觉不一样）。
#include <bits/stdc++.h>
using namespace std;

static const int L = 125;
static const int BASE = 10000;

int main() {
    int p;
    if (scanf("%d", &p) != 1) return 0;
    printf("%d\n", (int)(p * log10(2.0)) + 1);

    static int a[L];
    memset(a, 0, sizeof a);
    a[0] = 1;
    for (int t = 0; t < p; t++) {                 // 乘 P 次 2
        int carry = 0;
        for (int i = 0; i < L; i++) { int v = a[i] * 2 + carry; a[i] = v % BASE; carry = v / BASE; }
    }
    a[0] -= 1;

    char s[501];
    for (int i = 0; i < L; i++) { int v = a[i]; for (int d = 0; d < 4; d++) { s[500 - (i * 4 + d) - 1] = char('0' + v % 10); v /= 10; } }
    for (int r = 0; r < 10; r++) { fwrite(s + r * 50, 1, 50, stdout); putchar('\n'); }
    return 0;
}

// P1045 的度量程序：./p1045Count csv        （本页的数字都出自它）
//
//   line  : ★ 「不足 500 位」那条线精确在哪个 P（穷举求出）。
//   prec  : ★★ 位数用 double 算安不安全 —— 把 (1000, 3.1×10⁶) 里每个 P 的
//           frac(P·log₁₀2) 离整数最近的距离量出来，和 double 的误差比一比。
//   ops   : ★★★ 快速幂做几次「125×125 的乘法」↔ 连乘做几次「长度 125 的扫描」。
//   ms    : 顶格 P = 3.1×10⁶ 的秒表：正解 ／ 连乘 ／ 不截断（后者只能算，跑不动）。
//   len   : 不截断的话，2^P 有多少位、多少个 limb、一次乘法多少次基本运算。
//
// ⚠ 秒表用 steady_clock 在进程内量，每档 3 次取中位数。
#include <bits/stdc++.h>
#include <chrono>
using namespace std;
using namespace std::chrono;
typedef long long ll;

static const int L = 125, BASE = 10000;
static const int PMAX = 3100000;

struct Big { int a[L]; Big() { memset(a, 0, sizeof a); } };
static ll mulOps;
static Big mul(const Big& x, const Big& y) {
    static ll tmp[L]; memset(tmp, 0, sizeof tmp);
    for (int i = 0; i < L; i++) { if (!x.a[i]) continue; for (int j = 0; i + j < L; j++) { tmp[i + j] += (ll)x.a[i] * y.a[j]; mulOps++; } }
    Big z; ll c = 0;
    for (int i = 0; i < L; i++) { ll v = tmp[i] + c; z.a[i] = (int)(v % BASE); c = v / BASE; }
    return z;
}
static double med3(double a, double b, double c) { return max(min(a, b), min(max(a, b), c)); }

int main(int argc, char** argv) {
    bool csv = (argc > 1 && string(argv[1]) == "csv");
    const double LG2 = log10(2.0);

    /* ① 「不足 500 位」那条线 */
    int line = 0;
    for (int p = 1001; p <= 3000; p++) if ((int)(p * LG2) + 1 >= 500) { line = p; break; }

    /* ② 位数公式用 double 安不安全：frac(P·log₁₀2) 离 0 / 1 最近多少 */
    double minFrac = 1.0;
    int minP = 0;
    for (int p = 1001; p < PMAX; p++) {
        long double v = (long double)p * (long double)LG2;
        long double f = v - floorl(v);
        double d = (double)min(f, (long double)1.0 - f);
        if (d < minFrac) { minFrac = d; minP = p; }
    }

    /* ③ 两种写法的基本运算次数 */
    ll fastMuls = 0, fastOps;
    {
        mulOps = 0;
        Big res, base; res.a[0] = 1; base.a[0] = 2;
        for (int b = PMAX - 1; b; b >>= 1) { if (b & 1) { res = mul(res, base); fastMuls++; } base = mul(base, base); fastMuls++; }
        fastOps = mulOps;
    }
    ll bruteOps = (ll)(PMAX - 1) * L;

    /* ④ 秒表：顶格 P */
    static volatile int VP = PMAX - 1;
    double msFast, msBrute;
    {
        double t[3];
        for (int r = 0; r < 3; r++) {
            int p = VP;
            auto t0 = steady_clock::now();
            Big res, base; res.a[0] = 1; base.a[0] = 2;
            for (int b = p; b; b >>= 1) { if (b & 1) res = mul(res, base); base = mul(base, base); }
            volatile int sink = res.a[0]; (void)sink;
            t[r] = duration<double, milli>(steady_clock::now() - t0).count();
        }
        msFast = med3(t[0], t[1], t[2]);
    }
    {
        double t[3];
        for (int r = 0; r < 3; r++) {
            int p = VP;
            auto t0 = steady_clock::now();
            static int a[L]; memset(a, 0, sizeof a); a[0] = 1;
            for (int i = 0; i < p; i++) { int c = 0; for (int j = 0; j < L; j++) { int v = a[j] * 2 + c; a[j] = v % BASE; c = v / BASE; } }
            volatile int sink = a[0]; (void)sink;
            t[r] = duration<double, milli>(steady_clock::now() - t0).count();
        }
        msBrute = med3(t[0], t[1], t[2]);
    }

    /* ⑤ 不截断的话有多大 */
    ll topDigits = (ll)((PMAX - 1) * LG2) + 1;
    ll topLimbs = (topDigits + 3) / 4;
    double topOnceOps = (double)topLimbs * (double)topLimbs;

    if (csv) {
        printf("line,%d\n", line);
        printf("minFrac,%.3e\n", minFrac);
        printf("minP,%d\n", minP);
        printf("fastMuls,%lld\n", fastMuls);
        printf("fastOps,%lld\n", fastOps);
        printf("bruteOps,%lld\n", bruteOps);
        printf("opsRatio,%.0f\n", (double)bruteOps / (double)fastOps);
        printf("msFast,%.2f\n", msFast);
        printf("msBrute,%.0f\n", msBrute);
        printf("bruteOverFast,%.0f\n", msBrute / msFast);
        printf("topDigits,%lld\n", topDigits);
        printf("topLimbs,%lld\n", topLimbs);
        printf("topOnceOps,%.3e\n", topOnceOps);
        return 0;
    }
    printf("① 「2^P − 1 不足 500 位」那条线：第一个满 500 位的 P 是 **%d**\n", line);
    printf("\n② 位数用 double 算安不安全：P ∈ (1000, 3.1×10⁶) 里，P·log₁₀2 离整数最近的距离是 %.3e（P = %d）\n",
           minFrac, minP);
    printf("   ⇒ double 在这个量级上的误差约 1e−10 ⇒ **够用，而且余量还有三四个数量级**\n");
    printf("\n③ 基本运算次数（顶格 P = %d）\n", PMAX - 1);
    printf("   快速幂：%lld 次高精度乘法，一共 %lld 次 limb 乘法\n", fastMuls, fastOps);
    printf("   连乘  ：%d 次「乘 2」，一共 %lld 次 limb 运算 ⇒ 是快速幂的 %.0f 倍\n",
           PMAX - 1, bruteOps, (double)bruteOps / (double)fastOps);
    printf("\n④ 秒表（独占实测，3 次取中位数；时限 1000 ms）\n");
    printf("   ★ 快速幂 %.2f ms ／ 连乘 %.0f ms ⇒ 差 %.0f 倍\n", msFast, msBrute, msBrute / msFast);
    printf("   ⚠ 而**两者都过得去** —— 连乘赢在「乘 2 是 O(len) 而不是 O(len²)」\n");
    printf("\n⑤ 如果不截断：2^P 有 %lld 位、%lld 个 limb，**一次**乘法就要 %.3e 次运算 ⇒ 想都不用想\n",
           topDigits, topLimbs, topOnceOps);
    return 0;
}