• 【知识总结】多项式全家桶(二)(ln和exp)


    上一篇:【知识总结】多项式全家桶(一)(NTT、加减乘除和求逆)

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

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

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

    另外,(n)次多项式(A(x))可以看成关于(x)(n)次函数,可以对其求导。显然,(A(x)=sumlimits_{i=0}^{n-1}a_ix^i)的导数是(A'(x)=sumlimits_{i=0}^{n-2}a_{i+1}x^i(i+1)),积分是(int A(x)mathrm{d} x=sumlimits_{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))的复合函数)。然后代入上面的公式:

    [egin{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)的有效长度会增加一倍。

    下一篇:【知识总结】多项式全家桶(三)(任意模数NTT)

    代码(洛谷4726):

    #include <cstdio>
    #include <algorithm>
    #include <cctype>
    #include <cstring>
    #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 _ 0
    using namespace std;
    
    namespace zyt
    {
    	template<typename T>
    	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<typename T>
    	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();
    }
    
  • 相关阅读:
    mybatis报错invalid types () or values ()解决方法
    windows下新安装的mysql修改root password问题
    2分钟在eclipse下使用SpringBoot搭建Spring MVC的WEB项目
    Windows 10 下mysql 安装后无法启动问题
    【Head-First设计模式】C#版-学习笔记-开篇及文章目录
    【博客美化】03.分享按钮
    【博客美化】02.公告栏显示个性化时间
    【博客美化】文章目录
    MySql字段类型及字节
    数据库设计原则
  • 原文地址:https://www.cnblogs.com/zyt1253679098/p/10234040.html
Copyright © 2020-2023  润新知