complex.c

Go to the documentation of this file.
00001 /*                       C O M P L E X . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1985-2012 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * version 2.1 as published by the Free Software Foundation.
00010  *
00011  * This library is distributed in the hope that it will be useful, but
00012  * WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014  * Lesser General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU Lesser General Public
00017  * License along with this file; see the file named COPYING for more
00018  * information.
00019  */
00020 /** @addtogroup complex */
00021 /** @{ */
00022 /** @file libbn/complex.c
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     /* Note: classical formula may cause unnecessary overflow */
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;                      /* <= 1 */
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;                      /* < 1 */
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;           /* "INFINITY" */
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     /* special cases are not necessary; they are here for speed */
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             /* ip->im < 0.0 */
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             /* ip->re < 0.0 */
00090             op->im = sqrt(-re);
00091             op->re = 0.0;
00092         }
00093     } else {
00094         fastf_t ampl, temp;
00095 
00096         /* no shortcuts */
00097         ampl = bn_cx_ampl(ip);
00098         if ((temp = (ampl - re) * 0.5) < 0.0) {
00099             /* This case happens rather often, when the hypot() in
00100              * bn_cx_ampl() returns an ampl ever so slightly smaller
00101              * than ip->re.  This can easily happen when ip->re ~=
00102              * 10**20.  Just ignore the imaginary part.
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                 /* ip->im < 0.0 */
00116                 op->re = -sqrt(temp);
00117             }
00118         }
00119     }
00120 }
00121 
00122 /** @} */
00123 
00124 /*
00125  * Local Variables:
00126  * mode: C
00127  * tab-width: 8
00128  * indent-tabs-mode: t
00129  * c-file-style: "stroustrup"
00130  * End:
00131  * ex: shiftwidth=4 tabstop=8
00132  */
Generated on Tue Dec 11 13:14:27 2012 for LIBBN by  doxygen 1.6.3