Chromium Code Reviews
chromiumcodereview-hr@appspot.gserviceaccount.com (chromiumcodereview-hr) | Please choose your nickname with Settings | Help | Chromium Project | Gerrit Changes | Sign out
(997)

Side by Side Diff: fusl/src/complex/csqrt.c

Issue 1714623002: [fusl] clang-format fusl (Closed) Base URL: git@github.com:domokit/mojo.git@master
Patch Set: headers too Created 4 years, 10 months ago
Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.
Jump to:
View unified diff | Download patch
OLDNEW
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
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 }
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698