如何用Matlab进行多项式除法运算

Posted

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了如何用Matlab进行多项式除法运算相关的知识,希望对你有一定的参考价值。

求多项式x4+8x3-10除以多项式2x2-x+3的结果
貌似要用deconv这个指令,但不知道怎么用

1、鼠标右击打开桌面上的matlab程序,如下图所示,matlab运行需要一定的时间,不要着急。

2、除法运算的调用:调用格式如下:【k,r】=deconv(p,q)其中k返回的是多项式p初以q的商,r是余式。

3、输入程序:clear all关闭所有正在运行的程序,这个是很重要的。

4、下面我们以两组向量为例子作为演示程序:p1=[2 3 4 0 -2];p2=[0 0 8 -5 6];注意向量的表达的方法,注意其的格式。

5、输入程序:p=p1+p2;poly2sym(p)先令成一个向量p,做求和;poly2sym(p)函数我们之前说过,是一种创建向量的方法。

参考技术A [q,r]=deconv([1 8 0 0 -10],[2 -1 3])%q是商,r是余数。数组从后到前表示从〇次项到各高次项的系数本回答被提问者采纳 参考技术B 多项式带余除法:[q,r]=deconv(x,y) 多项式y被x除,q是商式,r是余式。
例如:
>> [q,r]=deconv(b,a) %带余除法
q =
1 2 1
r =
0 0 0 1 -2
>> q=poly2str(q,'x') %商式
q =
x^2 + 2 x + 1
>> r=poly2str(r,'x') %余式
r =
x - 2

多项式 - 除法与取模

一类问题:给定一个 (n) 次多项式 (F(x)) 和一个 (m) 次多项式 (G(x)),请求出多项式 (Q(x))(R(x)),满足以下条件:

  • (Q(x)) 次数为 (n?m)(R(x)) 次数小于 (m)
  • (F(x)=Q(x)?G(x)+R(x))

所有的运算在模 (998244353) 意义下进行。
 

多项式求逆

放一篇博客。直接做多项式的除法是无从下手的,或者说没有很优秀的方法可以优化它的复杂度。

然而逆元是个很方便的东西,乘法和除法是一对互逆运算,于是想到可以找到多项式 (G(x)) 的逆元 (G^{-1}(x)) 而转变成 (NTT) 解决的多项式乘问题。

我们设 (B(x))(A(x)) 在模 (x^n) 意义下的逆元,有 (A(x)B(x) equiv 1 (mod x^n)),下面列出几条规律。

(n=1) 时,(A(x)=c)(c) 是一个常数,则 (B(x)=c^{-1})

(n>1) 时,假设在 ((mod x^{lceil frac{n}{2} ceil)}) 意义下 (A(x)) 的逆元是 (B^′(x)),并且我们已经求出,那么:(A(x)B^′(x) equiv 1(mod x^{lceil frac{n}{2} ceil}))

而且在 (mod x^{lceil frac{n}{2} ceil}) 下也有:(A(x)B(x) equiv 1(mod x^{lceil frac{n}{2} ceil}))

然后就可以得到:(B(x)?B^′(x) equiv 0(mod x^{lceil frac{n}{2} ceil}))

两边平方 (B^2(x)?2B^′(x)B(x)+B^{′2}(x) equiv 0(mod x^n))

然后同时乘上 (A(x)),移项可以得到 (B(x) equiv 2B^′(x)?A(x)B^{′2}(x)(mod x^n)),由此式即可得到逆元。
 

(Code:)

#include <cmath>
#include <queue>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long u64;

const int mod = 998244353;
const int maxn = 600000 + 10;
int n, r[maxn], bit[maxn]; u64 f[maxn], inv[maxn], tmp[maxn];

inline int read() {
  register char ch = 0; register int w = 0, x = 0;
  while( !isdigit(ch) ) w |= (ch == ‘-‘), ch = getchar();
  while( isdigit(ch) ) x = (x * 10) + (ch ^ 48), ch = getchar();
  return w ? -x : x;
}

inline u64 Fast_pow(u64 a, int p) {
  u64 x = a, ans = 1ll;
  for( ; p; x = x * x % mod, p = p >> 1) if( p & 1 ) ans = x * ans % mod;
  return ans;
}

inline void Number_theory_transform(u64 *a, int limit, int type) {
  int rado = bit[limit];
  for(int i = 0; i < limit; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (rado - 1));
  for(int i = 0; i < limit; ++i) if( i < r[i] ) swap(a[i], a[r[i]]);
  for(int mid = 1; mid < limit; mid = mid << 1) {
    u64 Base_p = Fast_pow(3ll, (mod - 1) / (mid << 1));
    if( type == -1 ) Base_p = Fast_pow(Base_p, mod - 2);
    for(int l = 0, length = mid << 1; l < limit; l = l + length) {
      for(int k = 0, p = 1; k < mid; ++k, p = p * Base_p % mod) {
        u64 x = a[l + k], y = p * a[l + mid + k] % mod;
        a[l + k] = (x + y) % mod, a[l + mid + k] = (x - y + mod) % mod;
      }
    }
  }
  if( type == -1 ) for(int i = 0; i < limit; ++i) a[i] = a[i] * Fast_pow(limit, mod - 2) % mod;
}

inline void Inverse_transform(u64 *a, u64 *b, int limit) {
  if( limit == 1 ) { b[0] = Fast_pow(a[0], mod - 2); return ; }
  Inverse_transform(a, b, (limit + 1) >> 1);
  int len = limit << 1;
  for(int i = 0; i < limit; ++i) tmp[i] = a[i];
  for(int i = limit; i < len; ++i) tmp[i] = 0;
  Number_theory_transform(tmp, len, 1);
  Number_theory_transform(b, len, 1);
  for(int i = 0; i < len; ++i) b[i] = (2 * b[i] % mod - (b[i] * b[i] % mod) * tmp[i] % mod + mod) % mod;
  Number_theory_transform(b, len, -1);
  for(int i = limit; i < len; ++i) b[i] = 0;
}

int main(int argc, char const *argv[])
{
  freopen("..\nanjolno.in", "r", stdin);
  freopen("..\nanjolno.out", "w", stdout);

  scanf("%d", &n);
  for(int i = 0; i < n; ++i) f[i] = read();
  int limit = 1, rado = 0;
  while( limit <= n ) limit = limit << 1, ++rado;
  for(int i = 0; i < 20; ++i) bit[1 << i] = i;
  Inverse_transform(f, inv, limit);
  for(int i = 0; i < n; ++i) printf("%lld ", inv[i]);

  fclose(stdin), fclose(stdout);
  return 0;
}

 

多项式除法与取模

有了逆元,事情就简单了很多,但是依然只是个开始。

考虑最初的式子 (F(x)=Q(x)?G(x)+R(x)),要求得 (Q(x)) 最大的障碍就是 (R(x)),于是就想把 (R(x)) 通过一定的方法试试能不能消掉。

把上式的 (x) 全部用 (frac{1}{x}) 来替代,然后等式两边同时乘上 (x^n),这意味着对于一个多项式其系数会全部反转,等式变成了:

(F^R(x)=Q^R(x)G^R(x)+x^{n?m+1}R^R(x))

这样,(Q(x)) 是要求的元素之一,反转后次数仍然不高于 (n?m),而 (x^{n?m+1}R(x)) 这个部分的最低次项高于 (n?m),放到 (mod x^{n?m+1}) 意义下,(R(x)) 的影响被消除了。

于是 (F^R(x)=Q^R(x)G^R(x)(mod x^{n?m+1}))

这样就只需要求一个 (G^R(x)) 的逆元,就可以利用多项式常用的倍增方法得到 (Q^R(x)),那么显然,(R(x)=F(x)-G(x)Q(x))
 

(Code:)

inline int Limit(int len) { for(int i = 1; ; i = i << 1) if( i > len ) return i; }

inline void Multiply_transform(u64 *a, u64 *b, u64 *c, int len_1, int len_2) {
  int len = Limit(len_1 + len_2);
  for(int i = 0; i <= len_1; ++i) tmp_1[i] = a[i];
  for(int i = 0; i <= len_2; ++i) tmp_2[i] = b[i];
  Number_theory_transform(tmp_1, len, 1);
  Number_theory_transform(tmp_2, len, 1);
  for(int i = 0; i < len; ++i) c[i] = tmp_1[i] * tmp_2[i] % mod, tmp_1[i] = tmp_2[i] = 0;
  Number_theory_transform(c, len, -1);
}

int main(int argc, char const *argv[])
{
  for(int i = 0; i < 20; ++i) bit[1 << i] = i;
  scanf("%d%d", &n, &m);
  for(int i = 0; i <= n; ++i) f_1[i] = read();
  for(int i = 0; i <= m; ++i) f_2[i] = read();
  reverse(f_1, f_1 + n + 1), reverse(f_2, f_2 + m + 1);
  Inverse_transform(f_2, inv, Limit(n - m));
  Multiply_transform(f_1, inv, ans, n - m, n - m);
  reverse(ans, ans + n - m + 1);
  for(int i = 0; i <= n - m; ++i) printf("%lld ", ans[i]); printf("
");
  reverse(f_1, f_1 + n + 1), reverse(f_2, f_2 + m + 1);
  Multiply_transform(ans, f_2, mul, n - m, m);
  for(int i = 0; i < m; ++i) printf("%lld ", (f_1[i] - mul[i] + mod) % mod); printf("
");

  return 0;
}

 
                     与君相遇,乃思长生。





以上是关于如何用Matlab进行多项式除法运算的主要内容,如果未能解决你的问题,请参考以下文章

MATLAB编程及应用-第5章 多项式与数据分析

MATLAB编程与应用系列-第5章 多项式与数据分析(1)

如何用matlab线性回归分析?

在多元线性回归中,如何用matlab求得各个变量的T统计值及其p值?

如何用C语言实现一元多项式简单计算器的设计

多项式 - 除法与取模