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