UVa 766 Sum of Powers

题目描述

一位学生想要计算如下形式的和：

$$S_k(n)=\sum _{i=1}^n{i}^k$$

其中 $k$ 是一个固定的自然数，$n$ 是变化的自然数。他发现，如果对每个 $i$ 计算 $i^k$ 再求和，当 $n$ 很大时计算量太大（需要 $O(n)$ 次运算）。实际上，存在一种方法只需常数次运算（不随 $n$ 变化）。

可以证明，$S_k(n)$ 是一个关于 $n$ 的 $k+1$ 次多项式，其系数为有理数，即：

$$S_k(n)=\frac{1}{M}\left(a_{k+1}{n}^{k+1}+a_k{n}^k+\dots +a_{1}n+a_0\right)$$

其中 $M,a_{k+1},a_k,\dots ,a_1,a_0$ 都是整数。

我们要求 $M$ 为正整数且尽可能小。在此条件下，这些整数是唯一确定的。你需要编写一个程序，对于给定的 $k$，找到这样一组系数，以便学生能快速计算。

输入格式

输入文件包含多组测试数据，每组数据包含一个整数 $k$（$0\le k\le 20$）。
第一行是数据集的数量，接着是一个空行。每组数据之间也有一个空行。

输出格式

对于每组数据，按顺序输出整数 $M,a_{k+1},a_k,\dots ,a_1,a_0$，数字之间用一个空格隔开。确保 $M$ 为正且尽可能小。
每组输出之间用一个空行隔开。

示例输入

1
2

示例输出

6 2 3 1 0

题目分析

问题本质

我们需要找到多项式 $S_k(n)$ 的有理系数表示形式。已知：

1. $S_k(n)$ 是 $n$ 的 $k+1$ 次多项式
2. 可以表示为 $S_k(n)=\frac{1}{M}{\sum }_{j=0}^{k+1}{a}_j{n}^j$
3. 要求 $M$ 为正整数且最小
4. $a_j$ 为整数，且 $a_0=0$（因为 $S_k(0)=0$）

数学背景

这个问题涉及到幂和公式（$\text{Faulhaber’s formula}$）和伯努利数（$\text{Bernoulli numbers}$）。具体公式为：

$$S_k(n)=\frac{1}{k+1}\sum _{j=0}^k(-1{)}^j(\genfrac{}{}{0ex}{}{k+1}{j})B_j{n}^{k+1-j}=\frac{1}{k+1}\sum _{j=0}^k(\genfrac{}{}{0ex}{}{k+1}{j})B_j^{+}{n}^{k+1-j}$$

$$B_j^{+}=(-1{)}^j{B}_j$$

其中 $B_j$ 是伯努利数。

伯努利数定义

伯努利数 $B_j$ 可以通过以下递推关系定义：

• $B_0=1$
• 对于 $m\ge 1$：${\sum }_{j=0}^m\left(\genfrac{}{}{0ex}{}{m+1}{j}\right)B_j=0$

由此可得：

• $B_1=-\frac{1}{2}$
• $B_2=\frac{1}{6}$
• $B_3=0$
• $B_4=-\frac{1}{30}$
• …

但是，题目示例显示对于 $k=2$，输出是 $62310$，对应多项式：
$$S_2(n)=\frac{1}{6}(2n^3+3n^2+n)$$

标准公式给出的是：
$$S_2(n)=\frac{1}{3}{n}^3+\frac{1}{2}{n}^2+\frac{1}{6}n=\frac{1}{6}(2n^3+3n^2+n)$$

是一致的。

解题思路

方法选择

我们选择使用伯努利数方法，因为：

1. 效率高：时间复杂度为 $O(k^2)$，对于 $k\le 20$ 完全足够
2. 数值稳定：使用分数运算避免浮点数精度问题
3. 实现简单：代码相对简洁

算法步骤

1. 预处理组合数：使用杨辉三角计算 $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$，存储在数组中
2. 计算伯努利数 $B_j^{+}$：
   • 使用递推公式：${\sum }_{j=0}^m\left(\genfrac{}{}{0ex}{}{m+1}{j}\right)B_j=0$
   • 设置 $B_0=1$，$B_1=\frac{1}{2}$
   • 使用分数运算（分子/分母对）
3. 计算多项式系数：
   • 对于 $j=0,1,\dots ,k$：$a_{k+1-j}=\left(\genfrac{}{}{0ex}{}{k+1}{j}\right)B_j$
   • 所有系数乘以 $k+1$：得到 ${\tilde{a}}_j=(k+1)\cdot a_j$
4. 通分：
   • 找到所有 ${\tilde{a}}_j$ 分母的最小公倍数 $M$
   • 将所有系数乘以 $M$ 得到整数系数
5. 约分：
   • 找到所有系数和 $M$ 的最大公约数
   • 进行约分，确保 $M$ 最小
6. 符号处理：
   • 确保最高次系数 $a_{k+1}>0$
   • 确保 $M>0$
7. 处理特殊情况：
   • $k=0$ 时：$S_0(n)=n=1\cdot n+0$，直接输出 $110$

时间复杂度

• 计算组合数：$O(k^2)$
• 计算伯努利数：$O(k^2)$
• 总时间复杂度：$O(k^2)$，对于 $k\le 20$ 非常高效

空间复杂度

• 存储组合数：$O(k^2)$
• 存储伯努利数和系数：$O(k)$
• 总空间复杂度：$O(k^2)$

关键点

1. 伯努利数的符号约定：必须使用 $B_1^{+}=\frac{1}{2}$ 才能得到与题目示例一致的输出
2. 分数运算：使用分数表示避免精度问题，每次运算后立即约分
3. 边界情况：$k=0$ 需要特殊处理
4. 输出格式：注意组间空行

代码实现

// Sum of Powers
// UVa ID: 766
// Verdict: Accepted
// Submission Date: 2025-12-01
// UVa Run Time: 0.010s
//
// 版权所有（C）2025，邱秋。metaphysis # yeah dot net

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

const int MAXK = 21;

long long comb[MAXK + 2][MAXK + 2];

void initComb() {
    for (int i = 0; i <= MAXK + 1; ++i) {
        comb[i][0] = comb[i][i] = 1;
        for (int j = 1; j < i; ++j)
            comb[i][j] = comb[i - 1][j - 1] + comb[i - 1][j];
    }
}

// 使用递推计算伯努利数 B_j^+
vector<pair<long long, long long>> bernoulliPlus(int k) {
    vector<pair<long long, long long>> B(k + 1); // 存储分数 num/den
    
    B[0] = {1, 1}; // B0 = 1/1
    
    for (int m = 1; m <= k; ++m) {
        // 计算 sum_{j=0}^{m-1} C(m+1, j) * B_j
        pair<long long, long long> sum = {0, 1};
        for (int j = 0; j < m; ++j) {
            // 分数加法：a/b + c/d = (ad + bc)/(bd)
            long long a = comb[m + 1][j] * B[j].first;
            long long b = B[j].second;
            
            long long newDen = sum.second / __gcd(sum.second, b) * b;
            long long newNum = sum.first * (newDen / sum.second) + a * (newDen / b);
            
            long long g = __gcd(abs(newNum), newDen);
            sum = {newNum / g, newDen / g};
        }
        
        // B_m = -sum / (m+1)
        long long num = -sum.first;
        long long den = sum.second * (m + 1);
        long long g = __gcd(abs(num), den);
        B[m] = {num / g, den / g};
    }
    
    // 修正 B1^+ = +1/2（而不是 -1/2）
    B[1] = {1, 2};
    
    return B;
}

void solve(int k, vector<long long>& coeff, long long& M) {
    vector<pair<long long, long long>> B = bernoulliPlus(k);
    
    // 计算系数
    vector<pair<long long, long long>> polyCoeff(k + 2, {0, 1});
    
    for (int j = 0; j <= k; ++j) {
        long long binomial = comb[k + 1][j];
        
        // 分数乘法：(a/b) * c = (a*c)/b
        long long num = binomial * B[j].first;
        long long den = B[j].second;
        long long g = __gcd(abs(num), den);
        num /= g;
        den /= g;
        
        // 存储到对应位置
        polyCoeff[k + 1 - j] = {num, den};
    }
    
    // 所有系数乘以 (k+1) 得到 a_j/(k+1) 形式
    for (int j = 0; j <= k + 1; ++j) {
        polyCoeff[j].second *= (k + 1);
        long long g = __gcd(abs(polyCoeff[j].first), polyCoeff[j].second);
        if (g != 0) {
            polyCoeff[j].first /= g;
            polyCoeff[j].second /= g;
        }
    }
    
    // 找到所有分母的最小公倍数
    M = 1;
    for (int j = 0; j <= k + 1; ++j) {
        if (polyCoeff[j].first != 0) {
            M = M / __gcd(M, polyCoeff[j].second) * polyCoeff[j].second;
        }
    }
    
    // 转换为整数系数
    coeff.resize(k + 2);
    for (int j = 0; j <= k + 1; ++j) {
        coeff[j] = polyCoeff[j].first * (M / polyCoeff[j].second);
    }
    
    // 约分
    long long g = 0;
    for (int j = 0; j <= k + 1; ++j) {
        g = __gcd(g, abs(coeff[j]));
    }
    g = __gcd(g, M);
    
    if (g > 1) {
        M /= g;
        for (int j = 0; j <= k + 1; ++j) {
            coeff[j] /= g;
        }
    }
    
    // 确保最高次系数为正
    if (coeff[k + 1] < 0) {
        M = -M;
        for (int j = 0; j <= k + 1; ++j) {
            coeff[j] = -coeff[j];
        }
    }
    
    // 确保 M 为正
    if (M < 0) {
        M = -M;
        for (int j = 0; j <= k + 1; ++j) {
            coeff[j] = -coeff[j];
        }
    }
}

int main() {
    initComb();
    int t;
    cin >> t;
    
    for (int tc = 0; tc < t; ++tc) {
        if (tc > 0) cout << endl;
        
        int k;
        cin >> k;
        
        if (k == 0) {
            cout << "1 1 0\n";
            continue;
        }
        
        vector<long long> coeff;
        long long M;
        solve(k, coeff, M);
        
        cout << M;
        for (int j = k + 1; j >= 0; --j) {
            cout << " " << coeff[j];
        }
        cout << endl;
    }
    
    return 0;
}

总结

本题的关键在于理解幂和公式与伯努利数的关系，并注意伯努利数的符号约定。通过使用分数运算和递推计算伯努利数，我们可以高效地求解多项式系数。算法的时间复杂度和空间复杂度都在可接受范围内，能够正确处理 $k\le 20$ 的所有情况。

需要注意的细节包括：

1. 处理 $k=0$ 的特殊情况
2. 使用正确的伯努利数约定（$B_1^{+}=\frac{1}{2}$）
3. 确保输出格式正确（组间空行）
4. 进行通分和约分，使 $M$ 最小

这个解决方案展示了如何将数学理论转化为高效算法，是数学与计算机科学结合的典型例子。