Numbers : Faulhaber's Formula
|
$$ \sum_{k=1}^n k^p = 1^p+2^p+3^p+\cdots+n^p $$ $$ \sum n^m = \frac{1}{m+1}\left(B_0 n^{m+1}+{m+1 \choose 1}B_1^{+} n^m+{m+1 \choose 2}B_2 n^{m-1}+\cdots+{m+1 \choose m}B_m n \right) $$
Bernoulli number
|
$$ B_m^{+}=1- \sum_{k=0}^{m-1} {m \choose k} \frac{B_k^{+}}{m-k+1} $$
// -----------------------------------------------------------------------------
// Faulhaber Formula C#
// -----------------------------------------------------------------------------
static int[][] FiF(int nmax, int modulus)
{
int[][] fif = new int[2][];
fif[0] = new int[nmax + 1];
fif[0][0] = 1;
for (int i = 1; i <= nmax; i++)
fif[0][i] = (int)(((long)fif[0][i - 1] * i) % modulus);
long a = fif[0][nmax];
long b = modulus;
long p = 1;
long q = 0;
while (b > 0)
{
long c = a / b;
long d = a;
a = b;
b = d % b;
d = p;
p = q;
q = d - c * q;
}
fif[1] = new int[nmax + 1];
fif[1][nmax] = (int)(p < 0 ? p + modulus : p);
for (int i = nmax - 1; i >= 0; i--)
fif[1][i] = (int)(((long)fif[1][i + 1] * (i + 1)) % modulus);
return fif;
}
static long Choose(int n, int r, int[][] fif, int modulus)
{
if (n < 0 || r < 0 || r > n) return 0;
long factn = fif[0][n];
long invFactr = fif[1][r];
long invFactnr = fif[1][n - r];
long c1 = ((long)fif[0][n] * fif[1][r]) % modulus;
long c2 = (c1 * fif[1][n - r]) % modulus;
return ((((long)fif[0][n] * fif[1][r]) % modulus) * fif[1][n - r]) % modulus;
}
// -----------------------------------------------------------------------------
static long Mul(long a, long b, int modulus)
{
return (a * b) % modulus;
}
static long Div(long a, long b, int modulus)
{
return (a * Inv(b, modulus)) % modulus;
}
static long Pow(long a, long b, int modulus)
{
long res = 1;
while (b > 0)
{
if ((b & 1) == 1) res = (res * a) % modulus;
a = (a * a) % modulus;
b >>= 1;
}
return res;
}
static long Inv(long a, int modulus)
{
long b = modulus;
long p = 1, q = 0;
while (b > 0)
{
long c = a / b;
long d = a;
a = b;
b = d % b;
d = p;
p = q;
q = d - c * q;
}
return p < 0 ? p + modulus : p;
}
// -----------------------------------------------------------------------------
static int[] Bernoulli(int nmax, int[][] fif, int modulus)
{
int[] B = new int[nmax];
B[0] = 1;
for (int i = 1; i < nmax; i++)
{
long b = 0;
for (int j = 0; j < i; j++)
b = (b + Div(Mul(Choose(i, j, fif, modulus), B[j], modulus),
i - j + 1, modulus)) % modulus;
B[i] = (int)((1L + modulus - b) % modulus);
}
return B;
}
static int Faulhaber(int n, int p, int[] ber, int modulus)
{
long sum = 0;
if (n > 2)
{
long b = 1;
long n1 = (n - 1) % modulus;
for (int k = p; k >= 0; k--)
{
b = Mul(Div(Mul(b, k + 1, modulus), p - k + 1, modulus), n1, modulus);
sum = (sum + Mul(b, ber[k], modulus)) % modulus;
}
sum = Div(sum, p + 1, modulus);
}
return (int)sum;
}
// -----------------------------------------------------------------------------# ------------------------------------------------------------------------------
# Faulhaber's Formula Python
# ------------------------------------------------------------------------------
def FiF(nmax, modulus):
fif = [[0] * (nmax + 1) for i in range(2)]
fif[0][0] = 1
for i in range(1, nmax + 1):
fif[0][i] = int((fif[0][i - 1] * i) % modulus)
a, b = fif[0][nmax], modulus
p, q = 1, 0
while b > 0:
c, d = a // b, a
a, b = b, d % b
d = p
p, q = q, d - c * q
fif[1] = [0] * (nmax + 1)
fif[1][nmax] = int((p + modulus) if p < 0 else p)
for i in range(nmax - 1, -1, -1):
fif[1][i] = int((fif[1][i + 1] * (i + 1)) % modulus)
return fif
def Choose(n, r, fif, modulus):
if n < 0 or r < 0 or r > n: return 0
return int((((fif[0][n] * fif[1][r]) % modulus) * fif[1][n - r]) % modulus)
def Mul(a, b, modulus):
return (a * b) % modulus
def Div(a, b, modulus):
return (a * Inv(b, modulus)) % modulus
def Inv(a, modulus):
b = modulus
p, q = 1, 0
while b > 0:
c, d = a // b, a
a, b = b, d % b
d = p
p, q = q, d - c * q
return (p + modulus) if p < 0 else p
def Bernoulli(nmax, modulus):
FIF = FiF(nmax + 1, modulus)
B = [0] * nmax;
B[0] = 1
for i in range(1, nmax):
b = 0
for j in range(i):
b = (b + Div(Mul(Choose(i, j, FIF, modulus), B[j], modulus),
i - j + 1, modulus)) % modulus
B[i] = int((1 + modulus - b) % modulus)
return B
def Faulhaber(n, p, ber, modulus):
sum = 0
if n > 2:
b, n1 = 1, (n - 1) % modulus
for k in range(p, -1, -1):
b = Mul(Div(Mul(b, k + 1, modulus), p - k + 1, modulus), n1, modulus)
sum = (sum + Mul(b, ber[k], modulus)) % modulus
sum = Div(sum, p + 1, modulus)
return sum
# ------------------------------------------------------------------------------
if __name__ == "__main__":
pmax = 1010
modulus = 1000000009
ber = Bernoulli(pmax, modulus)
print(Faulhaber(10**18, 1000, ber, modulus))
# ------------------------------------------------------------------------------