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 |