多项式开根

最新推荐文章于 2021-12-28 18:29:47 发布
最新推荐文章于 2021-12-28 18:29:47 发布
阅读量555 收藏
点赞数
分类专栏： 多项式及生成函数
版权声明：本文为博主原创文章，遵循 CC 4.0 BY 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_45483201/article/details/113407372
版权

多项式开根

给定多项式 g ( x ) g(x) g(x)，求 f ( x ) f(x) f(x)，满足 f 2 ( x ) = g ( x ) f ^ 2(x) = g(x) f2(x)=g(x)。

假设我们已经得到了 g ( x ) g(x) g(x)，膜 x ⌈ n 2 ⌉ x ^{\lceil \frac{n}{2} \rceil} x⌈2n​⌉下的根 f 0 ( x ) f_0 (x) f0​(x)，要求膜 x n x ^ n xn下的根 f ( x ) f(x) f(x)

有 f 0 2 ( x ) ≡ g ( x ) ( m o d x ⌈ n 2 ⌉ ) f_0 ^2(x) \equiv g(x) \pmod {x ^{\lceil \frac{n}{2} \rceil}} f02​(x)≡g(x)(modx⌈2n​⌉)

移项再开方有 ( f 0 2 ( x ) − g ( x ) ) 2 ≡ 0 ( m o d x n ) \left(f_0 ^2(x) - g(x) \right) ^ 2 \equiv 0 \pmod {x ^ n} (f02​(x)−g(x))2≡0(modxn)

则， ( f 0 2 ( x ) + g ( x ) ) 2 ≡ 4 f 0 2 ( x ) g ( x ) ( m o d x n ) \left( f_0 ^ 2(x) + g(x) \right) ^ 2 \equiv 4 f_0 ^ 2(x) g(x) \pmod {x ^ n} (f02​(x)+g(x))2≡4f02​(x)g(x)(modxn)

g ( x ) ≡ ( f 0 2 ( x ) + g ( x ) 2 f 0 ( x ) ) 2 ( m o d x n ) g(x) \equiv \left(\frac{f_0 ^ 2(x) + g(x)}{2f_0 (x)} \right) ^ 2 \pmod {x ^ n} g(x)≡(2f0​(x)f02​(x)+g(x)​)2(modxn)

所以 f ( x ) ≡ f 0 2 ( x ) + g ( x ) 2 f 0 ( x ) ( m o d x n ) f(x) \equiv \frac{f_0 ^ 2(x) + g(x)}{2f_0(x)} \pmod {x ^ n} f(x)≡2f0​(x)f02​(x)+g(x)​(modxn)。

所以有 f ( x ) ≡ 2 − 1 f 0 ( x ) + 2 − 1 f 0 − 1 ( x ) g ( x ) ( m o d x n ) f(x) \equiv 2 ^{-1} f_0 (x) + 2 ^{-1} f_0 ^{-1}(x) g(x) \pmod {x ^ n} f(x)≡2−1f0​(x)+2−1f0−1​(x)g(x)(modxn)，

对于 g ( 0 ) = 1 g(0) = 1 g(0)=1的特殊情况

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 5e6 + 10, mod = 998244353, inv2 = mod + 1 >> 1;

int a[N], b[N], c[N], d[N], r[N];

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll *  a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv1(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = quick_pow(a[0], mod - 2);
    return ;
  }
  polyinv1(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * c[i] * b[i] % mod + mod) % mod;
    b[i] = 1ll * b[i] * cur % mod;
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

void polysqrt(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = 1;
    return ;
  }
  polysqrt(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  for (int i = 0; i < lim; i++) {
    d[i] = 0;
  }
  polyinv1(b, d, n);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  get_r(lim);
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  NTT(d, lim, 1);
  for (int i = 0; i < lim; i++) {
    b[i] = (1ll * inv2 * b[i] % mod + 1ll * inv2 * d[i] % mod * c[i] % mod) % mod; 
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n;
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  polysqrt(a, b, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}

二次剩余解一般情况

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int N = 5e6 + 10, mod = 998244353, inv2 = mod + 1 >> 1;

int a[N], b[N], c[N], d[N], r[N];

namespace Quadratic_residue {
  struct Complex {
    int r, i;

    Complex(int _r = 0, int _i = 0) : r(_r), i(_i) {}
  };

  int I2;

  Complex operator * (const Complex &a, Complex &b) {
    return Complex((1ll * a.r * b.r % mod  + 1ll * a.i * b.i % mod * I2 % mod) % mod, (1ll * a.r * b.i % mod + 1ll * a.i * b.r % mod) % mod);
  }

  Complex quick_pow(Complex a, int n) {
    Complex ans = Complex(1, 0);
    while (n) {
      if (n & 1) {
        ans = ans * a;
      }
      a = a * a;
      n >>= 1;
    }
    return ans;
  }

  int get_residue(int n) {
    mt19937 e(233);
    if (n == 0) {
      return 0;
    }
    if(quick_pow(n, (mod - 1) >> 1).r == mod - 1) {
      return -1;
    }
    uniform_int_distribution<int> r(0, mod - 1);
    int a = r(e);
    while(quick_pow((1ll * a * a % mod - n + mod) % mod, (mod - 1) >> 1).r == 1) {
      a = r(e);
    }
    I2 = (1ll * a * a % mod - n + mod) % mod;
    int x = quick_pow(Complex(a, 1), (mod + 1) >> 1).r, y = mod - x;
    if(x > y) swap(x, y);
    return x;
  }
}

int quick_pow(int a, int n) {
  int ans = 1;
  while (n) {
    if (n & 1) {
      ans = 1ll * ans * a % mod;
    }
    a = 1ll *  a * a % mod;
    n >>= 1;
  }
  return ans;
}

void get_r(int lim) {
  for (int i = 0; i < lim; i++) {
    r[i] = (i & 1) * (lim >> 1) + (r[i >> 1] >> 1);
  }
}

void NTT(int *f, int lim, int rev) {
  for (int i = 0; i < lim; i++) {
    if (i < r[i]) {
      swap(f[i], f[r[i]]);
    }
  }
  for (int mid = 1; mid < lim; mid <<= 1) {
    int wn = quick_pow(3, (mod - 1) / (mid << 1));
    for (int len = mid << 1, cur = 0; cur < lim; cur += len) {
      int w = 1;
      for (int k = 0; k < mid; k++, w = 1ll * w * wn % mod) {
        int x = f[cur + k], y = 1ll * w * f[cur + mid + k] % mod;
        f[cur + k] = (x + y) % mod, f[cur + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if (rev == -1) {
    int inv = quick_pow(lim, mod - 2);
    reverse(f + 1, f + lim);
    for (int i = 0; i < lim; i++) {
      f[i] = 1ll * f[i] * inv % mod;
    }
  }
}

void polyinv1(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = quick_pow(a[0], mod - 2);
    return ;
  }
  polyinv1(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  get_r(lim);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  for (int i = 0; i < lim; i++) {
    int cur = (2 - 1ll * c[i] * b[i] % mod + mod) % mod;
    b[i] = 1ll * b[i] * cur % mod;
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

void polysqrt(int *a, int *b, int n) {
  if (n == 1) {
    b[0] = Quadratic_residue::get_residue(a[0]);
    return ;
  }
  polysqrt(a, b, n + 1 >> 1);
  int lim = 1;
  while (lim < 2 * n) {
    lim <<= 1;
  }
  for (int i = 0; i < lim; i++) {
    d[i] = 0;
  }
  polyinv1(b, d, n);
  for (int i = 0; i < n; i++) {
    c[i] = a[i];
  }
  for (int i = n; i < lim; i++) {
    c[i] = 0;
  }
  get_r(lim);
  NTT(b, lim, 1);
  NTT(c, lim, 1);
  NTT(d, lim, 1);
  for (int i = 0; i < lim; i++) {
    b[i] = (1ll * inv2 * b[i] % mod + 1ll * inv2 * d[i] % mod * c[i] % mod) % mod; 
  }
  NTT(b, lim, -1);
  for (int i = n; i < lim; i++) {
    b[i] = 0;
  }
}

int main() {
  // freopen("in.txt", "r", stdin);
  // freopen("out.txt", "w", stdout);
  // ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
  int n;
  scanf("%d", &n);
  for (int i = 0; i < n; i++) {
    scanf("%d", &a[i]);
  }
  polysqrt(a, b, n);
  for (int i = 0; i < n; i++) {
    printf("%d%c", b[i], i + 1 == n ? '\n' : ' ');
  }
  return 0;
}
_lifehappy_
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
开根号的笔算算法图解_【learning】多项式开根详解+模板
weixin_39869593的博客
12-20 1689
概述多项式开跟是一个非常重要的知识点，许多多项式题目都要用到这一算法。用快速数论变换，多项式求逆元和倍增法可以在$O(n log n)$的时间复杂度下求出一个$n$次多项式的开根。前置技能快速数论变换(NTT)，多项式求逆元，二次剩余。多项式的开根给定一个多项式$A(x)$，其次数为$deg_A$，若存在一个多项式$B(x)$，使其满足$deg_B≤deg_A$，且$ B^2(x) \equiv ...
计算切比雪夫多项式的根。：计算任意次数的切比雪夫多项式的根。-matlab开发
06-01
U_n(x)的根位于开区间(-1, 1)内，除了n=0时，U_0(x) = 1，没有实根。 MATLAB提供了一个名为`chebroots`的函数，可以用来计算切比雪夫多项式的根。这个函数可以处理第一类或第二类切比雪夫多项式，并能选择是否进行...
【learning】多项式开根详解+模板
diaoyoutun2652的博客
07-26 200
概述 多项式开跟是一个非常重要的知识点，许多多项式题目都要用到这一算法。 用快速数论变换，多项式求逆元和倍增法可以在$O(n log n)$的时间复杂度下求出一个$n$次多项式的开根。 前置技能 快速数论变换（NTT），多项式求逆元，二次剩余。 多项式的开根 给定一个多项式$A(x)$，其次数为$deg_A$，若存在一个多项式$B(x)$，使其满足$deg_B≤d...
多项式开根小结
a_forever_dream的博客
05-20 292
多项式开根 跟 expexpexp 套路是一样的，套一个牛顿迭代就出来了。 给出多项式 A(x)A(x)A(x)，求出一个多项式 G(x)G(x)G(x) 满足 G2(x)≡A(x)(modxn)G^2(x)\equiv A(x)\pmod {x^n}G2(x)≡A(x)(modxn)。 设 F(G(x))=G2(x)−A(x)F(G(x))=G^2(x)-A(x)F(G(x))=G2(x)−A(x)，那么就是要求出 FFF 的零点。 考虑倍增，当 n=1n=1n=1 时，这题保证了 A(x)=1A(x)=
【luogu P5205】【模板】多项式开根（牛顿迭代）（多项式）
欢迎您的光顾awa
12-28 213
给你一个多项式 A(x)，要你找到一个多项式 B(x) 使得 B^2(x)=A(x) (mod x^n)。 如果有多个解选零次项系数最小的，然后保证 A(0)=1。
多项式牛顿迭代和多项式开根
litble的成(tui)长(fei)史
07-02 2746
orz picks 前置skill 基础求导 （全世界应该只有我不会求导了吧……=。=，但考虑到可能有小学生在看这篇文章，所以还是讲一下吧 导数，是微积分中的重要基础概念，关于它的具体内容就不讲了，在本文中，你只要知道，求导是用来拟合一个函数在各处的斜率的东西，并且有： 若f(x)=xaf(x)=xaf(x)=x^{a}，则f′(x)=axa−1f′(x)=axa−1f'(x)=ax^{...
多项式求逆与多项式开根
某蒟蒻的博客
02-13 1102
闲着没事干研究些黑科技（雾） 多项式求逆 求 A(x)*B(x)==1(mod x^n) 其中n为A(x)，B(x)的度的较大值 已知A(x)求B(x)，B(x)=A(x)^(-1)(mod x^n) 假设n=1，则B(x)=A(x)常数项在mod p 意义下的的逆元 假设n&gt;1 已知 A(x)*B’(x)==1(mod x^(n/2)) A(x)*B(x)==1(mod...
[学习笔记]多项式开根
weixin_30236595的博客
01-28 87
还是倍增思想 每一层 还要套一个多项式求逆 所以O(nlog^2n)，常数也不小 数组比较多 再trick一下 得到： $T=(T'+F*inv(T'))*inv2$ 可以只算一次多项式求逆，一次NTT void sqr(int *f,int *g,int n){ if(n==1){ g[0]=1;return; } s...
洛谷P5205 【模板】多项式开根
ACoder
05-29 284
牛顿迭代、多项式开根
多项式开根学习笔记
QQ82272760的博客
03-24 167
文章目录问题思路代码 问题 给定一个 n−1n-1n−1 次多项式 A(x)A(x)A(x) ，求一个 n−1n-1n−1 次多项式 B(x)B(x)B(x)，使得 B2(x)=A(x)B^2(x)=A(x)B2(x)=A(x)。 思路 考虑倍增。假设已经求出 n−1n-1n−1 次多项式 B0(x)B_0(x)B0​(x) 使得 B02(x)≡A(x) (mod xn)B_0^2(x)\equiv A(x)\ (\text{mod}\ x^n)B02​(x)≡A(x) (mo
beisaier.zip_图形描点_多项式绘制
09-14
通过以上步骤，我们可以清晰地看到2阶到7阶Legendre多项式的图形表现，从而理解其特性，例如奇偶性、单调性以及根的分布情况。这种方法不仅适用于Legendre多项式，也适用于其他类型的多项式函数，是一种有效的数值...
rational fractional polynomial_matlab_模态参数_polynomial_有理多项式_
10-01
这一步可能需要用到MATLAB的`roots`函数来寻找特征根。 4. 结果验证：计算拟合后的频率响应并与原始数据比较，检查拟合质量。此外，还可以通过比较模态参数与预期值或物理意义来验证结果的合理性。 在提供的压缩包...
Matlab多项式运算_matlab源码.rar
12-12
`roots`函数可以找到多项式的实根或复根，例如`r = roots(p1)`。 6. **多项式插值**： 当需要通过一系列数据点构造多项式时，可以使用`polyfit`函数。例如，`p = polyfit(x, y, n)`会找到最佳拟合的n次多项式。 ...
毕业设计MATLAB_有理分数多项式法.zip
最新发布
05-27
4. **稳定性和性能分析**：通过根轨迹分析或频域分析（如`bode`图，`nyquist`图），可以评估系统的稳定性和性能。MATLAB的控制系统工具箱提供了这些功能。 5. **仿真和验证**：最后，利用MATLAB的仿真功能，比如`...
快速傅里叶变换（完整推导过程 + 模板）
lifehappy
01-28 3317
快速傅里叶变换 多项式表示 系数表示法： 一个nnn次多项式可以用n+1n + 1n+1个系数表示出来：f(x)=a0+a1x+a2x2+⋯+an−1xn−1+anxnf(x) = a_0 + a_1 x + a_2 x ^ 2 + \dots + a_{n - 1} x ^{n- 1} + a_n x ^nf(x)=a0​+a1​x+a2​x2+⋯+an−1​xn−1+an​xn。 点值表示法： 通过线性代数，高斯消元我们可以知道，一个nnn次多项式可以通过n+1n + 1n+1个点联立方程组解得： f(
伯努利数（详解 + 例题 ：P3711 仓鼠的数学题）
lifehappy
03-05 1632
伯努利数 定义Sk(n)=∑i=0n−1ikS_k(n) = \sum\limits_{i = 0} ^{n - 1} i ^ kSk​(n)=i=0∑n−1​ik。 从二项式出发 (0+1)k+1=∑i=0kCk+1i0i+0k+1⋮(n−1+1)k+1=∑i=0kCk+1i(n−1)i+(n−1)k+1把次方k+1的移项，再整体相加，得nk+1=∑i=0kCk+1iSi(n)nk+1=∑i=0k−1Ck+1iSi(n)+(k+1)Sk(n)Sk(n)=1k+1(nk+1−∑i=0k−1Ck+1iSi(n
多项式对数函数ln f(x)
lifehappy
01-30 1581
多项式对数函数ln⁡f(x)\ln f(x)lnf(x) 如果存在解必然有[x0]f(x)=1[ x ^ 0]f(x) = 1[x0]f(x)=1， 对ln⁡f(x)\ln f(x)lnf(x)求导，有dln⁡f(x)dx≡f′(x)f(x)(modxn)\frac{d \ln f(x)}{dx} \equiv \frac{f'(x)}{f(x)} \pmod {x ^ n}dxdlnf(x)​≡f(x)f′(x)​(modxn)， dxdxdx乘到右边，再求积分有： ∫dln⁡f(x)≡∫f′(x))f(
多项式除法，多项式取模
lifehappy
02-05 1389
多项式除法 给定一个nnn次多项式F(x)F(x)F(x)和mmm次多项式G(x)G(x)G(x)，要求R(x),Q(x)R(x), Q(x)R(x),Q(x)，满足F(x)=R(x)G(x)+Q(x)F(x) = R(x)G(x) + Q(x)F(x)=R(x)G(x)+Q(x)。 R(x)R(x)R(x)是一个n−mn - mn−m阶多项式，Q(x)Q(x)Q(x)是一个小于mmm阶的多项式。 有F(x)≡R(x)G(x)+Q(x)(modxn+1)F(1x)≡R(1x)G(1x)+Q(1x)(modx
MATLAB实用教程：由根创建多项式
"MATLAB实用教程PPT，讲解了如何由根创建多项式，以及MATLAB的基本功能和使用环境" 在MATLAB中，函数`poly()`是一个非常实用的工具，它允许用户根据给定的根来创建对应的多项式。这个函数的核心思想是基于Vieta's ...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值