## File: complex.c

package info (click to toggle)
mathomatic 14.0.6-2
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304` ``````/* * Floating point complex number routines specifically for Mathomatic. * * Copyright (C) 1987-2008 George Gesslein II. */ #include "includes.h" /* * Convert doubles x and y from rectangular coordinates to polar coordinates. * * The amplitude is stored in *radiusp and the angle in radians is stored in *thetap. */ void rect_to_polar(x, y, radiusp, thetap) double x, y, *radiusp, *thetap; { *radiusp = sqrt(x * x + y * y); *thetap = atan2(y, x); } #if !LIBRARY /* * The roots command. */ int roots_cmd(cp) char *cp; { #define MAX_ROOT 10000.0 /* root limit needed because more roots become more inaccurate and take longer to calculate */ complexs c, c2, check; double d, k; double root; double radius, theta; double radius_root = 0.0; char buf[MAX_CMD_LEN]; if (*cp == '\0') { my_strlcpy(prompt_str, _("Enter root (positive integer): "), sizeof(prompt_str)); if ((cp = get_string(buf, sizeof(buf))) == NULL) return false; } root = strtod(cp, &cp); if ((*cp && !isspace(*cp)) || root < 0.0 || root > MAX_ROOT || fmod(root, 1.0) != 0.0) { printf(_("Root must be a positive integer less than or equal to %.0f.\n"), MAX_ROOT); return false; } cp = skip_space(cp); if (*cp == '\0') { my_strlcpy(prompt_str, _("Enter real part (X): "), sizeof(prompt_str)); if ((cp = get_string(buf, sizeof(buf))) == NULL) return false; } c.re = strtod(cp, &cp); if (*cp && !isspace(*cp)) { printf(_("Invalid number.\n")); return false; } cp = skip_space(cp); if (*cp == '\0') { my_strlcpy(prompt_str, _("Enter imaginary part (Y): "), sizeof(prompt_str)); if ((cp = get_string(buf, sizeof(buf))) == NULL) return false; } c.im = strtod(cp, &cp); if (*cp) { printf(_("Invalid number.\n")); return false; } if (c.re == 0.0 && c.im == 0.0) { return false; } /* convert to polar coordinates */ errno = 0; rect_to_polar(c.re, c.im, &radius, &theta); if (root) { radius_root = pow(radius, 1.0 / root); } check_err(); fprintf(gfp, _("\nThe polar coordinates before root taking are:\n%.12g amplitude and %.12g radians (%.12g degrees).\n\n"), radius, theta, theta * 180.0 / M_PI); if (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) { /* add constants to theta and convert back to rectangular coordinates */ c2.re = radius_root * cos((theta + 2.0 * k * M_PI) / root); c2.im = radius_root * sin((theta + 2.0 * k * M_PI) / root); complex_fixup(&c2); 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); } complex_fixup(&check); 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; } #endif /* * Approximate roots of complex numbers: * (complex^real) and (real^complex) and (complex^complex). * Only gives one root, when there may be many. * * Return true if expression was modified. */ int complex_root_simp(equation, np) token_type *equation; /* equation side pointer */ int *np; /* pointer to length of equation side */ { int i, j; int level; int len; complexs c, p; int modified = false; start_over: 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); if (!parse_complex(&equation[i+1], len, &p)) continue; for (j = i - 1; j >= 0 && equation[j].level >= level; j--) ; j++; if (!parse_complex(&equation[j], i - j, &c)) continue; if (c.im == 0.0 && p.im == 0.0) continue; i += len + 1; 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++; level++; equation[j].level = level; equation[j].kind = CONSTANT; equation[j].token.constant = c.im; j++; equation[j].level = level; equation[j].kind = OPERATOR; equation[j].token.operatr = TIMES; j++; equation[j].level = level; equation[j].kind = VARIABLE; equation[j].token.variable = IMAGINARY; modified = true; goto start_over; } if (modified) { warning(_("Complex number root approximated.")); } return modified; } /* * Get a constant, if the passed expression is a constant. * This code only handles the case of a single constant, where "n" == 1. * I don't know how to easily fix this, should at least handle "n" == 3, * which would allow 2 constants separated by any operator. * * Return true if successful, with number in "*dp". */ int get_constant(p1, n, dp) token_type *p1; /* expression pointer */ int n; /* length of expression */ double *dp; /* pointer to returned double */ { if (n != 1) return false; switch (p1[0].kind) { case CONSTANT: *dp = p1[0].token.constant; return true; case VARIABLE: if (var_is_const(p1[0].token.variable, dp)) { return true; } } return false; } /* * Parse a complex number expression. * * If successful return true with complex number in "*cp". */ int parse_complex(p1, n, cp) token_type *p1; /* expression pointer */ int n; /* length of expression */ complexs *cp; /* pointer to returned complex number */ { int j; int imag_cnt = 0, plus_cnt = 0, times_cnt = 0; complexs c; int level2; double junk; if (get_constant(p1, n, &c.re)) { c.im = 0.0; *cp = c; return true; } c.re = 0.0; c.im = 1.0; for (j = n - 1; j >= 0; j--) { switch (p1[j].kind) { case CONSTANT: break; case VARIABLE: if (var_is_const(p1[j].token.variable, &junk)) break; if (p1[j].token.variable != IMAGINARY) return false; imag_cnt++; break; case OPERATOR: level2 = p1[j].level; switch (p1[j].token.operatr) { case TIMES: if (++times_cnt > 1) return false; if (p1[j-1].level != level2 || p1[j+1].level != level2) return false; if (p1[j-1].kind == VARIABLE && p1[j-1].token.variable == IMAGINARY) { if (!get_constant(&p1[j+1], 1, &c.im)) return false; continue; } if (p1[j+1].kind == VARIABLE && p1[j+1].token.variable == IMAGINARY) { if (!get_constant(&p1[j-1], 1, &c.im)) return false; continue; } return false; case DIVIDE: if (++times_cnt > 1) return false; if (p1[j-1].level != level2 || p1[j+1].level != level2) return false; if (p1[j-1].kind == VARIABLE && p1[j-1].token.variable == IMAGINARY) { if (!get_constant(&p1[j+1], 1, &c.im)) return false; c.im = 1.0 / c.im; continue; } return false; case MINUS: if (imag_cnt) { c.im = -c.im; } case PLUS: if (++plus_cnt > 1) return false; if (p1[j-1].level == level2 && get_constant(&p1[j-1], 1, &c.re)) { continue; } if (p1[j+1].level == level2 && get_constant(&p1[j+1], 1, &c.re)) { if (p1[j].token.operatr == MINUS) c.re = -c.re; continue; } } default: return false; } } if (imag_cnt != 1) return false; *cp = c; return true; } ``````