complex.c
Go to the documentation of this file.
1 /* C O M P L E X . C
3  *
4  * Copyright (c) 1985-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20
21 #include "common.h"
22
23 #include <stdio.h>
24 #include <math.h>
25
26 #include "bu/log.h"
27 #include "vmath.h"
28 #include "bn/complex.h"
29
30
31 void
32 bn_cx_div(register bn_complex_t *ap, register const bn_complex_t *bp)
33 {
34  register fastf_t r, s;
35  register fastf_t ap__re;
36
37  /* Note: classical formula may cause unnecessary overflow */
38  ap__re = ap->re;
39  r = bp->re;
40  s = bp->im;
41  if (fabs(r) >= fabs(s)) {
42  if (ZERO(r))
43  goto err;
44  r = s / r; /* <= 1 */
45  s = 1.0 / (bp->re + r * s);
46  ap->re = (ap->re + ap->im * r) * s;
47  ap->im = (ap->im - ap__re * r) * s;
48  return;
49  } else {
50  if (ZERO(s))
51  goto err;
52  r = r / s; /* < 1 */
53  s = 1.0 / (s + r * bp->re);
54  ap->re = (ap->re * r + ap->im) * s;
55  ap->im = (ap->im * r - ap__re) * s;
56  return;
57  }
58 err:
59  bu_log("bn_cx_div: division by zero: %gR+%gI / %gR+%gI\n",
60  ap->re, ap->im, bp->re, bp->im);
61  ap->re = ap->im = 1.0e20; /* "INFINITY" */
62 }
63
64
65 void
67 {
68  const fastf_t re = ip->re;
69  const fastf_t im = ip->im;
70
71  /* special cases are not necessary; they are here for speed */
72  if (ZERO(re)) {
73  if (ZERO(im)) {
74  op->re = op->im = 0.0;
75  } else if (im > 0.0) {
76  op->re = op->im = sqrt(im * 0.5);
77  } else {
78  /* ip->im < 0.0 */
79  op->re = -(op->im = sqrt(im * -0.5));
80  }
81  } else if (ZERO(im)) {
82  if (re > 0.0) {
83  op->re = sqrt(re);
84  op->im = 0.0;
85  } else {
86  /* ip->re < 0.0 */
87  op->im = sqrt(-re);
88  op->re = 0.0;
89  }
90  } else {
91  fastf_t ampl, temp;
92
93  /* no shortcuts */
94  ampl = bn_cx_ampl(ip);
95  if ((temp = (ampl - re) * 0.5) < 0.0) {
96  /* This case happens rather often, when the hypot() in
97  * bn_cx_ampl() returns an ampl ever so slightly smaller
98  * than ip->re. This can easily happen when ip->re ~=
99  * 10**20. Just ignore the imaginary part.
100  */
101  op->im = 0.0;
102  } else {
103  op->im = sqrt(temp);
104  }
105
106  if ((temp = (ampl + re) * 0.5) < 0.0) {
107  op->re = 0.0;
108  } else {
109  if (im > 0.0) {
110  op->re = sqrt(temp);
111  } else {
112  /* ip->im < 0.0 */
113  op->re = -sqrt(temp);
114  }
115  }
116  }
117 }
118
119
120 /*
121  * Local Variables:
122  * mode: C
123  * tab-width: 8
124  * indent-tabs-mode: t
125  * c-file-style: "stroustrup"
126  * End:
127  * ex: shiftwidth=4 tabstop=8
128  */
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
if lu s
Definition: nmg_mod.c:3860
double im
imaginary part
Definition: complex.h:41
double re
real part
Definition: complex.h:40
unsigned char * bp
Definition: rot.c:56
#define bn_cx_ampl(cp)
Definition: complex.h:52
#define ZERO(val)
Definition: units.c:38
Complex numbers.
Definition: complex.h:39
void bn_cx_div(register bn_complex_t *ap, register const bn_complex_t *bp)
Definition: complex.c:32
double fastf_t
Definition: defines.h:300
void bn_cx_sqrt(bn_complex_t *op, const bn_complex_t *ip)
Compute square root of complex number.
Definition: complex.c:66