Codeforces D. The Sum of the k-th Powers(拉格朗日插值)

2023-10-25 12:50

本文主要是介绍Codeforces D. The Sum of the k-th Powers(拉格朗日插值),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

题目描述:

The Sum of the k-th Powers

time limit per test

2 seconds

memory limit per test

256 megabytes

input

standard input

output

standard output

There are well-known formulas: img, img, img. Also mathematicians found similar formulas for higher degrees.

Find the value of the sum img modulo 109 + 7 (so you should find the remainder after dividing the answer by the value 109 + 7).

Input

The only line contains two integers n, k (1 ≤ n ≤ 109, 0 ≤ k ≤ 106).

Output

Print the only integer a — the remainder after dividing the value of the sum by the value 109 + 7.

Examples

Input

Copy

4 1

Output

Copy

10

Input

Copy

4 2

Output

Copy

30

Input

Copy

4 3

Output

Copy

100

Input

Copy

4 0

Output

Copy

4

思路:

我从来没见过这样的题,今天是开了眼界了。听说是数学分析上的拉格朗日插值多项式(当时学的时候迷迷糊糊的)现在竟然就直接跳过理论来实际应用了。还是写代码,虽然知道是数值分析上的,本来就和计算机挂钩,但还是第一次直观感受到了数学在计算机的应用。以前学的时候老是疑惑我学了这些能干嘛?都说是应用广,但对我而言,应用只是一句空洞的话,现在看到了现实,还是有点小兴奋嗷~

在网上找了找资料,先摆一摆公式,以后弄明白了可能会回来补(可能吧)

拉格朗日插值公式

在数值分析中,拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律,而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测,在若干个不同的地方得到相应的观测值,拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日(插值)多项式。数学上来说,拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。——摘自百度百科

\[f(x)=\sum_{i=0}^n y_i\prod_{j!=i}\frac{x-x_j}{x_i-x_j}\].

简单的解释就是带入一个\(x_i\),对第\(i\)项来说后面的连乘是1,而对于非\(i\)项来说,后面的连乘是零。所以有\(f(x_i)=y_i\).

一种特殊的情况:\(x_i\)是连续的。

\(P_i=\prod_{j=0}^i(x-x_j)\)\(S_i=\prod_{j=i}^n(x-x_j)\).

而由于\(x_i\)是连续的正整数,连称号中的分母就是第一项\((1-2)(1-3)(1-4)...(1-n)\)和第二项\((2-1)(2-3)(2-4)...(2-n)等等\),可表示为\((-1^{n-i})(n-i)!(i-1)!\)

这是后面就变成了\[\prod_{i!=j} \frac{P_{i-1}S_{i+1}}{(-1)^{(n-i)}(n-i)!(i-1)!}=\frac{\prod_{j=0}^{n}{(x-j)!}}{(-1)^{n-i}(x-i)(n-i)!(i-1)!}\].

计算时用费马小定理求出\((x-i),(n-i)!,(i-1)!\)的逆元累加即可。

求逆元的复杂度是\(O(log_2n)\),累加是\(O(n)\),总的复杂度是\(O(nlong_2n)\).

求幂和

幂和可表示为\[S_k(n)=\dfrac{(n+1)^{k+1}-\sum\limits_{j=0}^{k-1}\binom{k+1}{j}S_j(n)-1}{k+1}\](具体推导见参考文章1)

可见为K+1次多项式。我们知道k+1次需要有k+2个点来确定,因此我们可以暴力计算出前k+2项的幂和,如果n小于等于k+2,我们通过这样的预处理可以直接得到答案。如果n>k+2,则使用上面提到的拉格朗日插值求得n的k次幂和。

注意最后得答案时加一些模在模一下保证是个正数。

代码:

#include <iostream>
#define max_n 1000006
const long long mod = 1000000007;
using namespace std;
int k;
int n;
long long fac[max_n];
long long y[max_n];
long long pi = 1;
long long q_mod(long long a,long long b,long long mod)
{long long res = 1;while(b){if(b&1){res = (res*a)%mod;}a = (a*a)%mod;b>>=1;}return res;
}int main()
{cin >> n >> k;y[0] = 0;for(int i = 1;i<=k+2;i++){y[i] = (y[i-1]+q_mod(i,k,mod))%mod;}if(n<=k+2){cout << y[n] << endl;return 0;}fac[0] = 1;for(int i = 1;i<=max_n;i++){fac[i] = (fac[i-1]*i)%mod;}for(int i = 1;i<=k+2;i++){pi = (pi*(n-i))%mod;}long long ans = 0;for(int i = 1;i<=k+2;i++){long long f = ((k+2-i)%2)?-1:1;long long inv = q_mod(fac[k+2-i]%mod*fac[i-1]%mod,mod-2,mod);long long inv2 = q_mod(n-i,mod-2,mod);ans = (ans+f*y[i]*pi%mod*inv%mod*inv2%mod)%mod;}cout << (ans+mod)%mod << endl;return 0;
}

参考文章:

Wolfycz,自然数幂和,https://www.cnblogs.com/Wolfycz/p/10622577.html,浅谈算法——拉格朗日插值,https://www.cnblogs.com/Wolfycz/p/10624678.html

qscqesze,Educational Codeforces Round 7 F. The Sum of the k-th Powers 拉格朗日插值法,https://www.cnblogs.com/qscqesze/p/5207132.html

sdfzchy,【CF622F】The Sum of the k-th Powers (拉格朗日插值法),https://blog.csdn.net/sdfzchy/article/details/78307534

另外还有一篇比较不错的差分求正整数的k次幂和的文章:

张永强.差分的应用及正整数的k 次方幂求和[J].唐山师范学院学报,2007,29(5) :140-141

转载于:https://www.cnblogs.com/zhanhonhao/p/11348661.html

这篇关于Codeforces D. The Sum of the k-th Powers(拉格朗日插值)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/282702

相关文章

最大流=最小割=最小点权覆盖集=sum-最大点权独立集

二分图最小点覆盖和最大独立集都可以转化为最大匹配求解。 在这个基础上,把每个点赋予一个非负的权值,这两个问题就转化为:二分图最小点权覆盖和二分图最大点权独立集。   二分图最小点权覆盖     从x或者y集合中选取一些点,使这些点覆盖所有的边,并且选出来的点的权值尽可能小。 建模:     原二分图中的边(u,v)替换为容量为INF的有向边(u,v),设立源点s和汇点t

Codeforces Round #240 (Div. 2) E分治算法探究1

Codeforces Round #240 (Div. 2) E  http://codeforces.com/contest/415/problem/E 2^n个数,每次操作将其分成2^q份,对于每一份内部的数进行翻转(逆序),每次操作完后输出操作后新序列的逆序对数。 图一:  划分子问题。 图二: 分而治之,=>  合并 。 图三: 回溯:

Codeforces Round #261 (Div. 2)小记

A  XX注意最后输出满足条件,我也不知道为什么写的这么长。 #define X first#define Y secondvector<pair<int , int> > a ;int can(pair<int , int> c){return -1000 <= c.X && c.X <= 1000&& -1000 <= c.Y && c.Y <= 1000 ;}int m

Codeforces Beta Round #47 C凸包 (最终写法)

题意慢慢看。 typedef long long LL ;int cmp(double x){if(fabs(x) < 1e-8) return 0 ;return x > 0 ? 1 : -1 ;}struct point{double x , y ;point(){}point(double _x , double _y):x(_x) , y(_y){}point op

Codeforces Round #113 (Div. 2) B 判断多边形是否在凸包内

题目点击打开链接 凸多边形A, 多边形B, 判断B是否严格在A内。  注意AB有重点 。  将A,B上的点合在一起求凸包,如果凸包上的点是B的某个点,则B肯定不在A内。 或者说B上的某点在凸包的边上则也说明B不严格在A里面。 这个处理有个巧妙的方法,只需在求凸包的时候, <=  改成< 也就是说凸包一条边上的所有点都重复点都记录在凸包里面了。 另外不能去重点。 int

Codeforces 482B 线段树

求是否存在这样的n个数; m次操作,每次操作就是三个数 l ,r,val          a[l] & a[l+1] &......&a[r] = val 就是区间l---r上的与的值为val 。 也就是意味着区间[L , R] 每个数要执行 | val 操作  最后判断  a[l] & a[l+1] &......&a[r] 是否= val import ja

如何导入sun.misc.BASE64Encoder和sum.misc.BASE64Decoder

右击项目名--->Build Path--->Configure Build Path...--->java Build Path--->Access rules:1 rule defined,added to all librar...   --->Edit --->Add...

Codeforces Round 971 (Div. 4) (A~G1)

A、B题太简单,不做解释 C 对于 x y 两个方向,每一个方向至少需要 x / k 向上取整的步数,取最大值。 由于 x 方向先移动,假如 x 方向需要的步数多于 y 方向的步数,那么最后 y 方向的那一步就不需要了,答案减 1 代码 #include <iostream>#include <algorithm>#include <vector>#include <string>

Codeforces#295(Div.2)A、B(模拟+BFS)

解题报告链接:点击打开链接 C. 题目链接:点击打开链接 解题思路: 对于给定的字符串,取出现次数最多的字母(可以同时有多个)。由这些字母组成长度为n的字符串,求有多少种组合。最后用数学知识即可。 完整代码: #include <algorithm>#include <iostream>#include <cstring>#include <climits>

Codeforces Round #281 (Div. 2)A(构造+暴力模拟)

题目链接:http://codeforces.com/problemset/problem/493/A 解题思路: 暴力的判断,分三种情况去判断即可。注意如果之前已经被罚下场后,那么在后面的罚下情况不应该算在输出结果内。 完整代码: #include <algorithm>#include <iostream>#include <cstring>#include <co