void poldiv(double u[], int n, double v[], int nv, double q[], double r[])
{
    for (int j = 1; j<=n; j++)
	{
        r[j] = u[j];
        q[j] = 0.0;
    }
    for (int k = n - nv; k>=0; k--)
	{
        q[k + 1] = r[nv + k] / v[nv];
        for (int j = nv + k - 1; j>=k + 1; j--)
		{
            r[j] = r[j] - q[k + 1] * v[j - k];
        }
    }
    r[nv] = 0.0;
}