Logo Search packages:      
Sourcecode: mathomatic version File versions

complex.c

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

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

/*
 * Multiply two complex numbers and return result.
 */
complexs
complex_mult(a, b)
complexs    a;
complexs    b;
{
      complexs    r;

      r.re = a.re * b.re - a.im * b.im;
      r.im = a.re * b.im + a.im * b.re;
      return(r);
}

/*
 * Take the natural logarithm of a complex number and return result.
 */
complexs
complex_log(a)
complexs    a;
{
      complexs    r;

      r.re = log(a.re * a.re + a.im * a.im) / 2.0;
      check_err();
      r.im = atan2(a.im, a.re);
      return(r);
}

/*
 * Raise the natural number (e) to the power of a complex number
 * and return result.
 */
complexs
complex_exp(a)
complexs    a;
{
      complexs    r;
      double            m;

      m = exp(a.re);
      check_err();
      r.re = m * cos(a.im);
      r.im = m * sin(a.im);
      return(r);
}

/*
 * Raise a complex number to the power of a complex number and return result.
 */
complexs
complex_pow(a, b)
complexs    a;
complexs    b;
{
      complexs    r;

      r = complex_log(a);
      r = complex_mult(r, b);
      r = complex_exp(r);
      return(r);
}

int
roots(cp)
char  *cp;
{
      complexs    c, c2, check;
      double            d, k;
      double            root;
      double            radius, theta;
      double            radius_root;
      char        *cp1;

      if (extra_garbage(cp))
            return false;
      for (;;) {
            strcpy(prompt_str, _("Enter root (positive integer): "));
            if ((cp1 = getstring((char *) scratch, n_tokens * sizeof(token_type))) == NULL)
                  return false;
            if (*cp1 == '\0') {
                  return false;
            }
            root = strtod(cp1, &cp1);
            if (*cp1 || root <= 0.0 || fmod(root, 1.0) != 0.0) {
                  printf(_("Root must be positive integer.\n"));
                  continue;
            }
            break;
      }
      do {
            strcpy(prompt_str, _("Enter real part: "));
            if ((cp1 = getstring((char *) scratch, n_tokens * sizeof(token_type))) == NULL)
                  return false;
            c.re = strtod(cp1, &cp1);
      } while (*cp1);
      do {
            strcpy(prompt_str, _("Enter imaginary part: "));
            if ((cp1 = getstring((char *) scratch, n_tokens * sizeof(token_type))) == NULL)
                  return false;
            c.im = strtod(cp1, &cp1);
      } while (*cp1);
/* convert to polar coordinates */
      radius = sqrt(c.re * c.re + c.im * c.im);
      theta = atan2(c.im, c.re);
      radius_root = pow(radius, 1.0 / root);
      if (c.im == 0.0) {
            fprintf(gfp, _("The %.12g roots of %.12g^(1/%.12g) are:\n\n"), root, c.re, root);
      } else {
            fprintf(gfp, _("The %.12g roots of (%.12g %+.12g*i#)^(1/%.12g) are:\n\n"), root, c.re, c.im, root);
      }
      for (k = 0.0; k < root; k += 1.0) {
            c2.re = radius_root * cos((theta + 2.0 * k * PI) / root);
            c2.im = radius_root * sin((theta + 2.0 * k * PI) / root);
            if (c2.im == 0.0) {
                  fprintf(gfp, "%.12g\n", c2.re);
            } else {
                  fprintf(gfp, "%.12g %+.12g*i#\n", c2.re, c2.im);
            }
            check = c2;
            for (d = 1.0; d < root; d += 1.0) {
                  check = complex_mult(check, c2);
            }
            if (check.im == 0.0) {
                  printf(_("Inverse Check: %.12g\n\n"), check.re);
            } else {
                  printf(_("Inverse Check: %.12g %+.12g*i#\n\n"), check.re, check.im);
            }
      }
      return true;
}

/*
 * Approximate roots of complex numbers (complex^real)
 * and real^complex and complex^complex.
 * Only gives one root, so this is not perfect.
 *
 * Returns true if expression was modified.
 */
int
complex_root_simp(equation, np)
token_type  *equation;
int         *np;
{
      int         i, j;
      int         imag_cnt, plus_cnt, times_cnt;
      int         level, level2;
      int         len;
      complexs    c, p;

      for (i = 1; i < *np; i += 2) {
            if (equation[i].token.operatr != POWER)
                  continue;
            level = equation[i].level;
            for (j = i + 2; j < *np && equation[j].level >= level; j += 2)
                  ;
            len = j - (i + 1);
            p.re = 0.0;
            p.im = 1.0;
            imag_cnt = 0;
            plus_cnt = 0;
            times_cnt = 0;
            for (j = i + len; j > i; j--) {
                  switch (equation[j].kind) {
                  case VARIABLE:
                        if (equation[j].token.variable != IMAGINARY)
                              goto end_loop;
                        imag_cnt++;
                        break;
                  case CONSTANT:
                        if (equation[j].level == level) {
                              p.re = equation[j].token.constant;
                              p.im = 0.0;
                              imag_cnt++;
                        }
                        break;
                  case OPERATOR:
                        switch (equation[j].token.operatr) {
                        case TIMES:
                              if (++times_cnt > 1)
                                    goto end_loop;
                              level2 = equation[j].level;
                              if (equation[j-1].level != level2
                                  || equation[j+1].level != level2)
                                    goto end_loop;
                              if (equation[j-1].kind == VARIABLE
                                  && equation[j-1].token.variable == IMAGINARY) {
                                    if (equation[j+1].kind != CONSTANT)
                                          goto end_loop;
                                    p.im = equation[j+1].token.constant;
                                    continue;
                              }
                              if (equation[j+1].kind == VARIABLE
                                  && equation[j+1].token.variable == IMAGINARY) {
                                    if (equation[j-1].kind != CONSTANT)
                                          goto end_loop;
                                    p.im = equation[j-1].token.constant;
                                    continue;
                              }
                              goto end_loop;
                        case MINUS:
                              if (imag_cnt) {
                                    p.im = -p.im;
                              }
                        case PLUS:
                              if (++plus_cnt > 1)
                                    goto end_loop;
                              level2 = equation[j].level;
                              if (equation[j-1].level == level2
                                  && equation[j-1].kind == CONSTANT) {
                                    p.re = equation[j-1].token.constant;
                                    continue;
                              }
                              if (equation[j+1].level == level2
                                  && equation[j+1].kind == CONSTANT) {
                                    p.re = equation[j+1].token.constant;
                                    if (equation[j].token.operatr == MINUS)
                                          p.re = -p.re;
                                    continue;
                              }
                        }
                  default:
                        goto end_loop;
                  }
            }
            if (imag_cnt != 1)
                  goto end_loop;
            c.re = 0.0;
            c.im = 1.0;
            imag_cnt = 0;
            plus_cnt = 0;
            times_cnt = 0;
            for (j = i - 1; j >= 0 && equation[j].level >= level; j--) {
                  switch (equation[j].kind) {
                  case VARIABLE:
                        if (equation[j].token.variable != IMAGINARY)
                              goto end_loop;
                        imag_cnt++;
                        break;
                  case CONSTANT:
                        if (equation[j].level == level) {
                              c.re = equation[j].token.constant;
                              c.im = 0.0;
                              imag_cnt++;
                        }
                        break;
                  case OPERATOR:
                        switch (equation[j].token.operatr) {
                        case TIMES:
                              if (++times_cnt > 1)
                                    goto end_loop;
                              level2 = equation[j].level;
                              if (equation[j-1].level != level2
                                  || equation[j+1].level != level2)
                                    goto end_loop;
                              if (equation[j-1].kind == VARIABLE
                                  && equation[j-1].token.variable == IMAGINARY) {
                                    if (equation[j+1].kind != CONSTANT)
                                          goto end_loop;
                                    c.im = equation[j+1].token.constant;
                                    continue;
                              }
                              if (equation[j+1].kind == VARIABLE
                                  && equation[j+1].token.variable == IMAGINARY) {
                                    if (equation[j-1].kind != CONSTANT)
                                          goto end_loop;
                                    c.im = equation[j-1].token.constant;
                                    continue;
                              }
                              goto end_loop;
                        case MINUS:
                              if (imag_cnt) {
                                    c.im = -c.im;
                              }
                        case PLUS:
                              if (++plus_cnt > 1)
                                    goto end_loop;
                              level2 = equation[j].level;
                              if (equation[j-1].level == level2
                                  && equation[j-1].kind == CONSTANT) {
                                    c.re = equation[j-1].token.constant;
                                    continue;
                              }
                              if (equation[j+1].level == level2
                                  && equation[j+1].kind == CONSTANT) {
                                    c.re = equation[j+1].token.constant;
                                    if (equation[j].token.operatr == MINUS)
                                          c.re = -c.re;
                                    continue;
                              }
                        }
                  default:
                        goto end_loop;
                  }
            }
            if (imag_cnt != 1)
                  goto end_loop;
            if (c.im == 0.0 && p.im == 0.0)
                  goto end_loop;
            j++;
            i += len + 1;
#if   !SILENT
            printf(_("Warning: Complex number root taken.\n"));
#endif
            c = complex_pow(c, p);
            if (*np + 5 - (i - j) > n_tokens) {
                  error_huge();
            }
            blt(&equation[j+5], &equation[i], (*np - i) * sizeof(token_type));
            *np += 5 - (i - j);
            equation[j].level = level;
            equation[j].kind = CONSTANT;
            equation[j].token.constant = c.re;
            j++;
            equation[j].level = level;
            equation[j].kind = OPERATOR;
            equation[j].token.operatr = PLUS;
            j++;
            equation[j].level = level + 1;
            equation[j].kind = CONSTANT;
            equation[j].token.constant = c.im;
            j++;
            equation[j].level = level + 1;
            equation[j].kind = OPERATOR;
            equation[j].token.operatr = TIMES;
            j++;
            equation[j].level = level + 1;
            equation[j].kind = VARIABLE;
            equation[j].token.variable = IMAGINARY;
            return true;
end_loop:
            ;
      }
      return false;
}

Generated by  Doxygen 1.6.0   Back to index