Project Euler 311

发表于 2026-01-07 分类于 Project Euler 阅读次数：本文字数： 2.7k 阅读时长 ≈ 2 分钟

Project Euler 311

题目

Biclinic Integral Quadrilaterals

\(ABCD\) is a convex, integer sided quadrilateral with \(1 \le AB < BC < CD < AD\).

\(BD\) has integer length. \(O\) is the midpoint of \(BD\). \(AO\) has integer length.

We’ll call \(ABCD\) a biclinic integral quadrilateral if \(AO = CO \le BO = DO\).

For example, the following quadrilateral is a biclinic integral quadrilateral:

\(AB = 19, BC = 29, CD = 37, AD = 43, BD = 48\) and \(AO = CO = 23\).

Let \(B(N)\) be the number of distinct biclinic integral quadrilaterals \(ABCD\) that satisfy \(AB^2+BC^2+CD^2+AD^2 \le N\).

We can verify that \(B(10 000) = 49\) and \(B(1 000 000) = 38239\).

Find \(B(10 000 000 000)\).

解决方案

把 \(O\) 放在原点，令对角线 \(BD\) 在 \(x\) 轴上：那么有\(B=(s,0), D=(-s,0)\)，其中\(BO=DO=s=\dfrac{BD}{2}.\)由于 \(AB,AD\) 都是整数，我们立刻得到）：

\[AB^2+AD^2 = 2(AO^2+BO^2)=2(r^2+s^2).\]

左边为整数，因此 \(2s^2\) 必须为整数。若 \(BD\) 为奇数，则 \(s\) 为半整数，从而 \(2s^2\) 为半整数，矛盾。故\(BD\) 必须为偶数，且 \(s\in\mathbb Z\)。

设\(r=AO=CO, s=BO=DO.\)用三角形中线定理可得：对任意点 \(P\)，\(PB^2+PD^2 = 2(PO^2+BO^2).\)

代入 \(P=A\)，得到\(AB^2+AD^2 = 2(AO^2+BO^2)=2(r^2+s^2).\)

同理代入 \(P=C\)，得到\(BC^2+CD^2 = 2(CO^2+DO^2)=2(r^2+s^2).\)

相加得\(AB^2+BC^2+CD^2+AD^2 = 4(r^2+s^2).\)

因此题目条件\(AB^2+BC^2+CD^2+AD^2 \le N\)，等价于\(r^2+s^2 \le \dfrac N4.\)令\(L=\left\lfloor\dfrac{N}{4}\right\rfloor.\)

对任意整数 \(x<y\)，令\(u=\dfrac{x+y}{2}, v=\dfrac{y-x}{2},\)则\(x=u-v, y=u+v,\)并且\(x^2+y^2 = (u-v)^2+(u+v)^2 = 2(u^2+v^2).\)

在本题中：

顶点 \(A\) 对应 \((AB,AD)\)
顶点 \(C\) 对应 \((BC,CD)\)
对角线中点体系对应 \((AO,BO)=(r,s)\)

于是同一个数\(k=r^2+s^2\)会需要 3 种表示：\(k = r^2+s^2,k = u_A^2+v_A^2,k = u_C^2+v_C^2.\)

令 \(R(k)\) 表示：将 \(k\) 写成 \(a^2+b^2\) 的方式数（只计正整数，且 \(1\le a\le b\)）。

那么每个满足条件的 biclinic 四边形（在给定严格边序约束下）等价于从 \(R(k)\) 个分解里选出 3 个互不相同的分解来对应“底”和“两条臂”的组合（细节会带来对称修正，但最终落在同一个简洁公式上）。最终可化为：

\[B(N)=\sum_{k\le L} \binom{R(k)}{3}.\]

使用两平方和定理：若\(k\)的质因子分解满足：\(\displaystyle{k=2^e\prod_{p\equiv1\pmod4}p^{\alpha_p}\prod_{q\equiv3\pmod4}q^{2\beta_q}}\)，也就是说，对于所有\(q\equiv 3 \pmod4\)的质因子\(q\)，其对应的指数都满足\(\beta_q\)整数（对于其它质因子则不做要求），则存在平方和表示，并且无符号且有序的表示数为\(S=\{(a,b):a,b\ge 0\}\)：

\[ r_2(k)=\prod_{p\equiv1\pmod4}(\alpha_p+1). \]

考虑两个需要注意的解：

若 \(k\) 是平方数：\((\sqrt k,0),(0,\sqrt k)\)在 \(S\)，
若 \(k=2t^2\)：\((t,t)\)，那么在 \(S\)。

因此，令示性函数\(\sigma_1(k)\)表示\(k\)是平方数，示性函数\(\sigma_2(k)\)表示\(k\)是满足\(k=2t^2\)，那么有

\[R(k)=R_{\ge0}(k)-\sigma_1(k)-\sigma_2(k)=\frac{r_2(k)-\sigma_1(k)-\sigma_2(k)}{2}.\]

令 \(L=N/4\)。我们对 \(k\) 的因子结构拆分为两部分：

\(m\)：由所有 \(p\equiv1\pmod4\) 的素因子（及其幂）构成
剩余部分：只允许\(2^e\)以及所有 \(q\equiv3\pmod4\) 以偶次幂出现（即 \(q^{2u}\)）

将“剩余部分可取的个数”预处理为一个前缀计数表 \(w[x]\)，表示：不超过 \(x\) 的数中，形如 \(4^a\cdot s^2\)（其中 \(s\) 只含 \(q\equiv3\pmod4\) 的素因子）的数量。注意，这里不是取\(2^e\)是希望通过两次查询（\(4^a\cdot s^2,4^a\cdot 2 \cdot s^2\)），防止占用空间过大。

然后 DFS 枚举 \(m\) 的素因子幂次组合，计算贡献即可：

主项：\(w\left[\left\lfloor \dfrac{L}{m}\right\rfloor\right]\cdot \dbinom{R(m)}{3}\)
\(2\) 奇次：\(w\left[\left\lfloor \dfrac{L}{2m}\right\rfloor\right]\cdot \dbinom{R(m)}{3}\)

为了让贡献尽可能尽早出现（也就是\(\dbinom{R(k)}{3}> 0\)），这就意味着\(R(k)\ge 3\)，最小满足贡献的\(M\)是\(5^2\cdot 13=325\)，因此剩余部分的积不超过\(\dfrac{L}{325}\)。

代码

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

// --------- 全局变量 ----------
// 题目要求计算 B(N)，这里默认 N = 1e10（可自行修改）
ll N = 1e10;

// L = N/4（推导中会出现 AB^2+BC^2+CD^2+AD^2 <= N  ⇔  某个 k <= N/4）
ll L;

// Kmax = L/325（只需要预处理到这个上界，原因是产生贡献的最小 m 为 325）
int Kmax;

// w_pref[x] = w[x] 的前缀计数（见 build_w_prefix 的定义）
vector<uint32_t> w_pref;

// primes1 = 所有 p ≡ 1 (mod 4) 的素数列表（用于 DFS 枚举 m 的 1(mod4) 部分）
vector<int> primes1;

// 最终答案累加器
ll ans;

// --------- 只筛奇数的素数表 ----------
vector<uint8_t> sieve_odd(int limit)
{
    // odd-only sieve up to limit (inclusive).
    // 返回 s，使得对奇数 n>=3，s[n//2] = 1 当且仅当 n 为素数。
    //
    // 说明：
    // - 只存奇数：n 映射到下标 n/2（整除）。
    // - 2 不在此表中（后续只关心奇素数）。
    if (limit < 2)
        return {};
    int size = limit / 2 + 1;
    vector<uint8_t> s(size, 1);
    s[0] = 0; // n=1 对应下标 0，不是素数

    int r = (int)std::sqrt((long double)limit);
    for (int p = 3; p <= r; p += 2)
    {
        if (!s[p / 2])
            continue;

        // 标记 p 的奇数倍合数：
        // 从 p^2 开始，步长为 2p（跳过偶数倍）
        ll step = 2LL * p;
        ll start = 1LL * p * p;
        for (ll x = start; x <= limit; x += step)
            s[(int)(x / 2)] = 0;
    }
    return s;
}

vector<int> primes_mod4_3_upto(int limit)
{
    // 生成 <=limit 的所有 p ≡ 3 (mod 4) 的素数（含 3）
    auto s = sieve_odd(limit);
    vector<int> res;
    for (int p = 3; p <= limit; p += 4)
    {
        if (s[p / 2])
            res.push_back(p);
    }
    return res;
}

// --------- 预处理 w[x] ----------
vector<uint32_t> build_w_prefix(int K)
{
    // w[x] = #{ t <= x : t = 4^a * s^2, 其中 s 的素因子全满足 q ≡ 3 (mod 4) } 的前缀计数
    //
    // 解释（对应数学推导中的“剩余部分”）：
    // - 我们把 k 的因子拆成 m * t，其中 m 只含 p≡1(mod4) 的素因子。
    // - 剩余部分 t 允许含 2^e 和所有 q≡3(mod4) 的偶次幂 q^{2u}。
    // - 将 e 的偶次部分写成 4^a，奇次部分额外乘一个 2，因此我们预处理的表只数 4^a*s^2；
    //   后续用两次查询（x 和 x/2）覆盖 t=4^a*s^2 与 t=2*4^a*s^2 两类情况。
    if (K <= 0)
        return vector<uint32_t>(1, 0);

    // 只需要枚举 s <= sqrt(K)，然后用 s^2 * 4^a 扩展到 <=K
    int sqrtK = (int)int_square(K);

    // primes3 = 所有 <=sqrtK 且 ≡3(mod4) 的素数
    vector<int> primes3 = primes_mod4_3_upto(sqrtK);

    // squares 存储所有可行的 s^2（s 只由 primes3 的幂构成，且 s<=sqrtK）
    vector<int> squares;
    squares.reserve(1 << 16);

    // DFS 枚举所有 s：
    // - curr 是当前 s
    // - idx 表示下一次可用的素数下标（保证素数严格递增，避免重复）
    function<void(int, int)> gen_s = [&](int curr, int idx)
    {
        squares.push_back(curr * curr);

        for (int j = idx; j < (int)primes3.size(); ++j)
        {
            int p = primes3[j];
            ll x = 1LL * curr * p;
            if (x > sqrtK)
                break;

            // 枚举 p 的不同指数：curr*p, curr*p^2, ...
            while (x <= sqrtK)
            {
                gen_s((int)x, j + 1);
                x *= p;
            }
        }
    };
    gen_s(1, 0);

    // mark[t] = 1 表示 t 可写成 4^a*s^2（s 的素因子全为 3(mod4)）
    vector<uint8_t> mark(K + 1, 0);
    for (int b : squares)
    { // b = s^2
        ll x = b;
        while (x <= K)
        {
            mark[(int)x] = 1;
            x *= 4; // 乘 4^a
        }
    }

    // 构造前缀和 w[x] = sum_{t<=x} mark[t]
    vector<uint32_t> w(K + 1, 0);
    uint32_t acc = 0;
    for (int i = 1; i <= K; ++i)
    {
        acc += mark[i];
        w[i] = acc;
    }
    return w;
}

// --------- 组合数 C(n,3) ----------
inline ll C3(ll n)
{
    // 若 n<3，此式自然为 0（但这里上游会确保 n>=0）
    return n * (n - 1) * (n - 2) / 6;
}

// --------- DFS 枚举 m（只含 p≡1(mod4) 的部分） ----------
void dfs(ll m, int idx, ll prod, int e)
{
    // m: 当前 p≡1(mod4) 部分的数值（即 m = ∏ p^{α_p}）
    // idx: 当前正处理 primes1[idx] 这个素数 p
    // prod: 已确定的之前素因子的 ∏(α+1)（不包含当前 p 的贡献）
    // e: 当前素数 p 的指数 α（>=1）
    //
    // 因此当前时刻：
    //   P = ∏(α_p+1) = prod * (e+1)
    // 该 P 决定了 R(k)（从两平方和定理推导而来，且对固定 m 与允许的 t 无关）
    ll P = prod * (e + 1); // P = ∏_{p≡1(mod4)} (α_p+1)

    // 只有当 P 足够大时，C(R,3) 才非 0；代码用 P>=5 作为门槛（P=5 也会得到 0）
    if (P >= 5)
    {
        // X = floor(L/m) ：t=4^a*s^2 这类剩余部分的上界
        ll X = L / m;

        // 对应 2 的指数为偶数：R = floor(P/2)
        ll R_even2 = P / 2;

        // 对应 2 的指数为奇数：
        // - 当 P 为奇数时，会导致取 ceil(P/2)（等价于对退化情形的一次“补偿”）
        // - 当 P 为偶数时，与偶数情形一致
        ll R_odd2 = (P & 1) ? ((P + 1) / 2) : R_even2;

        // t=4^a*s^2 的数量是 w[X]
        ans += (ll)w_pref[X] * C3(R_even2);

        // t=2*4^a*s^2 等价于 4^a*s^2 <= X/2，因此数量为 w[X/2]
        ans += (ll)w_pref[X / 2] * C3(R_odd2);
    }

    // 当前素数 p
    int p = primes1[idx];

    // 1) 尝试把当前 p 的指数再加 1：m *= p, e++
    ll mp = m * (ll)p;
    if (mp <= L)
        dfs(mp, idx, prod, e + 1);

    // 2) 固定当前 p 的指数 e，转入更大的素数 q（保持素数递增，避免重复）
    ll prod2 = prod * (e + 1);
    for (int j = idx + 1; j < (int)primes1.size(); ++j)
    {
        int q = primes1[j];
        ll mq = m * (ll)q;
        if (mq > L)
            break;
        dfs(mq, j, prod2, 1);
    }
}

// --------- 主函数：计算 B(N) ----------
ll B(ll N)
{
    // L = N/4
    L = N / 4;
    if (L <= 0)
        return 0;

    // 只有当 P=∏(α+1) 足够大时才会产生贡献；
    // 最小能产生贡献的 m 为 5^2*13=325，因此剩余部分上界为 L/325
    Kmax = (L >= 325) ? (int)(L / 325) : (int)L;

    // 预处理 w[x]（前缀计数）
    w_pref = build_w_prefix(Kmax);

    // primes1 只需要筛到 L/25：
    // 因为要让 m 可能产生贡献，至少得有 5^2（从而 m >= 25）
    int Pmax = (L >= 25) ? (int)(L / 25) : (int)L;
    auto odd_sieve = sieve_odd(Pmax);

    primes1.clear();
    for (int p = 5; p <= Pmax; p += 4)
    {
        // 只取 p≡1(mod4) 的素数
        if (odd_sieve[p / 2])
            primes1.push_back(p);
    }

    // 立方根界剪枝：
    // 这里沿用常见的枚举策略：从“至少有一个平方因子 p^2”起步，
    // 只需要起步素数 p <= cbrt(L)（否则 p^2 已经太大，无法再乘其它因子形成更复杂的 m）
    ll cb = (ll)cbrt((long double)L);
    while ((cb + 1) * (cb + 1) * (cb + 1) <= L)
        ++cb;
    while (cb * cb * cb > L)
        --cb;

    ans = 0;

    // 启动 DFS 的两种起点：
    //
    // (1) m = p^2（当前 p 的指数 e=2，prod=1）
    //     即 α_p = 2
    //
    // (2) m = p*q（p<q，两个素数指数都为 1），先固定 p 的贡献 (α_p+1)=2 进入 prod，
    //     然后从 q 作为当前素数开始递归
    for (int i = 0; i < (int)primes1.size(); ++i)
    {
        int p = primes1[i];
        if (p > cb)
            break;

        // 起点 (1)
        dfs(1LL * p * p, i, 1, 2);

        // 起点 (2)
        // 为保证还能再乘下去（形成至少一个平方项 q^2），通常限制 q <= sqrt(L/p)
        ll limit_q = int_square(L / p);
        for (int j = i + 1; j < (int)primes1.size(); ++j)
        {
            int q = primes1[j];
            if (q > limit_q)
                break;

            // 额外做一个与原 Python 同型的剪枝条件
            if (1LL * p * q * q <= L)
                dfs(1LL * p * q, j, 2, 1); // prod=2 表示已包含 (α_p+1)=2
        }
    }

    return ans;
}

int main()
{
    // 输出 B(1e10)
    cout << B(N) << "\n";
    return 0;
}