GRFreire/mathomatic
1
0
Fork 0
mirror of https://github.com/mfillpot/mathomatic.git synced 2026-01-08 04:29:39 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity
master
mathomatic/poly.c

2237 lines
61 KiB
C
Raw Permalink Blame History

/*
* Mathomatic simplifying and general polynomial routines.
* Includes polynomial and smart division, polynomial factoring, etc.
* Globals tlhs[] and trhs[] are used and wiped out by most of these routines.
*
* The polynomial division and GCD routines here are not recursive,
* due to the global static expression storage areas.
* This limits the polynomial GCD routines to mostly univariate operation.
* This also does not allow their use during solving.
* So far, these limitations have been a good thing,
* making Mathomatic faster and more stable and reliable.
* Making polynomial GCD calculations partially recursive
* with one level of recursion would enable multivariate operation for most cases, I think.
*
* Mathomatic has proved it is practical and efficient to do
* generalized polynomial operations. By generalized, I mean that
* the coefficients of polynomials are not specially treated or stored in any way.
* They are not separated out from the main variable of the polynomial.
* Maximum simplification of all expressions is not possible unless everything is generalized.
*
* Copyright (C) 1987-2012 George Gesslein II.
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
The chief copyright holder can be contacted at gesslein@mathomatic.org, or
George Gesslein II, P.O. Box 224, Lansing, NY 14882-0224 USA.
*/
#include "includes.h"
#define REMAINDER_IS_ZERO() (n_trhs == 1 && trhs[0].kind == CONSTANT && trhs[0].token.constant == 0.0)
/*
* The following static expression storage areas are of non-standard size
* and must only be used for temporary storage.
* Most Mathomatic expression manipulation and simplification routines should not be used
* on non-standard or constant size expression storage areas.
* Standard size expression storage areas that may be
* manipulated or simplified are the equation spaces, tlhs[], trhs[], and tes[] only.
*/
token_type divisor[DIVISOR_SIZE]; /* static expression storage areas for polynomial and smart division */
int n_divisor; /* length of expression in divisor[] */
token_type quotient[DIVISOR_SIZE];
int n_quotient; /* length of expression in quotient[] */
token_type gcd_divisor[DIVISOR_SIZE]; /* static expression storage area for polynomial GCD routine */
int len_d; /* length of expression in gcd_divisor[] */
static int pf_recurse(token_type *equation, int *np, int loc, int level, int do_repeat);
static int pf_sub(token_type *equation, int *np, int loc, int len, int level, int do_repeat);
static int save_factors(token_type *equation, int *np, int loc1, int len, int level);
static int do_gcd(long *vp);
static int mod_recurse(token_type *equation, int *np, int loc, int level);
static int polydiv_recurse(token_type *equation, int *np, int loc, int level);
static int pdiv_recurse(token_type *equation, int *np, int loc, int level, int code);
static int poly_div_sub(token_type *d1, int len1, token_type *d2, int len2, long *vp);
static int find_highest_count(token_type *p1, int n1, token_type *p2, int n2, long *vp1);
/*
* Compare function for qsort(3).
*/
static int
vcmp(p1, p2)
sort_type *p1, *p2;
{
if (p2->count == p1->count) {
if (p1->v < p2->v)
return -1;
if (p1->v == p2->v)
return 0;
return 1;
}
return(p2->count - p1->count);
}
/*
* Return true if passed expression is strictly a single polynomial term in variable v.
* The general form of a polynomial term is c*(v^d)
* The coefficient (c) and exponent (d) may be any expression,
* as long as it doesn't contain the variable v.
*/
int
poly_in_v_sub(p1, n, v, allow_divides)
token_type *p1; /* expression pointer */
int n; /* expression length */
long v; /* Mathomatic variable */
int allow_divides; /* if true, allow division by variable */
{
int i, k;
int level, vlevel;
int count;
level = min_level(p1, n);
for (i = 0, count = 0; i < n; i += 2) {
if (p1[i].kind == VARIABLE && p1[i].token.variable == v) {
count++;
if (count > 1)
return false;
vlevel = p1[i].level;
if (vlevel == level || vlevel == (level + 1)) {
for (k = 1; k < n; k += 2) {
if (p1[k].level == level) {
switch (p1[k].token.operatr) {
case DIVIDE:
if (!allow_divides && k == (i - 1))
return false;
case TIMES:
continue;
case POWER:
if (k == (i + 1))
continue;
default:
return false;
}
}
}
if (vlevel == (level + 1)) {
if ((i + 1) < n && p1[i+1].level == vlevel
&& p1[i+1].token.operatr == POWER) {
continue;
}
} else {
continue;
}
}
return false;
}
}
return true;
}
/*
* Return true if passed expression is a polynomial in the given variable v.
* The coefficients and exponents may be anything, as long as they don't contain the variable v.
* The passed expression should be fully unfactored, for a proper determination.
*/
int
poly_in_v(p1, n, v, allow_divides)
token_type *p1; /* expression pointer */
int n; /* expression length */
long v; /* Mathomatic variable */
int allow_divides; /* allow variable to be right of a divide (negative exponents) as a polynomial term */
{
int i, j;
for (i = 1, j = 0;; i += 2) {
if (i >= n || (p1[i].level == 1
&& (p1[i].token.operatr == PLUS || p1[i].token.operatr == MINUS))) {
if (!poly_in_v_sub(&p1[j], i - j, v, allow_divides)) {
return false;
}
j = i + 1;
}
if (i >= n)
break;
}
return true;
}
/*
* Factor polynomials by calling pf_sub() for every additive sub-expression.
* Factors repeated factor polynomials (like (x+1)^5) if do_repeat is true.
* And always factors multivariate polynomials with symbolic factors (like (x+a)*(x+b)).
*
* Does not factor polynomials with different numeric factors (like (x+1)*(x+2)).
*
* Because of accumulated floating point round-off error,
* this routine does not succeed at factoring most large polynomials.
*
* Return true if equation side was modified (factored).
*/
int
poly_factor(equation, np, do_repeat)
token_type *equation; /* pointer to the beginning of equation side */
int *np; /* pointer to length of equation side */
int do_repeat; /* factor repeated factors flag */
{
return pf_recurse(equation, np, 0, 1, do_repeat);
}
static int
pf_recurse(equation, np, loc, level, do_repeat)
token_type *equation;
int *np, loc, level;
int do_repeat;
{
int modified = false;
int i;
int count = 0, level_count = 0;
for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
switch (equation[i].token.operatr) {
case PLUS:
case MINUS:
count++;
if (equation[i].level == level) {
level_count++;
}
}
}
if (level_count && count > 1) { /* so we don't factor expressions with only one additive operator */
/* try to factor the sub-expression */
modified = pf_sub(equation, np, loc, i - loc, level, do_repeat);
}
for (i = loc; i < *np && equation[i].level >= level;) {
if (equation[i].level > level) {
modified |= pf_recurse(equation, np, i, level + 1, do_repeat);
i++;
for (; i < *np && equation[i].level > level; i += 2)
;
continue;
}
i++;
}
return modified;
}
/*
* Polynomial factoring subroutine.
* It can't factor everything, but it usually can factor polynomials
* if it will make the expression size smaller.
*
* Return true if equation side was modified (factored).
*/
static int
pf_sub(equation, np, loc, len, level, do_repeat)
token_type *equation; /* equation side holding the possible polynomial to factor */
int *np, /* pointer to length of equation side */
loc, /* index of start of polynomial in equation side */
len; /* length of polynomial */
int level; /* level of additive operators in polynomial */
int do_repeat; /* factor repeated factors flag */
{
token_type *p1;
int modified = false, symbolic_modified = false;
int i, j, k;
long v = 0, v1, last_v;
int len_first = 0;
int loc1, loc2, len2 = 0;
int loct, lent;
int count;
jmp_buf save_save;
int div_flag = 3;
int vc, cnt;
sort_type va[MAX_VARS];
double d;
int old_partial;
debug_string(3, "Entering pf_sub().");
old_partial = partial_flag;
loc2 = loc1 = loc;
find_greatest_power(&equation[loc1], len, &v, &d, &j, &k, &div_flag);
if (v == 0)
return false;
blt(save_save, jmp_save, sizeof(jmp_save));
if ((i = setjmp(jmp_save)) != 0) { /* trap errors */
partial_flag = old_partial;
blt(jmp_save, save_save, sizeof(jmp_save));
if (i == 13) { /* critical error code */
longjmp(jmp_save, i);
}
return(modified || symbolic_modified);
}
/* First factor polynomials with repeated factors */
/* using poly_gcd(polynomial, v * differentiate(polynomial, v)) to discover the factors: */
for (count = 1; do_repeat; count++) {
blt(trhs, &equation[loc1], len * sizeof(token_type));
n_trhs = len;
partial_flag = false;
uf_simp(trhs, &n_trhs);
partial_flag = old_partial;
if (level1_plus_count(trhs, n_trhs) < 2) {
/* must be at least 2 level 1 additive operators to be factorable */
goto skip_factor;
}
/* create a variable list with counts of the number of times each variable occurs: */
last_v = 0;
for (vc = 0; vc < ARR_CNT(va);) {
cnt = 0;
v1 = -1;
for (i = 0; i < n_trhs; i += 2) {
if (trhs[i].kind == VARIABLE && trhs[i].token.variable > last_v) {
if (v1 == -1 || trhs[i].token.variable < v1) {
v1 = trhs[i].token.variable;
cnt = 1;
} else if (trhs[i].token.variable == v1) {
cnt++;
}
}
}
if (v1 == -1)
break;
last_v = v1;
va[vc].v = v1;
va[vc].count = cnt;
vc++;
}
side_debug(3, &equation[loc1], len);
side_debug(3, trhs, n_trhs);
/* Find a valid polynomial base variable "v": */
cnt = -1;
if (v) {
if (vc > 1 && !poly_in_v(trhs, n_trhs, v, true)) {
v = 0;
}
}
for (i = 0; i < vc; i++) {
if ((va[i].v & VAR_MASK) <= SIGN) {
continue;
}
if (v == 0) {
if (poly_in_v(trhs, n_trhs, va[i].v, true)) {
v = va[i].v;
}
}
if (cnt < 0 || va[i].count < cnt) {
cnt = va[i].count;
}
}
if (cnt <= 1) /* A polynomial requires 2 or more instances of every normal variable to be factorable. */
goto skip_factor;
if (v == 0) {
#if 1
goto skip_factor;
#else
/* no polynomial variables found, try differentiating with respect to each normal variable until one works */
for (i = 0; i < vc; i++) {
if ((va[i].v & VAR_MASK) <= SIGN) {
continue;
}
v = va[i].v;
blt(tlhs, trhs, n_trhs * sizeof(token_type));
n_tlhs = n_trhs;
if (differentiate(tlhs, &n_tlhs, v)) {
break;
}
v = 0;
}
if (v == 0) {
goto skip_factor;
}
#endif
} else {
blt(tlhs, trhs, n_trhs * sizeof(token_type));
n_tlhs = n_trhs;
if (!differentiate(tlhs, &n_tlhs, v)) {
break;
}
}
#if !SILENT
if (debug_level >= 3) {
list_var(v, 0);
fprintf(gfp, _("Differentiation successful using variable %s.\n"), var_str);
}
#endif
simp_loop(tlhs, &n_tlhs);
if ((n_tlhs + 2) > min(DIVISOR_SIZE, n_tokens))
break;
for (i = 0; i < n_tlhs; i++)
tlhs[i].level++;
tlhs[n_tlhs].kind = OPERATOR;
tlhs[n_tlhs].level = 1;
tlhs[n_tlhs].token.operatr = TIMES;
n_tlhs++;
tlhs[n_tlhs].kind = VARIABLE;
tlhs[n_tlhs].level = 1;
tlhs[n_tlhs].token.variable = v;
n_tlhs++;
uf_simp(tlhs, &n_tlhs);
if (poly_gcd(&equation[loc1], len, tlhs, n_tlhs, v) <= 0)
break;
if (level1_plus_count(tlhs, n_tlhs) == 0)
break;
if (!save_factors(equation, np, loc1, len, level))
break;
loc1 += n_tlhs + 1;
len = n_trhs;
switch (count) {
case 1:
debug_string(1, "Polynomial with repeated factor factored.");
len_first = n_tlhs;
loc2 = loc1;
break;
case 2:
len2 = n_tlhs;
break;
}
modified = true;
}
/* Now factor polynomials with symbolic factors by grouping: */
if (!modified) {
last_v = 0;
next_v:
p1 = &equation[loc1];
blt(trhs, p1, len * sizeof(token_type));
n_trhs = len;
uf_simp_no_repeat(trhs, &n_trhs);
if (level1_plus_count(trhs, n_trhs) < 2) {
/* must be at least 2 level 1 additive operators to be factorable */
goto skip_factor;
}
for (;;) {
v = -1;
for (i = 0; i < len; i += 2) {
if (p1[i].kind == VARIABLE && p1[i].token.variable > last_v) {
if (v == -1 || p1[i].token.variable < v) {
v = p1[i].token.variable;
}
}
}
if (v == -1) {
break;
}
last_v = v;
/* make sure there is more than one "v" raised to the highest power: */
if (find_greatest_power(trhs, n_trhs, &v, &d, &j, &k, &div_flag) <= 1) {
continue;
}
blt(tlhs, trhs, n_trhs * sizeof(token_type));
n_tlhs = n_trhs;
/* do the grouping: */
while (factor_plus(tlhs, &n_tlhs, v, 0.0)) {
simp_loop(tlhs, &n_tlhs);
}
/* extract the highest power group: */
if (find_greatest_power(tlhs, n_tlhs, &v, &d, &j, &k, &div_flag) != 1) {
continue;
}
if (j) {
blt(tlhs, &tlhs[j], k * sizeof(token_type));
}
n_tlhs = k;
#if !SILENT
if (debug_level >= 3) {
fprintf(gfp, _("Trying factor: "));
list_proc(tlhs, n_tlhs, false);
fprintf(gfp, "\n");
}
#endif
if (poly_gcd(&equation[loc1], len, tlhs, n_tlhs, 0L) <= 0)
goto next_v;
if (level1_plus_count(tlhs, n_tlhs) == 0)
goto next_v;
if (!symbolic_modified) {
debug_string(1, "Symbolic polynomial factored.");
} else {
debug_string(1, "Found another symbolic factor.");
}
if (!save_factors(equation, np, loc1, len, level))
break;
len = n_tlhs;
symbolic_modified = true;
last_v = 0;
goto next_v;
}
}
skip_factor:
blt(jmp_save, save_save, sizeof(jmp_save));
if (modified) {
/* Repeated factor was factored out. */
/* See if we can factor out more of the repeated factor. */
if (len2) {
loct = loc2;
lent = len2;
} else {
loct = loc;
lent = len_first;
}
if (poly_gcd(&equation[loc1], len, &equation[loct], lent, v) > 0) {
if (save_factors(equation, np, loc1, len, level)) {
loc1 += n_tlhs + 1;
len = n_trhs;
}
}
if (len2) {
loc1 = loc2;
len = len2;
}
if (poly_gcd(&equation[loc], len_first, &equation[loc1], len, 0L) > 0) {
save_factors(equation, np, loc, len_first, level);
}
}
if (modified || symbolic_modified) {
for (i = loc; i < *np && equation[i].level >= level; i++)
;
#if DEBUG
if ((i & 1) != 1) {
error_bug("Error in result of pf_sub().");
}
#endif
debug_string(1, "Resulting factors of pf_sub():");
side_debug(1, &equation[loc], i - loc);
}
return(modified || symbolic_modified);
}
static int
save_factors(equation, np, loc1, len, level)
token_type *equation;
int *np, loc1, len, level;
{
int i, j;
i = n_tlhs + 1 + n_trhs;
if (i > (len * 3))
goto rejected;
if ((*np + (i - len)) > n_tokens)
goto rejected;
blt(&equation[loc1+i], &equation[loc1+len], (*np - (loc1 + len)) * sizeof(token_type));
*np += i - len;
blt(&equation[loc1], tlhs, n_tlhs * sizeof(token_type));
i = loc1 + n_tlhs;
equation[i].level = 0;
equation[i].kind = OPERATOR;
equation[i].token.operatr = TIMES;
i++;
blt(&equation[i], trhs, n_trhs * sizeof(token_type));
i += n_trhs;
for (j = loc1; j < i; j++)
equation[j].level += level;
return true;
rejected:
debug_string(1, "Polynomial factor rejected because too large.");
return false;
}
/*
* Remove level 1 trivial factors and divides from tlhs[].
* Makes it a simpler factor that looks a lot nicer and is easier to work with.
*
* Return true if result is a level 1 additive expression.
* If this returns false, nothing is removed.
*/
int
remove_factors(void)
{
int i, j, k;
int plus_flag = false, divide_flag = false;
int op;
debug_string(3, "Entering remove_factors() with: ");
side_debug(3, tlhs, n_tlhs);
do {
simp_ssub(tlhs, &n_tlhs, 0L, 1.0, false, true, 4);
} while (uf_power(tlhs, &n_tlhs));
for (i = 1, j = 0, k = 0;; i += 2) {
if (i >= n_tlhs) {
if (plus_flag && !divide_flag) {
if (k > 0)
j--;
blt(&scratch[k], &tlhs[j], (i - j) * sizeof(token_type));
k += i - j;
}
if (k <= 0) {
debug_string(3, "Leaving remove_factors() with false return and no change.");
return false;
}
blt(tlhs, scratch, k * sizeof(token_type));
n_tlhs = k;
debug_string(3, "Leaving remove_factors() with success and: ");
side_debug(3, tlhs, n_tlhs);
return true;
}
op = tlhs[i].token.operatr;
switch (tlhs[i].level) {
case 1:
switch (op) {
case PLUS:
case MINUS:
plus_flag = true;
continue;
case TIMES:
case DIVIDE:
break;
default:
debug_string(3, "Leaving remove_factors() with false return and no change.");
return false;
}
if (plus_flag && !divide_flag) {
if (k > 0)
j--;
blt(&scratch[k], &tlhs[j], (i - j) * sizeof(token_type));
k += i - j;
}
plus_flag = false;
divide_flag = (op == DIVIDE);
j = i + 1;
break;
case 2:
switch (op) {
case PLUS:
case MINUS:
plus_flag = true;
}
break;
}
}
}
/*
* This is the Euclidean GCD algorithm applied to polynomials.
* It needs to be made multivariate by making it recursive.
*
* Return the number of iterations (divisions), if successful,
* with the polynomial GCD result in gcd_divisor[] and len_d.
* Else return 0 or negative on failure. If 0, might work if operands switched,
* If a negative number, switching operands won't work either.
*/
static int
do_gcd(vp)
long *vp; /* polynomial base variable pointer */
{
int i;
int count;
for (count = 1; count < 50; count++) {
switch (poly_div(trhs, n_trhs, gcd_divisor, len_d, vp)) {
case 0:
/* divide failed */
return(1 - count);
case 2:
/* Total success! Remainder is zero. */
debug_string(2, "Found raw polynomial GCD:");
side_debug(2, gcd_divisor, len_d);
return count;
}
/* Do the Euclidean shuffle: swap trhs[] (remainder) and gcd_divisor[] */
if (len_d > n_tokens || n_trhs > DIVISOR_SIZE)
return 0;
blt(scratch, trhs, n_trhs * sizeof(token_type));
blt(trhs, gcd_divisor, len_d * sizeof(token_type));
blt(gcd_divisor, scratch, n_trhs * sizeof(token_type));
i = n_trhs;
n_trhs = len_d;
len_d = i;
}
return 0;
}
/*
* Compute the simplified and normalized polynomial Greatest Common Divisor
* of the expressions in "larger" and "smaller".
* This polynomial GCD routine is used for polynomial factoring, etc.
*
* Return a positive integer if successful.
* Return the GCD in trhs[].
* Return larger/GCD in tlhs[].
* The results are unfactored and simplified.
*/
int
poly_gcd(larger, llen, smaller, slen, v)
token_type *larger; /* larger polynomial */
int llen; /* larger polynomial length */
token_type *smaller; /* smaller polynomial */
int slen; /* smaller polynomial length */
long v; /* polynomial base variable */
{
int count;
debug_string(3, "Entering poly_gcd():");
side_debug(3, larger, llen);
side_debug(3, smaller, slen);
if (llen > n_tokens || slen > min(ARR_CNT(gcd_divisor), n_tokens))
return 0;
if (trhs != larger) {
blt(trhs, larger, llen * sizeof(token_type));
}
n_trhs = llen;
if (tlhs != smaller) {
blt(tlhs, smaller, slen * sizeof(token_type));
}
n_tlhs = slen;
if (!remove_factors())
return 0;
if (n_tlhs > ARR_CNT(gcd_divisor))
return 0;
blt(gcd_divisor, tlhs, n_tlhs * sizeof(token_type));
len_d = n_tlhs;
count = do_gcd(&v);
if (count <= 0)
return 0;
if (count > 1) {
if (len_d > n_tokens)
return 0;
blt(tlhs, gcd_divisor, len_d * sizeof(token_type));
n_tlhs = len_d;
if (!remove_factors())
return 0;
if (n_tlhs > ARR_CNT(gcd_divisor))
return 0;
blt(gcd_divisor, tlhs, n_tlhs * sizeof(token_type));
len_d = n_tlhs;
if (poly_div(larger, llen, gcd_divisor, len_d, &v) != 2) {
debug_string(1, "Polynomial GCD found, but larger divide failed in poly_gcd().");
return 0;
}
}
if (len_d > n_tokens)
return 0;
blt(trhs, gcd_divisor, len_d * sizeof(token_type));
n_trhs = len_d;
uf_simp(tlhs, &n_tlhs);
uf_simp(trhs, &n_trhs);
debug_string(3, "poly_gcd() successful.");
return(count);
}
/*
* Compute the polynomial Greatest Common Divisor of the expressions in "larger" and "smaller".
* This polynomial GCD routine is used by the division simplifiers.
*
* Return a positive integer if successful.
* Return larger/GCD in tlhs[].
* Return smaller/GCD in trhs[].
*/
int
poly2_gcd(larger, llen, smaller, slen, v, require_additive)
token_type *larger; /* larger polynomial */
int llen; /* larger polynomial length */
token_type *smaller; /* smaller polynomial */
int slen; /* smaller polynomial length */
long v; /* polynomial base variable */
int require_additive; /* require the GCD to contain addition or subtraction */
{
int i;
int count;
#if 0
jmp_buf save_save;
#endif
if (require_additive) {
/* Require an additive operator in both polynomials to continue. */
count = 0;
for (i = 1; i < llen; i += 2) {
if (larger[i].token.operatr == PLUS || larger[i].token.operatr == MINUS) {
count++;
break;
}
}
if (count == 0)
return 0;
count = 0;
for (i = 1; i < slen; i += 2) {
if (smaller[i].token.operatr == PLUS || smaller[i].token.operatr == MINUS) {
count++;
}
}
if (count == 0 /* || count > 200 */)
return 0;
}
debug_string(3, "Entering poly2_gcd():");
side_debug(3, larger, llen);
side_debug(3, smaller, slen);
if (llen > n_tokens || slen > min(ARR_CNT(gcd_divisor), n_tokens))
return 0;
blt(trhs, larger, llen * sizeof(token_type));
n_trhs = llen;
blt(tlhs, smaller, slen * sizeof(token_type));
n_tlhs = slen;
#if 0
if (require_additive) {
/* Require a level 1 additive operator in the divisor. */
blt(save_save, jmp_save, sizeof(jmp_save));
if ((i = setjmp(jmp_save)) != 0) { /* trap errors */
blt(jmp_save, save_save, sizeof(jmp_save));
if (i == 13) { /* critical error code */
longjmp(jmp_save, i);
}
return 0;
}
uf_simp(tlhs, &n_tlhs);
blt(jmp_save, save_save, sizeof(jmp_save));
if (level1_plus_count(tlhs, n_tlhs) == 0)
return 0;
}
#endif
if (n_tlhs > ARR_CNT(gcd_divisor))
return 0;
blt(gcd_divisor, tlhs, n_tlhs * sizeof(token_type));
len_d = n_tlhs;
count = do_gcd(&v);
if (count <= 0)
return count;
if (count > 1) {
if (require_additive && level1_plus_count(gcd_divisor, len_d) == 0)
return 0;
if (poly_div(smaller, slen, gcd_divisor, len_d, &v) != 2) {
debug_string(1, "Polynomial GCD found, but smaller divide failed in poly2_gcd().");
return 0;
}
blt(trhs, gcd_divisor, len_d * sizeof(token_type));
n_trhs = len_d;
if (n_tlhs > ARR_CNT(gcd_divisor))
return 0;
blt(gcd_divisor, tlhs, n_tlhs * sizeof(token_type));
len_d = n_tlhs;
blt(tlhs, trhs, n_trhs * sizeof(token_type));
n_tlhs = n_trhs;
if (poly_div(larger, llen, tlhs, n_tlhs, &v) != 2) {
debug_string(1, "Polynomial GCD found, but larger divide failed in poly2_gcd().");
return 0;
}
blt(trhs, gcd_divisor, len_d * sizeof(token_type));
n_trhs = len_d;
} else {
n_trhs = 1;
trhs[0] = one_token;
}
debug_string(3, "poly2_gcd() successful.");
return count;
}
/*
* This function returns true if the passed Mathomatic variable
* is of type integer.
*
* Integer variable names start with "integer".
*/
int
is_integer_var(v)
long v;
{
char *cp;
int (*strncmpfunc)();
if (case_sensitive_flag) {
strncmpfunc = strncmp;
} else {
strncmpfunc = strncasecmp;
}
cp = var_name(v);
if (cp && strncmpfunc(cp, V_INTEGER_PREFIX, strlen(V_INTEGER_PREFIX)) == 0)
return true;
else
return false;
}
/*
* This function is a strict test that
* returns true if passed expression is entirely integer and all variables
* are either "integer" or "sign".
*
* The result of evaluating the expression must be integer if this returns true.
* Should first be unfactored with uf_pplus() for a proper determination.
*/
int
is_integer_expr(p1, n)
token_type *p1; /* expression pointer */
int n; /* length of expression */
{
int i;
long v;
#if DEBUG
if (p1 == NULL || n < 1) {
error_bug("(p1 == NULL || n < 1) in is_integer_expr().");
}
#endif
for (i = 0; i < n; i++) {
switch (p1[i].kind) {
case OPERATOR:
if (p1[i].token.operatr == DIVIDE)
return false;
break;
case CONSTANT:
if (fmod(p1[i].token.constant, 1.0) != 0.0)
return false;
break;
case VARIABLE:
v = labs(p1[i].token.variable);
if (!is_integer_var(v) && (v & VAR_MASK) != SIGN)
return false;
break;
}
}
return true;
}
/*
* This routine is the modulus operator (%) simplifier for equation sides.
* Using "integer" variables will allow more simplification here.
* Globals tlhs[] and trhs[] are wiped out by the polynomial division.
*
* Return true if equation side was modified.
*/
int
mod_simp(equation, np)
token_type *equation; /* pointer to the beginning of equation side to simplify */
int *np; /* pointer to length of the equation side */
{
return mod_recurse(equation, np, 0, 1);
}
static int
mod_recurse(equation, np, loc, level)
token_type *equation;
int *np, loc, level;
{
int modified = false;
int i, j, k;
int i1, i2, i3, i4, i5;
int op, last_op2;
int len1, len2, len3;
int diff_sign;
for (i = loc; i < *np && equation[i].level >= level;) {
if (equation[i].level > level) {
modified |= mod_recurse(equation, np, i, level + 1);
i++;
for (; i < *np && equation[i].level > level; i += 2)
;
continue;
}
i++;
}
if (modified) /* make sure the deepest levels are completed, first */
return true;
for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
if (!(equation[i].level == level && equation[i].token.operatr == MODULUS))
continue;
for (k = i + 2;; k += 2) {
if (k >= *np || equation[k].level <= level)
break;
}
len1 = k - (i + 1);
last_op2 = 0;
for (j = loc; j < *np && equation[j].level >= level; j++) {
if (equation[j].level == level && equation[j].kind == OPERATOR) {
last_op2 = equation[j].token.operatr;
continue;
}
if (last_op2 == MODULUS) {
continue;
}
last_op2 = MODULUS;
op = 0;
for (i1 = k = j + 1; k < *np && equation[k].level > level; k += 2) {
if (equation[k].level == (level + 1)) {
op = equation[k].token.operatr;
i1 = k;
}
}
len2 = k - j;
switch (op) {
case MODULUS:
/* simplify (x%n)%n to x%n */
len3 = k - (i1 + 1);
if (se_compare(&equation[i+1], len1, &equation[i1+1], len3, &diff_sign)) {
blt(&equation[i1], &equation[k], (*np - k) * sizeof(token_type));
*np -= len3 + 1;
return true;
}
break;
case TIMES:
if (!is_integer_expr(&equation[j], len2))
break;
/* simplify (i%n*j)%n to (i*j)%n if j is integer */
for (i2 = i1 = j + 1;; i1 += 2) {
if (i1 >= k || equation[i1].level == (level + 1)) {
for (; i2 < i1; i2 += 2) {
if (equation[i2].level == (level + 2)
&& equation[i2].token.operatr == MODULUS) {
len3 = i1 - (i2 + 1);
if (se_compare(&equation[i+1], len1, &equation[i2+1], len3, &diff_sign)) {
blt(&equation[i2], &equation[i1], (*np - i1) * sizeof(token_type));
*np -= len3 + 1;
return true;
}
}
}
}
if (i1 >= k)
break;
}
break;
case PLUS:
case MINUS:
/* simplify (i%n+j)%n to (i+j)%n */
/* and ((i%n)*j+k)%n to (i*j+k)%n, when j is integer */
for (i2 = i1 = j + 1, i3 = j - 1;; i1 += 2) {
if (i1 >= k || equation[i1].level == (level + 1)) {
for (; i2 < i1; i2 += 2) {
if (equation[i2].level == (level + 2)) {
switch (equation[i2].token.operatr) {
case MODULUS:
len3 = i1 - (i2 + 1);
if (se_compare(&equation[i+1], len1, &equation[i2+1], len3, &diff_sign)) {
blt(&equation[i2], &equation[i1], (*np - i1) * sizeof(token_type));
*np -= len3 + 1;
return true;
}
break;
case TIMES:
i2 = i1 - 2;
if (!is_integer_expr(&equation[i3+1], i1 - (i3 + 1))) {
break;
}
for (i4 = i3 + 2; i4 < i1; i4 += 2) {
if (equation[i4].level == (level + 3)
&& equation[i4].token.operatr == MODULUS) {
for (i5 = i4 + 2; i5 < i1 && equation[i5].level > (level + 3); i5 += 2)
;
len3 = i5 - (i4 + 1);
if (se_compare(&equation[i+1], len1, &equation[i4+1], len3, &diff_sign)) {
blt(&equation[i4], &equation[i5], (*np - i5) * sizeof(token_type));
*np -= len3 + 1;
return true;
}
}
}
break;
}
}
}
i3 = i1;
}
if (i1 >= k)
break;
}
break;
}
#if 1
/* Remove integer*n multiples in x for x%n by doing */
/* polynomial division x/n and replacing with remainder%n. */
/* Globals tlhs[] and trhs[] are wiped out by the polynomial division here. */
if (poly_div(&equation[j], len2, &equation[i+1], len1, NULL)) {
uf_pplus(tlhs, &n_tlhs); /* so integer%(integer^integer) isn't simplified to 0 */
if (is_integer_expr(tlhs, n_tlhs)) {
if (n_trhs < len2 || REMAINDER_IS_ZERO()) { /* if it will make it smaller */
if ((*np + (n_trhs - len2)) > n_tokens)
error_huge();
for (k = 0; k < n_trhs; k++)
trhs[k].level += level;
blt(&equation[j+n_trhs], &equation[j+len2], (*np - (j + len2)) * sizeof(token_type));
*np += n_trhs - len2;
blt(&equation[j], trhs, n_trhs * sizeof(token_type));
debug_string(2, "Polynomial division successful in modulus simplification. The result is:");
side_debug(2, equation, *np);
return true;
}
}
}
#endif
}
}
return modified;
}
/*
* This routine is a division simplifier for equation sides.
* It reduces algebraic fractions.
*
* Return true if any polynomial Greatest Common Divisor (GCD)
* between a numerator and denominator was found and
* the expression was simplified by dividing the numerator
* and denominator by the GCD.
*/
int
poly_gcd_simp(equation, np)
token_type *equation;
int *np;
{
return polydiv_recurse(equation, np, 0, 1);
}
static int
polydiv_recurse(equation, np, loc, level)
token_type *equation;
int *np, loc, level;
{
int modified = false;
int i, j, k;
int last_op2;
int len1, len2;
int rv;
for (i = loc; i < *np && equation[i].level >= level;) {
if (equation[i].level > level) {
modified |= polydiv_recurse(equation, np, i, level + 1);
i++;
for (; i < *np && equation[i].level > level; i += 2)
;
continue;
}
i++;
}
start:
for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
#if DEBUG
if (equation[i].kind != OPERATOR)
error_bug("Bug in poly_gcd_simp().");
#endif
if (equation[i].level == level && equation[i].token.operatr == DIVIDE) {
for (k = i + 2;; k += 2) {
if (k >= *np || equation[k].level <= level)
break;
}
len1 = k - (i + 1);
last_op2 = 0;
for (j = loc; j < *np && equation[j].level >= level; j++) {
if (equation[j].level == level && equation[j].kind == OPERATOR) {
last_op2 = equation[j].token.operatr;
continue;
}
switch (last_op2) {
case DIVIDE:
continue;
case 0:
case TIMES:
break;
default:
error_bug("Expression is corrupt in poly_gcd_simp().");
}
last_op2 = DIVIDE;
for (k = j + 1;; k += 2) {
if (k >= *np || equation[k].level <= level)
break;
}
len2 = k - j;
if ((rv = poly2_gcd(&equation[i+1], len1, &equation[j], len2, 0L, true)) > 0) {
store_code:
for (k = 0; k < n_tlhs; k++)
tlhs[k].level += level;
for (k = 0; k < n_trhs; k++)
trhs[k].level += level;
if (((*np + (n_trhs - len2)) > n_tokens)
|| ((*np + (n_trhs - len2) + (n_tlhs - len1)) > n_tokens))
error_huge();
blt(&equation[j+n_trhs], &equation[j+len2], (*np - (j + len2)) * sizeof(token_type));
*np += n_trhs - len2;
if (i > j)
i += n_trhs - len2;
blt(&equation[j], trhs, n_trhs * sizeof(token_type));
blt(&equation[i+n_tlhs+1], &equation[i+1+len1], (*np - (i + 1 + len1)) * sizeof(token_type));
*np += n_tlhs - len1;
blt(&equation[i+1], tlhs, n_tlhs * sizeof(token_type));
debug_string(1, _("Division simplified with polynomial GCD."));
modified = true;
goto start;
}
if (rv == 0 && poly2_gcd(&equation[j], len2, &equation[i+1], len1, 0L, true) > 0) {
k = j - 1;
j = i + 1;
i = k;
k = len1;
len1 = len2;
len2 = k;
goto store_code;
}
}
}
}
return modified;
}
/*
* This routine is a division simplifier for equation sides.
* Check for divides and do polynomial and smart division.
* For example: (c - b + a)/(b - a) will be converted to (c/(b - a)) - 1.
* Will only make expressions with algebraic fractions smaller,
* otherwise it would go into an infinite loop.
* An expression with less operators is judged to be smaller.
* Complex fractions (fractions within fractions) are sometimes
* created by this, especially if quick_flag is false.
*
* Return true if expression was simplified.
*/
int
div_remainder(equation, np, poly_flag, quick_flag)
token_type *equation;
int *np;
int poly_flag; /* if true, try polynomial division first, then smart division */
int quick_flag; /* if true, keep algebraic fractions simpler */
{
int rv = false;
debug_string(3, "Entering div_remainder().");
if (quick_flag)
group_proc(equation, np);
rv = pdiv_recurse(equation, np, 0, 1, poly_flag);
if (quick_flag)
organize(equation, np);
debug_string(3, "Leaving div_remainder().");
return rv;
}
static int
pdiv_recurse(equation, np, loc, level, code)
token_type *equation;
int *np, loc, level, code;
{
int modified = false;
int i, j, k;
int op, op2, last_op2;
int len1, len2, real_len1;
int rv = 0;
int flag, power_flag, zero_remainder;
for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
if (!(equation[i].level == level && equation[i].token.operatr == DIVIDE))
continue;
for (k = i + 2;; k += 2) {
if (k >= *np || equation[k].level <= level)
break;
}
len1 = real_len1 = k - (i + 1);
last_op2 = 0;
for (j = loc; j < *np && equation[j].level >= level; j++) {
if (equation[j].level == level && equation[j].kind == OPERATOR) {
last_op2 = equation[j].token.operatr;
continue;
}
if (last_op2 == DIVIDE) {
continue;
}
last_op2 = DIVIDE;
op = 0;
for (k = j + 1; k < *np && equation[k].level > level; k += 2) {
if (equation[k].level == (level + 1)) {
op = equation[k].token.operatr;
}
}
if (op != PLUS && op != MINUS) {
continue;
}
len2 = k - j;
flag = code;
power_flag = false;
op = 0;
op2 = 0;
for (k = i + 2; k < *np && equation[k].level > level; k += 2) {
if (equation[k].level == level + 3) {
switch (equation[k].token.operatr) {
case PLUS:
case MINUS:
op2 = PLUS;
}
} else if (equation[k].level == level + 2) {
op = equation[k].token.operatr;
} else if (equation[k].level == level + 1) {
if (equation[k].token.operatr == POWER
&& (op == PLUS || op == MINUS || (op == TIMES && op2 == PLUS))) {
power_flag = true;
len1 = k - (i + 1);
}
break;
}
}
next_thingy:
if (!power_flag) {
len1 = real_len1;
}
if (flag || power_flag) {
rv = poly_div(&equation[j], len2, &equation[i+1], len1, NULL);
} else {
rv = smart_div(&equation[j], len2, &equation[i+1], len1);
}
zero_remainder = (rv > 0 && REMAINDER_IS_ZERO());
if (power_flag && !zero_remainder) {
rv = 0;
}
if (rv > 0) { /* if successful and the result is smaller than the original: */
if ((n_tlhs + 2 + n_trhs + len1) > n_tokens)
error_huge();
for (k = 0; k < n_tlhs; k++)
tlhs[k].level++;
tlhs[n_tlhs].level = 1;
tlhs[n_tlhs].kind = OPERATOR;
tlhs[n_tlhs].token.operatr = PLUS;
n_tlhs++;
for (k = 0; k < n_trhs; k++)
trhs[k].level += 2;
blt(&tlhs[n_tlhs], trhs, n_trhs * sizeof(token_type));
n_tlhs += n_trhs;
tlhs[n_tlhs].level = 2;
tlhs[n_tlhs].kind = OPERATOR;
tlhs[n_tlhs].token.operatr = DIVIDE;
n_tlhs++;
k = n_tlhs;
blt(&tlhs[n_tlhs], &equation[i+1], len1 * sizeof(token_type));
n_tlhs += len1;
for (; k < n_tlhs; k++)
tlhs[k].level += 2;
side_debug(3, &equation[j], len2);
side_debug(3, &equation[i+1], len1);
simpb_side(tlhs, &n_tlhs, false, true, 3); /* parameters should match what's used in simpa_side() */
side_debug(3, tlhs, n_tlhs);
if (power_flag) {
k = (var_count(tlhs, n_tlhs) <= var_count(&equation[j], len2));
} else {
k = ((var_count(tlhs, n_tlhs) + (n_tlhs >= len1 + 1 + len2)) <= (var_count(&equation[j], len2) + var_count(&equation[i+1], len1)));
}
if (k) {
for (k = 0; k < n_tlhs; k++)
tlhs[k].level += level;
if (power_flag) {
if ((*np - len2 + n_tlhs + 2) > n_tokens)
error_huge();
for (k = i + 2 + len1; k <= i + real_len1; k++) {
equation[k].level++;
}
blt(&equation[i+real_len1+3], &equation[k], (*np - k) * sizeof(token_type));
*np += 2;
equation[k].level = level + 2;
equation[k].kind = OPERATOR;
equation[k].token.operatr = MINUS;
k++;
equation[k].level = level + 2;
equation[k].kind = CONSTANT;
equation[k].token.constant = 1.0;
if (i < j) {
j += 2;
}
} else {
if ((*np - (len1 + 1 + len2) + n_tlhs) > n_tokens)
error_huge();
blt(&equation[i], &equation[i+1+len1], (*np - (i + 1 + len1)) * sizeof(token_type));
*np -= len1 + 1;
if (i < j) {
j -= len1 + 1;
}
}
blt(&equation[j+n_tlhs], &equation[j+len2], (*np - (j + len2)) * sizeof(token_type));
*np -= len2 - n_tlhs;
blt(&equation[j], tlhs, n_tlhs * sizeof(token_type));
if (flag || power_flag) {
debug_string(1, _("Polynomial division successful."));
} else {
debug_string(1, _("Smart division successful."));
}
side_debug(3, equation, *np);
return true;
}
}
if (power_flag) {
power_flag = false;
goto next_thingy;
}
if (flag == code) {
flag = !flag;
goto next_thingy;
}
}
}
for (i = loc; i < *np && equation[i].level >= level;) {
if (equation[i].level > level) {
modified |= pdiv_recurse(equation, np, i, level + 1, code);
i++;
for (; i < *np && equation[i].level > level; i += 2)
;
continue;
}
i++;
}
return modified;
}
/*
* Do generalized and accurate polynomial division.
* This is the main polynomial division function for Mathomatic.
* Works with any number of variables and any pair of expressions.
*
* Returns non-zero if successful: 2 if remainder is 0,
* 1 if result is smaller than the original two expressions,
* negative if result is larger.
* Quotient is returned in tlhs[] and remainder in trhs[].
*
* If *vp is 0, automatically select the best polynomial base variable and return it in *vp.
*/
int
poly_div(d1, len1, d2, len2, vp)
token_type *d1; /* pointer to dividend */
int len1; /* length of dividend */
token_type *d2; /* pointer to divisor */
int len2; /* length of divisor */
long *vp; /* pointer to polynomial base variable */
{
int i;
int rv;
int old_partial;
jmp_buf save_save;
old_partial = partial_flag;
partial_flag = false; /* We want full unfactoring during polynomial division. */
blt(save_save, jmp_save, sizeof(jmp_save));
if ((i = setjmp(jmp_save)) != 0) { /* Trap errors so we almost always return normally. */
blt(jmp_save, save_save, sizeof(jmp_save));
partial_flag = old_partial;
if (i == 13) { /* critical error code */
longjmp(jmp_save, i);
}
return false;
}
rv = poly_div_sub(d1, len1, d2, len2, vp);
blt(jmp_save, save_save, sizeof(jmp_save));
partial_flag = old_partial;
return rv;
}
#define MAX_GREATEST_POWER_TERMS 50 /* a limit to keep polynomial division from taking too long */
/*
* Do the actual polynomial division using a simple,
* generalized, polynomial long division algorithm.
*/
static int
poly_div_sub(d1, len1, d2, len2, vp)
token_type *d1;
int len1;
token_type *d2;
int len2;
long *vp;
{
int i;
int t1, len_t1;
int t2, len_t2;
int sign;
int divide_flag;
double last_power, divisor_power, d;
int sum_size;
int last_count, count; /* dividend number of terms raised to the highest power */
int divisor_count; /* divisor number of terms raised to the highest power */
long tmp_v = 0;
if (vp == NULL)
vp = &tmp_v;
if (len1 > n_tokens || len2 > n_tokens)
return false;
/* Copy the source polynomials to where we want them (tlhs and trhs). */
if (trhs != d1) {
blt(trhs, d1, len1 * sizeof(token_type));
}
n_trhs = len1;
if (tlhs != d2) {
blt(tlhs, d2, len2 * sizeof(token_type));
}
n_tlhs = len2;
/* Do the basic unfactoring and simplification of the dividend and divisor. */
uf_simp(trhs, &n_trhs);
uf_simp(tlhs, &n_tlhs);
if (*vp == 0) {
/* Select the best polynomial base variable. */
if (!find_highest_count(trhs, n_trhs, tlhs, n_tlhs, vp))
return false;
}
#if !SILENT
/* Display debugging info. */
if (debug_level >= 3) {
list_var(*vp, 0);
fprintf(gfp, _("poly_div() starts using base variable %s:\n"), var_str);
side_debug(3, trhs, n_trhs);
side_debug(3, tlhs, n_tlhs);
}
#endif
/* Determine divide_flag and if the polynomials can be divided. */
divide_flag = 2;
last_count = find_greatest_power(trhs, n_trhs, vp, &last_power, &t1, &len_t1, &divide_flag);
divisor_count = find_greatest_power(tlhs, n_tlhs, vp, &divisor_power, &t2, &len_t2, &divide_flag);
if (divisor_power <= 0 || last_power < divisor_power) {
divide_flag = !divide_flag;
last_count = find_greatest_power(trhs, n_trhs, vp, &last_power, &t1, &len_t1, &divide_flag);
divisor_count = find_greatest_power(tlhs, n_tlhs, vp, &divisor_power, &t2, &len_t2, &divide_flag);
if (divisor_power <= 0 || last_power < divisor_power) {
return false;
}
}
if (divisor_count > 1 || last_count > MAX_GREATEST_POWER_TERMS) /* Need recursion and factoring to do more than this. */
return false;
/* Initialize the quotient. */
n_quotient = 1;
quotient[0] = zero_token;
/* Store the divisor. */
if (n_tlhs > ARR_CNT(divisor))
return false;
blt(divisor, tlhs, n_tlhs * sizeof(token_type));
n_divisor = n_tlhs;
sum_size = n_trhs + n_quotient;
/* Loop until polynomial division is finished. */
for (;;) {
if (t1 > 0 && trhs[t1-1].token.operatr == MINUS)
sign = MINUS;
else
sign = PLUS;
if (t2 > 0 && divisor[t2-1].token.operatr == MINUS) {
if (sign == MINUS)
sign = PLUS;
else
sign = MINUS;
}
if ((len_t1 + len_t2 + 1) > n_tokens)
return false;
blt(tlhs, &trhs[t1], len_t1 * sizeof(token_type));
n_tlhs = len_t1;
for (i = 0; i < n_tlhs; i++)
tlhs[i].level++;
tlhs[n_tlhs].level = 1;
tlhs[n_tlhs].kind = OPERATOR;
tlhs[n_tlhs].token.operatr = DIVIDE;
n_tlhs++;
blt(&tlhs[n_tlhs], &divisor[t2], len_t2 * sizeof(token_type));
i = n_tlhs;
n_tlhs += len_t2;
for (; i < n_tlhs; i++)
tlhs[i].level++;
if (!simp_loop(tlhs, &n_tlhs))
return false;
if ((n_quotient + 1 + n_tlhs) > min(ARR_CNT(quotient), n_tokens))
return false;
for (i = 0; i < n_tlhs; i++)
tlhs[i].level++;
quotient[n_quotient].level = 1;
quotient[n_quotient].kind = OPERATOR;
quotient[n_quotient].token.operatr = sign;
n_quotient++;
blt(&quotient[n_quotient], tlhs, n_tlhs * sizeof(token_type));
n_quotient += n_tlhs;
if ((n_trhs + n_tlhs + n_divisor + 2) > n_tokens)
return false;
blt(&trhs[t1+1], &trhs[t1+len_t1], (n_trhs - (t1 + len_t1)) * sizeof(token_type));
n_trhs -= len_t1 - 1;
trhs[t1] = zero_token;
for (i = 0; i < n_trhs; i++)
trhs[i].level++;
trhs[n_trhs].level = 1;
trhs[n_trhs].kind = OPERATOR;
if (sign == PLUS)
trhs[n_trhs].token.operatr = MINUS;
else
trhs[n_trhs].token.operatr = PLUS;
n_trhs++;
blt(&trhs[n_trhs], tlhs, n_tlhs * sizeof(token_type));
i = n_trhs;
n_trhs += n_tlhs;
for (; i < n_trhs; i++)
trhs[i].level++;
trhs[n_trhs].level = 2;
trhs[n_trhs].kind = OPERATOR;
trhs[n_trhs].token.operatr = TIMES;
n_trhs++;
i = n_trhs;
blt(&trhs[n_trhs], divisor, t2 * sizeof(token_type));
n_trhs += t2;
trhs[n_trhs] = zero_token;
n_trhs++;
blt(&trhs[n_trhs], &divisor[t2+len_t2], (n_divisor - (t2 + len_t2)) * sizeof(token_type));
n_trhs += (n_divisor - (t2 + len_t2));
for (; i < n_trhs; i++)
trhs[i].level += 2;
side_debug(3, trhs, n_trhs);
uf_repeat(trhs, &n_trhs);
uf_tsimp(trhs, &n_trhs);
side_debug(4, trhs, n_trhs);
count = find_greatest_power(trhs, n_trhs, vp, &d, &t1, &len_t1, &divide_flag);
if (d < divisor_power) {
/* Success! Polynomial division ends here. */
debug_string(3, "Successful polynomial division!");
blt(tlhs, quotient, n_quotient * sizeof(token_type));
n_tlhs = n_quotient;
debug_string(3, "Quotient:");
side_debug(3, tlhs, n_tlhs);
debug_string(3, "Remainder:");
side_debug(3, trhs, n_trhs);
if (REMAINDER_IS_ZERO())
return 2;
if ((n_trhs + n_quotient) >= (sum_size /* - (sum_size / 10) */)) {
if ((n_trhs + 1) > sum_size && n_trhs > n_divisor)
return -2;
else
return -1;
}
return 1;
} else if (d < last_power) {
last_power = d;
last_count = count;
if (last_count > MAX_GREATEST_POWER_TERMS) /* Need factoring to do more than this. */
return false;
} else if (d > last_power) {
return false;
} else {
if (count >= last_count) {
return false;
}
last_count = count;
}
}
}
/*
* Do smart division.
*
* Smart division is heuristic division much like polynomial division,
* however instead of basing the division on the highest powers of a base variable,
* every term in the dividend is tried, and if a trial makes the
* expression smaller, we go with that.
*
* Only one term in the divisor is tried, to save time.
* Currently, the divisor term used is the one with the least number
* of variables, not including any terms with no variables.
*
* Returns true if successful.
* Quotient is returned in tlhs[] and remainder in trhs[].
*/
int
smart_div(d1, len1, d2, len2)
token_type *d1; /* pointer to dividend */
int len1; /* length of dividend */
token_type *d2; /* pointer to divisor */
int len2; /* length of divisor */
{
int i, j, k;
int t1, len_t1;
int t2 = 0, len_t2 = 0;
int sign;
int old_n_quotient;
int trhs_size;
int term_size = 0, term_count;
int term_pos = 0, skip_terms[100];
int skip_count;
token_type *qp;
int q_size;
int sum_size;
int count;
int dcount = 0; /* divisor term count */
int flag;
blt(trhs, d1, len1 * sizeof(token_type));
n_trhs = len1;
blt(tlhs, d2, len2 * sizeof(token_type));
n_tlhs = len2;
/* Do the basic unfactoring and simplification of the dividend and divisor. */
uf_simp_no_repeat(trhs, &n_trhs);
uf_simp_no_repeat(tlhs, &n_tlhs);
/* Display debugging info. */
debug_string(3, "smart_div() starts:");
side_debug(3, trhs, n_trhs);
side_debug(3, tlhs, n_tlhs);
/* Find which divisor term to use. */
for (i = 0, j = 0, k = 0, flag = false;; i++) {
if (i >= n_tlhs || (tlhs[i].kind == OPERATOR && tlhs[i].level == 1
&& (tlhs[i].token.operatr == PLUS || tlhs[i].token.operatr == MINUS))) {
dcount++;
if (flag) {
if (len_t2 == 0 || var_count(&tlhs[j], i - j) < k) {
len_t2 = i - j;
t2 = j;
k = var_count(&tlhs[t2], len_t2);
}
}
flag = false;
j = i + 1;
} else if (tlhs[i].kind == VARIABLE && tlhs[i].token.variable != IMAGINARY) {
flag = true;
}
if (i >= n_tlhs)
break;
}
if (len_t2 <= 0)
return false;
/* Initialize the quotient. */
n_quotient = 1;
quotient[0] = zero_token;
if (n_tlhs > ARR_CNT(divisor))
return false;
blt(divisor, tlhs, n_tlhs * sizeof(token_type));
n_divisor = n_tlhs;
try_one:
trhs_size = n_trhs;
for (skip_count = 0, count = 0;;) {
sum_size = n_trhs + n_quotient;
for (term_count = 1, q_size = 0;; term_count++) {
if (!get_term(trhs, n_trhs, term_count, &t1, &len_t1))
break;
flag = false;
for (i = 0; i < skip_count; i++) {
if (skip_terms[i] == t1) {
flag = true;
break;
}
}
if (flag)
continue;
if ((len_t1 + len_t2 + 1) > n_tokens)
return false;
blt(tlhs, &trhs[t1], len_t1 * sizeof(token_type));
n_tlhs = len_t1;
for (i = 0; i < n_tlhs; i++)
tlhs[i].level++;
tlhs[n_tlhs].level = 1;
tlhs[n_tlhs].kind = OPERATOR;
tlhs[n_tlhs].token.operatr = DIVIDE;
n_tlhs++;
blt(&tlhs[n_tlhs], &divisor[t2], len_t2 * sizeof(token_type));
i = n_tlhs;
n_tlhs += len_t2;
for (; i < n_tlhs; i++)
tlhs[i].level++;
if (!simp_loop(tlhs, &n_tlhs))
continue;
if (basic_size(tlhs, n_tlhs) <= basic_size(&trhs[t1], len_t1)) {
q_size = n_tlhs;
term_pos = t1;
term_size = len_t1;
break;
}
}
if (q_size <= 0) {
if (count <= 0) {
if (dcount > 1) {
dcount = 1;
t2 = 0;
len_t2 = n_divisor;
goto try_one; /* Try using the whole divisor as the divisor term. */
}
return false;
}
end_div:
if (dcount > 1) {
if (n_quotient + n_trhs >= trhs_size + 1) {
return false;
}
}
end_div2:
blt(tlhs, quotient, n_quotient * sizeof(token_type));
n_tlhs = n_quotient;
side_debug(3, tlhs, n_tlhs);
side_debug(3, trhs, n_trhs);
return true; /* Success! */
}
t1 = term_pos;
len_t1 = term_size;
if (t1 > 0 && trhs[t1-1].token.operatr == MINUS)
sign = MINUS;
else
sign = PLUS;
if (t2 > 0 && divisor[t2-1].token.operatr == MINUS) {
if (sign == MINUS)
sign = PLUS;
else
sign = MINUS;
}
if ((len_t1 + len_t2 + 1) > n_tokens)
return false;
blt(tlhs, &trhs[t1], len_t1 * sizeof(token_type));
n_tlhs = len_t1;
for (i = 0; i < n_tlhs; i++)
tlhs[i].level++;
tlhs[n_tlhs].level = 1;
tlhs[n_tlhs].kind = OPERATOR;
tlhs[n_tlhs].token.operatr = DIVIDE;
n_tlhs++;
blt(&tlhs[n_tlhs], &divisor[t2], len_t2 * sizeof(token_type));
i = n_tlhs;
n_tlhs += len_t2;
for (; i < n_tlhs; i++)
tlhs[i].level++;
simp_loop(tlhs, &n_tlhs);
if ((n_quotient + 1 + n_tlhs) > min(ARR_CNT(quotient), n_tokens))
return false;
for (i = 0; i < n_tlhs; i++)
tlhs[i].level++;
old_n_quotient = n_quotient;
quotient[n_quotient].level = 1;
quotient[n_quotient].kind = OPERATOR;
quotient[n_quotient].token.operatr = sign;
n_quotient++;
qp = &quotient[n_quotient];
q_size = n_tlhs;
blt(&quotient[n_quotient], tlhs, n_tlhs * sizeof(token_type));
n_quotient += n_tlhs;
if ((n_trhs + q_size + n_divisor + 2) > n_tokens)
return false;
blt(tlhs, trhs, n_trhs * sizeof(token_type));
n_tlhs = n_trhs;
blt(&trhs[t1+1], &trhs[t1+len_t1], (n_trhs - (t1 + len_t1)) * sizeof(token_type));
n_trhs -= len_t1 - 1;
trhs[t1] = zero_token;
for (i = 0; i < n_trhs; i++)
trhs[i].level++;
trhs[n_trhs].level = 1;
trhs[n_trhs].kind = OPERATOR;
if (sign == PLUS)
trhs[n_trhs].token.operatr = MINUS;
else
trhs[n_trhs].token.operatr = PLUS;
n_trhs++;
blt(&trhs[n_trhs], qp, q_size * sizeof(token_type));
i = n_trhs;
n_trhs += q_size;
for (; i < n_trhs; i++)
trhs[i].level++;
trhs[n_trhs].level = 2;
trhs[n_trhs].kind = OPERATOR;
trhs[n_trhs].token.operatr = TIMES;
n_trhs++;
i = n_trhs;
blt(&trhs[n_trhs], divisor, t2 * sizeof(token_type));
n_trhs += t2;
trhs[n_trhs] = zero_token;
n_trhs++;
blt(&trhs[n_trhs], &divisor[t2+len_t2], (n_divisor - (t2 + len_t2)) * sizeof(token_type));
n_trhs += (n_divisor - (t2 + len_t2));
for (; i < n_trhs; i++)
trhs[i].level += 2;
side_debug(3, trhs, n_trhs);
uf_tsimp(trhs, &n_trhs);
side_debug(4, trhs, n_trhs);
if (REMAINDER_IS_ZERO())
goto end_div2;
if (dcount > 1) {
if ((n_trhs + n_quotient) >= sum_size) {
if (skip_count >= ARR_CNT(skip_terms)) {
if (count == 0) {
return false;
} else {
n_quotient = old_n_quotient;
blt(trhs, tlhs, n_tlhs * sizeof(token_type));
n_trhs = n_tlhs;
goto end_div;
}
}
skip_terms[skip_count] = term_pos;
skip_count++;
n_quotient = old_n_quotient;
blt(trhs, tlhs, n_tlhs * sizeof(token_type));
n_trhs = n_tlhs;
debug_string(3, "Skipping last operation.");
continue;
}
}
if (n_trhs == 1 && trhs[0].kind == CONSTANT)
goto end_div;
skip_count = 0;
count++;
}
}
/*
* Return the size of a sub-expression,
* minus any constant multiplier.
*/
int
basic_size(p1, len)
token_type *p1;
int len;
{
int i, j;
int level;
int rv;
int constant_flag = true;
rv = len;
level = min_level(p1, len);
for (i = 0, j = -1; i < len; i++) {
if (p1[i].kind == OPERATOR) {
if (p1[i].level == level
&& (p1[i].token.operatr == TIMES || p1[i].token.operatr == DIVIDE)) {
if (constant_flag) {
rv -= i - j;
}
j = i;
constant_flag = true;
}
} else if (p1[i].kind != CONSTANT)
constant_flag = false;
}
if (constant_flag) {
rv -= i - j;
}
return rv;
}
int
get_term(p1, n1, count, tp1, lentp1)
token_type *p1;
int n1;
int count;
int *tp1, *lentp1;
{
int i, j;
int no;
for (no = 0, i = 1, j = 0;; i += 2) {
if (i >= n1 || (p1[i].level == 1
&& (p1[i].token.operatr == PLUS
|| p1[i].token.operatr == MINUS))) {
no++;
if (no >= count) {
*tp1 = j;
*lentp1 = i - j;
return true;
}
j = i + 1;
}
if (i >= n1)
return false;
}
}
/*
* This routine automatically finds the best variable to do polynomial division with.
*
* Returns 0 if nothing was found,
* otherwise the selected polynomial base variable is returned in *vp1,
* and the number of times it occurs in the dividend is returned.
*/
static int
find_highest_count(p1, n1, p2, n2, vp1)
token_type *p1; /* pointer to dividend expression */
int n1; /* length of dividend */
token_type *p2; /* pointer to divisor expression */
int n2; /* length of divisor */
long *vp1; /* variable pointer to return base variable with */
{
int i;
int vc, cnt;
long v1, last_v, cv; /* Mathomatic variables */
sort_type va[MAX_VARS];
int divide_flag;
int t1, len_t1;
int t2, len_t2;
double d1, d2;
int count1, count2;
last_v = 0;
for (vc = 0; vc < ARR_CNT(va);) {
cnt = 0;
v1 = -1;
for (i = 0; i < n1; i += 2) {
if (p1[i].kind == VARIABLE && p1[i].token.variable > last_v) {
if (v1 == -1 || p1[i].token.variable < v1) {
v1 = p1[i].token.variable;
cnt = 1;
} else if (p1[i].token.variable == v1) {
cnt++;
}
}
}
if (v1 == -1)
break;
last_v = v1;
va[vc].v = v1;
va[vc].count = cnt;
vc++;
}
if (vc <= 0)
return 0;
qsort((char *) va, vc, sizeof(*va), vcmp);
for (cv = IMAGINARY + 1; cv > 0; cv--) {
for (i = 0; i < vc; i++) {
if ((cv > IMAGINARY) ? ((va[i].v & VAR_MASK) <= SIGN) : (va[i].v != cv)) {
continue;
}
*vp1 = va[i].v;
divide_flag = 2;
count1 = find_greatest_power(p1, n1, vp1, &d1, &t1, &len_t1, &divide_flag);
count2 = find_greatest_power(p2, n2, vp1, &d2, &t2, &len_t2, &divide_flag);
if (d2 <= 0 || d1 < d2 || count2 > count1) {
divide_flag = !divide_flag;
count1 = find_greatest_power(p1, n1, vp1, &d1, &t1, &len_t1, &divide_flag);
count2 = find_greatest_power(p2, n2, vp1, &d2, &t2, &len_t2, &divide_flag);
if (d2 <= 0 || d1 < d2 || count2 > count1) {
continue;
}
}
return va[i].count;
}
#if 1
break;
#endif
}
return 0;
}
#define VALUE_CNT 3
void
term_value(dp, p1, n1, loc)
double *dp;
token_type *p1;
int n1, loc;
{
int i, j, k;
int divide_flag = false;
int level, div_level = 0;
double d, sub_count, sub_sum;
for (i = 0; i < VALUE_CNT; i++)
dp[i] = 0.0;
for (i = loc; i < n1; i++) {
level = p1[i].level;
if (p1[i].kind == VARIABLE) {
if (divide_flag) {
dp[0] -= 1.0;
dp[1] -= p1[i].token.variable;
dp[2] -= p1[i].token.variable;
} else {
dp[0] += 1.0;
dp[1] += p1[i].token.variable;
dp[2] += p1[i].token.variable;
}
} else if (p1[i].kind == OPERATOR) {
if (level == 1 && (p1[i].token.operatr == PLUS || p1[i].token.operatr == MINUS))
break;
if (p1[i].token.operatr == DIVIDE) {
if (divide_flag && level >= div_level)
continue;
div_level = level;
divide_flag = true;
} else if (divide_flag && level <= div_level) {
divide_flag = false;
}
}
}
divide_flag = false;
for (j = loc + 1; j < i; j += 2) {
level = p1[j].level;
if (p1[j].token.operatr == DIVIDE) {
if (divide_flag && level >= div_level)
continue;
div_level = level;
divide_flag = true;
} else if (divide_flag && level <= div_level) {
divide_flag = false;
}
if (p1[j].token.operatr == POWER && level == p1[j+1].level && p1[j+1].kind == CONSTANT) {
d = p1[j+1].token.constant - 1.0;
sub_count = 0.0;
sub_sum = 0.0;
for (k = j - 1; k >= loc && p1[k].level >= level; k--) {
if (p1[k].kind == VARIABLE) {
sub_count += 1;
sub_sum += p1[k].token.variable;
}
}
if (divide_flag) {
dp[0] -= d * sub_count;
dp[2] -= d * sub_sum;
} else {
dp[0] += d * sub_count;
dp[2] += d * sub_sum;
}
}
}
}
/*
* This routine finds the additive term raised the highest power,
* with ordering and much flexibility.
*
* The number of terms raised to the highest power is returned, unless *vp1 == 0,
* then 0 is returned and *vp1 is set to a valid polynomial base variable, if any.
*/
int
find_greatest_power(p1, n1, vp1, pp1, tp1, lentp1, dcodep)
token_type *p1; /* pointer to expression */
int n1; /* length of expression */
long *vp1; /* polynomial base variable */
double *pp1; /* returned power of the returned term, which will be the highest power */
int *tp1, *lentp1; /* the returned term index and length */
int *dcodep; /* divide flag pointer indicates if term is a denominator; */
/* for example: for x^5 it is false, for 1/x^5 it is true. */
/* If equal to 2, set it; if equal to 3, ignore it. */
{
int i, j, k, ii;
double d;
int flag, divide_flag = false;
int level, div_level = 0;
long v = 0;
int was_power = false;
double last_va[VALUE_CNT];
double va[VALUE_CNT];
int rv;
int count = 0;
CLEAR_ARRAY(last_va);
*pp1 = 0.0;
*tp1 = -1;
*lentp1 = 0;
rv = *dcodep;
for (j = 0, i = 1;; i += 2) {
if (i >= n1 || ((p1[i].token.operatr == PLUS || p1[i].token.operatr == MINUS) && p1[i].level == 1)) {
divide_flag = false;
if (!was_power && *pp1 <= 1.0) {
for (k = j; k < i; k++) {
if (p1[k].kind == VARIABLE) {
if (*dcodep <= 1 && *dcodep != divide_flag)
continue;
if (*vp1) {
if (p1[k].token.variable == *vp1) {
term_value(va, p1, n1, j);
flag = (*pp1 == 1.0 && (rv > divide_flag));
if (*pp1 == 1.0 && rv == divide_flag) {
if (*tp1 != j)
count++;
for (ii = 0; ii < ARR_CNT(va); ii++) {
if (va[ii] == last_va[ii])
continue;
if (va[ii] < last_va[ii])
flag = true;
break;
}
} else if (*pp1 < 1.0 || flag) {
count = 1;
}
if (*pp1 < 1.0 || flag) {
blt(last_va, va, sizeof(last_va));
*pp1 = 1.0;
*tp1 = j;
rv = divide_flag;
}
break;
}
} else if ((p1[k].token.variable & VAR_MASK) > SIGN) {
v = p1[k].token.variable;
*pp1 = 1.0;
*tp1 = j;
rv = divide_flag;
break;
}
} else if (p1[k].kind == OPERATOR) {
if (p1[k].token.operatr == DIVIDE) {
if (divide_flag && p1[k].level >= div_level)
continue;
div_level = p1[k].level;
divide_flag = true;
} else if (divide_flag && p1[k].level <= div_level) {
divide_flag = false;
}
if (p1[k].token.operatr == POWER) {
level = p1[k].level;
do {
k += 2;
} while (k < i && p1[k].level > level);
k--;
}
}
}
}
if (i >= n1)
break;
j = i + 1;
was_power = false;
divide_flag = false;
continue;
}
level = p1[i].level;
if (p1[i].token.operatr == DIVIDE) {
if (divide_flag && level >= div_level)
continue;
div_level = level;
divide_flag = true;
} else if (divide_flag && level <= div_level) {
divide_flag = false;
}
if (p1[i].token.operatr == POWER && p1[i+1].kind == CONSTANT
&& (*vp1 || level == p1[i+1].level)) {
if (*dcodep <= 1 && *dcodep != divide_flag)
continue;
d = p1[i+1].token.constant;
for (k = i;;) {
if (p1[k-1].kind == VARIABLE) {
if (*vp1) {
if (p1[k-1].token.variable == *vp1) {
was_power = true;
term_value(va, p1, n1, j);
flag = (d == *pp1 && (rv > divide_flag));
if (d == *pp1 && rv == divide_flag) {
if (*tp1 != j)
count++;
for (ii = 0; ii < ARR_CNT(va); ii++) {
if (va[ii] == last_va[ii])
continue;
if (va[ii] < last_va[ii])
flag = true;
break;
}
} else if (d > *pp1 || flag) {
count = 1;
}
if (d > *pp1 || flag) {
blt(last_va, va, sizeof(last_va));
*pp1 = d;
*tp1 = j;
rv = divide_flag;
}
break;
}
} else if ((p1[k-1].token.variable & VAR_MASK) > SIGN) {
was_power = true;
if (d > *pp1) {
v = p1[k-1].token.variable;
*pp1 = d;
*tp1 = j;
rv = divide_flag;
}
break;
}
}
k -= 2;
if (k <= j)
break;
if (p1[k].level <= level)
break;
}
}
}
if (*vp1 == 0) {
*vp1 = v;
}
if (*tp1 >= 0) {
for (i = *tp1 + 1; i < n1; i += 2) {
if ((p1[i].token.operatr == PLUS || p1[i].token.operatr == MINUS)
&& p1[i].level == 1) {
break;
}
}
*lentp1 = i - *tp1;
}
if (*dcodep == 2)
*dcodep = rv;
return count;
}
Reference in New Issue View Git Blame Copy Permalink