complex.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #include "common.h"
00028
00029 #include <stdio.h>
00030 #include <math.h>
00031 #include "bu.h"
00032 #include "vmath.h"
00033 #include "bn.h"
00034
00035 void
00036 bn_cx_div(register bn_complex_t *ap, register const bn_complex_t *bp)
00037 {
00038 register fastf_t r, s;
00039 register fastf_t ap__re;
00040
00041
00042 ap__re = ap->re;
00043 r = bp->re;
00044 s = bp->im;
00045 if (fabs(r) >= fabs(s)) {
00046 if (ZERO(r))
00047 goto err;
00048 r = s / r;
00049 s = 1.0 / (bp->re + r * s);
00050 ap->re = (ap->re + ap->im * r) * s;
00051 ap->im = (ap->im - ap__re * r) * s;
00052 return;
00053 } else {
00054 if (ZERO(s))
00055 goto err;
00056 r = r / s;
00057 s = 1.0 / (s + r * bp->re);
00058 ap->re = (ap->re * r + ap->im) * s;
00059 ap->im = (ap->im * r - ap__re) * s;
00060 return;
00061 }
00062 err:
00063 bu_log("bn_cx_div: division by zero: %gR+%gI / %gR+%gI\n",
00064 ap->re, ap->im, bp->re, bp->im);
00065 ap->re = ap->im = 1.0e20;
00066 }
00067
00068 void
00069 bn_cx_sqrt(bn_complex_t *op, const bn_complex_t *ip)
00070 {
00071 const fastf_t re = ip->re;
00072 const fastf_t im = ip->im;
00073
00074
00075 if (ZERO(re)) {
00076 if (ZERO(im)) {
00077 op->re = op->im = 0.0;
00078 } else if (im > 0.0) {
00079 op->re = op->im = sqrt(im * 0.5);
00080 } else {
00081
00082 op->re = -(op->im = sqrt(im * -0.5));
00083 }
00084 } else if (ZERO(im)) {
00085 if (re > 0.0) {
00086 op->re = sqrt(re);
00087 op->im = 0.0;
00088 } else {
00089
00090 op->im = sqrt(-re);
00091 op->re = 0.0;
00092 }
00093 } else {
00094 fastf_t ampl, temp;
00095
00096
00097 ampl = bn_cx_ampl(ip);
00098 if ((temp = (ampl - re) * 0.5) < 0.0) {
00099
00100
00101
00102
00103
00104 op->im = 0.0;
00105 } else {
00106 op->im = sqrt(temp);
00107 }
00108
00109 if ((temp = (ampl + re) * 0.5) < 0.0) {
00110 op->re = 0.0;
00111 } else {
00112 if (im > 0.0) {
00113 op->re = sqrt(temp);
00114 } else {
00115
00116 op->re = -sqrt(temp);
00117 }
00118 }
00119 }
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132