首页 > 其他分享 >闲话 23.1.7

闲话 23.1.7

时间:2023-01-07 20:13:44浏览次数:38  
标签:int 闲话 inv ret 23.1 s2 s1 mod

闲话

数学这里证明是个避不过去的东西。你总得面对他。

—— 巨佬

今日 ABC(
先写一点 慢慢补上

拉格朗日插值

熟知平面上 \(n + 1\) 个点对应着一个 \(n\) 次多项式的系数。可是这系数要怎么确定呢?我们已经有了高斯消元法,可是其复杂度是 \(O(n^3)\) 的,不优秀。拉格朗日提供了一种复杂度为 \(O(n^2)\) 的做法,其思想就是直接配凑。

考虑平面上的 \(n\) 个点,第 \(i\) 个点记作 \((x_i, y_i)\),并记 \(F(x)\) 为这 \(n\) 个点对应的次数最小的多项式。这里假设 \(x_i, y_i\) 分别互不相同。如果能对每个点 \(i\) 找到一个函数 \(f_i(x)\),使得 \(f_i(x_i) = 1\),并且 \(\forall j\neq i,\ f_i(x_j) = 0\),那我们直接对所有点值求和得到

\[F(x) = \sum_{i=1}^n y_i \times f_{i}(x) \]

这样 \(F(x)\) 肯定经过这 \(n\) 个点。随后我们就需要满足一下次数的最小性了。可以发现我们需要的 \(f_i(x)\) 应当是 \(n-1\) 次多项式,因此这启发我们通过 \(n-1\) 个单项式相乘的方法去限制 \(x = x_j\) 时 \(f_i(x)\) 为 \(0\),直接乘入一个单项式 \((x - x_j)\) 即可,这就得到了 \(f_i(x)\) 的部分构造。然后我们发现还差一个 \(f_i(x_i) = 1\),这如何满足呢?考虑一个 \(p / p = 1\),我们将目前构造出的多项式除以 \(\prod (x_i - x_j)\),这样在 \(x = x_i\) 时上下相同,取到 \(1\)。
很容易验证这样得出的 \(f_i(x)\) 是满足所有要求的。因此我们构造出了 \(F(x)\) 的一般形式。

总结一下能得到拉格朗日插值法的公式:

\[F(x) = \sum_{i=1}^n y_i \prod_{j\neq i}\frac{x - x_j}{x_i - x_j} \]

实现
inline int Lagrange(int n, int x[], int y[], int x_qry) {
	int ret = 0;
	rep(i,0,n) {
		int s1 = 1, s2 = 1;
		rep(j,0,n) {
			if (i != j) {
				s1 = 1ll * s1 * (x_qry - x[j]) % mod;
				s2 = 1ll * s2 * (x[i] - x[j]) % mod;
			} 
		} ret = (ret + 1ll * y[i] * s1 % mod * qp(s2, mod - 2) % mod) % mod;
	} return (ret + mod) % mod; 
} 

由这个式子可以得到一种 \(x_i\) 点值连续时的 \(O(n)\) 做法。不展开,留作练习。

实现
int s1[N], s2[N], inv[N];
inline int CLagrange(int n, int x[], int y[], int x_qry) {
	int ret = 0;
	s1[0] = (x_qry - x[0]) % mod, s2[n+1] = 1;
	rep(i,1,n) s1[i] = 1ll * s1[i-1] * (x_qry - x[i]) % mod;
	pre(i,n,0) s2[i] = 1ll * s2[i+1] * (x_qry - x[i]) % mod;
	inv[0] = inv[1] = 1;
	rep(i,2,n) inv[i] = -1ll * mod / i * inv[mod % i] % mod;
	rep(i,2,n) inv[i] = 1ll * inv[i-1] * inv[i] % mod;
	rep(i,0,n) ret = (ret + (1ll * y[i] * (i == 0 ? 1 : s1[i-1]) % mod * s2[i+1] % mod * inv[i] % mod * (((n-i) & 1) ? -1 : 1) * inv[n-i]) % mod) % mod;
	return (ret + mod) % mod;
} 

然后有一种动态做这个东西的方法,叫重心拉格朗日插值法,加入一个点的复杂度是 \(O(n)\) 的。

具体的式子不再推了,就是提取出一个 \(\prod_{i} (x - x_i)\),最后的部分除了形如 \((x_i - x_j)\) 的系数就只剩下 \(y_i / (x - x_i)\) 了。我们发现把系数组合得到 \(t_i = y_i / \prod_{j\neq i} (x_i - x_j)\) 后就能得到一个每个点计算贡献为 \(O(n)\) 的式子。

一个具体的应用是求

\[S_k(n) = \sum_{i=1}^n i^k \]

我们可以通过作差的方法得到 \(S_k(n)\) 是一个关于 \(n\) 的 \(k + 1\) 阶多项式。
对 \(S_k(n)\) 做差得到了 \(n^k\),这个东西是关于 \(n\) 的 \(k\) 阶多项式,因此 \(S_k(n)\) 是一个关于 \(n\) 的 \(k + 1\) 阶多项式。

随后可以直接通过拉格朗日差值得到系数和点值。

我们通过上面的讨论得知在点值连续时可以 \(O(k)\) 地得到这个多项式,这里只需要计算点值所以容易些。对于 \(i\in [1, k + 2]\) 计算 \(i^{k}\) 是容易 \(O(k)\) 的,于是我们就得到了一种 \(O(k)\) 计算 \(S_k(n)\) 的方法。

实现
int S(int n, int k) {
    vector<int> pref(k + 4), suff(k + 4), pw(k + 4), prime(k + 4); 
    vector<bool> __vis(k + 4, 0);
    pw[1] = 1;
    for (int i = 2, cnt = 0; i <= k + 2; ++ i) {
        if (!__vis[i]) prime[++ cnt] = i, pw[i] = qp(i, k);
        for (int j = 1; j <= cnt and i * prime[j] <= k + 2; ++ j) {
            __vis[i * prime[j]] = 1;
            pw[i * prime[j]] = mul(pw[i], pw[prime[j]]);
            if (i % prime[j] == 0) break;
        }
    } rep(i,2,k+2) pw[i] = add(pw[i], pw[i - 1]);
    pref[0] = suff[k + 3] = 1;
    rep(i, 1, k + 2) pref[i] = mul(pref[i - 1], (n - i + mod));
    pre(i, k + 2, 1) suff[i] = mul(suff[i + 1], (n - i + mod));
    int tmp = 0, ret = 0;
    rep(i, 1, k + 2)
        ret = add(ret, mul(((k - i) & 1) ? mod - 1 : 1, pw[i], pref[i - 1], suff[i + 1], ifc(i - 1), ifc(k + 2 - i)));
    return ret;
}

未完。

标签:int,闲话,inv,ret,23.1,s2,s1,mod
From: https://www.cnblogs.com/joke3579/p/chitchat230107.html

相关文章

  • 2023.1.7 DP 学习日志
    上午编辑距离(AcWing.899)思路:同最短编辑距离,只不过要做\(m\)次。code:#include<bits/stdc++.h>#definelllonglong#defineN1005usingnamespacestd;inlinel......
  • 2023.1.7学习笔记
    、经典类与新式类经典类:​不继承object的类或者其子类的类新式类:​继承object或者其之类的类​在python3中,只有新式类,所有类都默认继承object​在python2中,区......
  • 力扣每日一题2023.1.6---2180. 统计各位数字之和为偶数的整数个数
    给你一个正整数num,请你统计并返回小于或等于num且各位数字之和为偶数的正整数的数目。正整数的各位数字之和是其所有位上的对应数字相加的结果。示例1:输入:num=......
  • 2023.1.6 (Codeforces Round #842 (Div. 2))
    A.GreatestConvexLinkhttps://codeforces.com/contest/1768/problem/ADescription求出最大的\(x(1\leqx<k)\),使得\(x!+(x-1)!\)是\(k\)的倍数。Soluti......
  • 2023.1.6
    搜索了某些学校的心理学培养方案,并没有找到相应的课本和参考教材。有几个可以参考下:【资料整理】如何看到国内外各大学的课表、课程大纲,https://zhuanlan.zhihu.com/p/4......
  • 闲话 23.1.6
    闲话没自觉地发现活干不完菜菜菜今日推歌:《MarineBloomin’》八王子P×鹿乃昨天讲解是不是复读含量太高了?推荐阅读:活下去了byEntropyIncreaser人会约束自己。因......
  • the fourteenth——20223.1.6
    #include<stdio.h>intmain(){ 3,4,5;//这是一条语句 //把上面这条语句的值赋值给变量a inta=(3,4,5); printf("a=%d\n",a);}输出结果:a=5因为a的值是整......
  • 关系运算符、逻辑运算符——the thirteenth——2023.1.5
    关系运算符  在C语言中=是赋值的意思,而==才是等于的意思  逻辑运算符一共有三种:&&(并且)、||(或者)、!(非)年龄:取值16-50岁。身高:取值150cm-190cm。身材:1-火辣;2-......
  • 2023.1.6 DP 学习日志
    今天还是学DP,干了两道题1.最长上升子序列II(AcWing.896)数据加强版的最长上升子序列不能直接DP,还得二分(其实有点像贪心)比较简单思路就不写了。其实我WA了五次code:#inc......
  • 2023.1.06 java打印杨辉三角(二维数组)
    publicclassyanghui{publicstaticvoidmain(String[]args){int[][]yanghui=newint[10][];for(inti=0;i<yanghui.length;i++){......