博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
@总结 - 5@ 牛顿迭代法的应用——多项式开方,对数,指数,三角与幂函数
阅读量:4453 次
发布时间:2019-06-07

本文共 8644 字,大约阅读时间需要 28 分钟。

目录


@0 - 参考资料@

@0.5 - 多项式平方根@

已知一个多项式 \(A(x)\),求一个多项式 \(F(x)\),使得 \(F^2(x) = A(x)\)

我们仿照多项式求逆的方法,加上一个模数变成 \(F^2(x) = A(x) \mod x^n\)

然后考虑倍增,假如我们已知:
\[F_0^2(x) \equiv A(x) \mod x^{\lceil\frac{n}{2}\rceil}\]
套路移项,得到:
\[F_0^2(x)-A(x) \equiv 0 \mod x^{\lceil\frac{n}{2}\rceil}\]
套路平方,得到:
\[F_0^4(x)-2F_0^2(x)A(x)+A^2(x) \equiv 0 \mod x^{n}\]
你看右边空荡荡的,心里肯定有点儿不爽对吧:
\[F_0^4(x)+2F_0^2(x)A(x)+A^2(x) \equiv 4F_0^2(x)A(x) \mod x^{n}\]
你看右边已经有个 \(A(x)\) 了,我们就再把 \(A(x)\) 单独孤立出来吧。
\[\dfrac{F_0^4(x)+2F_0^2(x)A(x)+A^2(x)}{4F_0^2(x)} \equiv A(x) \mod x^{n}\]
你看左边已经是一个平方式了哦~
我们令 \(F(x) = \dfrac{F_0^2(x)+A(x)}{2F_0(x)}\),则 \(F^2(x)\equiv A(x) \mod x^{n}\)

然后就可以一路倍增上去了。边界情况下是一个模意义下的平方剩余。

这个多项式开方和多项式求逆之间有许多共通的思想,但有些地方是很具有构造性的。

但如果题目要我们求立方根呢?如果题目要我们解高次方程呢?我们这样推导肯定是不合适的。
我们是否可以找到一个通用的方法,对于这两个,甚至更多的问题都适用呢?

@1 - 牛顿迭代法@

@数学上的定义@

本节跟我们要讲的没有任何关系。

对于某一个函数 \(f(x)\) 与函数上某一个点 \((x_0, f(x_0))\),我们可以用过该点的函数切线解析式去近似地代替整个函数。

根据导函数的定义,过 \((x_0, f(x_0))\) 的函数切线解析式为 \(g(x) = f(x_0)+f'(x_0)*(x-x_0)\)。则我们可以用 \(g(x)\) 去近似替代 \(f(x)\)
假如 \(g(x)\) 的零点为 \(g(z)=0\)\(z = x_0 - \frac{f(x_0)}{f'(x_0)}\)),因为 \(g(x)\)\(f(x)\) 的近似,所以也许 \(f(z) = 0\)
但是这个概率非常的低。于是我们再选择 \((z,f(z))\) 进行迭代,求出 \((z,f(z))\) 对应的 \(g(x)\),重复十几二十百遍,应该就非常非常接近零点了。

然而……这个方法并不具有普适性,有些时候迭代会出现循环。但是对于多项式的相关操作,它的确是普适的。

根据泰勒展开

\[f(x)=f(x_0)+f'(x_0)(x-x_0)+\dfrac{f''(x_0)(x-x_0)^2}{2!}+\dots\]
我们只取前两项(因为它们是线性的),也可以得到上面的结果。

@对于多项式的定义@

总结多项式逆元与多项式开方,我们可以得到通用的一些东西。

对于一个自变量和因变量都是多项式的函数 \(G(F(x))\),我们要找到多项式 \(F_0(x)\),使得:

\[G(F_0(x))=0\]
\(G(F(x)) = A(x)*F(x) - 1\) 时相当于逆元,当 \(G(F(x)) = F^2(x)-A(x)\) 时相当于开方。

怎么弄呢?我们套路取模,转换为 \(G(F(x))\equiv 0\mod x^n\)

再套路倍增,假如我们已知:
\[G(F_0(x))\equiv 0\mod x^{\lceil\frac{n}{2}\rceil}\]
\(G(F(x))\)\(F_0(x)\) 处泰勒展开,得到:
\[G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))+\dfrac{G''(F_0(x))(F(x)-F_0(x))^2}{2!}+\dots\mod x^n\]
因为在模意义下 \(G(F(x)) \equiv 0\mod x^{\lceil\frac{n}{2}\rceil}\) 只有唯一解,所以:
\[F(x)\equiv F_0(x)\mod x^{\lceil\frac{n}{2}\rceil}\]
然后套路移项 + 套路平方:
\[(F(x)-F_0(x))^2\equiv 0\mod x^n\]
因为泰勒展开第三项及其以后的项都含 \((F(x)-F_0(x))^2\),所以我们可以直接忽略。
又因为根据定义 \(G(F(x))\equiv 0\mod x^n\),故:
\[0\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\mod x^n\]
变一下形:
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\mod x^n\]
好的我们就可以倍增了。

@2 - 牛顿迭代的应用@

@重新推导 - 多项式逆元@

对于 \(G(F(x)) = A(x)*F(x) - 1\) ,求 \(G(F(x))\equiv 0\mod x^n\)\(F(x)\)

根据迭代式:
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-\dfrac{A(x)*F_0(x) - 1}{A(x)}\mod x^n\]
根据严密的推导,我们可以得到:\(\dfrac{A(x)*F_0(x) - 1}{A(x)}\equiv F_0(x)(A(x)*F_0(x) - 1)\mod x^n\)
所以:
\[F(x)\equiv 2F_0(x)-A(x)*F_0^2(x)\mod x^n\]

时间复杂度 \(O(n\log n)\)

@重新推导 - 多项式平方根@

对于 \(G(F(x)) = F^2(x)-A(x)\) ,求 \(G(F(x))\equiv 0\mod x^n\)\(F(x)\)

根据迭代式:
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-\dfrac{F_0^2(x) - A(x)}{2F_0(x)}\mod x^n\]
稍微恒等变换一下,就可以得到我们上边那个公式。

时间复杂度 \(O(n\log n)\)

@多项式对数函数@

这个其实和牛顿迭代没啥关系。。。

已知 \(F(x) \equiv \ln A(x)\mod x^n\) ,求 \(F(x)\)

怎么去理解呢?其实它就是个泰勒展开:
\[\ln(1-P(x))=-P(x)-\frac{P(x)^2}{2}-\frac{P(x)^3}{3}-\dots\]

直接求不好做的,我们考虑两边求个导数:

\[F'(x) \equiv \frac{A'(x)}{A(x)}\mod x^n\]
这样我们就可以顺利求出 \(F'(x)\),再对 \(F'(x)\) 积分一下就可以了。
但是不定积分的常数项?一般来说题目会保证 \(A(x)\) 的常数项为 1,根据泰勒展开可得 \(F(x)\) 的常数项为 0。

时间复杂度 \(O(n\log n)\)

@多项式指数函数@

已知 \(F(x) \equiv e^{A(x)}\mod x^n\) ,求 \(F(x)\)

如果我们求导是没啥用的,毕竟 \(e^x\) 的导数还是它自己。
我们考虑两边求对数:
\[\ln F(x) \equiv A(x)\mod x^n\]
然后记 \(G(F(x)) = \ln F(x)-A(x)\),于是又可以愉快地进行牛顿迭代啦。
\[F(x)\equiv F_0(x)-\dfrac{G(F_0(x))}{G'(F_0(x))}\equiv F_0(x)-(\ln F_0- A(x))*F_0(x)\mod x^n\]

时间复杂度 \(O(n\log n)\)

@多项式幂函数@

\[F^k(x)=(e^{\ln F(x)})^k=e^{k\ln F(x)}\]

时间复杂度 \(O(n\log n)\)

@多项式三角函数@

\[e^{i*F(x)}=\cos F(x) + i*\sin F(x)\]

把运算改成支持复数的运算。

时间复杂度 \(O(n\log n)\)

@3 - 一些参考代码@

@对数函数@

本代码为 的 AC 代码。

#include
#include
using namespace std;const int G = 3;const int MAXN = 400000;const int MOD = 998244353;int pow_mod(int b, int p) { int ret = 1; while( p ) { if( p & 1 ) ret = 1LL*ret*b%MOD; b = 1LL*b*b%MOD; p >>= 1; } return ret;}int inv[MAXN + 5];void init() { inv[1] = 1; for(int i=2;i<=MAXN;i++) inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD;}struct Polynomial{ void poly_copy(int *A, int *B, int n) { for(int i=0;i
>1);(j^=l)
>=1); } for(int s=2;s<=n;s<<=1) { int t = (s>>1); int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s); for(int i=0;i
< (n<<1);len <<= 1); poly_inv(A, B, (n + 1) >> 1); poly_copy(tmp1, A, n); poly_clear(tmp1, n, len); ntt(tmp1, len, 1); ntt(B, len, 1); for(int i=0;i
< (n+m-1);len <<= 1); poly_copy(tmp2, A, n), poly_copy(tmp3, B, m); ntt(tmp2, len, 1), ntt(tmp3, len, 1); for(int i=0;i
=0;i--) A[i+1] = 1LL*B[i]*inv[i+1]%MOD; A[0] = 0; } void poly_deri(int *A, int *B, int n) { for(int i=1;i

@指数函数@

本代码为 的 AC 代码。

#include
#include
using namespace std;const int G = 3;const int MAXN = 400000;const int MOD = 998244353;int pow_mod(int b, int p) { int ret = 1; while( p ) { if( p & 1 ) ret = 1LL*ret*b%MOD; b = 1LL*b*b%MOD; p >>= 1; } return ret;}int inv[MAXN + 5];void init() { inv[1] = 1; for(int i=2;i<=MAXN;i++) inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD;}struct Polynomial{ void poly_copy(int *A, int *B, int n) { for(int i=0;i
>1);(j^=l)
>=1); } for(int s=2;s<=n;s<<=1) { int t = (s>>1); int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s); for(int i=0;i
< (n<<1);len <<= 1); poly_inv(A, B, (n + 1) >> 1); poly_copy(tmp1, A, n); ntt(tmp1, len, 1); ntt(B, len, 1); for(int i=0;i
< (n+m-1);len <<= 1); poly_copy(tmp2, A, n), poly_copy(tmp3, B, m); ntt(tmp2, len, 1), ntt(tmp3, len, 1); for(int i=0;i
=0;i--) A[i+1] = 1LL*B[i]*inv[i+1]%MOD; A[0] = 0; } void poly_deri(int *A, int *B, int n) { for(int i=1;i
< (n<<1);len <<= 1); poly_exp(A, B, (n + 1) >> 1); poly_copy(tmp5, A, n); poly_ln(B, tmp6, n); ntt(tmp5, len, 1), ntt(tmp6, len, 1), ntt(B, len, 1); for(int i=0;i

@4 - 算法应用@

例题:@COGS - 2189@ 帕秋莉的超级多项式

已知 \(F(x)\) 是个 n 次多项式,\(k\) 是一个整数,且已知:
\[G(x)\equiv(1+ \ln (1+\frac{1}{exp(\int \frac{1}{\sqrt{F(x)}})}))^k\mod x^n\]
\(G'(x)\),先取模再求导。

这道题……我暂时没有找到哪里可以交,但是找网上的程序对拍可以拍上几千组……应该是 AC 了的。

solution:按题意模拟 qwq。
主程序里面要把多项式的长度扩展成 2 的整数次幂才能 AC,但是依照我的理解是不用扩展的,不知道为什么。

#include
#include
#include
using namespace std;const int G = 3;const int MAXN = 1000000;const int MOD = 998244353;int pow_mod(int b, int p) { int ret = 1; while( p ) { if( p & 1 ) ret = 1LL*ret*b%MOD; b = 1LL*b*b%MOD; p >>= 1; } return ret;}int inv[MAXN + 5];void init() { inv[1] = 1; for(int i=2;i<=MAXN;i++) inv[i] = 1LL*(MOD - MOD/i)*inv[MOD%i]%MOD;}struct Polynomial{ void poly_copy(int *A, int *B, int n) { for(int i=0;i
>1);(j^=l)
>=1); } for(int s=2;s<=n;s<<=1) { int t = (s>>1); int u = (type == 1) ? pow_mod(G, (MOD-1)/s) : pow_mod(G, (MOD-1)-(MOD-1)/s); for(int i=0;i
< (n<<1);len <<= 1); poly_inv(A, B, (n + 1) >> 1); poly_copy(tmp1, A, n); ntt(tmp1, len, 1); ntt(B, len, 1); for(int i=0;i
< (n+m-1);len <<= 1); poly_copy(tmp2, A, n), poly_copy(tmp3, B, m); ntt(tmp2, len, 1), ntt(tmp3, len, 1); for(int i=0;i
=0;i--) A[i+1] = 1LL*B[i]*inv[i+1]%MOD; A[0] = 0; } void poly_deri(int *A, int *B, int n) { for(int i=1;i
< (n<<1);len <<= 1); poly_exp(A, B, (n + 1) >> 1); poly_copy(tmp5, A, n); poly_ln(B, tmp6, n); ntt(tmp5, len, 1), ntt(tmp6, len, 1), ntt(B, len, 1); for(int i=0;i
< (n<<1);len <<= 1); poly_sqrt(A, B, (n + 1) >> 1); poly_copy(tmp8, A, n); poly_inv(B, tmp9, n); ntt(tmp8, len, 1), ntt(tmp9, len, 1), ntt(B, len, 1); for(int i=0;i
<= n;m <<= 1); oper.poly_sqrt(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); oper.poly_inv(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); oper.poly_int(g, f, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); oper.poly_exp(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); oper.poly_inv(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); f[0]++; oper.poly_ln(f, g, m); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); f[0]++; oper.poly_pow(f, g, m, k); oper.poly_copy(f, g, m); oper.poly_clear(g, 0, m); oper.poly_deri(g, f, n); oper.poly_copy(f, g, n); oper.poly_clear(g, 0, n); for(int i=0;i

转载于:https://www.cnblogs.com/Tiw-Air-OAO/p/10200363.html

你可能感兴趣的文章
WOJ 46 完全背包
查看>>
spring框架学习笔记(十)
查看>>
排球计分程序(现场记分员)
查看>>
流媒体技术笔记(视频编码相关)
查看>>
神马16核的服务器你让我单线程跑ffmpeg
查看>>
block,inline,inline-block的区别
查看>>
html表单
查看>>
const关键字——读《嵌入式c语言进阶之道》整理
查看>>
libevent入门(1)
查看>>
CSS 样式显示为小手
查看>>
关联模型错误的蛋疼错误提示
查看>>
JS当心隐式的强制转换
查看>>
通过ros节点发布Twist Messages控制机器人--10
查看>>
STL--list
查看>>
maven 学习之路之二(1)
查看>>
爬虫第二课:解析网页中的元素
查看>>
域名注册必须实名认证 《互联网域名管理办法》11月1日实施
查看>>
Convert a byte[] array to readable string format. This makes the "hex" readable!
查看>>
记录一下学习过程中,比较好的文章
查看>>
python之Django学习笔记(一)---搭建Django开发环境和一些基本命令
查看>>