Polynomial Greatest Common Divisor C++

Question

Here is my attempt to implement a program that finds GCD of two polynomials. I realize there is a problem in division method. A while loop decrementing the degree of a resulting polynomial in division() goes into "infinity" in some cases, and I can't understand in which exactly.

Any clues on what goes wrong here?

#include<iostream>
#include<stdlib.h>



using namespace std;

//polynomial------------------------------
struct polynomial {
    float *coeff;
    int degree;
};

/*function declaration */
int get_data(polynomial *P); 
int display(polynomial *P);
int division(polynomial *P, polynomial *Q, polynomial *H, polynomial *R);
int polcpy(polynomial *p, polynomial *q);
void GCDPol(polynomial *P,polynomial *Q, polynomial *R);

//GCD of two polynomials------------------
void GCDpol(polynomial *P,polynomial *Q, polynomial *R) {
    polynomial *res, *v;
    res = (polynomial *) calloc(1,sizeof(polynomial));  
    v = (polynomial *) calloc(1,sizeof(polynomial));  
    while(1) {
        division(P, Q, v, R);
        if(R->degree==0 && R->coeff[0]==0)
            break;
        else
            polcpy(P,Q);
            polcpy(Q,R);
    }
    polcpy(R,Q);
    free(res);
    free(v);
}

//pol copy--------------------------------
int polcpy(polynomial *p, polynomial *q) {
    p->degree=q->degree;
    p->coeff = new float[p->degree + 1];
    for (int i=0; i<=p->degree; i++) 
        p->coeff[i]=q->coeff[i];
    return 0;
}

//division--------------------------------
int division(polynomial *P, polynomial *Q, polynomial *H, polynomial *R) {
    float u;
    int x;
    polynomial *nh, *nr;
    nh = (polynomial *) calloc(1,sizeof(polynomial));
    nr = (polynomial *) calloc(1,sizeof(polynomial));

    /*Euclidian Long Division*/
    polcpy(nr, P);  
    nh->degree = P->degree - Q->degree;
    nh->coeff = new float[nh->degree + 1];
    for (int i=nh->degree; i>=0; i--)  {
        nh->coeff[i] = nr->coeff[nr->degree] / Q->coeff[Q->degree];
        for (int j=i; j <= nr->degree; j++) {
            u = nh->coeff[i] * Q->coeff[j-i];
            nr->coeff[j] = nr->coeff[j] - u;
        }
        if (nr->degree > 0)  
            nr->degree--;
    }

    /*Quotient*/
    polcpy(H, nh);  
    /*Remainder*/   
    polcpy(R, nr);
    while(R->coeff[R->degree] == 0) {
        R->degree--;
    }   
    free(nh); 
    free(nr);
    return 0;
}

//display-------------------------------
int display(polynomial *P) {
    int i, j;
    for (i = P->degree; i >= 0; i--) {
        cout << P->coeff[i] << "x^" << i;
        if ((i - 1) != -1)
            cout << "+";
    }
    cout << "\n";
    return 0;
}

//get_data------------------------------
int get_data(polynomial *P) {
    cout << "Enter Degree Of Polynomial:";
    cin >> P->degree;
    P->coeff = new float[P->degree + 1];
    for (int i = P->degree; i >= 0; i--) {
        cout << "Enter coefficient of x^" << i << ":";
        cin >> P->coeff[i];
    }
    return 0;
}


int main() {
    polynomial *P, *Q, *R, *H;
    P = (polynomial *) calloc(1,sizeof(polynomial));
    Q = (polynomial *) calloc(1,sizeof(polynomial));
    R = (polynomial *) calloc(1,sizeof(polynomial));
    H = (polynomial *) calloc(1,sizeof(polynomial));
    cout<<"GDC\n";
    get_data(P);
    get_data(Q);
    cout << "Polynomial1:";
    display(P);
    cout << "Polynomial2:";
    display(Q);
    GCDpol(P,Q,R);
    display(R);
    free(R);
    free(P);
    free(Q);
    free(H);
    return 0;
}

Answer 1

This C++ library implements polynomial division: PolynomialDivHang_test_no_hang

Answer 2

It seems likely that the line

if(R->degree==0 && R->coeff[0]==0)
        break;

is what's broken. Your coefficients are floats. Because computers are (unfortunately) finite, there will be small errors in your float computation. The code only exits the while loop if the coefficient is exactly 0. It seems likely that on some inputs, although it should divide evenly, you get R->coeff[0] = 0.0000000001 or some other very small value that is not exactly 0.

Try checking if the absolute value is within some very small tolerance (something like 10^-10 or 10^-12.) If you actually want the exact value, you'll need to look into exact precision floats, which is its own can of worms.

(Looking at the code more closely, you do exact equality checks in other places too - these should all be changed to checking the absolute value is very small.)

Answer 3

while(R->coeff[R->degree] == 0) {
    R->degree--;
}

This goes badly wrong if R is the zero polynomial.