1 | /* e_fmodl.c -- long double version of e_fmod.c. |
2 | */ |
3 | /* |
4 | * ==================================================== |
5 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
6 | * |
7 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
8 | * Permission to use, copy, modify, and distribute this |
9 | * software is freely granted, provided that this notice |
10 | * is preserved. |
11 | * ==================================================== |
12 | */ |
13 | |
14 | /* __ieee754_remainderl(x,p) |
15 | * Return : |
16 | * returns x REM p = x - [x/p]*p as if in infinite |
17 | * precise arithmetic, where [x/p] is the (infinite bit) |
18 | * integer nearest x/p (in half way case choose the even one). |
19 | * Method : |
20 | * Based on fmodl() return x-[x/p]chopped*p exactlp. |
21 | */ |
22 | |
23 | #include <math.h> |
24 | #include <math_private.h> |
25 | #include <libm-alias-finite.h> |
26 | |
27 | static const long double zero = 0.0L; |
28 | |
29 | |
30 | long double |
31 | __ieee754_remainderl(long double x, long double p) |
32 | { |
33 | int64_t hx,hp; |
34 | uint64_t sx,lx,lp; |
35 | long double p_half; |
36 | double xhi, xlo, phi, plo; |
37 | |
38 | ldbl_unpack (x, &xhi, &xlo); |
39 | EXTRACT_WORDS64 (hx, xhi); |
40 | EXTRACT_WORDS64 (lx, xlo); |
41 | ldbl_unpack (p, &phi, &plo); |
42 | EXTRACT_WORDS64 (hp, phi); |
43 | EXTRACT_WORDS64 (lp, plo); |
44 | sx = hx&0x8000000000000000ULL; |
45 | lp ^= hp & 0x8000000000000000ULL; |
46 | hp &= 0x7fffffffffffffffLL; |
47 | lx ^= sx; |
48 | hx &= 0x7fffffffffffffffLL; |
49 | if (lp == 0x8000000000000000ULL) |
50 | lp = 0; |
51 | if (lx == 0x8000000000000000ULL) |
52 | lx = 0; |
53 | |
54 | /* purge off exception values */ |
55 | if(hp==0) return (x*p)/(x*p); /* p = 0 */ |
56 | if((hx>=0x7ff0000000000000LL)|| /* x not finite */ |
57 | (hp>0x7ff0000000000000LL)) /* p is NaN */ |
58 | return (x*p)/(x*p); |
59 | |
60 | |
61 | if (hp<=0x7fdfffffffffffffLL) x = __ieee754_fmodl(x,p+p); /* now x < 2p */ |
62 | if (((hx-hp)|(lx-lp))==0) return zero*x; |
63 | x = fabsl(x: x); |
64 | p = fabsl(x: p); |
65 | if (hp<0x0020000000000000LL) { |
66 | if(x+x>p) { |
67 | x-=p; |
68 | if(x+x>=p) x -= p; |
69 | } |
70 | } else { |
71 | p_half = 0.5L*p; |
72 | if(x>p_half) { |
73 | x-=p; |
74 | if(x>=p_half) x -= p; |
75 | } |
76 | } |
77 | if (sx) |
78 | x = -x; |
79 | return x; |
80 | } |
81 | libm_alias_finite (__ieee754_remainderl, __remainderl) |
82 | |