1 /* $OpenBSD: s_csqrtf.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 /* csqrtf() 19 * 20 * Complex square root 21 * 22 * 23 * 24 * SYNOPSIS: 25 * 26 * float complex csqrtf(); 27 * float complex z, w; 28 * 29 * w = csqrtf( 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 * 55 * Relative error: 56 * arithmetic domain # trials peak rms 57 * IEEE -10,+10 1,000,000 1.8e-7 3.5e-8 58 * 59 */ 60 61 #include <complex.h> 62 #include <math.h> 63 64 float complex 65 csqrtf(float complex z) 66 { 67 float complex w; 68 float x, y, r, t, scale; 69 70 x = creal(z); 71 y = cimag(z); 72 73 if(y == 0.0f) { 74 if (x < 0.0f) { 75 w = 0.0f + sqrtf(-x) * I; 76 return (w); 77 } 78 else if (x == 0.0f) { 79 return (0.0f + y * I); 80 } 81 else { 82 w = sqrtf(x) + y * I; 83 return (w); 84 } 85 } 86 87 if (x == 0.0f) { 88 r = fabsf(y); 89 r = sqrtf(0.5f*r); 90 if(y > 0) 91 w = r + r * I; 92 else 93 w = r - r * I; 94 return (w); 95 } 96 97 /* Rescale to avoid internal overflow or underflow. */ 98 if ((fabsf(x) > 4.0f) || (fabsf(y) > 4.0f)) { 99 x *= 0.25f; 100 y *= 0.25f; 101 scale = 2.0f; 102 } 103 else { 104 x *= 6.7108864e7f; /* 2^26 */ 105 y *= 6.7108864e7f; 106 scale = 1.220703125e-4f; /* 2^-13 */ 107 #if 0 108 x *= 4.0f; 109 y *= 4.0f; 110 scale = 0.5f; 111 #endif 112 } 113 w = x + y * I; 114 r = cabsf(w); 115 if (x > 0) { 116 t = sqrtf( 0.5f * r + 0.5f * x ); 117 r = scale * fabsf((0.5f * y) / t); 118 t *= scale; 119 } 120 else { 121 r = sqrtf(0.5f * r - 0.5f * x); 122 t = scale * fabsf((0.5f * y) / r); 123 r *= scale; 124 } 125 126 if (y < 0) 127 w = t - r * I; 128 else 129 w = t + r * I; 130 return (w); 131 } 132