mathomatic/super.c

541 lines
14 KiB
C

/*
* Group and combine algebraic fractions for Mathomatic.
*
* 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"
static int sf_recurse(token_type *equation, int *np, int loc, int level, int start_flag);
static int sf_sub(token_type *equation, int *np, int loc, int i1, int n1, int i2, int n2, int level, int start_flag);
static void
group_recurse(equation, np, loc, level)
token_type *equation; /* equation side pointer */
int *np; /* pointer to length of equation side */
int loc; /* starting location within equation side */
int level; /* current level of parentheses within the sub-expression starting at "loc" */
{
int i;
int len;
int di = -1, edi;
int group_flag = false;
int e1;
for (i = loc; i < *np && equation[i].level >= level;) {
if (equation[i].level > level) {
group_recurse(equation, np, i, level + 1);
i++;
for (; i < *np && equation[i].level > level; i += 2)
;
continue;
}
i++;
}
edi = e1 = i;
for (i = loc + 1; i < e1; i += 2) {
if (equation[i].level == level) {
if (equation[i].token.operatr == DIVIDE) {
if (di < 0) {
di = i;
continue;
}
group_flag = true;
for (len = i + 2; len < e1; len += 2) {
if (equation[len].level == level && equation[len].token.operatr != DIVIDE)
break;
}
len -= i;
if (edi == e1) {
i += len;
edi = i;
continue;
}
blt(scratch, &equation[i], len * sizeof(token_type));
blt(&equation[di+len], &equation[di], (i - di) * sizeof(token_type));
blt(&equation[di], scratch, len * sizeof(token_type));
edi += len;
i += len - 2;
} else {
if (di >= 0 && edi == e1)
edi = i;
}
}
}
if (group_flag) {
for (i = di + 1; i < edi; i++) {
if (equation[i].level == level && equation[i].kind == OPERATOR) {
#if DEBUG
if (equation[i].token.operatr != DIVIDE) {
error_bug("Bug in group_recurse().");
}
#endif
equation[i].token.operatr = TIMES;
}
equation[i].level++;
}
}
}
/*
* Group denominators of algebraic fractions together in an equation side.
* Grouping here means converting "a/b/c/d*e" to "a*e/(b*c*d)" or "a/(b*c*d)*e".
* Not guaranteed to put the grouped divisors last, reorder() puts divisors last.
*/
void
group_proc(equation, np)
token_type *equation; /* equation side pointer */
int *np; /* pointer to length of equation side */
{
group_recurse(equation, np, 0, 1);
}
/*
* Make equation side ready for display.
* Basically simplify, then convert non-integer constants in an equation side to fractions,
* when exactly equal to simple fractions.
* Also groups algebraic fraction denominators with group_proc() above.
*
* Return true if any fractions were created.
*/
int
fractions_and_group(equation, np)
token_type *equation;
int *np;
{
int rv = false;
elim_loop(equation, np);
if (fractions_display) {
rv = make_fractions(equation, np);
}
group_proc(equation, np);
return rv;
}
/*
* This function is the guts of the display command.
* Makes an equation space ready for display.
*
* Return true if any fractions were created.
*/
int
make_fractions_and_group(n)
int n; /* equation space number */
{
int rv = false;
if (empty_equation_space(n))
return false;
if (fractions_and_group(lhs[n], &n_lhs[n]))
rv = true;
if (n_rhs[n] > 0) {
if (fractions_and_group(rhs[n], &n_rhs[n]))
rv = true;
}
return rv;
}
/*
* Efficiently combine algebraic fractions added together
* by putting all terms over a common denominator.
* This means converting "(a/b)+(c/d)+f" directly to "(a*d+c*b+b*d*f)/b/d".
* The resulting expression is always equivalent to the original expression.
*
* If start_flag is 0, only combine denominators to convert complex fractions to simple fractions.
* Level one addition of fractions will be unchanged.
* Can make an expression over-complicated.
*
* If start_flag is 1, always combine denominators no matter what they are.
* This can easily make an expression large and complicated,
* but always results in a single, simple fraction.
* Used when solving for zero.
*
* If start_flag is 2, combine denominators
* and reduce the result by removing any polynomial GCD between them.
* Note that this wipes out globals tlhs[] and trhs[].
* This usually results in the simplest and most correct result.
*
* If start_flag is 3: same as start_flag = 2,
* except absolute value and imaginary denominators are combined, too.
* This simplifies all algebraic fractions into a single simple fraction.
* Removing the polynomial GCD usually prevents making a more complicated
* (or larger) algebraic fraction.
* Used by the fraction command.
*
* Return true if the equation side was modified.
*/
int
super_factor(equation, np, start_flag)
token_type *equation; /* pointer to the beginning of the equation side to process */
int *np; /* pointer to the length of the equation side */
int start_flag;
{
int rv;
group_proc(equation, np);
rv = sf_recurse(equation, np, 0, 1, start_flag);
organize(equation, np);
return rv;
}
static int
sf_recurse(equation, np, loc, level, start_flag)
token_type *equation;
int *np, loc, level, start_flag;
{
int modified = false;
int i, j, k;
int op;
int len1, len2;
if (!start_flag) {
for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
if (equation[i].level == level && equation[i].token.operatr == DIVIDE) {
start_flag = true;
break;
}
}
}
op = 0;
for (i = loc; i < *np && equation[i].level >= level;) {
if (equation[i].level > level) {
modified |= sf_recurse(equation, np, i, level + 1, start_flag);
i++;
for (; i < *np && equation[i].level > level; i += 2)
;
continue;
} else if (equation[i].kind == OPERATOR) {
op = equation[i].token.operatr;
}
i++;
}
if (modified || !start_flag)
return modified;
switch (op) {
case PLUS:
case MINUS:
break;
default:
return modified;
}
sf_again:
i = loc;
for (k = i + 1; k < *np && equation[k].level > level; k += 2)
;
len1 = k - i;
for (j = i + len1 + 1; j < *np && equation[j-1].level >= level; j += len2 + 1) {
for (k = j + 1; k < *np && equation[k].level > level; k += 2)
;
len2 = k - j;
#if 0
side_debug(0, &equation[i], len1);
side_debug(0, &equation[j], len2);
#endif
if (sf_sub(equation, np, loc, i, len1, j, len2, level + 1, start_flag)) {
#if 0
int junk;
printf("start_flag = %d\n", start_flag);
debug_string(0, "super_factor() successful:");
for (junk = 1; (loc + junk) < *np && equation[loc+junk].level > level; junk += 2)
;
side_debug(0, &equation[loc], junk);
#endif
modified = true;
goto sf_again;
}
}
return modified;
}
static int
sf_sub(equation, np, loc, i1, n1, i2, n2, level, start_flag)
token_type *equation;
int *np, loc, i1, n1, i2, n2, level, start_flag;
{
int i, j, k;
int b1, b2;
int len;
int e1, e2;
int op1, op2;
token_type *p1, *p2;
int np1, np2;
int div_flag1 = false, div_flag2 = false;
int rv;
e1 = i1 + n1;
e2 = i2 + n2;
op2 = equation[i2-1].token.operatr;
if (i1 <= loc) {
op1 = PLUS;
} else {
op1 = equation[i1-1].token.operatr;
}
for (i = i1 + 1; i < e1; i += 2) {
if (equation[i].level == level && equation[i].token.operatr == DIVIDE) {
div_flag1 = true;
break;
}
}
b1 = i + 1;
if (div_flag1) {
for (i += 2; i < e1; i += 2) {
if (equation[i].level <= level)
break;
}
}
for (j = i2 + 1; j < e2; j += 2) {
if (equation[j].level == level && equation[j].token.operatr == DIVIDE) {
div_flag2 = true;
break;
}
}
b2 = j + 1;
if (div_flag2) {
for (j += 2; j < e2; j += 2) {
if (equation[j].level <= level)
break;
}
}
if (!div_flag1 && !div_flag2)
return false;
#if 1
if (start_flag <= 2 && start_flag != 1) {
#if 0 /* Leave as 0; 1 will not simplify many imaginary number expressions, for example the tangent expansion. */
if (div_flag1 && found_var(&equation[b1], i - b1, IMAGINARY)) {
return false;
}
if (div_flag2 && found_var(&equation[b2], j - b2, IMAGINARY)) {
return false;
}
#endif
if (div_flag1) {
if (exp_is_absolute(&equation[b1], i - b1))
return false;
}
if (div_flag2) {
if (exp_is_absolute(&equation[b2], j - b2))
return false;
}
}
#endif
if (start_flag >= 2 && div_flag1 && div_flag2) {
#if DEBUG
debug_string(1, "Trying to find a polynomial GCD between 2 denominators in sf_sub():");
side_debug(1, &equation[b1], i - b1);
side_debug(1, &equation[b2], j - b2);
#endif
if ((rv = poly2_gcd(&equation[b1], i - b1, &equation[b2], j - b2, 0L, true)) > 0) {
p1 = tlhs;
np1 = n_tlhs;
p2 = trhs;
np2 = n_trhs;
goto do_gcd_super;
}
if (rv == 0 && poly2_gcd(&equation[b2], j - b2, &equation[b1], i - b1, 0L, true) > 0) {
p1 = trhs;
np1 = n_trhs;
p2 = tlhs;
np2 = n_tlhs;
goto do_gcd_super;
}
#if DEBUG
debug_string(1, "Done; polynomial GCD not found.");
#endif
}
if (n1 + n2 + (i - b1) + (j - b2) + 8 > n_tokens) {
error_huge();
}
if (!div_flag1) {
for (k = i1; k < e1; k++)
equation[k].level++;
}
if (!div_flag2) {
for (k = i2; k < e2; k++)
equation[k].level++;
}
len = (b1 - i1) - 1;
blt(scratch, &equation[i1], len * sizeof(token_type));
if (op1 == MINUS) {
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
scratch[len].level = level;
scratch[len].kind = CONSTANT;
scratch[len].token.constant = -1.0;
len++;
}
if (div_flag1) {
blt(&scratch[len], &equation[i], (e1 - i) * sizeof(token_type));
len += e1 - i;
}
if (div_flag2) {
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
blt(&scratch[len], &equation[b2], (j - b2) * sizeof(token_type));
len += j - b2;
}
for (k = 0; k < len; k++)
scratch[k].level += 2;
scratch[len].level = level + 1;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = op2;
len++;
k = len;
blt(&scratch[len], &equation[i2], (b2 - i2 - 1) * sizeof(token_type));
len += b2 - i2 - 1;
if (div_flag2) {
blt(&scratch[len], &equation[j], (e2 - j) * sizeof(token_type));
len += e2 - j;
}
if (div_flag1) {
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
blt(&scratch[len], &equation[b1], (i - b1) * sizeof(token_type));
len += i - b1;
}
for (; k < len; k++)
scratch[k].level += 2;
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = DIVIDE;
len++;
k = len;
if (div_flag1) {
blt(&scratch[len], &equation[b1], (i - b1) * sizeof(token_type));
len += i - b1;
}
if (div_flag1 && div_flag2) {
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
}
if (div_flag2) {
blt(&scratch[len], &equation[b2], (j - b2) * sizeof(token_type));
len += j - b2;
}
for (; k < len; k++)
scratch[k].level++;
end_mess:
if (*np + len - n1 - (n2 + 1) > n_tokens) {
error_huge();
}
if (op1 == MINUS) {
equation[i1-1].token.operatr = PLUS;
}
blt(&equation[i2-1], &equation[e2], (*np - e2) * sizeof(token_type));
*np -= n2 + 1;
blt(&equation[i1+len], &equation[e1], (*np - e1) * sizeof(token_type));
*np += len - n1;
blt(&equation[i1], scratch, len * sizeof(token_type));
return true;
do_gcd_super:
#if DEBUG
debug_string(1, "Found a polynomial GCD between 2 denominators in sf_sub().");
#endif
if (5 - i1 + e1 + (2*np2) + b2 - i2 + e2 - j + np1 > n_tokens) {
error_huge();
}
for (k = 0; k < np1; k++) {
p1[k].level += level;
}
for (k = 0; k < np2; k++) {
p2[k].level += level;
}
len = (b1 - i1) - 1;
blt(scratch, &equation[i1], len * sizeof(token_type));
if (op1 == MINUS) {
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
scratch[len].level = level;
scratch[len].kind = CONSTANT;
scratch[len].token.constant = -1.0;
len++;
}
/* if (div_flag1) { */
blt(&scratch[len], &equation[i], (e1 - i) * sizeof(token_type));
len += e1 - i;
/* } */
/* if (div_flag2) { */
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
blt(&scratch[len], p2, np2 * sizeof(token_type));
len += np2;
/* } */
for (k = 0; k < len; k++)
scratch[k].level += 2;
scratch[len].level = level + 1;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = op2;
len++;
k = len;
blt(&scratch[len], &equation[i2], (b2 - i2 - 1) * sizeof(token_type));
len += b2 - i2 - 1;
/* if (div_flag2) { */
blt(&scratch[len], &equation[j], (e2 - j) * sizeof(token_type));
len += e2 - j;
/* } */
/* if (div_flag1) { */
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
blt(&scratch[len], p1, np1 * sizeof(token_type));
len += np1;
/* } */
for (; k < len; k++)
scratch[k].level += 2;
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = DIVIDE;
len++;
k = len;
/* if (div_flag1) { */
blt(&scratch[len], &equation[b1], (i - b1) * sizeof(token_type));
len += (i - b1);
/* } */
/* if (div_flag1 && div_flag2) { */
scratch[len].level = level;
scratch[len].kind = OPERATOR;
scratch[len].token.operatr = TIMES;
len++;
/* } */
/* if (div_flag2) { */
blt(&scratch[len], p2, np2 * sizeof(token_type));
len += np2;
/* } */
for (; k < len; k++)
scratch[k].level++;
goto end_mess;
}