1 /*-
2 * Copyright (c) 2012 Brendan Fabeny
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
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 *
26 */
27
28 /*
29 * Stopgap replacements for missing standard libm functions. Slower than many
30 * other libm functions, but more accurate than many naive replacements. No
31 * exception handling, and no use of the ternary values (yet).
32 */
33
34 #pragma STDC FENV_ACCESS ON
35
36 #include <sys/param.h>
37
38 #include <complex.h>
39 #include <fenv.h>
40 #include <float.h>
41
42 #if DBL_MANT_DIG == LDBL_MANT_DIG && DBL_MAX_EXP == LDBL_MAX_EXP
43 #define LDBL_IS_DBL
44 #endif
45
46 #ifdef LDBL_IS_DBL
47 #include <math.h>
48 #else
49 #include <stdint.h>
50 #endif
51
52 #if 0
53 #define MPFR_USE_FILE
54 #define MPFR_USE_VA_LIST
55 #endif
56 #include "mpc.h"
57
58 #define FUNCP(lib, constprefix, tonearest, downward, upward, towardzero) \
59 lib ## _rnd_t \
60 getrnd_ ## lib(void) \
61 { \
62 lib ## _rnd_t rnd; \
63 int rndmode; \
64 \
65 rndmode=fegetround(); \
66 \
67 switch (rndmode) { \
68 case FE_TONEAREST: \
69 rnd = constprefix ## tonearest; \
70 case FE_DOWNWARD: \
71 rnd = constprefix ## downward; \
72 case FE_UPWARD: \
73 rnd = constprefix ## upward; \
74 case FE_TOWARDZERO: \
75 rnd = constprefix ## towardzero; \
76 } \
77 \
78 return(rnd); \
79 }
80
FUNCP(mpc,MPC_RND,NN,DD,UU,ZZ)81 FUNCP(mpc, MPC_RND, NN, DD, UU, ZZ)
82 FUNCP(mpfr, MPFR_RND, N, D, U, Z)
83
84 #define prec_flt FLT_MANT_DIG
85 #define prec_d DBL_MANT_DIG
86 #define prec_ld LDBL_MANT_DIG
87
88 /*
89 * Because mpc_get/set_fc are missing from the API, we just cast to doubles
90 * and back
91 */
92
93 #define prec_f prec_d
94 #define mpc_set_fc(mx, x, rnd) mpc_set_dc(mx, (double complex)x, rnd)
95 #define mpc_get_fc(mx, rnd) (float complex)mpc_get_dc(mx,rnd)
96
97 #define init_mpc(typesuffix, mx) \
98 mpc_init3(mx, prec_ ## typesuffix, prec_ ## typesuffix)
99
100 #define assign_mpc(typesuffix, mx, x, rnd) \
101 mpc_set_ ## typesuffix ## c(mx, x, rnd)
102
103 #define convert_mpc(typesuffix, mx, rnd) \
104 mpc_get_ ## typesuffix ## c(mx,rnd)
105
106 #define cleanup_mpc
107
108 #define init_mpfr(typesuffix, mx) \
109 mpfr_init2(mx, prec_ ## typesuffix)
110
111 #define assign_mpfr(typesuffix, mx, x, rnd) \
112 mpfr_set_ ## typesuffix(mx, x, rnd)
113
114 #define convert_mpfr(typesuffix, mx, rnd) \
115 mpfr_get_ ## typesuffix(mx,rnd)
116
117 #define cleanup_mpfr mpfr_free_cache()
118
119 /*
120 * 1 input and 1 output of similar types, using only one of the
121 * auxiliary libraries
122 */
123
124 #define FUNC1A(libmfunc, libmtype, typesuffix, auxlib, auxfunc) \
125 libmtype \
126 libmfunc(libmtype x) \
127 { \
128 int i, j; \
129 libmtype y; \
130 auxlib ## _rnd_t rnd; \
131 auxlib ## _t mx; \
132 \
133 rnd=getrnd_ ## auxlib(); \
134 init_ ## auxlib(typesuffix, mx); \
135 i=assign_ ## auxlib(typesuffix, mx, x, rnd); \
136 j=auxlib ## _ ## auxfunc(mx, mx, rnd); \
137 y=convert_ ## auxlib(typesuffix, mx, rnd); \
138 auxlib ## _clear(mx); \
139 cleanup_ ## auxlib; \
140 \
141 return(y); \
142 }
143
144 /* 1 complex input, 1 real or integer output */
145
146 #define FUNC1B(libmfunc, libmintype, libmouttype, mpctypesuffix, \
147 mpfrtypesuffix, mpcfunc) \
148 libmouttype \
149 libmfunc(libmintype x) \
150 { \
151 int i, j; \
152 libmouttype y; \
153 mpc_rnd_t mpcrnd; \
154 mpc_t mx; \
155 mpfr_rnd_t mpfrrnd; \
156 mpfr_t my; \
157 \
158 mpcrnd=getrnd_mpc(); \
159 mpfrrnd=MPC_RND_RE(mpcrnd); \
160 init_mpc(mpctypesuffix, mx); \
161 init_mpfr(mpfrtypesuffix, my); \
162 i=assign_mpc(mpctypesuffix, mx, x, mpcrnd); \
163 j=mpc ## _ ## mpcfunc(my, mx, mpfrrnd); \
164 y=convert_mpfr(mpfrtypesuffix, my, mpfrrnd); \
165 mpc_clear(mx); \
166 mpfr_clear(my); \
167 cleanup_mpfr; \
168 \
169 return(y); \
170 }
171
172 /*
173 * 1 input and 1 output of different types, using only one of the
174 * auxiliary libraries
175 */
176
177 #define FUNC1C(libmfunc, libmintype, libmouttype, insuffix, outsuffix, \
178 auxlib, auxfunc) \
179 libmouttype \
180 libmfunc(libmintype x) \
181 { \
182 int i, j; \
183 libmouttype y; \
184 auxlib ## _rnd_t rnd; \
185 auxlib ## _t mx; \
186 \
187 rnd=getrnd_ ## auxlib(); \
188 init_ ## auxlib(insuffix, mx); \
189 i=assign_ ## auxlib(insuffix, mx, x, rnd); \
190 j=auxlib ## _ ## auxfunc(mx, mx, rnd); \
191 y=convert_ ## auxlib(outsuffix, mx, rnd); \
192 auxlib ## _clear(mx); \
193 cleanup_ ## auxlib; \
194 \
195 return(y); \
196 }
197
198 /*
199 * 2 inputs and 1 output of similar types, using only one of the
200 * auxiliary libraries
201 */
202
203 #define FUNC2(libmfunc, libmtype, typesuffix, auxlib, auxfunc) \
204 libmtype \
205 libmfunc(libmtype x, libmtype y) \
206 { \
207 int i, j, k; \
208 libmtype z; \
209 auxlib ## _rnd_t rnd; \
210 auxlib ## _t mx, my; \
211 \
212 rnd=getrnd_ ## auxlib(); \
213 init_ ## auxlib(typesuffix, mx); \
214 i=assign_ ## auxlib(typesuffix, mx, x, rnd); \
215 init_ ## auxlib(typesuffix, my); \
216 j=assign_ ## auxlib(typesuffix, my, y, rnd); \
217 k=auxlib ## _ ## auxfunc(mx, mx, my, rnd); \
218 z=convert_ ## auxlib(typesuffix, mx, rnd); \
219 auxlib ## _clear(mx); \
220 auxlib ## _clear(my); \
221 cleanup_ ## auxlib; \
222 \
223 return(z); \
224 }
225
226 /*
227 * cast a single input from a longer to a shorter type, use the existing
228 * libm function for the shorter type to produce a single output of the
229 * shorter type, and then cast the output back to the longer type
230 */
231
232 #define CASTFUNC1A(libmfunc, outcast, auxfunc, incast) \
233 outcast \
234 libmfunc(outcast x) \
235 { \
236 outcast y; \
237 \
238 y=(outcast)auxfunc((incast)x); \
239 \
240 return(y); \
241 }
242
243 /*
244 * cast a single input from a longer to a shorter type, and then use the
245 * existing libm function for the shorter type to produce a single
246 * output of a third type
247 */
248
249 #define CASTFUNC1B(libmfunc, intype, outtype, auxfunc, incast) \
250 outtype \
251 libmfunc(intype x) \
252 { \
253 outtype y; \
254 \
255 y=auxfunc((incast)x); \
256 \
257 return(y); \
258 }
259
260 /*
261 * cast two inputs from the same longer type to the same shorter type,
262 * use the existing libm function for these inputs to produce a single
263 * output of the shorter type, and then cast the output back to the
264 * longer type
265 */
266
267 #define CASTFUNC2(libmfunc, outcast, auxfunc, incast) \
268 outcast \
269 libmfunc(outcast x, outcast y) \
270 { \
271 outcast z; \
272 \
273 z=(outcast)auxfunc((incast)x, (incast)y); \
274 \
275 return(z); \
276 }
277
278 #if __FreeBSD_version < 704102 || \
279 (__FreeBSD_version >= 800000 && __FreeBSD_version < 802502) || \
280 (__FreeBSD_version >= 900000 && __FreeBSD_version < 900027)
281
282 FUNC1A(log2, double, d, mpfr, log2)
283 FUNC1A(log2f, float, flt, mpfr, log2)
284
285 #endif
286 #if __FreeBSD_version < 800007
287
288 /* nan* will require explicit definitions:
289 * nan
290 * nanf
291 * nanl
292 */
293
294 #ifdef LDBL_IS_DBL
295
296 CASTFUNC1A(logbl, long double, logb, double)
297
298 #else
299
300 long double
301 logbl(long double x)
302 {
303 int i;
304 long *exp=0;
305 long double y;
306 mpfr_rnd_t rnd;
307 mpfr_t mx;
308
309 rnd=getrnd_mpfr();
310 mpfr_init2(mx, prec_ld);
311 i=mpfr_set_ld(mx, x, rnd);
312 y=mpfr_get_ld_2exp(exp, mx, rnd);
313 mpfr_clear(mx);
314 mpfr_free_cache();
315
316 return((long double)*exp);
317 }
318
319 #endif
320
321 FUNC1B(carg, double complex, double, d, d, arg)
322 FUNC1B(cargf, float complex, float, f, flt, arg)
323
324 FUNC1A(csqrt, double complex, d, mpc, sqrt)
325 FUNC1A(csqrtf, float complex, f, mpc, sqrt)
326
327 #endif
328 #if __FreeBSD_version < 800012
329 #ifdef LDBL_IS_DBL
330
331 CASTFUNC1B(llrintl, long double, long long, llrint, double)
332 CASTFUNC1B(lrintl, long double, long, lrint, double)
333
334 CASTFUNC1A(nearbyintl, long double, nearbyint, double)
335 CASTFUNC1A(rintl, long double, rint, double)
336 CASTFUNC1A(exp2l, long double, exp2, double)
337
338 #else
339
340 /*
341 * Since we cannot assume intmax_t is long long, and there is
342 * no long long analogue of mpfr_get_sj, use a cast for llrintl
343 */
344
345 FUNC1C(imrintl, long double, intmax_t, ld, sj, mpfr, rint)
346
347 long long
348 llrintl (long double x)
349 {
350 long long y;
351
352 y=(long long)imrintl(x);
353
354 return(y);
355 }
356
357 FUNC1C(lrintl, long double, long, ld, si, mpfr, rint)
358
359 /*
360 * rint* and nearbyint* are conflated: to be strictly
361 * correct, they should be distinguished by different handling
362 * of the inexact flag.
363 */
364
365 FUNC1A(nearbyintl, long double, ld, mpfr, rint)
366 FUNC1A(rintl, long double, ld, mpfr, rint)
367 FUNC1A(exp2l, long double, ld, mpfr, exp2)
368
369 #endif
370 #endif
371 #if __FreeBSD_version < 800022
372 #ifdef LDBL_IS_DBL
373
374 CASTFUNC1A(sinl, long double, sin, double)
375 CASTFUNC1A(cosl, long double, cos, double)
376 CASTFUNC1A(tanl, long double, tan, double)
377
378 #else
379
380 FUNC1A(sinl, long double, ld, mpfr, sin)
381 FUNC1A(cosl, long double, ld, mpfr, cos)
382 FUNC1A(tanl, long double, ld, mpfr, tan)
383
384 #endif
385
386 FUNC1A(tgammaf, float, flt, mpfr, gamma)
387
388 #endif
389 #if __FreeBSD_version < 800025
390 #ifdef LDBL_IS_DBL
391
392 CASTFUNC1A(sqrtl, long double, sqrt, double)
393
394 #else
395
396 FUNC1A(sqrtl, long double, ld, mpfr, sqrt)
397
398 #endif
399 #endif
400 #if __FreeBSD_version < 800030
401 #ifdef LDBL_IS_DBL
402
403 CASTFUNC2(hypotl, long double, hypot, double)
404 CASTFUNC2(remainderl, long double, remainder, double)
405
406 long double
407 remquol(long double x, long double y, int *quo)
408 {
409 long double z;
410
411 z=(long double)remquo((double)x,(double)y, quo);
412
413 return(z);
414 }
415
416 CASTFUNC1B(cabsl, long double complex, long double, cabs, double complex)
417
418 CASTFUNC1A(csqrtl, long double complex, csqrt, double complex)
419
420 #else
421
422 FUNC2(hypotl, long double, ld, mpfr, hypot)
423 FUNC2(remainderl, long double, ld, mpfr, remainder)
424
425 long double
426 remquol(long double x, long double y, int *quo)
427 {
428 int i, j, k;
429 long *q=0;
430 long double z;
431 mpfr_rnd_t rnd;
432 mpfr_t mx, my;
433
434 rnd=getrnd_mpfr();
435 mpfr_inits2(prec_ld, mx, my, (mpfr_ptr)0);
436 i=mpfr_set_ld(mx, x, rnd);
437 j=mpfr_set_ld(my, y, rnd);
438 k=mpfr_remquo(my, q, mx, my, rnd);
439 z=mpfr_get_ld(my, rnd);
440 *quo=(int)*q;
441 mpfr_clears(mx, my, (mpfr_ptr)0);
442 mpfr_free_cache();
443
444 return(z);
445 }
446
447 FUNC1B(cabsl, long double complex, long double, ld, ld, abs)
448
449 FUNC1A(csqrtl, long double complex, ld, mpc, sqrt)
450
451 #endif
452 #endif
453 #if __FreeBSD_version < 800040
454
455 #ifdef LDBL_IS_DBL
456
457 CASTFUNC2(fmodl, long double, fmod, double)
458
459 #else
460
461 FUNC2(fmodl, long double, ld, mpfr, fmod)
462
463 #endif
464 #endif
465 #if __FreeBSD_version < 800042
466 #ifdef LDBL_IS_DBL
467
468 CASTFUNC1A(acosl, long double, acos, double)
469 CASTFUNC1A(asinl, long double, asin, double)
470 CASTFUNC1A(atanl, long double, atan, double)
471
472 CASTFUNC2(atan2l, long double, atan2, double)
473
474 CASTFUNC1B(cargl, long double complex, long double, carg, double complex)
475
476 #else
477
478 FUNC1A(acosl, long double, ld, mpfr, acos)
479 FUNC1A(asinl, long double, ld, mpfr, asin)
480 FUNC1A(atanl, long double, ld, mpfr, atan)
481
482 FUNC2(atan2l, long double, ld, mpfr, atan2)
483
484 FUNC1B(cargl, long double complex, long double, ld, ld, arg)
485
486 #endif
487
488 FUNC1A(cproj, double complex, d, mpc, proj)
489 FUNC1A(cprojf, float complex, f, mpc, proj)
490 FUNC1A(cprojl, long double complex, ld, mpc, proj)
491
492 #endif
493 #if __FreeBSD_version < 900034
494
495 FUNC1A(cexp, double complex, d, mpc, exp)
496 FUNC1A(cexpf, float complex, f, mpc, exp)
497
498 #endif
499 #if __FreeBSD_version < 900035
500 #ifdef LDBL_IS_DBL
501
502 CASTFUNC1A(cbrtl, long double, cbrt, double)
503
504 #else
505
506 FUNC1A(cbrtl, long double, ld, mpfr, cbrt)
507
508 #endif
509 #endif
510 #if __FreeBSD_version < 1000001
511
512 FUNC1A(csin, double complex, d, mpc, sin)
513 FUNC1A(csinf, float complex, f, mpc, sin)
514 FUNC1A(csinh, double complex, d, mpc, sinh)
515 FUNC1A(csinhf, float complex, f, mpc, sinh)
516 FUNC1A(ccos, double complex, d, mpc, cos)
517 FUNC1A(ccosf, float complex, f, mpc, cos)
518 FUNC1A(ccosh, double complex, d, mpc, cosh)
519 FUNC1A(ccoshf, float complex, f, mpc, cosh)
520 FUNC1A(ctan, double complex, d, mpc, tan)
521 FUNC1A(ctanf, float complex, f, mpc, tan)
522 FUNC1A(ctanh, double complex, d, mpc, tanh)
523 FUNC1A(ctanhf, float complex, f, mpc, tanh)
524
525 #endif
526 #if __FreeBSD_version < 1000016
527 #ifdef LDBL_IS_DBL
528
529 CASTFUNC1A(expl, long double, exp, double)
530
531 #else
532
533 FUNC1A(expl, long double, ld, mpfr, exp)
534
535 #endif
536 #endif
537
538 #if __FreeBSD_version < 1000034
539
540 FUNC1A(cacos, double complex, d, mpc, acos)
541 FUNC1A(cacosf, float complex, f, mpc, acos)
542 FUNC1A(cacosh, double complex, d, mpc, acosh)
543 FUNC1A(cacoshf, float complex, f, mpc, acosh)
544 FUNC1A(casin, double complex, d, mpc, asin)
545 FUNC1A(casinf, float complex, f, mpc, asin)
546 FUNC1A(casinh, double complex, d, mpc, asinh)
547 FUNC1A(casinhf, float complex, f, mpc, asinh)
548 FUNC1A(catan, double complex, d, mpc, atan)
549 FUNC1A(catanf, float complex, f, mpc, atan)
550 FUNC1A(catanh, double complex, d, mpc, atanh)
551 FUNC1A(catanhf, float complex, f, mpc, atanh)
552
553 #ifdef LDBL_IS_DBL
554
555 CASTFUNC1A(log10l, long double, log10, double)
556 CASTFUNC1A(log1pl, long double, log1p, double)
557 CASTFUNC1A(log2l, long double, log2, double)
558 CASTFUNC1A(logl, long double, log, double)
559
560 #else
561
562 FUNC1A(log10l, long double, ld, mpfr, log10)
563 FUNC1A(log1pl, long double, ld, mpfr, log1p)
564 FUNC1A(log2l, long double, ld, mpfr, log2)
565 FUNC1A(logl, long double, ld, mpfr, log)
566
567 #endif
568 #endif
569
570 #if __FreeBSD_version < 1000035
571 #ifdef LDBL_IS_DBL
572
573 CASTFUNC1A(expm1l, long double, expm1, double)
574
575 #else
576
577 FUNC1A(expm1l, long double, ld, mpfr, expm1)
578
579 #endif
580 #endif
581
582 #if __FreeBSD_version < 1000037
583 #ifdef LDBL_IS_DBL
584
585 CASTFUNC1A(acoshl, long double, acosh, double)
586 CASTFUNC1A(asinhl, long double, asinh, double)
587 CASTFUNC1A(atanhl, long double, atanh, double)
588
589 #else
590
591 FUNC1A(acoshl, long double, ld, mpfr, acosh)
592 FUNC1A(asinhl, long double, ld, mpfr, asinh)
593 FUNC1A(atanhl, long double, ld, mpfr, atanh)
594
595 #endif
596 #endif
597
598 FUNC1A(clog, double complex, d, mpc, log)
599 FUNC1A(clogf, float complex, f, mpc, log)
600 FUNC1A(clogl, long double complex, ld, mpc, log)
601
602 FUNC2(cpow, double complex, d, mpc, pow)
603 FUNC2(cpowf, float complex, f, mpc, pow)
604 FUNC2(cpowl, long double complex, ld, mpc, pow)
605
606 #ifdef LDBL_IS_DBL
607
608 CASTFUNC1A(cacosl, long double complex, cacos, double complex)
609 CASTFUNC1A(cacoshl, long double complex, cacosh, double complex)
610 CASTFUNC1A(casinl, long double complex, casin, double complex)
611 CASTFUNC1A(casinhl, long double complex, casinh, double complex)
612 CASTFUNC1A(catanl, long double complex, catan, double complex)
613 CASTFUNC1A(catanhl, long double complex, catanh, double complex)
614
615 CASTFUNC1A(ccoshl, long double complex, ccosh, double complex)
616 CASTFUNC1A(ccosl, long double complex, ccos, double complex)
617 CASTFUNC1A(cexpl, long double complex, cexp, double complex)
618
619 CASTFUNC1A(coshl, long double, cosh, double)
620
621 CASTFUNC1A(csinhl, long double complex, csinh, double complex)
622 CASTFUNC1A(csinl, long double complex, csin, double complex)
623 CASTFUNC1A(ctanhl, long double complex, ctanh, double complex)
624 CASTFUNC1A(ctanl, long double complex, ctan, double complex)
625
626 CASTFUNC1A(erfcl, long double, erfc, double)
627 CASTFUNC1A(erfl, long double, erf, double)
628
629 CASTFUNC1A(lgammal, long double, lgamma, double)
630
631 CASTFUNC1A(sinhl, long double, sinh, double)
632 CASTFUNC1A(tanhl, long double, tanh, double)
633 CASTFUNC1A(tgammal, long double, tgamma, double)
634
635 CASTFUNC2(powl, long double, pow, double)
636
637 #else
638
639 FUNC1A(cacosl, long double complex, ld, mpc, acos)
640 FUNC1A(cacoshl, long double complex, ld, mpc, acosh)
641 FUNC1A(casinl, long double complex, ld, mpc, asin)
642 FUNC1A(casinhl, long double complex, ld, mpc, asinh)
643 FUNC1A(catanl, long double complex, ld, mpc, atan)
644 FUNC1A(catanhl, long double complex, ld, mpc, atanh)
645
646 FUNC1A(ccoshl, long double complex, ld, mpc, cosh)
647 FUNC1A(ccosl, long double complex, ld, mpc, cos)
648 FUNC1A(cexpl, long double complex, ld, mpc, exp)
649
650 FUNC1A(coshl, long double, ld, mpfr, cosh)
651
652 FUNC1A(csinhl, long double complex, ld, mpc, sinh)
653 FUNC1A(csinl, long double complex, ld, mpc, sin)
654 FUNC1A(ctanhl, long double complex, ld, mpc, tanh)
655 FUNC1A(ctanl, long double complex, ld, mpc, tan)
656
657 FUNC1A(erfcl, long double, ld, mpfr, erfc)
658 FUNC1A(erfl, long double, ld, mpfr, erf)
659
660 FUNC1A(lgammal, long double, ld, mpfr, lngamma)
661
662 FUNC1A(sinhl, long double, ld, mpfr, sinh)
663 FUNC1A(tanhl, long double, ld, mpfr, tanh)
664 FUNC1A(tgammal, long double, ld, mpfr, gamma)
665
666 FUNC2(powl, long double, ld, mpfr, pow)
667
668 #endif
669