| OLD | NEW |
| 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_csqrt.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_csqrt.c */ |
| 2 /*- | 2 /*- |
| 3 * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG> | 3 * Copyright (c) 2007 David Schultz <das@FreeBSD.ORG> |
| 4 * All rights reserved. | 4 * All rights reserved. |
| 5 * | 5 * |
| 6 * Redistribution and use in source and binary forms, with or without | 6 * Redistribution and use in source and binary forms, with or without |
| 7 * modification, are permitted provided that the following conditions | 7 * modification, are permitted provided that the following conditions |
| 8 * are met: | 8 * are met: |
| 9 * 1. Redistributions of source code must retain the above copyright | 9 * 1. Redistributions of source code must retain the above copyright |
| 10 * notice, this list of conditions and the following disclaimer. | 10 * notice, this list of conditions and the following disclaimer. |
| (...skipping 19 matching lines...) Expand all Loading... |
| 30 /* | 30 /* |
| 31 * gcc doesn't implement complex multiplication or division correctly, | 31 * gcc doesn't implement complex multiplication or division correctly, |
| 32 * so we need to handle infinities specially. We turn on this pragma to | 32 * so we need to handle infinities specially. We turn on this pragma to |
| 33 * notify conforming c99 compilers that the fast-but-incorrect code that | 33 * notify conforming c99 compilers that the fast-but-incorrect code that |
| 34 * gcc generates is acceptable, since the special cases have already been | 34 * gcc generates is acceptable, since the special cases have already been |
| 35 * handled. | 35 * handled. |
| 36 */ | 36 */ |
| 37 #pragma STDC CX_LIMITED_RANGE ON | 37 #pragma STDC CX_LIMITED_RANGE ON |
| 38 | 38 |
| 39 /* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */ | 39 /* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */ |
| 40 #define THRESH 0x1.a827999fcef32p+1022 | 40 #define THRESH 0x1.a827999fcef32p+1022 |
| 41 | 41 |
| 42 double complex csqrt(double complex z) | 42 double complex csqrt(double complex z) { |
| 43 { | 43 double complex result; |
| 44 » double complex result; | 44 double a, b; |
| 45 » double a, b; | 45 double t; |
| 46 » double t; | 46 int scale; |
| 47 » int scale; | |
| 48 | 47 |
| 49 » a = creal(z); | 48 a = creal(z); |
| 50 » b = cimag(z); | 49 b = cimag(z); |
| 51 | 50 |
| 52 » /* Handle special cases. */ | 51 /* Handle special cases. */ |
| 53 » if (z == 0) | 52 if (z == 0) |
| 54 » » return CMPLX(0, b); | 53 return CMPLX(0, b); |
| 55 » if (isinf(b)) | 54 if (isinf(b)) |
| 56 » » return CMPLX(INFINITY, b); | 55 return CMPLX(INFINITY, b); |
| 57 » if (isnan(a)) { | 56 if (isnan(a)) { |
| 58 » » t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ | 57 t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ |
| 59 » » return CMPLX(a, t); /* return NaN + NaN i */ | 58 return CMPLX(a, t); /* return NaN + NaN i */ |
| 60 » } | 59 } |
| 61 » if (isinf(a)) { | 60 if (isinf(a)) { |
| 62 » » /* | 61 /* |
| 63 » » * csqrt(inf + NaN i) = inf + NaN i | 62 * csqrt(inf + NaN i) = inf + NaN i |
| 64 » » * csqrt(inf + y i) = inf + 0 i | 63 * csqrt(inf + y i) = inf + 0 i |
| 65 » » * csqrt(-inf + NaN i) = NaN +- inf i | 64 * csqrt(-inf + NaN i) = NaN +- inf i |
| 66 » » * csqrt(-inf + y i) = 0 + inf i | 65 * csqrt(-inf + y i) = 0 + inf i |
| 67 » » */ | 66 */ |
| 68 » » if (signbit(a)) | 67 if (signbit(a)) |
| 69 » » » return CMPLX(fabs(b - b), copysign(a, b)); | 68 return CMPLX(fabs(b - b), copysign(a, b)); |
| 70 » » else | 69 else |
| 71 » » » return CMPLX(a, copysign(b - b, b)); | 70 return CMPLX(a, copysign(b - b, b)); |
| 72 » } | 71 } |
| 73 » /* | 72 /* |
| 74 » * The remaining special case (b is NaN) is handled just fine by | 73 * The remaining special case (b is NaN) is handled just fine by |
| 75 » * the normal code path below. | 74 * the normal code path below. |
| 76 » */ | 75 */ |
| 77 | 76 |
| 78 » /* Scale to avoid overflow. */ | 77 /* Scale to avoid overflow. */ |
| 79 » if (fabs(a) >= THRESH || fabs(b) >= THRESH) { | 78 if (fabs(a) >= THRESH || fabs(b) >= THRESH) { |
| 80 » » a *= 0.25; | 79 a *= 0.25; |
| 81 » » b *= 0.25; | 80 b *= 0.25; |
| 82 » » scale = 1; | 81 scale = 1; |
| 83 » } else { | 82 } else { |
| 84 » » scale = 0; | 83 scale = 0; |
| 85 » } | 84 } |
| 86 | 85 |
| 87 » /* Algorithm 312, CACM vol 10, Oct 1967. */ | 86 /* Algorithm 312, CACM vol 10, Oct 1967. */ |
| 88 » if (a >= 0) { | 87 if (a >= 0) { |
| 89 » » t = sqrt((a + hypot(a, b)) * 0.5); | 88 t = sqrt((a + hypot(a, b)) * 0.5); |
| 90 » » result = CMPLX(t, b / (2 * t)); | 89 result = CMPLX(t, b / (2 * t)); |
| 91 » } else { | 90 } else { |
| 92 » » t = sqrt((-a + hypot(a, b)) * 0.5); | 91 t = sqrt((-a + hypot(a, b)) * 0.5); |
| 93 » » result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); | 92 result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); |
| 94 » } | 93 } |
| 95 | 94 |
| 96 » /* Rescale. */ | 95 /* Rescale. */ |
| 97 » if (scale) | 96 if (scale) |
| 98 » » result *= 2; | 97 result *= 2; |
| 99 » return result; | 98 return result; |
| 100 } | 99 } |
| OLD | NEW |