博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
【知识总结】多项式全家桶(二)(ln和exp)
阅读量:4556 次
发布时间:2019-06-08

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

上一篇:

一、对数函数\(\ln(A)\)

求一个多项式\(B(x)\),满足\(B(x)=\ln(A(x))\)

这里需要一些最基本的微积分知识(不会?戳我(暂时戳不动):【知识总结】微积分初步挖坑待填)。

另外,\(n\)次多项式\(A(x)\)可以看成关于\(x\)\(n\)次函数,可以对其求导。显然,\(A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\)的导数是\(A'(x)=\sum\limits_{i=0}^{n-2}a_{i+1}x^i(i+1)\),积分是\(\int A(x)\mathrm{d} x=\sum\limits_{i=1}^{n}\frac{a_{i-1}}{i}x^i\)。可以写出如下代码(非常简单):

void derivative(const int *a, int *b, const int n){    for (int i = 1; i < n; i++)        b[i - 1] = (ll)a[i] * i % p;    b[n - 1] = 0;}void integral(const int *a, int *b, const int n){    for (int i = n - 1; i >= 0; i--)        b[i + 1] = (ll)a[i] * inv(i + 1) % p;    b[0] = 0;}

\(f(x)=\ln(x)\)的导数是\(f'(x)=\frac{1}{x}\)。回到原问题,对两边同时求导,得到(要用一下链式法则\(g(f(x))\)的导数是\(g'(f(x))f'(x)\)):

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

求个\(A(x)\)的逆元(多项式求逆)然后乘上\(A'(x)\),最后把\(B(x)\)积分回去就好了。

至于代码……下面算多项式指数函数的时候要算对数函数,所以暂时省略。

二、指数函数\(\exp(x)\)

求多项式\(B(x)\)满足\(B(x)=e^{A(x)}\)

首先,这个式子相当于求\(\ln B(x)=A(x)\)\(\ln B(x)-A(x)=0\)

设关于多项式的函数\(F(B(x))=\ln B(x)-A(x)\),那么问题就是求这个函数的零点(\(A(x)\)是给定的,视作常数)。

求函数零点的方法之一是牛顿迭代,公式如下(\(i\)是迭代次数,\(x\)是自变量,\(F(x)\)是要求零点的函数,\(F'(x_0)\)\(F(x)\)\(x_0\)处的导数):

\[x_{i+1}=x_i-\frac{F(x_i)}{F'(x_i)}\]

\(F(B(x))=\ln B(x)-A(x)\)求导,得到\(F'(B(x))=\frac{1}{B(x)}\)(注意自变量是\(B(x)\)不是\(x\)。这不是一个\(F(x)\)\(B(x)\)的复合函数)。然后代入上面的公式:

\[ \begin{aligned} B_{i+1}(x)&=B_i(x)-\frac{\ln B_i(x)-A(x)}{\frac{1}{B_i(x)}}\\ &=B_i(x)-B_i(\ln B_i(x)-A(x))\\ &=B_i(x)(1-\ln B_i(x)-A(x)) \end{aligned} \]

由于多项式乘法的存在,每迭代一次\(B\)的有效长度会增加一倍。

下一篇:

代码():

#include 
#include
#include
#include
#undef i#undef j#undef k#undef true#undef false#undef min#undef max#undef swap#undef sort#undef if#undef for#undef while#undef printf#undef scanf#undef putchar#undef getchar#define _ 0using namespace std;namespace zyt{ template
inline bool read(T &x) { char c; bool f = false; x = 0; do c = getchar(); while (c != EOF && c != '-' && !isdigit(c)); if (c == EOF) return false; if (c == '-') f = true, c = getchar(); do x = x * 10 + c - '0', c = getchar(); while (isdigit(c)); if (f) x = -x; return true; } template
inline void write(T x) { static char buf[20]; char *pos = buf; if (x < 0) putchar('-'), x = -x; do *pos++ = x % 10 + '0'; while (x /= 10); while (pos > buf) putchar(*--pos); } typedef long long ll; const int N = 1e5 + 10, LEN = (N << 2), p = 998244353, g = 3; namespace Polynomial { inline int power(int a, int b) { int ans = 1; while (b) { if (b & 1) ans = (ll)ans * a % p; a = (ll)a * a % p; b >>= 1; } return ans; } inline int inv(const int a) { return power(a, p - 2); } int omega[LEN], winv[LEN], rev[LEN]; void init(const int n, const int lg2) { int w = power(g, (p - 1) / n), wi = inv(w); omega[0] = winv[0] = 1; for (int i = 1; i < n; i++) { omega[i] = (ll)omega[i - 1] * w % p; winv[i] = (ll)winv[i - 1] * wi % p; } for (int i = 0; i < n; i++) rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (lg2 - 1))); } void ntt(int *a, const int *w, const int n) { for (int i = 0; i < n; i++) if (i < rev[i]) swap(a[i], a[rev[i]]); for (int l = 1; l < n; l <<= 1) for (int i = 0; i < n; i += (l << 1)) for (int k = 0; k < l; k++) { int tmp = (a[i + k] - (ll)w[n / (l << 1) * k] * a[i + l + k] % p + p) % p; a[i + k] = (a[i + k] + (ll)w[n / (l << 1) * k] * a[i + l + k] % p) % p; a[i + l + k] = tmp; } } void mul(const int *a, const int *b, int *c, const int n) { static int x[LEN], y[LEN]; int m = 1, lg2 = 0; while (m < (n << 1) - 1) m <<= 1, ++lg2; init(m, lg2); memcpy(x, a, sizeof(int[n])); memset(x + n, 0, sizeof(int[m - n])); memcpy(y, b, sizeof(int[n])); memset(y + n, 0, sizeof(int[m - n])); ntt(x, omega, m), ntt(y, omega, m); for (int i = 0; i < m; i++) x[i] = (ll)x[i] * y[i] % p; ntt(x, winv, m); int invm = inv(m); for (int i = 0; i < m; i++) x[i] = (ll)x[i] * invm % p; memcpy(c, x, sizeof(int[n])); } void _inv(const int *a, int *b, const int n) { if (n == 1) b[0] = inv(a[0]); else { static int tmp[LEN]; _inv(a, b, (n + 1) >> 1); int m = 1, lg2 = 0; while (m < (n << 1) + 1) m <<= 1, ++lg2; init(m, lg2); memcpy(tmp, a, sizeof(int[n])); memset(tmp + n, 0, sizeof(int[m - n])); memset(b + ((n + 1) >> 1), 0, sizeof(int[m - ((n + 1) >> 1)])); ntt(tmp, omega, m); ntt(b, omega, m); for (int i = 0; i < m; i++) b[i] = (b[i] * 2LL % p - (ll)tmp[i] * b[i] % p * b[i] % p + p) % p; ntt(b, winv, m); int invm = inv(m); for (int i = 0; i < m; i++) b[i] = (ll)b[i] * invm % p; memset(b + n, 0, sizeof(int[m - n])); } } void inv(const int *a, int *b, const int n) { static int tmp[LEN]; memcpy(tmp, a, sizeof(int[n])); _inv(tmp, b, n); } void derivative(const int *a, int *b, const int n) { for (int i = 1; i < n; i++) b[i - 1] = (ll)a[i] * i % p; b[n - 1] = 0; } void integral(const int *a, int *b, const int n) { for (int i = n - 1; i >= 0; i--) b[i + 1] = (ll)a[i] * inv(i + 1) % p; b[0] = 0; } void ln(const int *a, int *b, const int n) { static int tmp[LEN], inva[LEN]; derivative(a, tmp, n); inv(a, inva, n - 1); mul(inva, tmp, b, n - 1); integral(b, b, n - 1); } void _exp(const int *a, int *b, const int n) { if (n == 1) b[0] = 1; else { static int tmp[LEN]; _exp(a, b, (n + 1) >> 1); ln(b, tmp, n); for (int i = 0; i < n; i++) tmp[i] = (-tmp[i] + a[i] + p) % p; tmp[0] = (tmp[0] + 1) % p; mul(b, tmp, b, n); } } void exp(const int *a, int *b, const int n) { static int tmp[LEN]; memcpy(tmp, a, sizeof(int[n])); _exp(tmp, b, n); } } int work() { static int a[LEN]; int n; read(n); for (int i = 0; i < n; i++) read(a[i]); Polynomial::exp(a, a, n); for (int i = 0; i < n; i++) write(a[i]), putchar(' '); return (0^_^0); }}int main(){ return zyt::work();}

转载于:https://www.cnblogs.com/zyt1253679098/p/10234040.html

你可能感兴趣的文章
Wp7中,使你的程序在程序列表中的名称自适应手机语言设置
查看>>
error: .repo/manifests/: contains uncommitted changes
查看>>
逆元(inv)
查看>>
CentOS Docker 安装
查看>>
debian(kali Linux) 安装net Core
查看>>
centos 7防火墙设置
查看>>
自定义进度条(圆形、横向进度条)
查看>>
spark-streaming-kafka采坑
查看>>
9.Mongodb与python交互
查看>>
18-[JavaScript]-函数,Object对象,定时器,正则表达式
查看>>
读取短信回执
查看>>
EF 数据初始化
查看>>
PreparedStatement与Statement
查看>>
WebService -- Java 实现之 CXF ( 使用CXF工具生成client 程序)
查看>>
Factorial
查看>>
Android开发中,9-patch 图片设置背景带来的问题
查看>>
剑指Offer-反转链表
查看>>
[LeetCode]Two Sum
查看>>
java 常用集合list与Set、Map区别及适用场景总结
查看>>
HTML与javascript语法
查看>>