首页 技术 正文

牛顿迭代，多项式求逆，除法，开方，exp，ln，求幂

技术 2022年11月9日
0 收藏 952 点赞 2,174 浏览 8211 个字

牛顿迭代

若

\[G(F_0(x))\equiv 0(mod\ x^{2^t})
\]

牛顿迭代

\[F(x)\equiv F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}(mod\ x^{2^{t+1}})
\]

以下多数都可以牛顿迭代公式一步得到

多项式求逆

给定\(A(x)\)求满足\(A(x)*B(x)=1\)的\(B(x)\)

写成

\[A(x)*B(x)=1(mod \ x^n)
\]

我们会求$$A(x)*B(x)=1(mod \ x^1)$$

然后我们考虑求$$A(x)*B(x)=1(mod \ x^t)$$

\[(A(x)*B(x)-1)^2=0(mod \ x^{2t})
\]

\[A(x)(2B(x)-A(x)*B^2(x))=1(mod \ x^{2t})
\]

把\(2B(x)-A(x)*B^2(x)\)当作新的\(B\)倍增算

从\(mod \ x^1\)倍增到大于等于\(n\)

# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;
const int _(4e5 + 5);
const int Zsy(998244353);
const int Phi(998244352);
const int G(3);IL int Input(){
RG int x = 0, z = 1; RG char c = getchar();
for(; c < '0' || c > '9'; c = getchar()) z = c == '-' ? -1 : 1;
for(; c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
return x * z;
}IL int Pow(RG ll x, RG ll y){
RG ll ret = 1;
for(; y; x = x * x % Zsy, y >>= 1)
if(y & 1) ret = ret * x % Zsy;
return ret;
}int N, r[_], l, A[_], B[_];IL void NTT(RG int *P, RG int opt){
for(RG int i = 0; i < N; ++i) if(i < r[i]) swap(P[i], P[r[i]]);
for(RG int i = 1; i < N; i <<= 1){
RG int wn = Pow(G, Phi / (i << 1));
if(opt == -1) wn = Pow(wn, Zsy - 2);
for(RG int j = 0, p = i << 1; j < N; j += p)
for(RG int w = 1, k = 0; k < i; w = 1LL * w * wn % Zsy, ++k){
RG int X = P[k + j], Y = 1LL * w * P[k + j + i] % Zsy;
P[k + j] = (X + Y) % Zsy, P[k + j + i] = (X - Y + Zsy) % Zsy;
}
}
}IL void Mul(RG int *a, RG int *b, RG int len){
for(l = 0, N = 1, len += len - 1; N <= len; N <<= 1) ++l;
for(RG int i = 0; i < N; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
for(RG int i = 0; i < N; ++i) A[i] = B[i] = 0;
for(RG int i = 0; i <= len >> 1; ++i) A[i] = a[i], B[i] = b[i];
NTT(A, 1), NTT(B, 1);
for(RG int i = 0; i < N; ++i) A[i] = 1LL * A[i] * B[i] % Zsy * B[i] % Zsy;
NTT(A, -1);
RG int inv = Pow(N, Zsy - 2);
for(RG int i = 0; i <= len; ++i) A[i] = 1LL * A[i] * inv % Zsy;
}int n, a[_], b[_], c[_], m;int main(RG int argc, RG char *argv[]){
n = Input();
for(RG int i = 0; i < n; ++i) a[i] = Input() % Zsy;
for(m = 1; m < n; m <<= 1);
b[0] = Pow(a[0], Zsy - 2);
for(RG int t = 2; t <= m; t <<= 1){
for(RG int i = 0; i < t; ++i) c[i] = b[i];
Mul(a, b, t);
for(RG int i = 0; i < t; ++i)
b[i] = ((c[i] + c[i]) % Zsy - A[i] + Zsy) % Zsy;
}
for(RG int i = 0; i < n; ++i) printf("%d ", b[i]);
return puts(""), 0;
}

多项式除法

已知\(n\)次多项式\(A(x)\)和\(m\)次多项式\(B(x)\)，\(n>m\)

求\(A(x)\)除以\(B(x)\)的商\(C(x)\)及余式\(D(x)\)

也就是

\[A(x)=B(x)C(x)+D(x)
\]

\[A(\frac{1}{x})=B(\frac{1}{x})C(\frac{1}{x})+D(\frac{1}{x})
\]

\[x^nA(\frac{1}{x})=(x^mB(\frac{1}{x}))(x^{n-m}C(\frac{1}{x}))+x^{n-m+1}(x^{n-1}D(\frac{1}{x}))
\]

对\(x^{n-m+1}\)取模

\[A^R(x)=B^R(x)C^R(x)(mod \ x^{n-m+1})
\]

\(R\)表示把系数倒置，即\(swap\)前后

那么可以多项式求逆，再相乘得到\(C^R(x)\)

然后去掉\(R\)带入原式相乘再相减得到\(D(x)\)

# include <bits/stdc++.h>
# define IL inline
# define RG register
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;IL int Input(){
RG int x = 0, z = 1; RG char c = getchar();
for(; c < '0' || c > '9'; c = getchar()) z = c == '-' ? -1 : 1;
for(; c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
return x * z;
}const int mod(998244353);
const int phi(998244352);
const int maxn(8e5 + 5);
const int g(3);IL int Pow(RG ll x, RG ll y){
RG ll ret = 1;
for(; y; y >>= 1, x = x * x % mod)
if(y & 1) ret = ret * x % mod;
return ret;
}int x[maxn], y[maxn], z[maxn], s[maxn], c[maxn], n, m, a[maxn], b[maxn], r[maxn], len;IL void NTT(RG int *p, RG int m, RG int opt){
RG int l = 0, n = 1;
for(; n < m; n <<= 1) ++l;
for(RG int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
for(RG int i = 0; i < n; ++i) if(r[i] < i) swap(p[r[i]], p[i]);
for(RG int i = 1; i < n; i <<= 1){
RG int wn = Pow(g, phi / (i << 1));
if(opt == -1) wn = Pow(wn, mod - 2);
for(RG int j = 0, t = i << 1; j < n; j += t)
for(RG int k = 0, w = 1; k < i; ++k, w = 1LL * w * wn % mod){
RG int x = p[j + k], y = 1LL * w * p[j + k + i] % mod;
p[j + k] = (x + y) % mod, p[j + k + i] = (x - y + mod) % mod;
}
}
if(opt == -1){
RG int inv = Pow(n, mod - 2);
for(RG int i = 0; i < n; ++i) p[i] = 1LL * p[i] * inv % mod;
}
}IL void Inv(RG int *p, RG int *q, RG int num){
if(num == 1){
q[0] = Pow(p[0], mod - 2);
return;
}
Inv(p, q, num >> 1);
for(RG int i = 0; i < num; ++i) a[i] = p[i], b[i] = q[i];
RG int l = num << 1;
NTT(a, l, 1), NTT(b, l, 1);
for(RG int i = 0; i < l; ++i) a[i] = 1LL * a[i] * b[i] % mod * b[i] % mod;
NTT(a, l, -1);
for(RG int i = 0; i < num; ++i) q[i] = ((2 * q[i] - a[i]) % mod + mod) % mod;
for(RG int i = 0; i < l; ++i) a[i] = b[i] = 0;
}int main(){
n = Input(), m = Input();
for(RG int i = 0; i <= n; ++i) y[i] = Input();
for(RG int i = 0; i <= m; ++i) x[i] = Input();
reverse(x, x + m + 1), reverse(y, y + n + 1);
for(RG int i = 0; i <= n - m; ++i) s[i] = x[i];
for(len = 1; len <= n - m; len <<= 1);
Inv(s, z, len);
for(len = 1; len <= (n - m) << 1; len <<= 1);
for(RG int i = 0; i <= n - m; ++i) a[i] = y[i];
for(RG int i = 0; i <= n - m; ++i) b[i] = z[i];
len <<= 1;
NTT(a, len, 1), NTT(b, len, 1);
for(RG int i = 0; i < len; ++i) b[i] = 1LL * b[i] * a[i] % mod;
NTT(b, len, -1), reverse(b, b + n - m + 1);
for(RG int i = 0; i <= n - m; ++i) c[i] = b[i], printf("%d ", c[i]);
puts(""), reverse(y, y + n + 1), reverse(x, x + m + 1);
for(RG int i = 0; i < len; ++i) a[i] = b[i] = 0;
for(len = 1; len <= n; len <<= 1);
for(RG int i = 0; i <= n; ++i) a[i] = x[i];
for(RG int i = 0; i <= n - m; ++i) b[i] = c[i];
NTT(a, len, 1), NTT(b, len, 1);
for(RG int i = 0; i < len; ++i) a[i] = 1LL * a[i] * b[i] % mod;
NTT(a, len, -1);
for(RG int i = 0; i < m; ++i) y[i] = (y[i] - a[i] + mod) % mod, printf("%d ", y[i]);
return 0;
}

多项式开方

给定\(B(x)\)，求\(A(x)\)使\(A^2(x)=B(x)\)

写成

\[A^2(x)=B(x) (mod \ x^n)
\]

然后

\[A^2(x)=B(x)(mod \ x^1)
\]

是可以求的，好像是什么二次剩余

留坑在这里以后补

考虑求

\[A^2(x)=B(x)(mod \ x^{2t})
\]

令

\[C^2(x)=B(x)(mod \ x^t)
\]

那么

\[(A^2(x)-C^2(x))^2=0(mod \ x^{2t})
\]

\[A^4(x)-2A^2(x)C^2(x)+C^$(x)=0(mod \ x^{2t})
\]

\[A^2(x)=\frac{A^4(x)+C^4(x)}{2C^2(x)}(mod \ x^{2t})
\]

\[2A^2(x)=\frac{A^4(x)+2C^2(x)A^2(x)+C^4(x)}{2C^2(x)}(mod \ x^{2t})
\]

\[A^2(x)=\frac{(A^2(x)+C^2(x))^2}{(2C(x))^2}(mod \ x^{2t})
\]

\[A(x)=\frac{B(x)+C^2(x)}{2C(x)}(mod \ x^{2t})
\]

同样的还是倍增求

多项式ln

设\(B(x)=ln(A(x))\)

同时求导

就是

\[B'(x)=\frac{A'(x)}{A(x)}
\]

那么多项式求导，然后多项式求逆，然后多项式积分就好了

多项式exp

设\(B(x)=e^{A(x)}\)

那么

\[ln(B(x))=A(x)
\]

牛顿迭代

得

\[B(x)(mod\ 2^{t+1})=B(x)(mod\ 2^t)(1-ln(B(x)(mod\ 2^t))+A(x))
\]

多项式求幂

设\(B(x)=A^k(x)\)

那么\(ln(B(x))=kln(A(x))\)

然后就是求\(ln\)然后\(exp\)就好了

模板！！！

COGS2189. [HZOI 2015] 帕秋莉的超级多项式

# include <bits/stdc++.h>
# define RG register
# define IL inline
# define Fill(a, b) memset(a, b, sizeof(a))
using namespace std;
typedef long long ll;IL int Input(){
    RG int x = 0, z = 1; RG char c = getchar();
    for(; c < '0' || c > '9'; c = getchar()) z = c == '-' ? -1 : 1;
    for(; c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
    return x * z;
}const int maxn(8e5 + 5);
const int mod(998244353);
const int phi(998244352);
const int inv2(499122177);
const int g(3);int a[maxn], b[maxn], c[maxn], d[maxn], r[maxn];IL int Pow(RG ll x, RG ll y){
RG ll ret = 1;
for(; y; y >>= 1, x = x * x % mod)
if(y & 1) ret = ret * x % mod;
return ret;
}IL void NTT(RG int *p, RG int m, RG int opt){
RG int n = 1, l = 0;
for(; n < m; n <<= 1) ++l;
for(RG int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
for(RG int i = 0; i < n; ++i) if(r[i] < i) swap(p[i], p[r[i]]);
for(RG int i = 1; i < n; i <<= 1){
RG int t = i << 1, wn = Pow(g, phi / t);
if(opt == -1) wn = Pow(wn, mod - 2);
for(RG int j = 0; j < n; j += t)
for(RG int k = 0, w = 1; k < i; ++k, w = 1LL * w * wn % mod){
RG int x = p[k + j], y = 1LL * w * p[k + j + i] % mod;
p[k + j] = (x + y) % mod, p[k + j + i] = (x - y + mod) % mod;
}
}
if(opt == -1){
RG int inv = Pow(n, mod - 2);
for(RG int i = 0; i < n; ++i) p[i] = 1LL * p[i] * inv % mod;
}
}IL void Inv(RG int *p, RG int *q, RG int len){
if(len == 1){
q[0] = Pow(p[0], mod - 2);
return;
}
Inv(p, q, len >> 1);
for(RG int i = 0; i < len; ++i) a[i] = p[i], b[i] = q[i];
RG int tmp = len << 1;
NTT(a, tmp, 1), NTT(b, tmp, 1);
for(RG int i = 0; i < tmp; ++i) a[i] = 1LL * a[i] * b[i] % mod * b[i] % mod;
NTT(a, tmp, -1);
for(RG int i = 0; i < len; ++i) q[i] = ((2 * q[i] - a[i]) % mod + mod) % mod;
for(RG int i = 0; i < tmp; ++i) a[i] = b[i] = 0;
}IL void Sqrt(RG int *p, RG int *q, RG int len){
if(len == 1){
q[0] = sqrt(p[0]); //???
return;
}
Sqrt(p, q, len >> 1), Inv(q, c, len);
RG int tmp = len << 1;
for(RG int i = 0; i < len; ++i) a[i] = p[i];
NTT(a, tmp, 1), NTT(c, tmp, 1);
for(RG int i = 0; i < tmp; ++i) a[i] = 1LL * a[i] * c[i] % mod;
NTT(a, tmp, -1);
for(RG int i = 0; i < len; ++i) q[i] = 1LL * (q[i] + a[i]) % mod * inv2 % mod;
for(RG int i = 0; i < tmp; ++i) a[i] = c[i] = 0;
}IL void ICalc(RG int *p, RG int *q, RG int len){
q[len - 1] = 0;
for(RG int i = 1; i < len; ++i) q[i - 1] = 1LL * p[i] * i % mod;
}IL void Calc(RG int *p, RG int *q, RG int len){
q[0] = 0;
for(RG int i = 1; i < len; ++i) q[i] = 1LL * Pow(i, mod - 2) * p[i - 1] % mod;
}IL void Ln(RG int *p, RG int *q, RG int len){
Inv(p, c, len), ICalc(p, a, len);
RG int tmp = len << 1;
NTT(c, tmp, 1), NTT(a, tmp, 1);
for(RG int i = 0; i < tmp; ++i) c[i] = 1LL * c[i] * a[i] % mod;
NTT(c, tmp, -1), Calc(c, q, len);
for(RG int i = 0; i < tmp; ++i) a[i] = c[i] = 0;
}IL void Exp(RG int *p, RG int *q, RG int len){
if(len == 1){
q[0] = 1;
return;
}
Exp(p, q, len >> 1), Ln(q, b, len);
for(RG int i = 0; i < len; ++i) b[i] = (mod - b[i] + p[i]) % mod, c[i] = q[i];
b[0] = (b[0] + 1) % mod;
RG int tmp = len << 1;
NTT(b, tmp, 1), NTT(c, tmp, 1);
for(RG int i = 0; i < tmp; ++i) b[i] = 1LL * b[i] * c[i] % mod;
NTT(b, tmp, -1);
for(RG int i = 0; i < len; ++i) q[i] = b[i];
for(RG int i = 0; i < tmp; ++i) b[i] = c[i] = 0;
}IL void CalcPow(RG int *p, RG int *q, RG int len, RG int y){
Ln(p, d, len);
for(RG int i = 0; i < len; ++i) d[i] = 1LL * d[i] * y % mod;
Exp(d, q, len);
for(RG int i = 0; i < len; ++i) d[i] = 0;
}int f[maxn], h[maxn], n, k, len;IL void Record(){
for(RG int i = 0; i < len; ++i) f[i] = h[i], h[i] = 0;
}int main(){
freopen("polynomial.in", "r", stdin);
freopen("polynomial.out", "w", stdout);
n = Input() - 1, k = Input();
for(RG int i = 0; i <= n; ++i) f[i] = Input();
for(len = 1; len <= n; len <<= 1);
Sqrt(f, h, len), Record();
Inv(f, h, len), Record();
Calc(f, h, n + 1), Record();
Exp(f, h, len), Record();
Inv(f, h, len), Record();
f[0] = (f[0] + 1) % mod;
Ln(f, h, len), Record();
f[0] = (f[0] + 1) % mod;
CalcPow(f, h, len, k), Record();
ICalc(f, h, n + 1);
for(RG int i = 0; i <= n; ++i) printf("%d ", h[i]);
    return 0;
}
点赞 952

贡献者

上一篇： IntelliJ IDEA16 热部署，解决每次修改java文件就得重启tomcat的问题
下一篇： 汉诺塔matlab实现
相关推荐
python开发_常用的python模块及安装方法
python开发_常用的python模块及安装方法
adodb：我们领导推荐的数据库连接组件bsddb3：BerkeleyDB的连接组件Cheetah-1.0：我比较喜欢这个版本的cheeta…
技术
日期：2022-11-24 点赞：878 阅读：9,024
Educational Codeforces Round 11 C. Hard Process 二分
Educational Codeforces Round 11 C. Hard Process 二分
C. Hard Process题目连接：http://www.codeforces.com/contest/660/problem/CDes…
技术
日期：2022-11-24 点赞：807 阅读：5,514
下载Ubuntn 17.04 内核源代码
下载Ubuntn 17.04 内核源代码
zengkefu@server1:/usr/src$ uname -aLinux server1 4.10.0-19-generic #21…
技术
日期：2022-11-24 点赞：569 阅读：6,362
可用Active Desktop Calendar V7.86 注册码序列号
可用Active Desktop Calendar V7.86 注册码序列号
可用Active Desktop Calendar V7.86 注册码序列号Name: www.greendown.cn Code: &nb…
技术
日期：2022-11-24 点赞：733 阅读：6,143
Android调用系统相机、自定义相机、处理大图片
Android调用系统相机、自定义相机、处理大图片
Android调用系统相机和自定义相机实例本博文主要是介绍了android上使用相机进行拍照并显示的两种方式，并且由于涉及到要把拍到的照片显…
技术
日期：2022-11-24 点赞：512 阅读：7,776
Struts的使用
Struts的使用
一、Struts2的获取　　Struts的官方网站为：http://struts.apache.org/　　下载完Struts2的jar包,…
技术
日期：2022-11-24 点赞：671 阅读：4,854
Lomu

贡献者

这家伙很懒，只想把你留下。

文章
55
收藏
0
人气
172,844
粉丝
进主页
文章展示
  • [ 技术 ] python开发_常用的python模块及安装方法
  • [ 技术 ] Educational Codeforces Round 11 C. Hard Process 二分
  • [ 技术 ] 下载Ubuntn 17.04 内核源代码
  • [ 技术 ] 可用Active Desktop Calendar V7.86 注册码序列号
  • [ 技术 ] Android调用系统相机、自定义相机、处理大图片
  • [ 技术 ] Struts的使用
    • 个人收藏笔记记录
    我的笔记
    首页 首页 菜单 首页 首页