1 /* $NetBSD: catrigl.c,v 1.3 2022/04/19 20:32:16 rillig Exp $ */
2 /*-
3 * Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
13 * documentation and/or other materials provided with the distribution.
14 *
15 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
16 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25 * SUCH DAMAGE.
26 */
27
28 /*
29 * The algorithm is very close to that in "Implementing the complex arcsine
30 * and arccosine functions using exception handling" by T. E. Hull, Thomas F.
31 * Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
32 * Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
33 * http://dl.acm.org/citation.cfm?id=275324.
34 *
35 * The code for catrig.c contains complete comments.
36 */
37 #include <sys/cdefs.h>
38 __RCSID("$NetBSD: catrigl.c,v 1.3 2022/04/19 20:32:16 rillig Exp $");
39
40 #include "namespace.h"
41 #ifdef __weak_alias
42 __weak_alias(casinl, _casinl)
43 #endif
44 #ifdef __weak_alias
45 __weak_alias(catanl, _catanl)
46 #endif
47
48
49 #include <sys/param.h>
50 #include <complex.h>
51 #include <float.h>
52 #include <math.h>
53 #ifdef notyet // missing log1pl __HAVE_LONG_DOUBLE
54
55 #include "math_private.h"
56
57 #undef isinf
58 #define isinf(x) (fabsl(x) == INFINITY)
59 #undef isnan
60 #define isnan(x) ((x) != (x))
61 #define raise_inexact() do { volatile float junk __unused = /*LINTED*/1 + tiny; } while (0)
62 #undef signbit
63 #define signbit(x) (__builtin_signbitl(x))
64
65 #if __HAVE_LONG_DOUBLE + 0 == 128
66 // Ok
67 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
68 // XXX: Byte order
69 #define EXT_EXPBITS 15
70 struct ieee_ext {
71 uint64_t ext_frac;
72 uint16_t ext_exp:EXT_EXPBITS;
73 uint16_t ext_sign:1;
74 uint16_t ext_pad;
75 };
76 #define extu_exp extu_ext.ext_exp
77 #define extu_sign extu_ext.ext_sign
78 #define extu_frac extu_ext.ext_frac
79 union ieee_ext_u {
80 long double extu_ld;
81 struct ieee_ext extu_ext;
82 };
83 #else
84 #error "unsupported long double format"
85 #endif
86
87 #define GET_LDBL_EXPSIGN(r, s) \
88 do { \
89 union ieee_ext_u u; \
90 u.extu_ld = s; \
91 r = u.extu_sign; \
92 r >>= EXT_EXPBITS - 1; \
93 } while (0)
94 #define SET_LDBL_EXPSIGN(s, r) \
95 do { \
96 union ieee_ext_u u; \
97 u.extu_ld = s; \
98 u.extu_exp &= __BITS(0, EXT_EXPBITS - 1); \
99 u.extu_exp |= (r) << (EXT_EXPBITS - 1); \
100 s = u.extu_ld; \
101 } while (0)
102
103 static const long double
104 A_crossover = 10,
105 B_crossover = 0.6417,
106 FOUR_SQRT_MIN = 0x1p-8189L,
107 QUARTER_SQRT_MAX = 0x1p8189L,
108 RECIP_EPSILON = 1/LDBL_EPSILON,
109 SQRT_MIN = 0x1p-8191L;
110
111 static const long double
112 m_e = 2.71828182845904523536028747135266250e0L, /* 0x15bf0a8b1457695355fb8ac404e7a.0p-111 */
113 m_ln2 = 6.93147180559945309417232121458176568e-1L, /* 0x162e42fefa39ef35793c7673007e6.0p-113 */
114 pio2_hi = 1.5707963267948966192313216916397514L, /* pi/2 */
115 SQRT_3_EPSILON = 2.40370335797945490975336727199878124e-17L, /* 0x1bb67ae8584caa73b25742d7078b8.0p-168 */
116 SQRT_6_EPSILON = 3.39934988877629587239082586223300391e-17L; /* 0x13988e1409212e7d0321914321a55.0p-167 */
117
118 static const volatile double
119 pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */
120 static const volatile float
121 tiny = 0x1p-100;
122
123 static long double complex clog_for_large_values(long double complex z);
124
125 inline static long double
f(long double a,long double b,long double hypot_a_b)126 f(long double a, long double b, long double hypot_a_b)
127 {
128 if (b < 0)
129 return ((hypot_a_b - b) / 2);
130 if (b == 0)
131 return (a / 2);
132 return (a * a / (hypot_a_b + b) / 2);
133 }
134
135 inline static void
do_hard_work(long double x,long double y,long double * rx,int * B_is_usable,long double * B,long double * sqrt_A2my2,long double * new_y)136 do_hard_work(long double x, long double y, long double *rx, int *B_is_usable, long double *B, long double *sqrt_A2my2, long double *new_y)
137 {
138 long double R, S, A;
139 long double Am1, Amy;
140
141 R = hypotl(x, y+1);
142 S = hypotl(x, y-1);
143
144 A = (R + S) / 2;
145 if (A < 1)
146 A = 1;
147
148 if (A < A_crossover) {
149 if (y == 1 && x < LDBL_EPSILON*LDBL_EPSILON/128) {
150 *rx = sqrtl(x);
151 } else if (x >= LDBL_EPSILON * fabsl(y-1)) {
152 Am1 = f(x, 1+y, R) + f(x, 1-y, S);
153 *rx = log1pl(Am1 + sqrtl(Am1*(A+1)));
154 } else if (y < 1) {
155 *rx = x/sqrtl((1-y)*(1+y));
156 } else {
157 *rx = log1pl((y-1) + sqrtl((y-1)*(y+1)));
158 }
159 } else
160 *rx = logl(A + sqrtl(A*A-1));
161
162 *new_y = y;
163
164 if (y < FOUR_SQRT_MIN) {
165 *B_is_usable = 0;
166 *sqrt_A2my2 = A * (2 / LDBL_EPSILON);
167 *new_y= y * (2 / LDBL_EPSILON);
168 return;
169 }
170
171 *B = y/A;
172 *B_is_usable = 1;
173
174 if (*B > B_crossover) {
175 *B_is_usable = 0;
176 if (y == 1 && x < LDBL_EPSILON/128) {
177 *sqrt_A2my2 = sqrtl(x)*sqrtl((A+y)/2);
178 } else if (x >= LDBL_EPSILON * fabsl(y-1)) {
179 Amy = f(x, y+1, R) + f(x, y-1, S);
180 *sqrt_A2my2 = sqrtl(Amy*(A+y));
181 } else if (y > 1) {
182 *sqrt_A2my2 = x * (4/LDBL_EPSILON/LDBL_EPSILON) * y /
183 sqrtl((y+1)*(y-1));
184 *new_y = y * (4/LDBL_EPSILON/LDBL_EPSILON);
185 } else {
186 *sqrt_A2my2 = sqrtl((1-y)*(1+y));
187 }
188 }
189 }
190
191 long double complex
casinhl(long double complex z)192 casinhl(long double complex z)
193 {
194 long double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
195 int B_is_usable;
196 long double complex w;
197
198 x = creall(z);
199 y = cimagl(z);
200 ax = fabsl(x);
201 ay = fabsl(y);
202
203 if (isnan(x) || isnan(y)) {
204 if (isinf(x))
205 return (CMPLXL(x, y+y));
206 if (isinf(y))
207 return (CMPLXL(y, x+x));
208 if (y == 0) return (CMPLXL(x+x, y));
209 return (CMPLXL(x+0.0L+(y+0), x+0.0L+(y+0)));
210 }
211
212 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
213 if (signbit(x) == 0)
214 w = clog_for_large_values(z) + m_ln2;
215 else
216 w = clog_for_large_values(-z) + m_ln2;
217 return (CMPLXL(copysignl(creall(w), x), copysignl(cimagl(w), y)));
218 }
219
220 if (x == 0 && y == 0)
221 return (z);
222
223 raise_inexact();
224
225 if (ax < SQRT_6_EPSILON/4 && ay < SQRT_6_EPSILON/4)
226 return (z);
227
228 do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
229 if (B_is_usable)
230 ry = asinl(B);
231 else
232 ry = atan2l(new_y, sqrt_A2my2);
233 return (CMPLXL(copysignl(rx, x), copysignl(ry, y)));
234 }
235
236 long double complex
casinl(long double complex z)237 casinl(long double complex z)
238 {
239 long double complex w = casinhl(CMPLXL(cimagl(z), creall(z)));
240 return (CMPLXL(cimagl(w), creall(w)));
241 }
242
243 long double complex
cacosl(long double complex z)244 cacosl(long double complex z)
245 {
246 long double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
247 int sx, sy;
248 int B_is_usable;
249 long double complex w;
250
251 x = creall(z);
252 y = cimagl(z);
253 sx = signbit(x);
254 sy = signbit(y);
255 ax = fabsl(x);
256 ay = fabsl(y);
257
258 if (isnan(x) || isnan(y)) {
259 if (isinf(x))
260 return (CMPLXL(y+y, -INFINITY));
261 if (isinf(y))
262 return (CMPLXL(x+x, -y));
263 if (x == 0) return (CMPLXL(pio2_hi + pio2_lo, y+y));
264 return (CMPLXL(x+0.0L+(y+0), x+0.0L+(y+0)));
265 }
266
267 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
268 w = clog_for_large_values(z);
269 rx = fabsl(cimagl(w));
270 ry = creall(w) + m_ln2;
271 if (sy == 0)
272 ry = -ry;
273 return (CMPLXL(rx, ry));
274 }
275
276 if (x == 1 && y == 0)
277 return (CMPLXL(0, -y));
278
279 raise_inexact();
280
281 if (ax < SQRT_6_EPSILON/4 && ay < SQRT_6_EPSILON/4)
282 return (CMPLXL(pio2_hi - (x - pio2_lo), -y));
283
284 do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
285 if (B_is_usable) {
286 if (sx==0)
287 rx = acosl(B);
288 else
289 rx = acosl(-B);
290 } else {
291 if (sx==0)
292 rx = atan2l(sqrt_A2mx2, new_x);
293 else
294 rx = atan2l(sqrt_A2mx2, -new_x);
295 }
296 if (sy==0)
297 ry = -ry;
298 return (CMPLXL(rx, ry));
299 }
300
301 long double complex
cacoshl(long double complex z)302 cacoshl(long double complex z)
303 {
304 long double complex w;
305 long double rx, ry;
306
307 w = cacosl(z);
308 rx = creall(w);
309 ry = cimagl(w);
310 if (isnan(rx) && isnan(ry))
311 return (CMPLXL(ry, rx));
312 if (isnan(rx))
313 return (CMPLXL(fabsl(ry), rx));
314 if (isnan(ry))
315 return (CMPLXL(ry, ry));
316 return (CMPLXL(fabsl(ry), copysignl(rx, cimagl(z))));
317 }
318
319 static long double complex
clog_for_large_values(long double complex z)320 clog_for_large_values(long double complex z)
321 {
322 long double x, y;
323 long double ax, ay, t;
324
325 x = creall(z);
326 y = cimagl(z);
327 ax = fabsl(x);
328 ay = fabsl(y);
329 if (ax < ay) {
330 t = ax;
331 ax = ay;
332 ay = t;
333 }
334
335 if (ax > LDBL_MAX / 2)
336 return (CMPLXL(logl(hypotl(x / m_e, y / m_e)) + 1, atan2l(y, x)));
337
338 if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
339 return (CMPLXL(logl(hypotl(x, y)), atan2l(y, x)));
340
341 return (CMPLXL(logl(ax*ax + ay*ay) / 2, atan2l(y, x)));
342 }
343
344 inline static long double
sum_squares(long double x,long double y)345 sum_squares(long double x, long double y)
346 {
347 if (y < SQRT_MIN)
348 return (x*x);
349
350 return (x*x + y*y);
351 }
352
353 inline static long double
real_part_reciprocal(long double x,long double y)354 real_part_reciprocal(long double x, long double y)
355 {
356 long double scale;
357 uint16_t hx, hy;
358 int16_t ix, iy;
359
360 GET_LDBL_EXPSIGN(hx, x);
361 ix = hx & 0x7fff;
362 GET_LDBL_EXPSIGN(hy, y);
363 iy = hy & 0x7fff;
364 #define BIAS (LDBL_MAX_EXP - 1)
365 #define CUTOFF (LDBL_MANT_DIG / 2 + 1)
366 if (ix - iy >= CUTOFF || isinf(x))
367 return (1/x);
368 if (iy - ix >= CUTOFF)
369 return (x/y/y);
370 if (ix <= BIAS + LDBL_MAX_EXP / 2 - CUTOFF)
371 return (x/(x*x + y*y));
372 scale = 1;
373 SET_LDBL_EXPSIGN(scale, 0x7fff - ix);
374 x *= scale;
375 y *= scale;
376 return (x/(x*x + y*y) * scale);
377 }
378
379 long double complex
catanhl(long double complex z)380 catanhl(long double complex z)
381 {
382 long double x, y, ax, ay, rx, ry;
383
384 x = creall(z);
385 y = cimagl(z);
386 ax = fabsl(x);
387 ay = fabsl(y);
388
389 if (y == 0 && ax <= 1)
390 return (CMPLXL(atanhl(x), y)); /* XXX need atanhl() */
391
392 if (x == 0)
393 return (CMPLXL(x, atanl(y)));
394
395 if (isnan(x) || isnan(y)) {
396 if (isinf(x))
397 return (CMPLXL(copysignl(0, x), y+y));
398 if (isinf(y))
399 return (CMPLXL(copysignl(0, x), copysignl(pio2_hi + pio2_lo, y)));
400 return (CMPLXL(x+0.0L+(y+0), x+0.0L+(y+0)));
401 }
402
403 if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
404 return (CMPLXL(real_part_reciprocal(x, y), copysignl(pio2_hi + pio2_lo, y)));
405
406 if (ax < SQRT_3_EPSILON/2 && ay < SQRT_3_EPSILON/2) {
407 raise_inexact();
408 return (z);
409 }
410
411 if (ax == 1 && ay < LDBL_EPSILON) {
412 #if 0
413 if (ay > 2*LDBL_MIN)
414 rx = - logl(ay/2) / 2;
415 else
416 #endif
417 rx = - (logl(ay) - m_ln2) / 2;
418 } else
419 rx = log1pl(4*ax / sum_squares(ax-1, ay)) / 4;
420
421 if (ax == 1)
422 ry = atan2l(2, -ay) / 2;
423 else if (ay < LDBL_EPSILON)
424 ry = atan2l(2*ay, (1-ax)*(1+ax)) / 2;
425 else
426 ry = atan2l(2*ay, (1-ax)*(1+ax) - ay*ay) / 2;
427
428 return (CMPLXL(copysignl(rx, x), copysignl(ry, y)));
429 }
430
431 long double complex
catanl(long double complex z)432 catanl(long double complex z)
433 {
434 long double complex w = catanhl(CMPLXL(cimagl(z), creall(z)));
435 return (CMPLXL(cimagl(w), creall(w)));
436 }
437
438 #else
439 __strong_alias(_casinl, casin)
440 __strong_alias(_catanl, catan)
441 __strong_alias(cacoshl, cacosh)
442 __strong_alias(cacosl, cacos)
443 __strong_alias(casinhl, casinh)
444 __strong_alias(catanhl, catanh)
445 #endif
446