- // Trial division with polynomial interpolation
- int i, k;
-
- // Compute values at evaluation points 0..adeg
- vector<numeric> alpha; alpha.reserve(adeg + 1);
- exvector u; u.reserve(adeg + 1);
- numeric point = _num0();
- ex c;
- for (i=0; i<=adeg; i++) {
- ex bs = b.subs(*x == point);
- while (bs.is_zero()) {
- point += _num1();
- bs = b.subs(*x == point);
- }
- if (!divide_in_z(a.subs(*x == point), bs, c, var+1))
- return false;
- alpha.push_back(point);
- u.push_back(c);
- point += _num1();
- }
-
- // Compute inverses
- vector<numeric> rcp; rcp.reserve(adeg + 1);
- rcp.push_back(_num0());
- for (k=1; k<=adeg; k++) {
- numeric product = alpha[k] - alpha[0];
- for (i=1; i<k; i++)
- product *= alpha[k] - alpha[i];
- rcp.push_back(product.inverse());
- }
-
- // Compute Newton coefficients
- exvector v; v.reserve(adeg + 1);
- v.push_back(u[0]);
- for (k=1; k<=adeg; k++) {
- ex temp = v[k - 1];
- for (i=k-2; i>=0; i--)
- temp = temp * (alpha[k] - alpha[i]) + v[i];
- v.push_back((u[k] - temp) * rcp[k]);
- }
-
- // Convert from Newton form to standard form
- c = v[adeg];
- for (k=adeg-1; k>=0; k--)
- c = c * (*x - alpha[k]) + v[k];
-
- if (c.degree(*x) == (adeg - bdeg)) {
- q = c.expand();
- return true;
- } else
- return false;
+ // Trial division with polynomial interpolation
+ int i, k;
+
+ // Compute values at evaluation points 0..adeg
+ vector<numeric> alpha; alpha.reserve(adeg + 1);
+ exvector u; u.reserve(adeg + 1);
+ numeric point = _num0();
+ ex c;
+ for (i=0; i<=adeg; i++) {
+ ex bs = b.subs(*x == point);
+ while (bs.is_zero()) {
+ point += _num1();
+ bs = b.subs(*x == point);
+ }
+ if (!divide_in_z(a.subs(*x == point), bs, c, var+1))
+ return false;
+ alpha.push_back(point);
+ u.push_back(c);
+ point += _num1();
+ }
+
+ // Compute inverses
+ vector<numeric> rcp; rcp.reserve(adeg + 1);
+ rcp.push_back(_num0());
+ for (k=1; k<=adeg; k++) {
+ numeric product = alpha[k] - alpha[0];
+ for (i=1; i<k; i++)
+ product *= alpha[k] - alpha[i];
+ rcp.push_back(product.inverse());
+ }
+
+ // Compute Newton coefficients
+ exvector v; v.reserve(adeg + 1);
+ v.push_back(u[0]);
+ for (k=1; k<=adeg; k++) {
+ ex temp = v[k - 1];
+ for (i=k-2; i>=0; i--)
+ temp = temp * (alpha[k] - alpha[i]) + v[i];
+ v.push_back((u[k] - temp) * rcp[k]);
+ }
+
+ // Convert from Newton form to standard form
+ c = v[adeg];
+ for (k=adeg-1; k>=0; k--)
+ c = c * (*x - alpha[k]) + v[k];
+
+ if (c.degree(*x) == (adeg - bdeg)) {
+ q = c.expand();
+ return true;
+ } else
+ return false;