Symbolic Integration (Manuel Bronstein)

Well, now that i have a new book, Symbolic Integration I of Manuel Bronstein, i will implement some of the code examples on the initial chapter: polydivide, euclidean algorithm, square free factorization, all of them using GiNaC library.

PolyDivide

This one find the unique q and r in a field K[z] that holds a = b*q + r

lst PolyDivide(ex A, ex B, symbol x)
{
    ex Q = 0;
    ex R = A;

    int d = R.degree(x) - B.degree(x);

    while ((not R.is_zero()) and (d >=0))
    {
        ex T = R.lcoeff(x) / B.lcoeff(x) * pow(x,d);
        Q = expand(Q + T);
        R = expand(R - B*T);
        d = R.degree(x) - B.degree(x);
    }
    return lst(Q,R);
}

PolyPseudoDivide
Same that last, but in a integral domain, like Z.

lst PolyPseudoDivide(ex A, ex B, symbol x)
{
    ex b = B.lcoeff(x);
    int N = A.degree(x) - B.degree(x) + 1;
    ex Q = 0;
    ex R = A;

    int d = R.degree(x) - B.degree(x);

    while ((not R.is_zero()) and (d >=0))
    {
        ex T = R.lcoeff(x) * pow(x,d);
        N -= 1;
        Q = expand(b * Q + T);
        R = expand(b * R - T * B);
        d = R.degree(x) - B.degree(x);
    }

    b = pow(b, N);
    return lst( b * Q, b * R);
}

Well, testing both algorithm against the quo and rem that do the same computation in GiNaC, we have this:



The test is a division between two polynomials of the same n degree randomly generated, in x-axis is n and in y-axis is ellapsed time of the compute code in sec.

int main()
{
    symbol x("x");
    int i = 0;
    ex a=0, b=0;

    clock_t s, e;
    double t1,t2, t3, t4, t5;

    for(int n = 1; n < 100000 ; n += 1000)
    {
        srand ( time(NULL) );

        for(; i < n; i++)
        {
            a += (rand() % 200 - 100) * pow(x,i);
            b += (rand() % 200 - 100) * pow(x,i);
        }

        s = clock();
        PolyDivide(a,b,x);
        e = clock();
        t1 = (e-s);

        s = clock();
        PolyPseudoDivide(a,b,x);
        e = clock();
        t2 = (e-s);

        s = clock();
        quo(a,b,x); rem(a,b,x);
        e = clock();
        t3 = (e-s);

        ex q;
        s = clock();
        q = quo(a,b,x);a - b*q;
        e = clock();
        t4 = (e-s);

        ex r;
        s = clock();
        r = rem(a,b,x);(a-r)/b;
        e = clock();
        t5 = (e-s);


        t1 /= CLOCKS_PER_SEC;
        t2 /= CLOCKS_PER_SEC;
        t3 /= CLOCKS_PER_SEC;
        t4 /= CLOCKS_PER_SEC;
        t5 /= CLOCKS_PER_SEC;

        cout << n << " " << t1 << " " << t2 << " " << t3 << " " << t4 << " " << t5 << endl;
    }

    return 0;
}