123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391 |
- // Factor a polynomial
- #include "stdafx.h"
- #include "defs.h"
- extern "C" {
- #include "console.h"
- }
- static void rationalize_coefficients(int);
- static int get_factor(void);
- static void evalpoly(void);
- static void yydivpoly(void);
- static int expo;
- static U **polycoeff;
- #define POLY p1
- #define X p2
- #define Z p3
- #define A p4
- #define B p5
- #define Q p6
- #define RESULT p7
- #define FACTOR p8
- void
- factorpoly(void)
- {
- save();
- p2 = pop();
- p1 = pop();
- if (!find(p1, p2)) {
- push(p1);
- restore();
- return;
- }
- if (!ispoly(p1, p2)) {
- push(p1);
- restore();
- return;
- }
- if (!issymbol(p2)) {
- push(p1);
- restore();
- return;
- }
- push(p1);
- push(p2);
- yyfactorpoly();
- restore();
- }
- //-----------------------------------------------------------------------------
- //
- // Input: tos-2 true polynomial
- //
- // tos-1 free variable
- //
- // Output: factored polynomial on stack
- //
- //-----------------------------------------------------------------------------
- void
- yyfactorpoly(void)
- {
- int h, i;
- save();
- X = pop();
- POLY = pop();
- h = tos;
- if (isfloating(POLY))
- stop("floating point numbers in polynomial");
- polycoeff = stack + tos;
- push(POLY);
- push(X);
- expo = coeff() - 1;
- rationalize_coefficients(h);
- // for univariate polynomials we could do expo > 1
- while (expo > 0) {
- if (iszero(polycoeff[0])) {
- push_integer(1);
- A = pop();
- push_integer(0);
- B = pop();
- } else if (get_factor() == 0) {
- if (verbosing)
- printf("no factor found");
- break;
- }
- push(A);
- push(X);
- multiply();
- push(B);
- add();
- FACTOR = pop();
- if (verbosing) {
- printf("successFACTOR=");
- print(FACTOR);
- }
- // factor out negative sign (not req'd because A > 1)
- #if 0
- if (isnegativeterm(A)) {
- push(FACTOR);
- negate();
- FACTOR = pop();
- push(RESULT);
- negate_noexpand();
- RESULT = pop();
- }
- #endif
- push(RESULT);
- push(FACTOR);
- multiply_noexpand();
- RESULT = pop();
- yydivpoly();
- while (expo && iszero(polycoeff[expo]))
- expo--;
- }
- // unfactored polynomial
- push(zero);
- for (i = 0; i <= expo; i++) {
- push(polycoeff[i]);
- push(X);
- push_integer(i);
- power();
- multiply();
- add();
- }
- POLY = pop();
- if (verbosing) {
- printf("POLY=");
- print(POLY);
- }
- // factor out negative sign
- if (expo > 0 && isnegativeterm(polycoeff[expo])) {
- push(POLY);
- negate();
- POLY = pop();
- push(RESULT);
- negate_noexpand();
- RESULT = pop();
- }
- push(RESULT);
- push(POLY);
- multiply_noexpand();
- RESULT = pop();
- if (verbosing) {
- printf("RESULT=");
- print(RESULT);
- }
- stack[h] = RESULT;
- tos = h + 1;
- restore();
- }
- static void
- rationalize_coefficients(int h)
- {
- int i;
- // LCM of all polynomial coefficients
- RESULT = one;
- for (i = h; i < tos; i++) {
- push(stack[i]);
- denominator();
- push(RESULT);
- lcm();
- RESULT = pop();
- }
- // multiply each coefficient by RESULT
- for (i = h; i < tos; i++) {
- push(RESULT);
- push(stack[i]);
- multiply();
- stack[i] = pop();
- }
- // reciprocate RESULT
- push(RESULT);
- reciprocate();
- RESULT = pop();
- }
- static int
- get_factor(void)
- {
- int i, j, h;
- int a0, an, na0, nan;
- if (verbosing) {
- push(zero);
- for (i = 0; i <= expo; i++) {
- push(polycoeff[i]);
- push(X);
- push_integer(i);
- power();
- multiply();
- add();
- }
- POLY = pop();
- printf("POLY=");
- print(POLY);
- }
- h = tos;
- an = tos;
- push(polycoeff[expo]);
- divisors_onstack();
- nan = tos - an;
- a0 = tos;
- push(polycoeff[0]);
- divisors_onstack();
- na0 = tos - a0;
- if (verbosing) {
- printf("divisors of base term");
- for (i = 0; i < na0; i++) {
- printf(", ");
- print(stack[a0 + i]);
- }
- printf(" ");
- printf("divisors of leading term");
- for (i = 0; i < nan; i++) {
- printf(", ");
- print(stack[an + i]);
- }
- }
- // try roots
- for (i = 0; i < nan; i++) {
- for (j = 0; j < na0; j++) {
- A = stack[an + i];
- B = stack[a0 + j];
- push(B);
- push(A);
- divide();
- negate();
- Z = pop();
- evalpoly();
- if (verbosing) {
- printf("try A=");
- print(A);
- printf(", B=");
- print(B);
- printf(", root ");
- print(X);
- printf("=-B/A=");
- print(Z);
- printf(", POLY(");
- print(Z);
- printf(")=");
- print(Q);
- }
- if (iszero(Q)) {
- tos = h;
- return 1;
- }
- push(B);
- negate();
- B = pop();
- push(Z);
- negate();
- Z = pop();
- evalpoly();
- if (verbosing) {
- printf("try A=");
- print(A);
- printf(", B=");
- print(B);
- printf(", root ");
- print(X);
- printf("=-B/A=");
- print(Z);
- printf(", POLY(");
- print(Z);
- printf(")=");
- print(Q);
- }
- if (iszero(Q)) {
- tos = h;
- return 1;
- }
- }
- }
- tos = h;
- return 0;
- }
- //-----------------------------------------------------------------------------
- //
- // Divide a polynomial by Ax+B
- //
- // Input: polycoeff Dividend coefficients
- //
- // expo Degree of dividend
- //
- // A As above
- //
- // B As above
- //
- // Output: polycoeff Contains quotient coefficients
- //
- //-----------------------------------------------------------------------------
- static void
- yydivpoly(void)
- {
- int i;
- Q = zero;
- for (i = expo; i > 0; i--) {
- push(polycoeff[i]);
- polycoeff[i] = Q;
- push(A);
- divide();
- Q = pop();
- push(polycoeff[i - 1]);
- push(Q);
- push(B);
- multiply();
- subtract();
- polycoeff[i - 1] = pop();
- }
- polycoeff[0] = Q;
- }
- static void
- evalpoly(void)
- {
- int i;
- push(zero);
- for (i = expo; i >= 0; i--) {
- push(Z);
- multiply();
- push(polycoeff[i]);
- add();
- }
- Q = pop();
- }
|