zoukankan      html  css  js  c++  java
  • 「学习笔记」拉格朗日插值

    一般形式的构造过程

    ((n+1)) 个横坐标不同的点可以唯一确定一个 (n) 次多项式,假如知道答案是个关于某个值的多项式以及其项数和若干个点,那么就可以拉格朗日插值出多项式,算出答案。

    简单来讲,我们需要构造一个 (n) 次多项式,使得它满足在已知的点 (x_i) 处取值为 (y_i)

    构造 (f_i(x)),使得其仅在 (x_i) 处取值为 (1),其他 (x_j,j eq i) 处取值为 (0)

    那么 (y_if_i(x)) 仅在 (x_i) 处取值为 (y_i),其他 (x_j,j eq i) 处取值为 (0)

    (y_if_i(x)) 求和即得到了我们要求的多项式:在 (x_i) 处取值为 (y_i)(n) 次多项式。

    如何构造 (f_i(x))

    • 如果要在 (x_j,j eq i) 处取值为 (0),可以让 ((x-x_j),j eq i) 乘起来,即为 (prod_{j eq i}(x-x_j)),为描述方便记作 (g(x)),基于这个式子,不管怎么对其乘除,在 (x_j,j eq i) 处取值都为 (0)

    • 如果要在 (x_i) 处取值为 (1),则让 (f_i(x)=frac{g(x)}{g(x_i)}),这样当 (x=x_i) 时,(f_i(x_i)=frac{g_{x_i}}{g_{x_i}}=1),注意到 (x_i eq j,i eq j),所以不会出现除 (0) 的情况。

    综上所述,我们构造出了 (f_i(x)=prod_{j eq i}frac{x-x_j}{x_i-x_j}),题目中所求的 (n) 次多项式即为:

    [egin{aligned} &sum _{i=1}^{n}y_if_i(x) \ &sum_{i=1}^ny_iprod_{j eq i}frac{x-x_j}{x_i-x_j} end{aligned} ]

    (x_i) 是连续的 (n) 个整数的 (mathcal{O}(n)) 做法

    代入值到 (x) 中直接对这个式子计算,每一次枚举 (i) 后,最后一块算分母的逆元,只会算 (n) 次逆元,时间复杂度为 (mathcal{O}(n^2))

    特别地,若 (x_i) 是连续的 (n) 个整数,则可以 (mathcal{O}(n)) 计算。

    (pre_i=prod_{j=1}^i x-j,suf_i=prod_{j=i+1}^nx-j),也就是分子的前缀/后缀积,(fac_i=i!)

    [egin{aligned} f(x)&=sum_{i=1}^ny_iprod_{j eq i}frac{x-x_j}{x_i-x_j} \ &=sum_{i=1}^ny_ifrac{pre_{i-1}cdot suf_{i+1}}{fac_{i-1}cdot fac_{n-i}cdot (-1)^{n-i}} end{aligned} ]

    可以线性递推求逆元,(mathcal{O}(n+log mod)) 求逆元也是可以接受的,故复杂度为 (mathcal{O}(n))

    重心拉格朗日插值法

    如果每加一次点就要询问,如何解决?推一下式子试试。

    [egin{aligned} f(x)&=sum_{i=1}^ny_iprod_{j eq i}frac{x-x_j}{x_i-x_j} \ &=sum_{i=1}^ny_ifrac{prod_{j eq i}(x-x_j)}{prod_{j eq i}(x_i-x_j)} \ &=sum_{i=1}^ny_ifrac{prod_{j=1}^n(x-x_j)}{left(prod_{j=1}^n(x_i-x_j) ight )cdot (x-x_i)} \ &=prod_{j=1}^n(x-x_j)sum_{i=1}^nfrac{y_i}{left(prod_{j=1}^n(x_i-x_j) ight )cdot (x-x_i)} end{aligned} ]

    (g=prod_{j=1}^n(x-x_j),w(i)=prod_{j=1}^n(x_i-x_j))

    则原式 (=gsum_{i=1}^nfrac{y_i}{w(i)cdot (x-x_i)}).

    每加入一个点,(mathcal{O}(n)) 计算出 (g,w(i))(mathcal{O}(n)) 计算 (f(x)) 即可。

    栗题一

    ABC208F

    给定 (n,m,k),计算 (f(n,m)) 的值,模 (10^9+7)

    [egin{aligned} displaystyle f(n, m)& = egin{cases} 0 & (n = 0) ewline \ n^K & (n gt 0, m = 0) ewline \ f(n-1, m) + f(n, m-1) & (n gt 0, m gt 0) end{cases} end{aligned} ]

    (0 leq N leq 10^{18},0 leq M leq 30,1 leq K leq 2.5 imes 10^6)

    易发现答案为 (1^k,2^k,3^k,cdots ,n^k)(m) 阶前缀和的第 (n) 项值。

    考虑 (i^k) 对答案的贡献次数,是对 (f_0=1,f_k=0(k>0)) 的数组作 (m) 阶前缀和后的第 ((n-i)) 项,也就是 (f_{n-i})。其值为 (inom{n-i+m-1}{m-1})

    从组合意义上考虑,作一阶前缀和为 (1,1,1,cdots),设后面第 (i) 阶前缀和的 (f_j=g_{i,j}),有递推式 (g_{i,j}=g_{i-1,j}+g_{i,j-1}),恰为只能向右、下走的格点计数。

    故贡献的次数为 ((1,0) o (m,n-i)) 的路径数,即为 (inom{n-i+m-1}{m-1})

    所以答案为:

    [sum_{i=0}^n i^kinom{n-i+m-1}{m-1} ]

    考虑 (i^kinom{n-i+m-1}{m-1})(i^k) 为关于 (i)(k) 次单项式。

    (inom{n-i+m-1}{m-1}=frac{(n-i+m-1)!}{(m-1)!(n-i)!}),其中 ((n-i+m-1)!) 为关于 (i)((n-i+m-1)) 次多项式,((n-i)!) 为关于 (i)((m-1)) 次多项式,((m-1)!) 是常数。

    (i^kinom{n-i+m-1}{m-1}) 是关于 (i)((m+k-1)) 次多项式。

    对于 (iin [0,n]),对这个多项式求和即为答案。

    其可以写成 (ans=a_0s_0+a_1s_1+a_2s_2+cdots+a_{m+k-1}s_{m+k-1}) 的形式,其中 (s_i=sumlimits_{j=0}^ij^k)

    由于 (s_i) 是关于 (i)((i+1)) 次多项式,故答案为关于 (n)((m+k)) 次多项式。

    (l=m+k+1),选出前 (l) 个连续的整数,算出其作 (m) 阶前缀和的答案,拉格朗日插值即可。

    时间复杂度 (mathcal{O}(mk))

    #include<iostream>
    #include<cstdio>
    #include<algorithm>
    typedef long long ll;
    const ll mod = 1000000007;
    template <typename T> T Max(T x, T y) { return x > y ? x : y; }
    template <typename T> T Min(T x, T y) { return x < y ? x : y; }
    template <typename T> T Add(T x, T y) { return (x + y >= mod) ? (x + y - mod) : (x + y); }
    template <typename T> T Mod(T x) { return (x >= mod) ? (x - mod) : (x < 0 ? (x + mod) : x); }
    template <typename T> T Mul(T x, T y) { return x * y % mod; }
    template <typename T>
    T &read(T &r) {
    	r = 0; bool w = 0; char ch = getchar();
    	while(ch < '0' || ch > '9') w = ch == '-' ? 1 : 0, ch = getchar();
    	while(ch >= '0' && ch <= '9') r = (r << 3) + (r <<1) + (ch ^ 48), ch = getchar();
    	return r = w ? -r : r;
    }
    ll qpow(ll x, ll y) { ll sumq = 1; while(y) { if(y & 1) sumq = sumq * x % mod; x = x * x % mod; y >>= 1; } return sumq; }
    const int N = 2600100;
    ll n;
    int m, k;
    ll a[N], fac[N], inv[N], pre[N], suf[N], ans;
    signed main() {
    	read(n); read(m); read(k); int l = m+k+1;
    	if(!n) { puts("0");	return 0; }
    	for(int i = 1; i <= l; ++i) a[i] = qpow(i, k);
    	for(int j = 1; j <= m; ++j)
    		for(int i = 1; i <= l; ++i)
    			a[i] = Add(a[i], a[i-1]);
    	inv[0] = fac[0] = 1; for(int i = 1; i <= l; ++i) fac[i] = fac[i-1] * i % mod;
    	inv[l] = qpow(fac[l], mod-2); n %= mod;
    	for(int i = l-1; i; --i) inv[i] = inv[i+1] * (i+1) % mod;
    	pre[0] = 1; for(int i = 1; i <= l; ++i) pre[i] = pre[i-1] * (n - i) % mod;
    	suf[l+1] = 1; for(int i = l; i; --i) suf[i] = suf[i+1] * (n - i) % mod;
    	for(int i = 1; i <= l; ++i)
    		ans = Add(ans, Mod(a[i] * pre[i-1] % mod * suf[i+1] % mod * inv[i-1] % mod * inv[l-i] % mod * (((l-i)&1) ? -1 : 1)));
    	printf("%lld
    ", ans);
    	return 0;
    }
    
  • 相关阅读:
    浮动清除
    解剖JavaScript中的null和undefined【转】
    关于innerHTML以及html2dom
    javascript 作用域
    4390. 【GDOI2016模拟3.16】图计数 (Standard IO)
    5049. 【GDOI2017模拟一试4.11】腐女的生日
    4273_NOIP2015模拟10.28B组_圣章-精灵使的魔法语
    jzoj_5631_(NOI2018模拟4.5)_A
    jzoj_1001_最难的问题_Floyd
    jzoj_3385_黑魔法师之门
  • 原文地址:https://www.cnblogs.com/do-while-true/p/14971131.html
Copyright © 2011-2022 走看看