Deleted Added
sdiff udiff text old ( 251404 ) new ( 275819 )
full compact
catrigf.c (251404) catrigf.c (275819)
1/*-
2 * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright

--- 25 unchanged lines hidden (view full) ---

34 * See catrig.c for complete comments.
35 *
36 * XXX comments were removed automatically, and even short ones on the right
37 * of statements were removed (all of them), contrary to normal style. Only
38 * a few comments on the right of declarations remain.
39 */
40
41#include <sys/cdefs.h>
1/*-
2 * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright

--- 25 unchanged lines hidden (view full) ---

34 * See catrig.c for complete comments.
35 *
36 * XXX comments were removed automatically, and even short ones on the right
37 * of statements were removed (all of them), contrary to normal style. Only
38 * a few comments on the right of declarations remain.
39 */
40
41#include <sys/cdefs.h>
42__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 251404 2013-06-05 05:33:01Z das $");
42__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 275819 2014-12-16 09:21:56Z ed $");
43
44#include <complex.h>
45#include <float.h>
46
47#include "math.h"
48#include "math_private.h"
49
50#undef isinf

--- 100 unchanged lines hidden (view full) ---

151
152 x = crealf(z);
153 y = cimagf(z);
154 ax = fabsf(x);
155 ay = fabsf(y);
156
157 if (isnan(x) || isnan(y)) {
158 if (isinf(x))
43
44#include <complex.h>
45#include <float.h>
46
47#include "math.h"
48#include "math_private.h"
49
50#undef isinf

--- 100 unchanged lines hidden (view full) ---

151
152 x = crealf(z);
153 y = cimagf(z);
154 ax = fabsf(x);
155 ay = fabsf(y);
156
157 if (isnan(x) || isnan(y)) {
158 if (isinf(x))
159 return (cpackf(x, y + y));
159 return (CMPLXF(x, y + y));
160 if (isinf(y))
160 if (isinf(y))
161 return (cpackf(y, x + x));
161 return (CMPLXF(y, x + x));
162 if (y == 0)
162 if (y == 0)
163 return (cpackf(x + x, y));
164 return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
163 return (CMPLXF(x + x, y));
164 return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
165 }
166
167 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
168 if (signbit(x) == 0)
169 w = clog_for_large_values(z) + m_ln2;
170 else
171 w = clog_for_large_values(-z) + m_ln2;
165 }
166
167 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
168 if (signbit(x) == 0)
169 w = clog_for_large_values(z) + m_ln2;
170 else
171 w = clog_for_large_values(-z) + m_ln2;
172 return (cpackf(copysignf(crealf(w), x),
172 return (CMPLXF(copysignf(crealf(w), x),
173 copysignf(cimagf(w), y)));
174 }
175
176 if (x == 0 && y == 0)
177 return (z);
178
179 raise_inexact();
180
181 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
182 return (z);
183
184 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
185 if (B_is_usable)
186 ry = asinf(B);
187 else
188 ry = atan2f(new_y, sqrt_A2my2);
173 copysignf(cimagf(w), y)));
174 }
175
176 if (x == 0 && y == 0)
177 return (z);
178
179 raise_inexact();
180
181 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
182 return (z);
183
184 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
185 if (B_is_usable)
186 ry = asinf(B);
187 else
188 ry = atan2f(new_y, sqrt_A2my2);
189 return (cpackf(copysignf(rx, x), copysignf(ry, y)));
189 return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
190}
191
192float complex
193casinf(float complex z)
194{
190}
191
192float complex
193casinf(float complex z)
194{
195 float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
195 float complex w = casinhf(CMPLXF(cimagf(z), crealf(z)));
196
196
197 return (cpackf(cimagf(w), crealf(w)));
197 return (CMPLXF(cimagf(w), crealf(w)));
198}
199
200float complex
201cacosf(float complex z)
202{
203 float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
204 int sx, sy;
205 int B_is_usable;
206 float complex w;
207
208 x = crealf(z);
209 y = cimagf(z);
210 sx = signbit(x);
211 sy = signbit(y);
212 ax = fabsf(x);
213 ay = fabsf(y);
214
215 if (isnan(x) || isnan(y)) {
216 if (isinf(x))
198}
199
200float complex
201cacosf(float complex z)
202{
203 float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
204 int sx, sy;
205 int B_is_usable;
206 float complex w;
207
208 x = crealf(z);
209 y = cimagf(z);
210 sx = signbit(x);
211 sy = signbit(y);
212 ax = fabsf(x);
213 ay = fabsf(y);
214
215 if (isnan(x) || isnan(y)) {
216 if (isinf(x))
217 return (cpackf(y + y, -INFINITY));
217 return (CMPLXF(y + y, -INFINITY));
218 if (isinf(y))
218 if (isinf(y))
219 return (cpackf(x + x, -y));
219 return (CMPLXF(x + x, -y));
220 if (x == 0)
220 if (x == 0)
221 return (cpackf(pio2_hi + pio2_lo, y + y));
222 return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
221 return (CMPLXF(pio2_hi + pio2_lo, y + y));
222 return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
223 }
224
225 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
226 w = clog_for_large_values(z);
227 rx = fabsf(cimagf(w));
228 ry = crealf(w) + m_ln2;
229 if (sy == 0)
230 ry = -ry;
223 }
224
225 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
226 w = clog_for_large_values(z);
227 rx = fabsf(cimagf(w));
228 ry = crealf(w) + m_ln2;
229 if (sy == 0)
230 ry = -ry;
231 return (cpackf(rx, ry));
231 return (CMPLXF(rx, ry));
232 }
233
234 if (x == 1 && y == 0)
232 }
233
234 if (x == 1 && y == 0)
235 return (cpackf(0, -y));
235 return (CMPLXF(0, -y));
236
237 raise_inexact();
238
239 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
236
237 raise_inexact();
238
239 if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
240 return (cpackf(pio2_hi - (x - pio2_lo), -y));
240 return (CMPLXF(pio2_hi - (x - pio2_lo), -y));
241
242 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
243 if (B_is_usable) {
244 if (sx == 0)
245 rx = acosf(B);
246 else
247 rx = acosf(-B);
248 } else {
249 if (sx == 0)
250 rx = atan2f(sqrt_A2mx2, new_x);
251 else
252 rx = atan2f(sqrt_A2mx2, -new_x);
253 }
254 if (sy == 0)
255 ry = -ry;
241
242 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
243 if (B_is_usable) {
244 if (sx == 0)
245 rx = acosf(B);
246 else
247 rx = acosf(-B);
248 } else {
249 if (sx == 0)
250 rx = atan2f(sqrt_A2mx2, new_x);
251 else
252 rx = atan2f(sqrt_A2mx2, -new_x);
253 }
254 if (sy == 0)
255 ry = -ry;
256 return (cpackf(rx, ry));
256 return (CMPLXF(rx, ry));
257}
258
259float complex
260cacoshf(float complex z)
261{
262 float complex w;
263 float rx, ry;
264
265 w = cacosf(z);
266 rx = crealf(w);
267 ry = cimagf(w);
268 if (isnan(rx) && isnan(ry))
257}
258
259float complex
260cacoshf(float complex z)
261{
262 float complex w;
263 float rx, ry;
264
265 w = cacosf(z);
266 rx = crealf(w);
267 ry = cimagf(w);
268 if (isnan(rx) && isnan(ry))
269 return (cpackf(ry, rx));
269 return (CMPLXF(ry, rx));
270 if (isnan(rx))
270 if (isnan(rx))
271 return (cpackf(fabsf(ry), rx));
271 return (CMPLXF(fabsf(ry), rx));
272 if (isnan(ry))
272 if (isnan(ry))
273 return (cpackf(ry, ry));
274 return (cpackf(fabsf(ry), copysignf(rx, cimagf(z))));
273 return (CMPLXF(ry, ry));
274 return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z))));
275}
276
277static float complex
278clog_for_large_values(float complex z)
279{
280 float x, y;
281 float ax, ay, t;
282
283 x = crealf(z);
284 y = cimagf(z);
285 ax = fabsf(x);
286 ay = fabsf(y);
287 if (ax < ay) {
288 t = ax;
289 ax = ay;
290 ay = t;
291 }
292
293 if (ax > FLT_MAX / 2)
275}
276
277static float complex
278clog_for_large_values(float complex z)
279{
280 float x, y;
281 float ax, ay, t;
282
283 x = crealf(z);
284 y = cimagf(z);
285 ax = fabsf(x);
286 ay = fabsf(y);
287 if (ax < ay) {
288 t = ax;
289 ax = ay;
290 ay = t;
291 }
292
293 if (ax > FLT_MAX / 2)
294 return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
294 return (CMPLXF(logf(hypotf(x / m_e, y / m_e)) + 1,
295 atan2f(y, x)));
296
297 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
295 atan2f(y, x)));
296
297 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
298 return (cpackf(logf(hypotf(x, y)), atan2f(y, x)));
298 return (CMPLXF(logf(hypotf(x, y)), atan2f(y, x)));
299
299
300 return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
300 return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
301}
302
303static inline float
304sum_squares(float x, float y)
305{
306
307 if (y < SQRT_MIN)
308 return (x * x);

--- 32 unchanged lines hidden (view full) ---

341 float x, y, ax, ay, rx, ry;
342
343 x = crealf(z);
344 y = cimagf(z);
345 ax = fabsf(x);
346 ay = fabsf(y);
347
348 if (y == 0 && ax <= 1)
301}
302
303static inline float
304sum_squares(float x, float y)
305{
306
307 if (y < SQRT_MIN)
308 return (x * x);

--- 32 unchanged lines hidden (view full) ---

341 float x, y, ax, ay, rx, ry;
342
343 x = crealf(z);
344 y = cimagf(z);
345 ax = fabsf(x);
346 ay = fabsf(y);
347
348 if (y == 0 && ax <= 1)
349 return (cpackf(atanhf(x), y));
349 return (CMPLXF(atanhf(x), y));
350
351 if (x == 0)
350
351 if (x == 0)
352 return (cpackf(x, atanf(y)));
352 return (CMPLXF(x, atanf(y)));
353
354 if (isnan(x) || isnan(y)) {
355 if (isinf(x))
353
354 if (isnan(x) || isnan(y)) {
355 if (isinf(x))
356 return (cpackf(copysignf(0, x), y + y));
356 return (CMPLXF(copysignf(0, x), y + y));
357 if (isinf(y))
357 if (isinf(y))
358 return (cpackf(copysignf(0, x),
358 return (CMPLXF(copysignf(0, x),
359 copysignf(pio2_hi + pio2_lo, y)));
359 copysignf(pio2_hi + pio2_lo, y)));
360 return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
360 return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
361 }
362
363 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
361 }
362
363 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
364 return (cpackf(real_part_reciprocal(x, y),
364 return (CMPLXF(real_part_reciprocal(x, y),
365 copysignf(pio2_hi + pio2_lo, y)));
366
367 if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
368 raise_inexact();
369 return (z);
370 }
371
372 if (ax == 1 && ay < FLT_EPSILON)
373 rx = (m_ln2 - logf(ay)) / 2;
374 else
375 rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
376
377 if (ax == 1)
378 ry = atan2f(2, -ay) / 2;
379 else if (ay < FLT_EPSILON)
380 ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2;
381 else
382 ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
383
365 copysignf(pio2_hi + pio2_lo, y)));
366
367 if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
368 raise_inexact();
369 return (z);
370 }
371
372 if (ax == 1 && ay < FLT_EPSILON)
373 rx = (m_ln2 - logf(ay)) / 2;
374 else
375 rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
376
377 if (ax == 1)
378 ry = atan2f(2, -ay) / 2;
379 else if (ay < FLT_EPSILON)
380 ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2;
381 else
382 ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
383
384 return (cpackf(copysignf(rx, x), copysignf(ry, y)));
384 return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
385}
386
387float complex
388catanf(float complex z)
389{
385}
386
387float complex
388catanf(float complex z)
389{
390 float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
390 float complex w = catanhf(CMPLXF(cimagf(z), crealf(z)));
391
391
392 return (cpackf(cimagf(w), crealf(w)));
392 return (CMPLXF(cimagf(w), crealf(w)));
393}
393}