Logo Search packages:      
Sourcecode: mathomatic version File versions

integrate.c

/*
 * Algebraic manipulator integration routines.
 *
 * Copyright (c) 1996 George Gesslein II.
 */

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

static int  poly_in_v_sub();
static int  int_sub();
static int  laplace_sub();
static int  inv_laplace_sub();

/*
 * Return true if passed expression is a polynomial in
 * variable v.
 */
int
poly_in_v(equation, n, v, allow_divides)
token_type  *equation;
int         n;
long        v;
int         allow_divides;
{
      int   i, j, k;
      int   rv;

      rv = true;
      if (level1_plus(equation, n)) {
            for (i = 1, j = 0;; i += 2) {
                  if (i >= n || equation[i].level == 1) {
                        if (!poly_in_v_sub(&equation[j], i - j, v, allow_divides)) {
                              rv = false;
                              break;
                        }
                        j = i + 1;
                  }
                  if (i >= n)
                        break;
            }
      } else {
            if (!poly_in_v_sub(equation, n, v, allow_divides))
                  rv = false;
      }
      return rv;
}

/*
 * Return true if passed expression is strictly a single
 * polynomial term in v.
 */
static int
poly_in_v_sub(equation, n, v, allow_divides)
token_type  *equation;
int         n;
long        v;
int         allow_divides;
{
      int   i, j, k;
      int   level;
      int   count;

      level = min_level(equation, n);
      for (i = 0, count = 0; i < n; i++) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
                  if (count > 1)
                        return false;
                  if (equation[i].level == level || equation[i].level == (level + 1)) {
                        for (k = 1; k < n; k += 2) {
                              if (equation[k].level == level) {
                                    switch (equation[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 (equation[i].level == (level + 1)) {
                              if ((i + 1) < n && equation[i+1].level == equation[i].level
                                  && equation[i+1].token.operatr == POWER) {
                                    continue;
                              }
                        } else {
                              continue;
                        }
                  }
                  return false;
            }
      }
      return true;
}

/*
 * Make variable v always raised to a power.
 */
make_powers(equation, np, v)
token_type  *equation;
int         *np;
long        v;
{
      int   i;
      int   level;

      for (i = 0; i < *np; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  level = equation[i].level;
                  if ((i + 1) >= *np || equation[i+1].token.operatr != POWER) {
                        if (i > 0 && equation[i-1].token.operatr == POWER)
                              continue;
                        if (*np + 2 > n_tokens) {
                              error_huge();
                        }
                        level++;
                        equation[i].level = level;
                        i++;
                        blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
                        *np += 2;
                        equation[i].level = level;
                        equation[i].kind = OPERATOR;
                        equation[i].token.operatr = POWER;
                        i++;
                        equation[i].level = level;
                        equation[i].kind = CONSTANT;
                        equation[i].token.constant = 1.0;
                  }
            }
      }
}

/*
 * Integrate a polynomial expression with respect to v.
 *
 * Return true if successful.
 */
int
integrate(equation, np, v)
token_type  *equation;
int         *np;
long        v;
{
      int   i, j;

      make_powers(equation, np, v);
      if (level1_plus(equation, *np)) {
            for (j = 0, i = 1;; i += 2) {
                  if (i >= *np || equation[i].level == 1) {
                        if (!int_sub(equation, np, j, i, v))
                              return false;
                        for (i = j + 1;; i += 2) {
                              if (i >= *np) {
                                    return true;
                              }
                              if (equation[i].level == 1) {
                                    j = i + 1;
                                    break;
                              }
                        }
                  }
            }
      } else {
            if (!int_sub(equation, np, 0, *np, v))
                  return false;
      }
      return true;
}

/*
 * Do the actual integration.
 *
 * Return true if successful.
 */
static int
int_sub(equation, np, loc, eloc, v)
token_type  *equation;
int         *np;
int         loc, eloc;
long        v;
{
      int   i, j, k;
      int   len;
      int   level;
      int   mlevel;
      int   count;
      int   exp_pos;
      complexs    c;
      int   div_flag;

      exp_pos = -1;
      level = min_level(&equation[loc], eloc - loc);
      for (i = loc, count = 0; i < eloc; i++) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
                  if (count > 1)
                        return false;
                  if (i > loc && equation[i].level == equation[i-1].level
                      && equation[i-1].token.operatr == POWER) {
                        if (equation[i-2].level != equation[i].level) {
                              return false;
                        }
                        if (equation[i].level == level) {
                              exp_pos = i - 2;
                              continue;
                        } else if (equation[i].level == level + 1) {
                              for (k = loc + 1; k < eloc; k += 2) {
                                    if (equation[k].level == level) {
                                          switch (equation[k].token.operatr) {
                                          case TIMES:
                                          case DIVIDE:
                                                break;
                                          default:
                                                return false;
                                          }
                                    }
                              }
                              exp_pos = i - 2;
                              continue;
                        }
                        return false;
                  }
                  if (equation[i].level == level || equation[i].level == (level + 1)) {
                        for (k = loc + 1; k < eloc; k += 2) {
                              if (equation[k].level == level) {
                                    switch (equation[k].token.operatr) {
                                    case DIVIDE:
                                    case TIMES:
                                          continue;
                                    case POWER:
                                          if (k == (i + 1))
                                                continue;
                                    default:
                                          return false;
                                    }
                              }
                        }
                        if (equation[i].level == (level + 1)) {
                              if ((i + 1) < eloc && equation[i+1].level == equation[i].level
                                  && equation[i+1].token.operatr == POWER) {
                                    continue;
                              }
                        } else {
                              continue;
                        }
                  }
                  return false;
            }
      }
      if (exp_pos >= 0) {
            c.im = 0.0;
            switch (equation[exp_pos].kind) {
            case CONSTANT:
                  c.re = equation[exp_pos].token.constant;
                  break;
            case VARIABLE:
                  if (var_is_const(equation[exp_pos].token.variable, &c.re))
                        break;
            default:
                  return false;
            }
            c = complex_log(c);
            if (exp_pos > loc && equation[exp_pos-1].token.operatr == DIVIDE) {
                  c.re = -c.re;
                  c.im = -c.im;
            }
            if (*np + 6 > n_tokens) {
                  error_huge();
            }
            for (j = loc; j < eloc; j++)
                  equation[j].level++;
            blt(&equation[eloc+6], &equation[eloc], (*np - eloc) * sizeof(token_type));
            *np += 6;
            k = eloc;
            equation[k].level = level;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = DIVIDE;
            k++;
            equation[k].level = level + 1;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.re;
            k++;
            equation[k].level = level + 1;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = PLUS;
            k++;
            equation[k].level = level + 2;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.im;
            k++;
            equation[k].level = level + 2;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = TIMES;
            k++;
            equation[k].level = level + 2;
            equation[k].kind = VARIABLE;
            equation[k].token.variable = IMAGINARY;
            return true;
      }
      mlevel = level + 1;
      for (j = loc; j < eloc; j++)
            equation[j].level += 2;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  div_flag = (i > loc && equation[i-1].token.operatr == DIVIDE);
                  i++;
                  if (i >= eloc || equation[i].token.operatr != POWER)
                        return false;
                  level = equation[i].level;
                  i++;
                  if (div_flag) {
                        if (equation[i].level == level
                            && equation[i].kind == CONSTANT
                            && equation[i].token.constant == 1.0)
                              return false;
                        if (*np + 2 > n_tokens)
                              error_huge();
                        for (j = i; j < eloc && equation[j].level >= level; j++)
                              equation[j].level++;
                        equation[i-3].token.operatr = TIMES;
                        blt(&equation[i+2], &equation[i], (*np - i) * sizeof(token_type));
                        *np += 2;
                        eloc += 2;
                        equation[i].level = level + 1;
                        equation[i].kind = CONSTANT;
                        equation[i].token.constant = -1.0;
                        equation[i+1].level = level + 1;
                        equation[i+1].kind = OPERATOR;
                        equation[i+1].token.operatr = TIMES;
                  }
                  for (j = i; j < eloc && equation[j].level >= level; j++)
                        equation[j].level++;
                  len = j - i;
                  if (*np + len + 5 > n_tokens)
                        error_huge();
                  blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type));
                  *np += 2;
                  eloc += 2;
                  len += 2;
                  level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = PLUS;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = 1.0;
                  j++;
                  blt(&equation[eloc+len+1], &equation[eloc], (*np - eloc) * sizeof(token_type));
                  *np += len + 1;
                  k = eloc;
                  equation[k].level = mlevel;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = DIVIDE;
                  k++;
                  blt(&equation[k], &equation[i], len * sizeof(token_type));
                  return true;
            }
      }
      if (*np + 2 > n_tokens) {
            error_huge();
      }
      blt(&equation[eloc+2], &equation[eloc], (*np - eloc) * sizeof(token_type));
      *np += 2;
      k = eloc;
      equation[k].level = mlevel;
      equation[k].kind = OPERATOR;
      equation[k].token.operatr = TIMES;
      k++;
      equation[k].level = mlevel;
      equation[k].kind = VARIABLE;
      equation[k].token.variable = v;
      return true;
}

/*
 * The integrate command.
 */
int
integrate_cmd(cp)
char  *cp;
{
      long  v;
      int   i, j;
      long  l;
      int   er;
      int   order;

      v = 0;
      order = 1;
      if (notdefined(cur_equation)) {
            return false;
      }
      if (*cp) {
            cp = parse_var2(&v, cp);
            if (cp == NULL) {
                  return false;
            }
            if (*cp) {
                  order = atoi(cp);
                  cp = skip_num(cp);
            }
      }
      if (extra_garbage(cp)) {
            return false;
      }
      if (order < 1) {
            printf(_("Order out of range.\n"));
            return false;
      }
      if (v == 0) {
            if (!prompt_var(&v)) {
                  return false;
            }
      }
      i = next_espace();
      er = !solved_equation(cur_equation);
      if (!er) {
            er = ((lhs[cur_equation][0].token.variable & VAR_MASK) <= SIGN);
      }
      partial_flag = false;
      uf_simp(rhs[cur_equation], &n_rhs[cur_equation]);
      partial_flag = true;
      factorv(rhs[cur_equation], &n_rhs[cur_equation], v);
      blt(rhs[i], rhs[cur_equation], n_rhs[cur_equation] * sizeof(token_type));
      n_rhs[i] = n_rhs[cur_equation];
      for (j = 0; j < order; j++) {
            if (!integrate(rhs[i], &n_rhs[i], v)) {
                  printf(_("Integration failed.\n"));
                  return false;
            }
            simp_loop(rhs[i], &n_rhs[i]);
      }
      blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
      n_lhs[i] = n_lhs[cur_equation];
      if (!er) {
/* decrement the number of primes in the LHS variable */
            l = lhs[i][0].token.variable;
            for (j = 1; j <= order; j++) {
                  l -= PRIME_INCREMENT;
                  if (l < 0) {
                        l += PRIME_INCREMENT;
                        break;
                  }
            }
            lhs[i][0].token.variable = l;
      }
      cur_equation = i;
      list_sub(cur_equation);
      return true;
}

/*
 * Compute the Laplace transform of a polynomial expression with respect to v.
 *
 * Return true if successful.
 */
int
laplace(equation, np, v)
token_type  *equation;
int         *np;
long        v;
{
      int   i, j;

      make_powers(equation, np, v);
      if (level1_plus(equation, *np)) {
            for (j = 0, i = 1;; i += 2) {
                  if (i >= *np || equation[i].level == 1) {
                        if (!laplace_sub(equation, np, j, i, v))
                              return false;
                        for (i = j + 1;; i += 2) {
                              if (i >= *np) {
                                    return true;
                              }
                              if (equation[i].level == 1) {
                                    j = i + 1;
                                    break;
                              }
                        }
                  }
            }
      } else {
            if (!laplace_sub(equation, np, 0, *np, v))
                  return false;
      }
      return true;
}

static int
laplace_sub(equation, np, loc, eloc, v)
token_type  *equation;
int         *np;
int         loc, eloc;
long        v;
{
      int   i, j, k;
      int   len;
      int   level, level2;
      int   mlevel;
      int   count;
      int   exp_pos;
      complexs    c;

      exp_pos = -1;
      level = min_level(&equation[loc], eloc - loc);
      for (i = loc, count = 0; i < eloc; i++) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
                  if (count > 1)
                        return false;
                  if (i > loc && equation[i].level == equation[i-1].level
                      && equation[i-1].token.operatr == POWER) {
                        if (equation[i-2].level != equation[i].level) {
                              return false;
                        }
                        if (equation[i].level == level) {
                              exp_pos = i - 2;
                              continue;
                        } else if (equation[i].level == level + 1) {
                              for (k = loc + 1; k < eloc; k += 2) {
                                    if (equation[k].level == level) {
                                          switch (equation[k].token.operatr) {
                                          case TIMES:
                                          case DIVIDE:
                                                break;
                                          default:
                                                return false;
                                          }
                                    }
                              }
                              exp_pos = i - 2;
                              continue;
                        }
                        return false;
                  }
                  if (equation[i].level == level || equation[i].level == (level + 1)) {
                        for (k = loc + 1; k < eloc; k += 2) {
                              if (equation[k].level == level) {
                                    switch (equation[k].token.operatr) {
                                    case DIVIDE:
                                          if (k == (i - 1))
                                                return false;
                                    case TIMES:
                                          continue;
                                    case POWER:
                                          if (k == (i + 1))
                                                continue;
                                    default:
                                          return false;
                                    }
                              }
                        }
                        if (equation[i].level == (level + 1)) {
                              if ((i + 1) < eloc && equation[i+1].level == equation[i].level
                                  && equation[i+1].token.operatr == POWER) {
                                    continue;
                              }
                        } else {
                              continue;
                        }
                  }
                  return false;
            }
      }
      if (exp_pos >= 0) {
            if (*np + 6 > n_tokens) {
                  error_huge();
            }
            c.im = 0.0;
            switch (equation[exp_pos].kind) {
            case CONSTANT:
                  c.re = equation[exp_pos].token.constant;
                  break;
            case VARIABLE:
                  if (var_is_const(equation[exp_pos].token.variable, &c.re))
                        break;
            default:
                  return false;
            }
            c = complex_log(c);
            if (exp_pos > loc && equation[exp_pos-1].token.operatr == DIVIDE) {
                  equation[exp_pos-1].token.operatr = TIMES;
                  c.re = -c.re;
                  c.im = -c.im;
            }
            equation[exp_pos].kind = CONSTANT;
            equation[exp_pos].token.constant = 1.0;
            exp_pos++;
            equation[exp_pos].token.operatr = DIVIDE;
            exp_pos++;
            equation[exp_pos].level++;
            level2 = equation[exp_pos].level;
            exp_pos++;
            blt(&equation[exp_pos+6], &equation[exp_pos], (*np - exp_pos) * sizeof(token_type));
            *np += 6;
            k = exp_pos;
            equation[k].level = level2;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = MINUS;
            k++;
            equation[k].level = level2 + 1;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.re;
            k++;
            equation[k].level = level2 + 1;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = PLUS;
            k++;
            equation[k].level = level2 + 2;
            equation[k].kind = CONSTANT;
            equation[k].token.constant = c.im;
            k++;
            equation[k].level = level2 + 2;
            equation[k].kind = OPERATOR;
            equation[k].token.operatr = TIMES;
            k++;
            equation[k].level = level2 + 2;
            equation[k].kind = VARIABLE;
            equation[k].token.variable = IMAGINARY;
            return true;
      }
      mlevel = level + 1;
      for (j = loc; j < eloc; j++)
            equation[j].level += 2;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  i++;
                  if (i >= eloc || equation[i].token.operatr != POWER)
                        return false;
                  level = equation[i].level;
                  i++;
                  for (j = i; j < eloc && equation[j].level >= level; j++)
                        equation[j].level++;
                  len = j - i;
                  if (*np + len + 7 > n_tokens)
                        error_huge();
                  blt(&equation[j+4], &equation[j], (*np - j) * sizeof(token_type));
                  *np += 4;
                  eloc += 4;
                  level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = PLUS;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = 1.0;
                  j++;
                  for (k = i; k < j; k++)
                        equation[k].level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = TIMES;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = -1.0;
                  j++;
                  blt(&equation[eloc+len+3], &equation[eloc], (*np - eloc) * sizeof(token_type));
                  *np += len + 3;
                  k = eloc;
                  equation[k].level = mlevel;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = TIMES;
                  k++;
                  blt(&equation[k], &equation[i], len * sizeof(token_type));
                  k += len;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = FACTORIAL;
                  k++;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = CONSTANT;
                  equation[k].token.constant = 0.0;
                  k++;
                  return true;
            }
      }
      if (*np + 2 > n_tokens) {
            error_huge();
      }
      blt(&equation[eloc+2], &equation[eloc], (*np - eloc) * sizeof(token_type));
      *np += 2;
      k = eloc;
      equation[k].level = mlevel;
      equation[k].kind = OPERATOR;
      equation[k].token.operatr = DIVIDE;
      k++;
      equation[k].level = mlevel;
      equation[k].kind = VARIABLE;
      equation[k].token.variable = v;
      return true;
}

/*
 * Compute the inverse Laplace transform of a polynomial expression with respect to v.
 *
 * Return true if successful.
 */
int
inv_laplace(equation, np, v)
token_type  *equation;
int         *np;
long        v;
{
      int   i, j;

      make_powers(equation, np, v);
      if (level1_plus(equation, *np)) {
            for (j = 0, i = 1;; i += 2) {
                  if (i >= *np || equation[i].level == 1) {
                        if (!inv_laplace_sub(equation, np, j, i, v))
                              return false;
                        for (i = j + 1;; i += 2) {
                              if (i >= *np) {
                                    return true;
                              }
                              if (equation[i].level == 1) {
                                    j = i + 1;
                                    break;
                              }
                        }
                  }
            }
      } else {
            if (!inv_laplace_sub(equation, np, 0, *np, v))
                  return false;
      }
      return true;
}

static int
inv_laplace_sub(equation, np, loc, eloc, v)
token_type  *equation;
int         *np;
int         loc, eloc;
long        v;
{
      int   i, j, k;
      int   len;
      int   level;
      int   mlevel;
      int   count;

      count = 0;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  count++;
            }
      }
      if (count > 1)
            return false;
      mlevel = min_level(&equation[loc], eloc - loc) + 1;
      for (j = loc; j < eloc; j++)
            equation[j].level += 2;
      for (i = loc; i < eloc; i += 2) {
            if (equation[i].kind == VARIABLE && equation[i].token.variable == v) {
                  i++;
                  if (i >= eloc || equation[i].token.operatr != POWER)
                        return false;
                  if ((i - 2) <= loc || equation[i-2].token.operatr != DIVIDE)
                        return false;
                  level = equation[i].level;
                  i++;
                  for (j = i; j < eloc && equation[j].level >= level; j++)
                        equation[j].level++;
                  len = j - i;
                  if (*np + len + 7 > n_tokens)
                        error_huge();
                  equation[i-3].token.operatr = TIMES;
                  blt(&equation[j+2], &equation[j], (*np - j) * sizeof(token_type));
                  *np += 2;
                  eloc += 2;
                  len += 2;
                  level++;
                  equation[j].level = level;
                  equation[j].kind = OPERATOR;
                  equation[j].token.operatr = MINUS;
                  j++;
                  equation[j].level = level;
                  equation[j].kind = CONSTANT;
                  equation[j].token.constant = 1.0;
                  j++;
                  blt(&equation[eloc+len+3], &equation[eloc], (*np - eloc) * sizeof(token_type));
                  *np += len + 3;
                  k = eloc;
                  equation[k].level = mlevel;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = DIVIDE;
                  k++;
                  blt(&equation[k], &equation[i], len * sizeof(token_type));
                  k += len;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = OPERATOR;
                  equation[k].token.operatr = FACTORIAL;
                  k++;
                  equation[k].level = mlevel + 1;
                  equation[k].kind = CONSTANT;
                  equation[k].token.constant = 0.0;
                  k++;
                  return true;
            }
      }
      return false;
}

/*
 * The laplace command.
 */
int
laplace_cmd(cp)
char  *cp;
{
      long  v;
      int   i, j;
      int   inverse_flag;

      v = 0;
      if (notdefined(cur_equation)) {
            return false;
      }
      inverse_flag = (strcmp_tospace(cp, "inverse") == 0);
      if (inverse_flag) {
            cp = skip_param(cp);
      }
      if (*cp) {
            cp = parse_var2(&v, cp);
            if (cp == NULL) {
                  return false;
            }
      }
      if (extra_garbage(cp)) {
            return false;
      }
      if (v == 0) {
            if (!prompt_var(&v)) {
                  return false;
            }
      }
      i = next_espace();
      partial_flag = false;
      uf_simp(rhs[cur_equation], &n_rhs[cur_equation]);
      partial_flag = true;
      factorv(rhs[cur_equation], &n_rhs[cur_equation], v);
      blt(rhs[i], rhs[cur_equation], n_rhs[cur_equation] * sizeof(token_type));
      n_rhs[i] = n_rhs[cur_equation];
      if (inverse_flag) {
            if (!poly_in_v(rhs[i], n_rhs[i], v, true)) {
                  printf(_("Inverse Laplace failed.\n"));
                  return false;
            }
            if (!inv_laplace(rhs[i], &n_rhs[i], v)) {
                  printf(_("Inverse Laplace failed.\n"));
                  return false;
            }
      } else {
            if (!laplace(rhs[i], &n_rhs[i], v)) {
                  printf(_("Laplace failed.\n"));
                  return false;
            }
      }
      simp_loop(rhs[i], &n_rhs[i]);
      blt(lhs[i], lhs[cur_equation], n_lhs[cur_equation] * sizeof(token_type));
      n_lhs[i] = n_lhs[cur_equation];
      cur_equation = i;
      list_sub(cur_equation);
      return true;
}

Generated by  Doxygen 1.6.0   Back to index