Project Euler 566

题目

Cake Icing Puzzle

Adam plays the following game with his birthday cake.

He cuts a piece forming a circular sector of \(60\) degrees and flips the piece upside down, with the icing on the bottom.

He then rotates the cake by 60 degrees counterclockwise, cuts an adjacent \(60\) degree piece and flips it upside down.

He keeps repeating this, until after a total of twelve steps, all the icing is back on top.

Amazingly, this works for any piece size, even if the cutting angle is an irrational number: all the icing will be back on top after a finite number of steps.

Now, Adam tries something different: he alternates cutting pieces of size \(x=\dfrac{360}{9}\) degrees, \(y=\dfrac{360}{10}\) degrees and \(z=\dfrac{360 }{\sqrt{11}}\) degrees. The first piece he cuts has size \(x\) and he flips it. The second has size \(y\) and he flips it. The third has size \(z\) and he flips it. He repeats this with pieces of size \(x, y\) and \(z\) in that order until all the icing is back on top, and discovers he needs \(60\) flips altogether.

Let \(F(a, b, c)\) be the minimum number of piece flips needed to get all the icing back on top for pieces of size \(x=\dfrac{360}{a}\) degrees, \(y=\dfrac{360}{b}\) degrees and \(z=\dfrac{360}{\sqrt{c}}\) degrees.

Let \(G(n) = \sum_{9 \le a < b < c \le n} F(a,b,c)\), for integers \(a, b\) and \(c\).

You are given that \(F(9,10,11)=60, F(10,14,16)=506, F(15,16,17)=785232\).

You are also given \(G(11)=60, G(14)=58020\) and \(G(17)=1269260\).

Find \(G(53)\).

解决方案

把整块蛋糕的圆周长度归一化为 \(1\)（对应 \(360^\circ\)）。每一步做两件事：从当前起点沿圆周切下一个指定长度的扇形片（相当于取一段连续圆弧）。将这段片翻面（会让其覆盖的区域取反），然后把起点沿圆周前移这段长度，继续切下一段相邻的片。

题目的关键并不在这是不是糖霜一面，而在切口最终会形成怎样的有限分割；一旦分割固定，后续每 \(3\) 次翻片的作用就是一个确定的置换 + 翻转重复作用于这些小段，因此问题转化为求最小周期。

对固定 \(a,b,c\)，把圆周长度放大成一个方便的整尺度 \(L\)：

• 若 \(c=s^2\) 是完全平方数，则 \(z=1/s\) 为有理数，取\(L=\operatorname{lcm}(a,b,s),x=\dfrac{L}{a},y=\dfrac{L}{b},z=\dfrac{L}{s}.\)全部是整数。
• 若 \(c\) 不是平方数，取\(L=\operatorname{lcm}(a,b),x=\dfrac{L}{a},y=\dfrac{L}{b},z=\dfrac{L}{\sqrt c}.\)此时长度属于\(\mathbb Z\left[\dfrac1{\sqrt c}\right]=\left\{u+\dfrac{v}{\sqrt c}:u,v\in\mathbb Z\right\}.\)

我们用二元整数对 \((u,v)\) 表示数 \(u+v/\sqrt c\)。比较两个数大小等价于判断差\(d=u+\dfrac{v}{\sqrt c}\)的符号。因为 \(\sqrt c>0\)，若 \(u,v\) 同号则 \(d\) 同号；若异号，则比较 \(|u|\) 与 \(|v|/\sqrt c\)。用平方消去根号：\(u+\dfrac{v}{\sqrt c} > 0 \iff(u>0 \land u^2c>v^2)\lor(u<0 \land u^2c<v^2).\)

把所有出现过的切口位置收集起来，圆周被切成许多小弧段。把每一段小弧段称为一个 元素。一次切一片并翻面的效果对元素来说只有两点：该片覆盖的若干连续元素顺序会被反转；这些元素的糖霜朝向全部取反。

把当前起点作为 \(0\)，把圆周展开为直线长度 \(L\)。考虑连续三次翻片（依次取长度 \(x,y,z\) 的相邻三片）并且把起点移动到第三片末端（也就是下一轮的起点）。设这三片对应的直线分段分别为字符串（元素序列）\(X,Y,Z,\)剩余部分为 \(R\)，则展开后原始顺序为\(X,Y,Z,R.\)翻片等价于对各片内部反转，因此三次翻片并把起点移到第三片末端后，整体元素序列变为\(R,\overline{X},\overline{Y},\overline{Z},\)其中 \(\overline{W}\) 表示字符串反转。

这说明：每三次翻片对元素的作用是固定的置换，并且在 \(X,Y,Z\) 内的元素会翻转一次，其余不翻。接下来只需知道最终切口稳定后，总共有多少元素 \(N\)，并且在块的起点处，长度 \(x,y,z\) 各自覆盖多少个元素：记为 \(e_x,e_y,e_z\)。一旦得到 \((N,e_x,e_y,e_z)\)，三步合成的置换完全由这三个整数决定，与实际角度数值无关。

令\(X=x, Y=y, Z=z, S=X+Y+Z.\)把块的起点设为 \(0\)，则三片边界切口在位置 \(X,X+Y,S\)。圆周长度为 \(L\)。把 \(L\) 写成\(L = mS + E, 0<E\le S,\) 其中 \(m=\lfloor L/S\rfloor\)，\(E\) 是余长。这里的 \(\lfloor\cdot\rfloor\) 在 \(\mathbb Z[1/\sqrt c]\) 的有序结构里是良定义的（用上面的精确比较找最大 \(m\) 使 \((m+1)S\le L\)）。

切口稳定后，长度为 \(S\) 的每个完整块内切口结构相同，整个圆周由 \(m\) 个完整块加一个长度 \(E\) 的尾巴组成。因此：一个完整块内元素数为 \(e_x+e_y+e_z\)；全圆周元素数\(N = m(e_x+e_y+e_z) + e_E,\)其中 \(e_E\) 是尾巴部分包含的元素数（等价于切口位置 \(\le E\)的个数）。

于是只需找出第一块 \([0,S]\) 内所有稳定切口的位置（只需知道它们落在 \(X/Y/Z\) 的哪一段，以及是否 \(\le E\)），就能得到所有计数。

取第一块中任一切口位置 \(t\in(0,S]\)。在展开串\(X,Y,Z,R \to R,\overline X,\overline Y,\overline Z\)中跟踪这个点的位置：

• 若 \(t\in(0,X]\)（在 \(X\) 内），它在 \(\overline X\) 内距 \(\overline X\) 起点的位置为 \(X-t\)，\(\overline X\) 在新串中的起点距离为 \(|R|=L-S\)，因此新坐标为\((L-S)+(X-t)=L-(Y+Z)-t.\)
• 若 \(t\in(X,X+Y]\)，类似得到\(L+X-Z-t.\)
• 若 \(t\in(X+Y,S]\)，得到\(L+X+Y-t.\)

下一轮我们关心的是它在下一块起点下的相对位置，即对 \(S\) 取模并落回 \((0,S]\)：

\[ t' \equiv \begin{cases} L-(Y+Z)-t & t\in X,\\ L+X-Z-t & t\in Y,\\ L+X+Y-t & t\in Z. \end{cases} \pmod S \]

切口传播过程是确定且可逆的，因此从边界切口 \(X,X+Y,S\) 出发沿此映射迭代，只会在有限集合中循环。对题目范围内的所有三元组，该迭代都会回到某个边界切口，从而得到有限切口集合。

因此，从起始点 \(X,X+Y,S\) 分别出发，不断算 \(t\mapsto t'\)：若 \(t'\) 落在 \(X/Y/Z\) 哪段，就把对应 \(e_x/e_y/e_z\) 加一；若 \(t'\le E\)，则把 \(e_E\) 加一；一旦 \(t'\) 等于 \(X,X+Y,S\) 之一就停止该轨道。初始化时 \(e_x=e_y=e_z=1\)（每段至少有一个元素），并把 \(X,X+Y\) 是否 \(\le E\) 计入 \(e_E\)。这样即可在完全精确的整数运算下得到 \((N,e_x,e_y,e_z)\)，而不需要用高级数据结构去维护全部切口位置。

接下来，我们通过三步合成变成置换和翻转，并且周期可分解到环。

把元素按当前起点顺序编号 \(0,1,\dots,N-1\)。一次翻长度覆盖前 \(p\) 个元素的片并把起点移到片末端，等价于：前缀 \([0,p)\) 被反转并移到末尾；这 \(p\) 个元素翻转一次；其余元素向前平移 \(p\)。因此索引映射为\(i\mapsto\begin{cases}N-1-i & i<p,\\i-p & i\ge p.\end{cases}\)并且当 \(i<p\) 时该元素糖霜朝向取反；当 \(i\ge p\) 时糖霜朝向保持不变。

把该操作记为一个带翻转标记的置换，那么依次取 \(p=e_x,e_y,e_z\) 的三次合成，得到总置换 \(P\) 及每个位置的一步翻转标记。

把 \(P\) 分解为不相交的环。对某个环，按环上的顺序把元素写成\(v_0\to v_1\to \cdots \to v_{n-1}\to v_0,\)并定义二进制序列\(g_j\)表示一步从\(v_j\)走到\(v_{j+1}\)时是否翻转。

初始时所有元素糖霜朝上（\(0\)）。做 \(k\) 次三步合成后，元素 \(v_j\) 被翻转的奇偶为\(p_k(j)=g_j\oplus g_{j+1}\oplus\cdots\oplus g_{j+k-1},\)（下标按 \(n\) 取模）。要求全朝上等价于对该环所有 \(j\)，都有 \(p_k(j)=0\)。

设 \(g\) 的最小周期为 \(p\)（把 \(g\) 看作长度 \(n\) 的循环词，等价于在线性串中求最小重复块长度）。因为 \(g\) 是 \(p\)-周期的，长度为 \(p\) 的任意连续窗口包含同一段周期词，因此\(\displaystyle p_p(j)=\bigoplus_{t=0}^{p-1} g_t\)与 \(j\) 无关：做 \(p\) 次后三步合成，环上所有元素要么全不翻（若该异或为 \(0\)），要么全翻一次（若为 \(1\)）。若 \(\displaystyle \bigoplus_{t=0}^{p-1} g_t=0\)，则 \(k=p\) 就使该环全回到朝上，且更小的 \(k\) 不可能（否则会推出更小周期）。若该异或为 \(1\)，则 \(k=p\) 时全翻一次，还需要再做 \(p\) 次才能翻回，最小 \(k=2p\)。

因此每个环的回归次数是 \(p\) 或 \(2p\)，整体（要求所有环同时回归）所需的三步合成次数是这些环回归次数的最小公倍数 \(T\)，对应翻片次数为\(F_{\equiv 0}(a,b,c)=3T.\)

不过，翻片总次数可能不是 \(3\) 的倍数。若答案为 \(3k+1\)，等价于先单独做一次翻 \(X\)；然后重复做三步块，但块的顺序变为 \((Y,Z,X)\)，共做 \(k\) 次。

这时进入块循环时，初始糖霜状态不是全 \(0\)：恰好被那一次 \(X\) 涉及的元素翻过一次。在编号 \(0,\dots,N-1\) 的线性表示下，做完一次 \(X\) 后，未被翻的那段长度是 \(N-e_x\)，被翻过的那段是末尾的 \(e_x\) 个，因此初始目标向量为\(t(i)=\mathbf{1}\{i\ge N-e_x\}.\)对块置换 \(Q_{(Y,Z,X)}\) 的每个环：

• 若该环上的 \(t\) 全为 \(0\)，则该环起始时所有元素仍然朝上，因此它只要求块迭代次数 \(k\) 必须是该环回到全朝上的最小回归次数的倍数。设该回归次数为 \(T\)（上文已说明 \(T\) 只能是 \(p\) 或 \(2p\)），则约束为 \(k\equiv 0\pmod T\)，可以并入全局模数。
• 若该环上 \(t\) 混杂（环跨过了分界点），就不能用“倍数整除”直接给出 \(k\)，但环的周期仍然是上一节的 \(p\) 或 \(2p\)。此时可以在该环内顺着复合 \(Q\) 一步步迭代，最多试到周期长度，找最小 \(k\) 使得“累计翻转向量 \(= t\)”。若找不到则该同余类无解。

不同环给出同余条件 \(k\equiv r_i\pmod{m_i}\)。用带兼容性检查的中国剩余定理逐步合并，一旦冲突立刻判定无解。若合并成功得到最小正解 \(k\)，则\(F_{\equiv 1}(a,b,c)=3k+1.\)

同理，若答案为 \(3k+2\)，就是先做 \(X\) 再做 \(Y\)，然后块顺序为 \((Z,X,Y)\)，初始翻转段长度为 \(e_x+e_y\)，即可得到\(F_{\equiv 2}(a,b,c)=3k+2.\)

最终\(F(a,b,c)\)是\(F_{\equiv 0}(a,b,c),F_{\equiv 1}(a,b,c),F_{\equiv 2}(a,b,c)\)中的最小值（如果它们存在）。

大量不同的 \((a,b,c)\) 会得到相同的剖分四元组 \((N,e_x,e_y,e_z)\)，从而 \(F\) 也相同。用哈希表对该四元组记忆化即可显著加速。

代码

#include <bits/stdc++.h>
#include "tools.h"
using namespace std;

const int N = 53;
typedef long long ll;
static ll lcmll(ll a, ll b)
{
    return a / __gcd(a, b) * b;
}

struct Quad
{
    ll u = 0, v = 0;

    Quad() = default;
    Quad(ll uu, ll vv) : u(uu), v(vv) {}

    static int sgn(ll x) { return (x > 0) - (x < 0); }

    static int signum(ll du, ll dv, ll c)
    {
        if (du == 0)
            return sgn(dv);
        if (dv == 0)
            return sgn(du);
        int su = sgn(du), sv = sgn(dv);
        if (su == sv)
            return su;
        ll uu = du * du;
        ll vv = dv * dv;
        ll cmp = uu * c - vv;
        int sc = (cmp > 0) - (cmp < 0);
        return su * sc;
    }

    int cmp(const Quad &o, ll c) const
    {
        return signum(u - o.u, v - o.v, c);
    }

    Quad norm_mod(const Quad &mod, ll c) const
    {
        Quad x = *this;
        Quad zero;
        while (x.cmp(zero, c) <= 0)
            x = x + mod;
        while (x.cmp(mod, c) > 0)
            x = x - mod;
        return x;
    }

    friend Quad operator+(const Quad &a, const Quad &b) { return {a.u + b.u, a.v + b.v}; }
    friend Quad operator-(const Quad &a, const Quad &b) { return {a.u - b.u, a.v - b.v}; }
    friend Quad operator*(const Quad &a, ll k)
    {
        ll uu = a.u * k;
        ll vv = a.v * k;
        return {(ll)uu, (ll)vv};
    }
    friend Quad operator*(ll k, const Quad &a) { return a * k; }
    friend bool operator==(const Quad &a, const Quad &b) { return a.u == b.u && a.v == b.v; }
    friend bool operator!=(const Quad &a, const Quad &b) { return !(a == b); }
};

struct Dissection
{
    int N, ex, ey, ez;
};

static Dissection dissect(int a, int b, int c)
{
    ll s = int_square(c);
    ll L;
    Quad X, Y, Z;
    if (s * s == c)
    {
        L = lcmll(lcmll(a, b), (int)s);
        X = {(ll)(L / a), 0};
        Y = {(ll)(L / b), 0};
        Z = {(ll)(L / s), 0};
    }
    else
    {
        L = lcmll(a, b);
        X = {(ll)(L / a), 0};
        Y = {(ll)(L / b), 0};
        Z = {0, (ll)L};
    }

    Quad S = X + Y + Z;
    Quad total{(ll)L, 0};

    ll m = 0;
    while ((S * (m + 1)).cmp(total, c) <= 0)
        ++m;

    Quad extra = total - S * m;
    Quad pieceX = X;
    Quad pieceY = X + Y;
    Quad pieceZ = S;

    ll elemX = 1, elemY = 1, elemZ = 1, elemE = 0;
    if (pieceX.cmp(extra, c) <= 0)
        ++elemE;
    if (pieceY.cmp(extra, c) <= 0)
        ++elemE;

    auto step = [&](const Quad &src) -> Quad
    {
        Quad dst;
        if (src.cmp(pieceX, c) <= 0)
        {
            dst = total - Y - Z - src;
        }
        else if (src.cmp(pieceY, c) <= 0)
        {
            dst = (total + X - Z) - src;
        }
        else
        {
            dst = (total + X + Y) - src;
        }
        dst = dst - S * m;
        dst = dst.norm_mod(S, c);
        return dst;
    };

    vector<Quad> starts = {pieceX, pieceY, pieceZ};
    for (const Quad &st : starts)
    {
        Quad src = st;
        for (;;)
        {
            Quad dst = step(src);
            if (dst == pieceX || dst == pieceY || dst == pieceZ)
                break;
            if (dst.cmp(pieceX, c) < 0)
                ++elemX;
            else if (dst.cmp(pieceY, c) < 0)
                ++elemY;
            else
                ++elemZ;
            if (dst.cmp(extra, c) <= 0)
                ++elemE;
            src = dst;
        }
    }

    ll N = m * (elemX + elemY + elemZ) + elemE;
    return {(int)N, (int)elemX, (int)elemY, (int)elemZ};
}

struct Cycle
{
    vector<int> nxt;
    vector<uint8_t> fl;
    int flipped = 0;

    Cycle() = default;

    Cycle(int n, int piece)
    {
        nxt.resize(n);
        fl.resize(n);
        flipped = 0;
        for (int i = 0; i < piece; i++)
        {
            fl[i] = 1;
            nxt[i] = n - 1 - i;
            flipped++;
        }
        for (int i = piece; i < n; i++)
        {
            fl[i] = 0;
            nxt[i] = i - piece;
        }
    }

    static Cycle fromOrbit(const vector<uint8_t> &flip)
    {
        int n = (int)flip.size();
        Cycle c;
        c.nxt.resize(n);
        c.fl = flip;
        c.flipped = 0;
        for (int i = 0; i < n; i++)
        {
            c.nxt[i] = (i + 1 == n ? 0 : i + 1);
            if (c.fl[i])
                c.flipped++;
        }
        return c;
    }

    Cycle operator*(const Cycle &other) const
    {
        int n = (int)nxt.size();
        Cycle r;
        r.nxt.resize(n);
        r.fl.resize(n);
        r.flipped = 0;
        for (int i = 0; i < n; i++)
        {
            int mid = nxt[i];
            int to = other.nxt[mid];
            uint8_t f = fl[i] ^ other.fl[mid];
            r.nxt[i] = to;
            r.fl[i] = f;
            if (f)
                r.flipped++;
        }
        return r;
    }
};

static int min_period_binary(const vector<uint8_t> &s)
{
    int n = (int)s.size();
    vector<int> l(n);
    for (int i = 1; i < n; i++)
    {
        int j = l[i - 1];
        while (j > 0 && s[i] != s[j])
            j = l[j - 1];
        if (s[i] == s[j])
            j++;
        l[i] = j;
    }
    int p = n - l[n - 1];
    if (p > 0 && n % p == 0)
        return p;
    return n;
}

static ll period_allup(const vector<uint8_t> &flip)
{
    int n = (int)flip.size();
    if (n == 0)
        return 1;
    int p = min_period_binary(flip);
    int parity = 0;
    for (int i = 0; i < p; i++)
        parity ^= (int)flip[i];
    return parity ? 2LL * p : 1LL * p;
}

struct OrbitInfo
{
    vector<uint8_t> flip;
    vector<uint8_t> target;
    ll per;
};

static ll hybrid_k(const Cycle &perm, int head)
{
    int N = (int)perm.nxt.size();
    vector<uint8_t> vis(N, 0);
    ll res = 0, mod = 1;
    vector<OrbitInfo> hard;

    for (int i = 0; i < N; i++)
        if (!vis[i])
        {
            vector<int> idx;
            int cur = i;
            while (!vis[cur])
            {
                vis[cur] = 1;
                idx.push_back(cur);
                cur = perm.nxt[cur];
            }
            int n = (int)idx.size();
            vector<uint8_t> f(n), t(n);
            bool any = false;
            for (int j = 0; j < n; j++)
            {
                int id = idx[j];
                f[j] = perm.fl[id];
                t[j] = (id >= head);
                any |= t[j];
            }
            ll per = period_allup(f);
            if (!any)
            {
                ll g = __gcd(mod, per);
                mod = mod / g * per;
            }
            else
            {
                hard.push_back({move(f), move(t), per});
            }
        }

    sort(hard.begin(), hard.end(), [](const OrbitInfo &a, const OrbitInfo &b)
         { return a.flip.size() < b.flip.size(); });

    for (auto &o : hard)
    {
        Cycle base = Cycle::fromOrbit(o.flip);
        Cycle cur = base;
        ll early = -1;
        for (ll k = 1; k <= o.per; k++)
        {
            bool ok = true;
            for (size_t i = 0; i < o.flip.size(); i++)
                if (cur.fl[i] != o.target[i])
                {
                    ok = false;
                    break;
                }
            if (ok)
            {
                early = k;
                break;
            }
            cur = cur * base;
        }
        if (early < 0)
            return 0;
        ll merged = CRT2(res, mod, early % o.per, o.per);
        if (merged < 0)
            return 0;
        ll g = __gcd(mod, o.per);
        mod = mod / g * o.per;
        res = merged;
    }

    res %= mod;
    if (res == 0)
        res = mod;
    return res;
}

struct Key
{
    int N, ex, ey, ez;
};
struct KeyHash
{
    size_t operator()(const Key &k) const noexcept
    {
        size_t h = (uint32_t)k.N;
        h = h * 1315423911u + (uint32_t)k.ex;
        h = h * 1315423911u + (uint32_t)k.ey;
        h = h * 1315423911u + (uint32_t)k.ez;
        return h;
    }
};
struct KeyEq
{
    bool operator()(const Key &a, const Key &b) const noexcept
    {
        return a.N == b.N && a.ex == b.ex && a.ey == b.ey && a.ez == b.ez;
    }
};

static ll F(int a, int b, int c, unordered_map<Key, ll, KeyHash, KeyEq> &cache)
{
    Dissection d = dissect(a, b, c);
    Key key{d.N, d.ex, d.ey, d.ez};
    auto it = cache.find(key);
    if (it != cache.end())
        return it->second;

    int N = d.N;
    Cycle cx(N, d.ex);
    Cycle cy(N, d.ey);
    Cycle cz(N, d.ez);

    Cycle perm0 = cx * cy * cz;

    vector<uint8_t> vis(N, 0);
    ll l = 1;
    for (int i = 0; i < N; i++)
        if (!vis[i])
        {
            vector<uint8_t> f;
            int cur = i;
            while (!vis[cur])
            {
                vis[cur] = 1;
                f.push_back(perm0.fl[cur]);
                cur = perm0.nxt[cur];
            }
            ll per = period_allup(f);
            ll g = __gcd(l, per);
            l = l / g * per;
        }

    ll ans = 3 * l;

    Cycle perm1 = cy * cz * cx;
    ll k1 = hybrid_k(perm1, N - d.ex);
    if (k1)
        ans = min(ans, 3 * k1 + 1);

    Cycle perm2 = cz * cx * cy;
    ll k2 = hybrid_k(perm2, N - d.ex - d.ey);
    if (k2)
        ans = min(ans, 3 * k2 + 2);

    cache.emplace(key, ans);
    return ans;
}

static ll G(int n)
{
    unordered_map<Key, ll, KeyHash, KeyEq> cache;
    cache.reserve(20000);
    ll sum = 0;
    for (int c = 11; c <= n; c++)
        for (int b = 10; b < c; b++)
            for (int a = 9; a < b; a++)
                sum += F(a, b, c, cache);
    return sum;
}

int main()
{
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout << G(N) << "\n";
    return 0;
}