OLD | NEW |
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_csqrtf.c */ | 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_csqrtf.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 18 matching lines...) Expand all Loading... |
29 | 29 |
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 float complex csqrtf(float complex z) | 39 float complex csqrtf(float complex z) { |
40 { | 40 float a = crealf(z), b = cimagf(z); |
41 » float a = crealf(z), b = cimagf(z); | 41 double t; |
42 » double t; | |
43 | 42 |
44 » /* Handle special cases. */ | 43 /* Handle special cases. */ |
45 » if (z == 0) | 44 if (z == 0) |
46 » » return CMPLXF(0, b); | 45 return CMPLXF(0, b); |
47 » if (isinf(b)) | 46 if (isinf(b)) |
48 » » return CMPLXF(INFINITY, b); | 47 return CMPLXF(INFINITY, b); |
49 » if (isnan(a)) { | 48 if (isnan(a)) { |
50 » » t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ | 49 t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ |
51 » » return CMPLXF(a, t); /* return NaN + NaN i */ | 50 return CMPLXF(a, t); /* return NaN + NaN i */ |
52 » } | 51 } |
53 » if (isinf(a)) { | 52 if (isinf(a)) { |
54 » » /* | 53 /* |
55 » » * csqrtf(inf + NaN i) = inf + NaN i | 54 * csqrtf(inf + NaN i) = inf + NaN i |
56 » » * csqrtf(inf + y i) = inf + 0 i | 55 * csqrtf(inf + y i) = inf + 0 i |
57 » » * csqrtf(-inf + NaN i) = NaN +- inf i | 56 * csqrtf(-inf + NaN i) = NaN +- inf i |
58 » » * csqrtf(-inf + y i) = 0 + inf i | 57 * csqrtf(-inf + y i) = 0 + inf i |
59 » » */ | 58 */ |
60 » » if (signbit(a)) | 59 if (signbit(a)) |
61 » » » return CMPLXF(fabsf(b - b), copysignf(a, b)); | 60 return CMPLXF(fabsf(b - b), copysignf(a, b)); |
62 » » else | 61 else |
63 » » » return CMPLXF(a, copysignf(b - b, b)); | 62 return CMPLXF(a, copysignf(b - b, b)); |
64 » } | 63 } |
65 » /* | 64 /* |
66 » * The remaining special case (b is NaN) is handled just fine by | 65 * The remaining special case (b is NaN) is handled just fine by |
67 » * the normal code path below. | 66 * the normal code path below. |
68 » */ | 67 */ |
69 | 68 |
70 » /* | 69 /* |
71 » * We compute t in double precision to avoid overflow and to | 70 * We compute t in double precision to avoid overflow and to |
72 » * provide correct rounding in nearly all cases. | 71 * provide correct rounding in nearly all cases. |
73 » * This is Algorithm 312, CACM vol 10, Oct 1967. | 72 * This is Algorithm 312, CACM vol 10, Oct 1967. |
74 » */ | 73 */ |
75 » if (a >= 0) { | 74 if (a >= 0) { |
76 » » t = sqrt((a + hypot(a, b)) * 0.5); | 75 t = sqrt((a + hypot(a, b)) * 0.5); |
77 » » return CMPLXF(t, b / (2.0 * t)); | 76 return CMPLXF(t, b / (2.0 * t)); |
78 » } else { | 77 } else { |
79 » » t = sqrt((-a + hypot(a, b)) * 0.5); | 78 t = sqrt((-a + hypot(a, b)) * 0.5); |
80 » » return CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)); | 79 return CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)); |
81 » } | 80 } |
82 } | 81 } |
OLD | NEW |