1*bf3b7e6dSmaya /* $NetBSD: csqrt.c,v 1.4 2017/01/01 19:32:14 maya Exp $ */
29d8b5fa7Sdrochner
39d8b5fa7Sdrochner /*-
49d8b5fa7Sdrochner * Copyright (c) 2007 The NetBSD Foundation, Inc.
59d8b5fa7Sdrochner * All rights reserved.
69d8b5fa7Sdrochner *
79d8b5fa7Sdrochner * This code is derived from software written by Stephen L. Moshier.
89d8b5fa7Sdrochner * It is redistributed by the NetBSD Foundation by permission of the author.
99d8b5fa7Sdrochner *
109d8b5fa7Sdrochner * Redistribution and use in source and binary forms, with or without
119d8b5fa7Sdrochner * modification, are permitted provided that the following conditions
129d8b5fa7Sdrochner * are met:
139d8b5fa7Sdrochner * 1. Redistributions of source code must retain the above copyright
149d8b5fa7Sdrochner * notice, this list of conditions and the following disclaimer.
159d8b5fa7Sdrochner * 2. Redistributions in binary form must reproduce the above copyright
169d8b5fa7Sdrochner * notice, this list of conditions and the following disclaimer in the
179d8b5fa7Sdrochner * documentation and/or other materials provided with the distribution.
189d8b5fa7Sdrochner *
199d8b5fa7Sdrochner * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
209d8b5fa7Sdrochner * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
219d8b5fa7Sdrochner * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
229d8b5fa7Sdrochner * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
239d8b5fa7Sdrochner * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
249d8b5fa7Sdrochner * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
259d8b5fa7Sdrochner * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
269d8b5fa7Sdrochner * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
279d8b5fa7Sdrochner * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
289d8b5fa7Sdrochner * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
299d8b5fa7Sdrochner * POSSIBILITY OF SUCH DAMAGE.
309d8b5fa7Sdrochner */
319d8b5fa7Sdrochner
329d8b5fa7Sdrochner #include <complex.h>
339d8b5fa7Sdrochner #include <math.h>
349d8b5fa7Sdrochner
359d8b5fa7Sdrochner double complex
csqrt(double complex z)369d8b5fa7Sdrochner csqrt(double complex z)
379d8b5fa7Sdrochner {
389d8b5fa7Sdrochner double complex w;
399d8b5fa7Sdrochner double x, y, r, t, scale;
409d8b5fa7Sdrochner
419d8b5fa7Sdrochner x = creal (z);
429d8b5fa7Sdrochner y = cimag (z);
439d8b5fa7Sdrochner
44*bf3b7e6dSmaya /*
45*bf3b7e6dSmaya * input is a real number and imaginary part isn't -0.0.
46*bf3b7e6dSmaya * negative zero is on the branch cut.
47*bf3b7e6dSmaya */
4800e942c2Smaya if ((y == 0.0) && !signbit(y)) {
499d8b5fa7Sdrochner if (x == 0.0) {
509d8b5fa7Sdrochner w = 0.0 + y * I;
519d8b5fa7Sdrochner } else {
529d8b5fa7Sdrochner if (x < 0.0) {
531102a777Smaya r = sqrt(-x);
549d8b5fa7Sdrochner w = 0.0 + r * I;
559d8b5fa7Sdrochner } else {
561102a777Smaya r = sqrt(x);
571102a777Smaya w = r;
589d8b5fa7Sdrochner }
599d8b5fa7Sdrochner }
609d8b5fa7Sdrochner return w;
619d8b5fa7Sdrochner }
629d8b5fa7Sdrochner if (x == 0.0) {
631102a777Smaya if (y > 0) {
641102a777Smaya r = sqrt(0.5 * y);
659d8b5fa7Sdrochner w = r + r * I;
661102a777Smaya } else {
671102a777Smaya r = sqrt(-0.5 * y);
689d8b5fa7Sdrochner w = r - r * I;
691102a777Smaya }
709d8b5fa7Sdrochner return w;
719d8b5fa7Sdrochner }
729d8b5fa7Sdrochner /* Rescale to avoid internal overflow or underflow. */
739d8b5fa7Sdrochner if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) {
749d8b5fa7Sdrochner x *= 0.25;
759d8b5fa7Sdrochner y *= 0.25;
769d8b5fa7Sdrochner scale = 2.0;
779d8b5fa7Sdrochner } else {
789d8b5fa7Sdrochner #if 1
799d8b5fa7Sdrochner x *= 1.8014398509481984e16; /* 2^54 */
809d8b5fa7Sdrochner y *= 1.8014398509481984e16;
819d8b5fa7Sdrochner scale = 7.450580596923828125e-9; /* 2^-27 */
829d8b5fa7Sdrochner #else
839d8b5fa7Sdrochner x *= 4.0;
849d8b5fa7Sdrochner y *= 4.0;
859d8b5fa7Sdrochner scale = 0.5;
869d8b5fa7Sdrochner #endif
879d8b5fa7Sdrochner }
889d8b5fa7Sdrochner w = x + y * I;
899d8b5fa7Sdrochner r = cabs(w);
909d8b5fa7Sdrochner if (x > 0) {
919d8b5fa7Sdrochner t = sqrt(0.5 * r + 0.5 * x);
929d8b5fa7Sdrochner r = scale * fabs((0.5 * y) / t );
939d8b5fa7Sdrochner t *= scale;
949d8b5fa7Sdrochner } else {
959d8b5fa7Sdrochner r = sqrt(0.5 * r - 0.5 * x);
969d8b5fa7Sdrochner t = scale * fabs((0.5 * y) / r);
979d8b5fa7Sdrochner r *= scale;
989d8b5fa7Sdrochner }
99*bf3b7e6dSmaya if (y > 0)
1009d8b5fa7Sdrochner w = t + r * I;
101*bf3b7e6dSmaya else
102*bf3b7e6dSmaya w = t - r * I;
1039d8b5fa7Sdrochner return w;
1049d8b5fa7Sdrochner }
105