扩展欧几里得求两多项式最大公因式

#include <iostream>#include <string.h>#include <stdio.h>#include <math.h>using namespace std;typedef long long LL;const double eps = 1e-8;const int MOD = 999983;const int N = 55;struct Poly{    int n;    LL a[N];};Poly p[25];LL gcd(LL a,LL b){    return b? gcd(b,a%b):a;}Poly delete_zero(Poly x){    int i,j;    Poly tmp;    for(i=0;i<x.n && x.a[i] == 0;i++);    for(j=0;i<x.n;i++,j++) tmp.a[j] = x.a[i];    tmp.n = j;    return tmp;}Poly poly_gcd(Poly x,Poly y){    x = delete_zero(x);    y = delete_zero(y);    Poly yy = y,tmp;    tmp.n = 0;    int i=0;    if(x.n == y.n)    {        double k = x.a[0] / y.a[0];        for(i=1;i<x.n;i++)            if(fabs(x.a[i]*1.0 - k*y.a[i]) > eps) break;        if(i == x.n) return y;    }    LL g = gcd(x.a[0],y.a[0]);    LL tx = y.a[0] / g;    LL ty = x.a[0] / g;    for(i=0;i<x.n;i++)    {        x.a[i] *= tx;        x.a[i] %= MOD;    }    for(i=0;i<y.n;i++)    {        y.a[i] *= ty;        y.a[i] %= MOD;    }    if(x.n < y.n) swap(x,y);    for(i=1;i<y.n;i++)        tmp.a[i-1] = ((x.a[i] - y.a[i])%MOD + MOD)%MOD;    for(;i<x.n;i++)        tmp.a[i-1] = x.a[i];    tmp.n = x.n - 1;    tmp = delete_zero(tmp);    if(tmp.n == 0) return yy;    return poly_gcd(y,tmp);}