Logo Search packages:      
Sourcecode: mathomatic version File versions

factor.c

/*
 * Algebraic manipulator factorizing routines.
 *
 * Copyright (c) 1996 George Gesslein II.
 *
 * There are proper mathematical names for many algebraic rules.
 * I obviously don't know what they are.
 */

#include "am.h"
#include "externs.h"

static void try_factor();
static int  fplus_recurse();
static int  fplus_sub();
static int  big_fplus();
static int  ftimes_recurse();
static int  ftimes_sub();
static int  fpower_recurse();
static int  fpower_sub();
static int  fc_recurse();

static double nn, vv;
static double sq[] = {
      10, 2, 4, 2, 4, 6, 2, 6,
       4, 2, 4, 6, 6, 2, 6, 4,
       2, 6, 4, 6, 8, 4, 2, 4,
       2, 4, 8, 6, 4, 6, 2, 4,
       6, 2, 6, 6, 4, 2, 4, 6,
       2, 6, 4, 2, 4, 2,10, 2
};

/*
 * Factor the integer in "start".
 * Store unique factors in the "unique" array.
 * Return true if successful.
 */
int
factor_one(start)
double      start;
{
      int         i;
      double            ii;

      uno = 0;
      nn = start;
      if (nn == 0.0)
            return false;
      if (fabs(nn) >= MAX_K_INTEGER) {
            /* too large to factor */
            return false;
      }
      if (fmod(nn, 1.0) != 0.0) {
            /* not an integer */
            return false;
      }
      vv = 1.0 + sqrt(fabs(nn));
      try_factor(2.0);
      try_factor(3.0);
      try_factor(5.0);
      try_factor(7.0);
      ii = 1.0;
      while (ii <= vv) {
            for (i = 0; i < 48; i++) {
                  ii += sq[i];
                  try_factor(ii);
            }
      }
      if (fmod(nn, 1.0) != 0.0) {
            printf(_("Floating point error in factor_one()!\n"));
            return false;
      }
      if (nn != 1.0) {
            try_factor(nn);
      }
      return true;
}

/*
 * See if "arg" is one or more factors of "nn".
 * If so, save it and remove it from "nn".
 */
static void
try_factor(arg)
double      arg;
{
      for (;;) {
            if (fmod(nn, arg) == 0.0) {
                  if (uno > 0 && unique[uno-1] == arg) {
                        ucnt[uno-1]++;
                  } else {
                        unique[uno] = arg;
                        ucnt[uno] = 1;
                        uno++;
                  }
                  nn /= arg;
                  vv = 1.0 + sqrt(fabs(nn));
                  if (nn <= 1.0 && nn >= -1.0)
                        break;
            } else {
                  break;
            }
      }
}

/*
 * Convert unique[] back into an integer.
 */
double
multiply_out_unique()
{
      int   i, j;
      double      d;

      d = 1.0;
      for (i = 0; i < uno; i++) {
            for (j = 0; j < ucnt[i]; j++) {
                  d *= unique[i];
            }
      }
      return d;
}

/*
 * Factor integers in an equation side.
 */
factor_int(equation, np)
token_type  *equation;
int         *np;
{
      int   i, j;
      int   xsize;
      int   level;

      for (i = 0; i < *np; i += 2) {
            if (equation[i].kind == CONSTANT) {
                  if (factor_one(equation[i].token.constant) && uno > 0) {
                        level = equation[i].level;
                        if (uno > 1 && *np > 1)
                              level++;
                        xsize = -2;
                        for (j = 0; j < uno; j++) {
                              if (ucnt[j] > 1)
                                    xsize += 4;
                              else
                                    xsize += 2;
                        }
                        if (*np + xsize > n_tokens)
                              error_huge();
                        for (j = 0; j < uno; j++) {
                              if (ucnt[j] > 1)
                                    xsize = 4;
                              else
                                    xsize = 2;
                              if (j == 0)
                                    xsize -= 2;
                              if (xsize > 0) {
                                    blt(&equation[i+xsize], &equation[i], (*np - i) * sizeof(token_type));
                                    *np += xsize;
                                    if (j > 0) {
                                          i++;
                                          equation[i].kind = OPERATOR;
                                          equation[i].level = level;
                                          equation[i].token.operatr = TIMES;
                                          i++;
                                    }
                              }
                              equation[i].kind = CONSTANT;
                              equation[i].level = level;
                              equation[i].token.constant = unique[j];
                              if (ucnt[j] > 1) {
                                    equation[i].level = level + 1;
                                    i++;
                                    equation[i].kind = OPERATOR;
                                    equation[i].level = level + 1;
                                    equation[i].token.operatr = POWER;
                                    i++;
                                    equation[i].level = level + 1;
                                    equation[i].kind = CONSTANT;
                                    equation[i].token.constant = ucnt[j];
                              }
                        }
                  }
            }
      }
}

/*
 * Factor divides only.
 * (a/c + b/c) -> ((a+b)/c).
 *
 * Return true if equation side was modified.
 */
int
factor_divide(equation, np, v, d)
token_type  *equation;
int         *np;
long        v;
double            d;
{
      return fplus_recurse(equation, np, 0, 1, v, d, 2);
}

/*
 * Take care of subtraction and addition of the same expression
 * multiplied by constants.
 * (2*a + 3*a - a) -> (4*a).
 *
 * Return true if equation side was modified.
 */
int
subtract_itself(equation, np)
token_type  *equation;
int         *np;
{
      return fplus_recurse(equation, np, 0, 1, 0L, 0.0, 1);
}

/*
 * Factor equation side.
 * (a*c + b*c) -> (c*(a + b)).
 *
 * If "v" equals 0L, factor anything,
 * including identical bases raised to powers:
 * (x^2 + x) -> (x*(x + 1)).
 * If "d" equals 1.0, only factor identical bases raised to the power of a constant.
 *
 * If "v" is a variable, only factor expressions containing that variable.
 * If "d" is not equal to 0.0 or 1.0, factor only expressions containing "v" raised
 * to the power of "d".
 *
 * Return true if equation side was modified.
 */
int
factor_plus(equation, np, v, d)
token_type  *equation;
int         *np;
long        v;
double            d;
{
      return fplus_recurse(equation, np, 0, 1, v, d, 0);
}

static int
fplus_recurse(equation, np, loc, level, v, d, factor_code)
token_type  *equation;
int         *np, loc, level;
long        v;
double            d;
int         factor_code;
{
      int   modified;
      int   i, j, k;
      int   op;
      int   len1, len2;

      modified = false;
      op = 0;
      for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
            if (equation[i].level == level) {
                  op = equation[i].token.operatr;
                  break;
            }
      }
      switch (op) {
      case PLUS:
      case MINUS:
            break;
      default:
            goto skip_plus;
      }
      for (i = loc;;) {
f_again:
            for (k = i + 1;; k += 2) {
                  if (k >= *np || equation[k].level <= level)
                        break;
            }
            len1 = k - i;
            for (j = i + len1 + 1;; j += len2 + 1) {
                  if (j >= *np || equation[j-1].level < level)
                        break;
                  for (k = j + 1;; k += 2) {
                        if (k >= *np || equation[k].level <= level)
                              break;
                  }
                  len2 = k - j;
                  if (fplus_sub(equation, np, loc, i, len1, j, len2, level + 1, factor_code == 1, factor_code == 2, v, d)) {
                        modified = true;
                        goto f_again;
                  }
            }
            i += len1 + 1;
            if (i >= *np || equation[i-1].level < level)
                  break;
      }
      if (modified)
            return true;
skip_plus:
      for (i = loc;;) {
            if (i >= *np || equation[i].level < level)
                  break;
            if (equation[i].level > level) {
                  modified |= fplus_recurse(equation, np, i, level + 1, v, d, factor_code);
                  i++;
                  for (; i < *np && equation[i].level > level; i += 2)
                        ;
                  continue;
            }
            i++;
      }
      return modified;
}

static int
fplus_sub(equation, np, loc, i1, n1, i2, n2, level, whole_flag, div_only, v, d)
token_type  *equation;
int         *np, loc, i1, n1, i2, n2;
int         level, whole_flag, div_only;
long        v;
double            d;
{
      int   e1, e2;
      int   op1, op2;
      int   i, j, k, l, m;
      int   b1, b2;
      int   len;
      int   sop1;
      int   diff_sign;
      int   ai, aj;
      int   flag1, flag2;
      int   same_flag;
      double      save_k1, save_k2;
      double      save_d1, save_d2;
      double      power;
      double      d1, d2;

      e1 = i1 + n1;
      e2 = i2 + n2;
      op2 = equation[i2-1].token.operatr;
      if (i1 <= loc) {
            op1 = PLUS;
      } else
            op1 = equation[i1-1].token.operatr;
      i = i1 - 1;
f_outer:
      b1 = i + 1;
      if (b1 >= e1)
            return false;
      if (whole_flag) {
            i = e1;
            if (n1 > 1 && equation[b1].kind == CONSTANT
                && equation[b1+1].level == level
                && (equation[b1+1].token.operatr == TIMES
                || equation[b1+1].token.operatr == DIVIDE)) {
                  b1 += 2;
            }
      } else {
            i = b1 + 1;
            for (; i < e1; i += 2) {
                  if (equation[i].level == level
                      && (equation[i].token.operatr == TIMES
                      || equation[i].token.operatr == DIVIDE)) {
                        break;
                  }
            }
      }
      if (b1 <= i1)
            sop1 = TIMES;
      else
            sop1 = equation[b1-1].token.operatr;
      if (div_only && sop1 != DIVIDE)
            goto f_outer;
      if (i - b1 == 1 && equation[b1].kind == CONSTANT
          && fabs(equation[b1].token.constant) == 1.0) {
            goto f_outer;
      }
      if (!whole_flag && v && (v != MATCH_ANY)) {
            if (d == 0.0 || d == 1.0) {
                  for (k = b1;; k += 2) {
                        if (k >= i)
                              goto f_outer;
                        if (equation[k].kind == VARIABLE
                            && equation[k].token.variable == v) {
                              break;
                        }
                  }
            } else {
                  for (k = b1 + 1;; k += 2) {
                        if (k >= i)
                              goto f_outer;
                        if (equation[k].token.operatr == POWER
                            && equation[k].level == equation[k+1].level
                            && equation[k+1].kind == CONSTANT
                            && equation[k+1].token.constant == d) {
                              for (l = k - 1; l >= 0; l--) {
                                    if (equation[l].level < equation[k].level)
                                          break;
                                    if (equation[l].kind == VARIABLE
                                        && equation[l].token.variable == v)
                                          goto factor_this;
                              }
                        }
                  }
            }
      }
factor_this:
      j = i2 - 1;
f_inner:
      b2 = j + 1;
      if (b2 >= e2)
            goto f_outer;
      if (whole_flag) {
            j = e2;
            if (n2 > 1 && equation[b2].kind == CONSTANT
                && equation[b2+1].level == level
                && (equation[b2+1].token.operatr == TIMES
                || equation[b2+1].token.operatr == DIVIDE)) {
                  b2 += 2;
            }
      } else {
            j = b2 + 1;
            for (; j < e2; j += 2) {
                  if (equation[j].level == level
                      && (equation[j].token.operatr == TIMES
                      || equation[j].token.operatr == DIVIDE)) {
                        break;
                  }
            }
            if (b2 <= i2) {
                  if (sop1 == DIVIDE)
                        goto f_inner;
            } else {
                  if (equation[b2-1].token.operatr != sop1)
                        goto f_inner;
            }
      }
      if (j - b2 == 1 && equation[b2].kind == CONSTANT
          && fabs(equation[b2].token.constant) == 1.0) {
            goto f_inner;
      }
      ai = i;
      aj = j;
      if (whole_flag) {
            flag1 = false;
            flag2 = false;
            if (b1 > i1) {
                  b1 = i1;
                  save_k1 = equation[b1].token.constant;
                  equation[b1].token.constant = 1.0;
                  flag1 = true;
            }
            if (b2 > i2) {
                  b2 = i2;
                  save_k2 = equation[b2].token.constant;
                  equation[b2].token.constant = 1.0;
                  flag2 = true;
            }
      }
      same_flag = se_compare(&equation[b1], i - b1, &equation[b2], j - b2, &diff_sign);
      if (whole_flag) {
            if (flag1) {
                  equation[i1].token.constant = save_k1;
                  b1 += 2;
            }
            if (flag2) {
                  equation[i2].token.constant = save_k2;
                  b2 += 2;
            }
      }
      if (same_flag) {
            power = 1.0;
more_power:
            if (sop1 == DIVIDE) {
                  scratch[0].level = level;
                  scratch[0].kind = CONSTANT;
                  scratch[0].token.constant = 1.0;
                  scratch[1].level = level;
                  scratch[1].kind = OPERATOR;
                  scratch[1].token.operatr = DIVIDE;
                  len = 2;
            } else
                  len = 0;
            k = len;
            blt(&scratch[len], &equation[b1], (ai - b1) * sizeof(*equation));
            len += ai - b1;
            if (power != 1.0) {
                  for (; k < len; k++)
                        scratch[k].level += 2;
                  scratch[len].level = level + 1;
                  scratch[len].kind = OPERATOR;
                  scratch[len].token.operatr = POWER;
                  len++;
                  scratch[len].level = level + 1;
                  scratch[len].kind = CONSTANT;
                  scratch[len].token.constant = power;
                  len++;
                  if (always_positive(power))
                        diff_sign = false;
            } else if (b1 == i1 && ai == e1) {
                  for (; k < len; k++)
                        scratch[k].level++;
            }
            scratch[len].level = level;
            scratch[len].kind = OPERATOR;
            scratch[len].token.operatr = TIMES;
            len++;
            k = len;
            blt(&scratch[len], &equation[i1], (b1 - i1) * sizeof(*equation));
            len += b1 - i1;
            if (ai != i) {
                  l = len;
                  m = len + ai - b1;
                  blt(&scratch[len], &equation[b1], (i - b1) * sizeof(*equation));
                  len += i - b1;
                  if (b1 == i1 && i == e1)
                        for (; l < len; l++)
                              scratch[l].level++;
                  l = m;
                  m++;
                  for (; m < len; m++)
                        scratch[m].level++;
                  scratch[len].level = scratch[l].level + 1;
                  scratch[len].kind = OPERATOR;
                  scratch[len].token.operatr = MINUS;
                  len++;
                  scratch[len].level = scratch[l].level + 1;
                  scratch[len].kind = CONSTANT;
                  scratch[len].token.constant = power;
                  len++;
                  scratch[len].level = level;
                  scratch[len].kind = OPERATOR;
                  scratch[len].token.operatr = TIMES;
                  len++;
            }
            scratch[len].level = level;
            scratch[len].kind = CONSTANT;
            if (op1 == MINUS) {
                  scratch[len].token.constant = -1.0;
            } else
                  scratch[len].token.constant = 1.0;
            len++;
            blt(&scratch[len], &equation[i], (e1 - i) * sizeof(*equation));
            len += e1 - i;
            for (; k < len; k++)
                  scratch[k].level += 2;
            scratch[len].level = level + 1;
            scratch[len].kind = OPERATOR;
            diff_sign ^= (op2 == MINUS);
            if (diff_sign)
                  scratch[len].token.operatr = MINUS;
            else
                  scratch[len].token.operatr = PLUS;
            len++;
            k = len;
            if (aj != j) {
                  if (len + n2 + 2 > n_tokens) {
                        error_huge();
                  }
            } else {
                  if (len + (b2 - i2) + (e2 - j) + 1 > n_tokens) {
                        error_huge();
                  }
            }
            blt(&scratch[len], &equation[i2], (b2 - i2) * sizeof(*equation));
            len += b2 - i2;
            if (aj != j) {
                  m = len + aj - b2;
                  blt(&scratch[len], &equation[b2], (j - b2) * sizeof(*equation));
                  len += j - b2;
                  l = m;
                  m++;
                  for (; m < len; m++)
                        scratch[m].level++;
                  scratch[len].level = scratch[l].level + 1;
                  scratch[len].kind = OPERATOR;
                  scratch[len].token.operatr = MINUS;
                  len++;
                  scratch[len].level = scratch[l].level + 1;
                  scratch[len].kind = CONSTANT;
                  scratch[len].token.constant = power;
                  len++;
            } else {
                  scratch[len].level = level;
                  scratch[len].kind = CONSTANT;
                  scratch[len].token.constant = 1.0;
                  len++;
            }
            blt(&scratch[len], &equation[j], (e2 - j) * sizeof(*equation));
            len += e2 - j;
            for (; k < len; k++)
                  scratch[k].level += 2;
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(*equation));
            *np -= n2 + 1;
            blt(&equation[i1+len], &equation[e1], (*np - e1) * sizeof(*equation));
            *np += len - n1;
            blt(&equation[i1], scratch, len * sizeof(*equation));
            return true;
      }
      if (whole_flag)
            return false;
      if (v || div_only)
            goto f_inner;
      if (b1 == i1 && i == e1)
            k = level;
      else
            k = level + 1;
      save_d1 = 1.0;
      for (l = b1 + 1;; l += 2) {
            if (l >= i) {
                  break;
            }
            if (equation[l].level == k && equation[l].token.operatr == POWER) {
                  if (equation[l+1].level == k && equation[l+1].kind == CONSTANT) {
                        save_d1 = equation[l+1].token.constant;
                        if (save_d1 <= 0.0)
                              goto f_inner;
                  } else
                        save_d1 = -1.0;
                  ai = l;
                  break;
            }
      }
      if (b2 == i2 && j == e2)
            k = level;
      else
            k = level + 1;
      save_d2 = 1.0;
      for (l = b2 + 1;; l += 2) {
            if (l >= j) {
                  break;
            }
            if (equation[l].level == k && equation[l].token.operatr == POWER) {
                  if (equation[l+1].level == k && equation[l+1].kind == CONSTANT) {
                        save_d2 = equation[l+1].token.constant;
                        if (save_d2 <= 0.0)
                              goto f_inner;
                  } else
                        save_d2 = -1.0;
                  aj = l;
                  break;
            }
      }
      if (ai == i && aj == j)
            goto f_inner;
      if (ai - b1 == 1 && equation[b1].kind == CONSTANT)
            goto f_inner;
      if (d == 1.0 && (save_d1 < 0.0 || save_d2 < 0.0))
            goto f_inner;
      if (se_compare(&equation[b1], ai - b1, &equation[b2], aj - b2, &diff_sign)) {
            if (save_d1 > 0.0 || save_d2 > 0.0) {
                  if (save_d1 < 0.0)
                        power = save_d2;
                  else if (save_d2 < 0.0)
                        power = save_d1;
                  else {
                        power = min(save_d1, save_d2);
                        if (!diff_sign && (fmod(power, 1.0) != 0.0)) {
                              if (fmod(max(save_d1, save_d2) - power, 1.0) == 0.0) {
                                    goto more_power;
                              }
                        }
                  }
                  if (power < 1.0)
                        goto f_inner;
                  modf(power, &power);
                  goto more_power;
            }
            d1 = i - ai;
            d2 = j - aj;
            if (d1 == d2) {
                  d1 = 1.0;
                  d2 = 1.0;
                  if ((ai + 2) < i) {
                        k = equation[ai].level;
                        if (equation[ai+1].level == (k + 1)
                            && equation[ai+2].level == (k + 1)
                            && equation[ai+1].kind == CONSTANT
                            && (equation[ai+2].token.operatr == TIMES
                            || equation[ai+2].token.operatr == DIVIDE)) {
                              d1 = fabs(equation[ai+1].token.constant);
                        }
                  }
                  if ((aj + 2) < j) {
                        k = equation[aj].level;
                        if (equation[aj+1].level == (k + 1)
                            && equation[aj+2].level == (k + 1)
                            && equation[aj+1].kind == CONSTANT
                            && (equation[aj+2].token.operatr == TIMES
                            || equation[aj+2].token.operatr == DIVIDE)) {
                              d2 = fabs(equation[aj+1].token.constant);
                        }
                  }
            }
            if (d1 <= d2) {
                  len = big_fplus(equation, level, diff_sign, sop1,
                      op1, op2, i1, i2, b1, b2, ai, aj, i, j, e1, e2);
            } else {
                  len = big_fplus(equation, level, diff_sign, sop1,
                      op2, op1, i2, i1, b2, b1, aj, ai, j, i, e2, e1);
            }
            goto end_mess;
      }
      goto f_inner;
}

static int
big_fplus(equation, level, diff_sign, sop1, op1, op2, i1, i2, b1, b2, ai, aj, i, j, e1, e2)
token_type  *equation;
int         level;
int         diff_sign;
int         sop1;
int         op1, op2;
int         i1, i2;
int         b1, b2;
int         ai, aj;
int         i, j;
int         e1, e2;
{
      int   k, l, m, n, o;
      int   len;

      if (sop1 == DIVIDE) {
            scratch[0].level = level;
            scratch[0].kind = CONSTANT;
            scratch[0].token.constant = 1.0;
            scratch[1].level = level;
            scratch[1].kind = OPERATOR;
            scratch[1].token.operatr = DIVIDE;
            len = 2;
      } else
            len = 0;
      k = len;
      o = len + ai - b1;
      blt(&scratch[len], &equation[b1], (i - b1) * sizeof(*equation));
      len += i - b1;
      if (b1 == i1 && i == e1) {
            for (; k < len; k++)
                  scratch[k].level++;
      }
      scratch[len].level = level;
      scratch[len].kind = OPERATOR;
      scratch[len].token.operatr = TIMES;
      len++;
      k = len;
      blt(&scratch[len], &equation[i1], (b1 - i1) * sizeof(*equation));
      len += b1 - i1;
      scratch[len].level = level;
      scratch[len].kind = CONSTANT;
      if (op1 == MINUS) {
            scratch[len].token.constant = -1.0;
      } else
            scratch[len].token.constant = 1.0;
      len++;
      blt(&scratch[len], &equation[i], (e1 - i) * sizeof(*equation));
      len += e1 - i;
      for (; k < len; k++)
            scratch[k].level += 2;
      scratch[len].level = level + 1;
      scratch[len].kind = OPERATOR;
      if (op2 == MINUS)
            scratch[len].token.operatr = MINUS;
      else
            scratch[len].token.operatr = PLUS;
      len++;
      k = len;
      blt(&scratch[len], &equation[i2], (b2 - i2) * sizeof(*equation));
      len += b2 - i2;
      if (len + (e2 - b2) + 2 * (i - ai) + 2 > n_tokens) {
            error_huge();
      }
      n = len;
      m = len + aj - b2;
      blt(&scratch[len], &equation[b2], (j - b2) * sizeof(*equation));
      len += j - b2;
      l = m;
      m++;
      for (; m < len; m++)
            scratch[m].level++;
      if (diff_sign && b2 == i2 && j == e2) {
            for (; n < len; n++)
                  scratch[n].level++;
      }
      scratch[len].level = scratch[l].level + 1;
      scratch[len].kind = OPERATOR;
      scratch[len].token.operatr = MINUS;
      len++;
      m = len;
      blt(&scratch[len], &equation[ai+1], (i - (ai + 1)) * sizeof(*equation));
      len += i - (ai + 1);
      n = min_level(&equation[ai+1], i - (ai + 1));
      n = scratch[l].level + 2 - n;
      for (; m < len; m++)
            scratch[m].level += n;
      if (diff_sign) {
            scratch[len].level = level;
            scratch[len].kind = OPERATOR;
            if (sop1 == DIVIDE)
                  scratch[len].token.operatr = TIMES;
            else
                  scratch[len].token.operatr = DIVIDE;
            len++;
            scratch[len].level = level + 1;
            scratch[len].kind = CONSTANT;
            scratch[len].token.constant = -1.0;
            len++;
            blt(&scratch[len], &scratch[o], (i - ai) * sizeof(*equation));
            len += i - ai;
      }
      blt(&scratch[len], &equation[j], (e2 - j) * sizeof(*equation));
      len += e2 - j;
      for (; k < len; k++)
            scratch[k].level += 2;
      return len;
}

/*
 * Factor equation side.
 * a^b * a^c -> a^(b + c).
 * Return true if equation side was modified.
 */
int
factor_times(equation, np)
token_type  *equation;
int         *np;
{
      return ftimes_recurse(equation, np, 0, 1);
}

static int
ftimes_recurse(equation, np, loc, level)
token_type  *equation;
int         *np, loc, level;
{
      int   modified;
      int   i, j, k;
      int   op;
      int   len1, len2;

      modified = false;
      op = 0;
      for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
            if (equation[i].level == level) {
                  op = equation[i].token.operatr;
                  break;
            }
      }
      switch (op) {
      case TIMES:
      case DIVIDE:
            break;
      default:
            goto skip_times;
      }
      for (i = loc;;) {
f_again:
            for (k = i + 1;; k += 2) {
                  if (k >= *np || equation[k].level <= level)
                        break;
            }
            len1 = k - i;
            for (j = i + len1 + 1;; j += len2 + 1) {
                  if (j >= *np || equation[j-1].level < level)
                        break;
                  for (k = j + 1;; k += 2) {
                        if (k >= *np || equation[k].level <= level)
                              break;
                  }
                  len2 = k - j;
                  if (ftimes_sub(equation, np, loc, i, len1, j, len2, level + 1)) {
                        modified = true;
                        goto f_again;
                  }
            }
            i += len1 + 1;
            if (i >= *np || equation[i-1].level < level)
                  break;
      }
      if (modified)
            return true;
skip_times:
      for (i = loc;;) {
            if (i >= *np || equation[i].level < level)
                  break;
            if (equation[i].level > level) {
                  modified |= ftimes_recurse(equation, np, i, level + 1);
                  i++;
                  for (; i < *np && equation[i].level > level; i += 2)
                        ;
                  continue;
            }
            i++;
      }
      return modified;
}

static int
ftimes_sub(equation, np, loc, i1, n1, i2, n2, level)
token_type  *equation;
int         *np, loc, i1, n1, i2, n2, level;
{
      int   e1, e2;
      int   op1, op2;
      int   i, j, k;
      int   len, rlen1, len2;
      int   diff_sign;
      int   both_divide;

      op2 = equation[i2-1].token.operatr;
      e1 = i1 + n1;
      e2 = i2 + n2;
      if (i1 <= loc) {
            op1 = TIMES;
      } else {
            op1 = equation[i1-1].token.operatr;
      }
      if ((n1 == 1 && equation[i1].kind == CONSTANT)
          || (n2 == 1 && equation[i2].kind == CONSTANT)) {
            return false;
      }
      both_divide = (op1 == DIVIDE && op2 == DIVIDE);
#if   true
      if (se_compare(&equation[i1], e1 - i1, &equation[i2], e2 - i2, &diff_sign)) {
            i = e1;
            j = e2;
            goto common_base;
      }
#endif
      for (i = i1 + 1; i < e1; i += 2) {
            if (equation[i].level == level
                && equation[i].token.operatr == POWER) {
                  break;
            }
      }
      for (j = i2 + 1; j < e2; j += 2) {
            if (equation[j].level == level
                && equation[j].token.operatr == POWER) {
                  break;
            }
      }
      if (se_compare(&equation[i1], i - i1, &equation[i2], j - i2, &diff_sign)) {
            goto common_base;
      }
#if   false
      if (i != e1 && se_compare(&equation[i1], e1 - i1, &equation[i2], j - i2, &diff_sign)) {
            i = e1;
            goto common_base;
      }
      if (j != e2 && se_compare(&equation[i1], i - i1, &equation[i2], e2 - i2, &diff_sign)) {
            j = e2;
            goto common_base;
      }
#endif
      return false;

common_base:
      rlen1 = ((i == e1) ? 1 : (e1 - i - 1));
      len = (i - i1) + 1 + ((op1 == DIVIDE && !both_divide) ? 2 : 0)
          + rlen1 + 1
          + ((j == e2) ? 1 : (e2 - j - 1));
      len -= n1;
      if (j - i2 == 1 && equation[i2].kind == CONSTANT && equation[i2].token.constant == -1.0)
            return false;
      if (diff_sign) {
            if (j - i2 == 1 && equation[i2].kind == CONSTANT)
                  return false;
            len2 = 2 + e2 - j;
            if (*np + len2 + len > n_tokens) {
                  error_huge();
            }
            blt(&equation[e2+len2], &equation[e2], (*np - e2) * sizeof(token_type));
            *np += len2;
            equation[e2].level = level - 1;
            equation[e2].kind = OPERATOR;
            equation[e2].token.operatr = op2;
            equation[e2+1].level = level;
            equation[e2+1].kind = CONSTANT;
            equation[e2+1].token.constant = -1.0;
            blt(&equation[e2+2], &equation[j], (e2 - j) * sizeof(token_type));
      }
      if (*np + len > n_tokens) {
            error_huge();
      }
      blt(&equation[e1+len], &equation[e1], (*np - e1) * sizeof(*equation));
      *np += len;
      if (i == e1) {
            for (k = i1; k < e1; k++)
                  equation[k].level++;
            equation[i].level = level;
            equation[i].kind = OPERATOR;
            equation[i].token.operatr = POWER;
            equation[i+1].level = level;
            equation[i+1].kind = CONSTANT;
            equation[i+1].token.constant = 1.0;
      }
      if (op1 == DIVIDE && !both_divide) {
            equation[i1-1].token.operatr = TIMES;
            blt(&equation[i+3], &equation[i+1], rlen1 * sizeof(*equation));
            i++;
            equation[i].level = level;
            equation[i].kind = CONSTANT;
            equation[i].token.constant = -1.0;
            i++;
            equation[i].level = level;
            equation[i].kind = OPERATOR;
            equation[i].token.operatr = TIMES;
            binary_parenthesize(equation, i + 1 + rlen1, i);
      }
      i += rlen1 + 1;
      equation[i].level = level;
      equation[i].kind = OPERATOR;
      if (op2 == DIVIDE && !both_divide)
            equation[i].token.operatr = MINUS;
      else
            equation[i].token.operatr = PLUS;
      if (j == e2) {
            equation[i+1].level = level;
            equation[i+1].kind = CONSTANT;
            equation[i+1].token.constant = 1.0;
      } else {
            blt(&equation[i+1], &equation[j+len+1], (e2 - j - 1) * sizeof(*equation));
      }
      binary_parenthesize(equation, i + e2 - j, i);
      blt(&equation[i2+len-1], &equation[e2+len], (*np - (e2 + len)) * sizeof(*equation));
      *np -= n2 + 1;
      return true;
}

/*
 * Factor equation side.
 * a^c * b^c -> (a * b)^c.
 * Return true if equation side was modified.
 */
int
factor_power(equation, np)
token_type  *equation;
int         *np;
{
      return fpower_recurse(equation, np, 0, 1);
}

static int
fpower_recurse(equation, np, loc, level)
token_type  *equation;
int         *np, loc, level;
{
      int   modified;
      int   i, j, k;
      int   op;
      int   len1, len2;

      modified = false;
      op = 0;
      for (i = loc + 1; i < *np && equation[i].level >= level; i += 2) {
            if (equation[i].level == level) {
                  op = equation[i].token.operatr;
                  break;
            }
      }
      switch (op) {
      case TIMES:
      case DIVIDE:
            break;
      default:
            goto skip_power;
      }
      for (i = loc;;) {
f_again:
            for (k = i + 1;; k += 2) {
                  if (k >= *np || equation[k].level <= level)
                        break;
            }
            len1 = k - i;
            for (j = i + len1 + 1;; j += len2 + 1) {
                  if (j >= *np || equation[j-1].level < level)
                        break;
                  for (k = j + 1;; k += 2) {
                        if (k >= *np || equation[k].level <= level)
                              break;
                  }
                  len2 = k - j;
                  if (fpower_sub(equation, np, loc, i, len1, j, len2, level + 1)) {
                        modified = true;
                        goto f_again;
                  }
            }
            i += len1 + 1;
            if (i >= *np || equation[i-1].level < level)
                  break;
      }
skip_power:
      for (i = loc;;) {
            if (i >= *np || equation[i].level < level)
                  break;
            if (equation[i].level > level) {
                  modified |= fpower_recurse(equation, np, i, level + 1);
                  i++;
                  for (; i < *np && equation[i].level > level; i += 2)
                        ;
                  continue;
            }
            i++;
      }
      return modified;
}

static int
fpower_sub(equation, np, loc, i1, n1, i2, n2, level)
token_type  *equation;
int         *np, loc, i1, n1, i2, n2, level;
{
      int   e1, e2;
      int   op1, op2;
      int   pop1;
      int   start2;
      int   i, j, k;
      int   b1, b2;
      int   len;
      int   diff_sign;
      int   all_divide;
      token_type  temp;

      e1 = i1 + n1;
      e2 = i2 + n2;
      op2 = equation[i2-1].token.operatr;
      if (i1 <= loc) {
            op1 = TIMES;
      } else
            op1 = equation[i1-1].token.operatr;
      for (i = i1; i < e1; i++) {
            if (equation[i].level == level
                && equation[i].kind == OPERATOR
                && equation[i].token.operatr == POWER) {
                  break;
            }
      }
      for (j = i2; j < e2; j++) {
            if (equation[j].level == level
                && equation[j].kind == OPERATOR
                && equation[j].token.operatr == POWER) {
                  break;
            }
      }
      if (i >= e1 || j >= e2)
            return false;
      start2 = j;
fp_outer:
      pop1 = equation[i].token.operatr;
      if (pop1 == POWER)
            pop1 = TIMES;
      b1 = i + 1;
      if (b1 >= e1)
            return false;
#if   false
      i = e1;
#else
      i = b1 + 1;
      for (; i < e1; i += 2) {
            if (equation[i].level == (level + 1)
                && (equation[i].token.operatr == TIMES
                || equation[i].token.operatr == DIVIDE)) {
                  break;
            }
      }
#endif
      temp.level = 1;
      temp.kind = CONSTANT;
      temp.token.constant = 1.0;
      if (se_compare(&equation[b1], i - b1, &temp, 1, &diff_sign)) {
            goto fp_outer;
      }
      j = start2;
fp_inner:
      b2 = j + 1;
      if (b2 >= e2)
            goto fp_outer;
#if   false
      j = e2;
#else
      j = b2 + 1;
      for (; j < e2; j += 2) {
            if (equation[j].level == (level + 1)
                && (equation[j].token.operatr == TIMES
                || equation[j].token.operatr == DIVIDE)) {
                  break;
            }
      }
#endif
      if (equation[b2-1].token.operatr == POWER) {
            if (pop1 != TIMES)
                  goto fp_inner;
      } else if (equation[b2-1].token.operatr != pop1) {
            goto fp_inner;
      }
      if (se_compare(&equation[b1], i - b1, &equation[b2], j - b2, &diff_sign)) {
            if (op2 == DIVIDE)
                  diff_sign = !diff_sign;
            all_divide = (op1 == DIVIDE && diff_sign);
            len = 0;
            blt(&scratch[len], &equation[i1], (b1 - i1) * sizeof(*equation));
            len += b1 - i1;
            scratch[len].level = level + 1;
            scratch[len].kind = CONSTANT;
            if (!all_divide && op1 == DIVIDE) {
                  scratch[len].token.constant = -1.0;
            } else {
                  scratch[len].token.constant = 1.0;
            }
            len++;
            blt(&scratch[len], &equation[i], (e1 - i) * sizeof(*equation));
            len += e1 - i;
            for (k = 0; k < len; k++)
                  scratch[k].level += 2;
            scratch[len].level = level + 1;
            scratch[len].kind = OPERATOR;
            scratch[len].token.operatr = TIMES;
            len++;
            k = len;
            blt(&scratch[len], &equation[i2], (b2 - i2) * sizeof(*equation));
            len += b2 - i2;
            scratch[len].level = level + 1;
            scratch[len].kind = CONSTANT;
            if (!all_divide && diff_sign) {
                  scratch[len].token.constant = -1.0;
            } else {
                  scratch[len].token.constant = 1.0;
            }
            len++;
            blt(&scratch[len], &equation[j], (e2 - j) * sizeof(*equation));
            len += e2 - j;
            for (; k < len; k++)
                  scratch[k].level += 2;
            scratch[len].level = level;
            scratch[len].kind = OPERATOR;
            scratch[len].token.operatr = POWER;
            len++;
            if (pop1 == DIVIDE) {
                  scratch[len].level = level + 1;
                  scratch[len].kind = CONSTANT;
                  scratch[len].token.constant = 1.0;
                  len++;
                  scratch[len].level = level + 1;
                  scratch[len].kind = OPERATOR;
                  scratch[len].token.operatr = DIVIDE;
                  len++;
            }
            k = len;
            blt(&scratch[len], &equation[b1], (i - b1) * sizeof(*equation));
            len += i - b1;
            for (; k < len; k++)
                  scratch[k].level++;
            if (*np + len - n1 - (n2 + 1) > n_tokens) {
                  error_huge();
            }
            if (!all_divide && op1 == DIVIDE) {
                  equation[i1-1].token.operatr = TIMES;
            }
            blt(&equation[i2-1], &equation[e2], (*np - e2) * sizeof(*equation));
            *np -= n2 + 1;
            blt(&equation[i1+len], &equation[e1], (*np - e1) * sizeof(*equation));
            *np += len - n1;
            blt(&equation[i1], scratch, len * sizeof(*equation));
            return true;
      }
      goto fp_inner;
}

/*
 * Factor constants in equation side.
 *
 * This routine is often necessary because the expression compare (se_compare())
 * does not return a multiplier (except for +/-1.0).
 * This routine is not used during polynomial operations.
 * It is required for simplification of algebraic fractions, etc.
 *
 * If "level_code" is 0, all additive expressions are normalized
 * by making at least one coefficient unity by factoring out
 * the smallest constant.
 *
 * If "level_code" is 1, any level 1 additive expression is factored
 * nicely for readability, while all deeper levels are normalized.
 *
 * If "level_code" is 2, nothing is normalized unless it increases
 * readability.
 *
 * If "level_code" is 3, nothing is done.
 *
 * Return true if equation side was modified.
 */
int
factor_constants(equation, np, level_code)
token_type  *equation;
int         *np;
int         level_code;
{
      if (level_code > 2)
            return false;
      return fc_recurse(equation, np, 0, 1, level_code);
}

static int
fc_recurse(equation, np, loc, level, level_code)
token_type  *equation;
int         *np, loc, level;
int         level_code;
{
      int   modified;
      int   i, j, k;
      int   op;
      int   neg_flag;
      double      minimum;
      double      d;
      int   first;
      int   count;

      modified = false;
      for (i = loc;;) {
            if (i >= *np || equation[i].level < level)
                  break;
            if (equation[i].level > level) {
                  modified |= fc_recurse(equation, np, i, level + 1, level_code);
                  i++;
                  for (; i < *np && equation[i].level > level; i += 2)
                        ;
                  continue;
            }
            i++;
      }
      if (modified)
            return true;
      minimum = 1.0;
      neg_flag = true;
      first = true;
      count = 0;
      for (i = loc;;) {
break_cont:
            if (i >= *np || equation[i].level < level)
                  break;
            if (equation[i].level == level) {
                  switch (equation[i].kind) {
                  case CONSTANT:
                        if (i == loc && equation[i].token.constant >= 0.0)
                              neg_flag = false;
                        d = fabs(equation[i].token.constant);
                        if (first) {
                              minimum = d;
                              first = false;
                        } else if (d < minimum) {
                              minimum = d;
                        }
                        break;
                  case OPERATOR:
                        count++;
                        switch (equation[i].token.operatr) {
                        case PLUS:
                              neg_flag = false;
                        case MINUS:
                              break;
                        default:
                              return modified;
                        }
                        break;
                  default:
                        if (i == loc)
                              neg_flag = false;
                        if (first) {
                              minimum = 1.0;
                              first = false;
                        } else if (1.0 < minimum)
                              minimum = 1.0;
                        break;
                  }
            } else {
                  op = 0;
                  for (j = i + 1;; j += 2) {
                        if (j >= *np || equation[j].level <= level)
                              break;
                        if (equation[j].level == level + 1) {
                              op = equation[j].token.operatr;
                        }
                  }
                  if (op == TIMES || op == DIVIDE) {
                        for (k = i; k < j; k++) {
                              if (equation[k].level == (level + 1)
                                  && equation[k].kind == CONSTANT) {
                                    if (i == loc && equation[k].token.constant >= 0.0)
                                          neg_flag = false;
                                    d = fabs(equation[k].token.constant);
                                    if (first) {
                                          minimum = d;
                                          first = false;
                                    } else if (d < minimum) {
                                          minimum = d;
                                    }
                                    i = j;
                                    goto break_cont;
                              }
                        }
                  }
                  if (i == loc)
                        neg_flag = false;
                  if (first) {
                        minimum = 1.0;
                        first = false;
                  } else if (1.0 < minimum)
                        minimum = 1.0;
                  i = j;
                  continue;
            }
            i++;
      }
      if (first || count == 0 || (!neg_flag && minimum == 1.0))
            return modified;
      if (!isfinite(minimum))
            return modified;
      if (level_code > 1 || (level_code && (level == 1))) {
            for (i = loc;;) {
                  d = 1.0;
                  if (equation[i].kind == CONSTANT) {
                        if (equation[i].level == level) {
                              d = equation[i].token.constant;
                        } else if ((i + 1) < *np
                            && equation[i+1].level == (level + 1)
                            && (equation[i+1].token.operatr == TIMES
                            || equation[i+1].token.operatr == DIVIDE)) {
                              d = equation[i].token.constant;
                        }
                  }
                  if ((minimum < 1.0 && fmod(d, 1.0) == 0.0)
                      || (fmod(d / minimum, 1.0) != 0.0)) {
                        if (neg_flag) {
                              minimum = 1.0;
                              break;
                        }
                        return modified;
                  }
                  i++;
                  for (; i < *np && equation[i].level > level; i += 2)
                        ;
                  if (i >= *np || equation[i].level < level)
                        break;
                  i++;
            }
      }
      if (neg_flag)
            minimum = -minimum;
      if (*np + ((count + 2) * 2) > n_tokens) {
            error_huge();
      }
      for (i = loc;;) {
            if (i >= *np || equation[i].level < level)
                  break;
            if (equation[i].kind != OPERATOR) {
                  for (j = i;;) {
                        equation[j].level++;
                        j++;
                        if (j >= *np || equation[j].level <= level)
                              break;
                  }
                  blt(&equation[j+2], &equation[j], (*np - j) * sizeof(*equation));
                  *np += 2;
                  equation[j].level = level + 1;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = DIVIDE;
                  j++;
                  equation[j].level = level + 1;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = minimum;
                  i = j;
            }
            i++;
      }
      for (i = loc;; i++) {
            if (i >= *np || equation[i].level < level)
                  break;
            equation[i].level++;
      }
      blt(&equation[i+2], &equation[i], (*np - i) * sizeof(*equation));
      *np += 2;
      equation[i].level = level;
      equation[i].kind = OPERATOR;
      equation[i].token.operatr = TIMES;
      i++;
      equation[i].level = level;
      equation[i].kind = CONSTANT;
      equation[i].token.constant = minimum;
      return true;
}

Generated by  Doxygen 1.6.0   Back to index