月份：2019年8月

题目链接：HDU-6561

给不超过 10^4 个阻值为 1 \Omega 的电阻，通过串联和并联使得最后形成的电路有效阻值为 r。

这道题可以考虑有理逼近中的连分数。

a_0 + \dfrac{1}{a_1 + \dfrac{1}{a_2 + \dfrac{1}{\ddots + \dfrac{1}{a_n}}}}

我们通常将它记为 [a_0, a_1, a_2, \cdots, a_n]。

记 \frac{p_k}{q_k} 为 [a_0, a_1, \cdots, a_k]，我们称 \frac{p_k}{q_k} 为 [a_0, a_1, a_2, \cdots, a_n] 的第 k 个渐进分数。

可以证明，所有有理数都可以表示为这样的一个有限连分数。

考虑斐波那契的 \phi = \frac{1 + \sqrt5}{2}，我们可以令连分数序列为无限的 [1, 1, 1, \cdots, 1, \cdots]。而且我们可以发现，\frac{p_k}{q_k} 正好等于 \frac{F_{k+1}}{F_k}。数论中还会研究，任何无理数也都可以用连分数来逼近。

可以看出来，将这样的一个有理数摆进来，首先将它的整数部分去掉，然后再求倒数，就变成了 [a_1, a_2, \cdots, a_n]，问题规模得以缩小。当我们迭代到没有小数部分的时候，这个有理数已经被精确表示了。类似于一个辗转相除法的例子，对不对？时间复杂度为 O(\log q)。

我们需要(cai)知(chu)道(lai)的是，渐进级数越精确，误差越小。如果我们用某个连分数来表示这道题的答案呢？

一定有 a_0 = 0。那么我们不妨设定一个数列 \lbrace b_i \rbrace，其中

b_i = \begin{cases} \frac{1}{a_i + \frac{1}{b_{i+1}}}, & i = 2k+1 \\ a_i + b_{i+1}, & i=2k \end{cases}

那么可以认为 b_i 是由某一些电阻构成的复合电阻。当 i 为奇数时，我们将 a_i 个 1 \Omega 的电阻与 b_{i+1} 并联起来；当 i 为偶数时，我们将 a_i 个 1 \Omega 的电阻与 b_{i+1} 串联起来。注意当超过连分数长度 n 时，b_{i+1} 为 0 或 \infty，取决于奇偶性。

首先题目给定的 r 输入进来以后一定是一个有理数。那我们不妨先求出 [ r – \epsilon, r + \epsilon ] 中分母最小的分数 \frac{p_i}{q_i}。

我们将 \frac{p_i}{q_i} 的连分数表示出来。但是由于题目要求所用电阻数量不能超过 10^4，我们可以将超过的部分略去，留下一堆多余的电阻，这样比略去的逼近效果略好。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll EPS = 10000000;
const int MAXN = 10000;
int kase;

ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; }

void euclid(ll a, ll b, ll x, ll y, ll &p, ll &q)
{
    ll K = a / b; a -= K * b, x -= K * y;
    if (x > y) p = q = 1; else euclid(y, x, b, a, q, p);
    p += K * q;
}

void solve()
{
    static char buf[1000];
    static ll lfs[1000];
    scanf("%s", buf);
    ll fz = 0, fm = 1, fzm, fmm, fzM, fmM;

    for (int i = 2; buf[i]; i++)
    {
        fz = fz * 10 + buf[i] - '0';
        fm *= 10;
    }

    if (fm < EPS) fz *= EPS / fm, fm = EPS;
    fzm = fz - fm / EPS, fmm = gcd(fzm, fm), fzm /= fmm, fmm = fm / fmm;
    fzM = fz + fm / EPS, fmM = gcd(fzM, fm), fzM /= fmM, fmM = fm / fmM;
    euclid(fzm, fmm, fzM, fmM, fz, fm);
    if (fmM < fm) fz = fzM, fm = fmM;
    if (fmm < fm) fz = fzm, fm = fmm;
    int tot = 1, sum = 0, idx = 0, idg, id1, id2, idq = -1, op;

    while (fz && sum < MAXN)
    {
        swap(fz, fm);
        sum += (lfs[tot++] = min(fz / fm, MAXN - sum + 0ll));
        fz %= fm;
    }

    // printf("[ 0"); for (int i = 1; i < tot; i++) printf(", %lld", lfs[i]); printf(" ]\n");
    printf("Case %d:\n%d %d\n", ++kase, sum, sum-1);

    idg = sum-1;
    for (int j = tot-1; j > 0; j--)
    {
        op = j & 1;
        id1 = idx++;

        for (int i = 1; i < lfs[j]; i++)
        {
            id2 = idx++;
            printf("%d %d %d\n", op, id1, id2);
            id1 = ++idg;
        }

        if (~idq)
        {
            id2 = idq;
            printf("%d %d %d\n", op, id1, id2);
            ++idg;
        }

        idq = idg;
    }
}

int main()
{
    int T; scanf("%d", &T);
    while (T--) solve();
    return 0;
}

这道题总而言之感觉还是有点怪怪的，例如 0.000001 一定是无解的，但是竟然也能 AC。就这样吧，应该是良心出题人没有出卡人数据。

题目链接：HDU 6632

已知数字 a,b,p，并且 p=2^a3^b+1。想求解离散对数 a^x \equiv b \mod p。

需要一个东西叫做 Pohlig-Hellman 算法。以下是来自维基百科内容的翻译，将使用群论语言。

p^e 阶循环群

这个算法用于计算 p^e 阶循环群 G=\langle g \rangle 的离散对数 g^x = h。

1. 初始化 x_0 := 0

2. 计算 \gamma := g^{p^{e-1}}，由拉格朗日配集分解计数定理知，ord~\gamma = p。

3. 对于 k \in \lbrace 0,\cdots,e-1 \rbrace，分别做

   1. 计算 h_k := (g^{-x_k}h)^{p^{e-1-k}}。由此构造可知，此元素的阶一定整除 p，即阶为 1 或 p，因此 h_k \in \langle \gamma \rangle。

   2. 使用 BSGS 算法，计算 d_k \in \lbrace 0, \cdots, p-1 \rbrace，其中 \gamma^{d_k} = h_k。这个操作的复杂度为 O(\sqrt{p})

   3. 令 x_{k+1} := x_k + p^k d_k

4. 返回 x_e

当 e \ll p 时，这个算法优于 BSGS，复杂度为 O(e\sqrt{p})。

一般合数阶循环群

首先我们将 G = \langle g \rangle 的阶数刻画为 n = \prod_{i=1}^r {p_i}^{e_i}。现在我们解方程 g^x = h。

1. 对循环群阶数做因数分解

2. 对于 i \in \lbrace 1, \cdots, r \rbrace，分别做

   1. 计算 g_i := g^{n / {p_i}^{e_i}}，由拉格朗日定理知此元素的阶为 {p_i}^{e_i}。

   2. 计算 h_i := h^{n/{p_i}^{e_i}}，由此构造知 h_i \in \langle g_i \rangle

   3. 利用上面算法计算 x_i 使得 {g_i}^{x_i} = h_i

3. 利用中国剩余定理解方程组 x \equiv x_i \mod {p_i}^{e_i}

4. 返回 x

因此总的时间复杂度为

T(n) = O\left( \sum_{i=1}^r e_i (\log n + \sqrt{p_i}) \right)

回到本题

显然因为 p_i = 2, 3，我们不用 BSGS……好的。

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;

lld fpow(lld _a, lld k, lld MOD) {
    __int128 b = 1, a = _a % MOD;
    while (k) {
        if (k & 1) b = b * a % MOD;
        a = a * a % MOD;
        k >>= 1;
    }
    return b;
}

__int128 extgcd(__int128 a, __int128 b, __int128 &x, __int128 &y) {
    if (!b) return x = 1, y = 0, a;
    __int128 d = extgcd(b, a % b, y, x);
    return y -= (a / b) * x, d;
}

int T, e2, e3, G;
lld MOD, n, p2, p3;

lld solvePE(lld g, lld h, lld p, lld pe, int e)
{
    lld x = 0, pp = 1;
    lld y = fpow(g, pe/p, MOD);

    for (int k = 0; k < e; k++)
    {
        __int128 hk = fpow(g, pe-x, MOD);
        hk = hk * h % MOD;
        hk = fpow(hk, pe/p/pp, MOD);
        __int128 d, yy;
        for (d = 0, yy = 1; d < p; d++, yy = yy * y % MOD)
            if (yy == hk) break;
        x = (x + pp * d) % MOD, pp *= p;
    }

    return x;
}

__int128 solveN(lld h)
{
    if (e3 == 0) return solvePE(G, h, 2, p2, e2);
    else if (e2 == 0) return solvePE(G, h, 3, p3, e3);
    __int128 x2 = solvePE(fpow(G, p3, MOD), fpow(h, p3, MOD), 2, p2, e2);
    __int128 x3 = solvePE(fpow(G, p2, MOD), fpow(h, p2, MOD), 3, p3, e3);
    __int128 invp2, invp3; assert(extgcd(p2, p3, invp2, invp3) == 1);
    return (x2 * p3 * invp3 + x3 * p2 * invp2) % (MOD-1);
}

int main()
{
    cin >> T;

    while (T--)
    {
        long long a, b;
        cin >> MOD >> a >> b;
        n = MOD-1, e2 = e3 = 0, p2 = p3 = 1;
        while (n % 2 == 0) n /= 2, p2 *= 2, e2++;
        while (n % 3 == 0) n /= 3, p3 *= 3, e3++;
        for (G = 2; (e2 > 0 && fpow(G, (MOD-1)/2, MOD) == 1)
                 || (e3 > 0 && fpow(G, (MOD-1)/3, MOD) == 1); G++);
        __int128 inda = solveN(a), indb = solveN(b), x, y;
        __int128 d = extgcd(inda, MOD-1, x, y);
        lld ans = indb % d ? -1
            : lld((indb / d * x % (MOD-1) + (MOD-1)) % ((MOD-1)/d));
        cout << ans << endl;
    }

    return 0;
}

我竟然都忘了 extgcd 是怎么求多解的，我好菜啊QAQ