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

#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);
}