1 /* $OpenBSD: s_csqrt.c,v 1.1 2008/09/07 20:36:09 martynas Exp $ */ 2 /* 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 4 * 5 * Permission to use, copy, modify, and distribute this software for any 6 * purpose with or without fee is hereby granted, provided that the above 7 * copyright notice and this permission notice appear in all copies. 8 * 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 11 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 12 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 13 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 14 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 15 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 16 */ 17 18 /* csqrt() 19 * 20 * Complex square root 21 * 22 * 23 * 24 * SYNOPSIS: 25 * 26 * double complex csqrt(); 27 * double complex z, w; 28 * 29 * w = csqrt (z); 30 * 31 * 32 * 33 * DESCRIPTION: 34 * 35 * 36 * If z = x + iy, r = |z|, then 37 * 38 * 1/2 39 * Re w = [ (r + x)/2 ] , 40 * 41 * 1/2 42 * Im w = [ (r - x)/2 ] . 43 * 44 * Cancellation error in r-x or r+x is avoided by using the 45 * identity 2 Re w Im w = y. 46 * 47 * Note that -w is also a square root of z. The root chosen 48 * is always in the right half plane and Im w has the same sign as y. 49 * 50 * 51 * 52 * ACCURACY: 53 * 54 * Relative error: 55 * arithmetic domain # trials peak rms 56 * DEC -10,+10 25000 3.2e-17 9.6e-18 57 * IEEE -10,+10 1,000,000 2.9e-16 6.1e-17 58 * 59 */ 60 61 #include <complex.h> 62 #include <math.h> 63 64 double complex 65 csqrt(double complex z) 66 { 67 double complex w; 68 double x, y, r, t, scale; 69 70 x = creal (z); 71 y = cimag (z); 72 73 if (y == 0.0) { 74 if (x == 0.0) { 75 w = 0.0 + y * I; 76 } 77 else { 78 r = fabs (x); 79 r = sqrt (r); 80 if (x < 0.0) { 81 w = 0.0 + r * I; 82 } 83 else { 84 w = r + y * I; 85 } 86 } 87 return (w); 88 } 89 if (x == 0.0) { 90 r = fabs (y); 91 r = sqrt (0.5*r); 92 if (y > 0) 93 w = r + r * I; 94 else 95 w = r - r * I; 96 return (w); 97 } 98 /* Rescale to avoid internal overflow or underflow. */ 99 if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) { 100 x *= 0.25; 101 y *= 0.25; 102 scale = 2.0; 103 } 104 else { 105 x *= 1.8014398509481984e16; /* 2^54 */ 106 y *= 1.8014398509481984e16; 107 scale = 7.450580596923828125e-9; /* 2^-27 */ 108 #if 0 109 x *= 4.0; 110 y *= 4.0; 111 scale = 0.5; 112 #endif 113 } 114 w = x + y * I; 115 r = cabs(w); 116 if (x > 0) { 117 t = sqrt(0.5 * r + 0.5 * x); 118 r = scale * fabs((0.5 * y) / t); 119 t *= scale; 120 } 121 else { 122 r = sqrt( 0.5 * r - 0.5 * x ); 123 t = scale * fabs( (0.5 * y) / r ); 124 r *= scale; 125 } 126 if (y < 0) 127 w = t - r * I; 128 else 129 w = t + r * I; 130 return (w); 131 } 132