POK
e_remainder.c
1 /*
2  * POK header
3  *
4  * The following file is a part of the POK project. Any modification should
5  * made according to the POK licence. You CANNOT use this file or a part of
6  * this file is this part of a file for your own project
7  *
8  * For more information on the POK licence, please see our LICENCE FILE
9  *
10  * Please follow the coding guidelines described in doc/CODING_GUIDELINES
11  *
12  * Copyright (c) 2007-2009 POK team
13  *
14  * Created by julien on Fri Jan 30 14:41:34 2009
15  */
16 
17 /* @(#)e_remainder.c 5.1 93/09/24 */
18 /*
19  * ====================================================
20  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
21  *
22  * Developed at SunPro, a Sun Microsystems, Inc. business.
23  * Permission to use, copy, modify, and distribute this
24  * software is freely granted, provided that this notice
25  * is preserved.
26  * ====================================================
27  */
28 
29 
30 /* __ieee754_remainder(x,p)
31  * Return :
32  * returns x REM p = x - [x/p]*p as if in infinite
33  * precise arithmetic, where [x/p] is the (infinite bit)
34  * integer nearest x/p (in half way case choose the even one).
35  * Method :
36  * Based on fmod() return x-[x/p]chopped*p exactlp.
37  */
38 
39 #ifdef POK_NEEDS_LIBMATH
40 #include <libm.h>
41 #include "math_private.h"
42 
43 static const double zero = 0.0;
44 
45 
46 double
47 __ieee754_remainder(double x, double p)
48 {
49  int32_t hx,hp;
50  uint32_t sx,lx,lp;
51  double p_half;
52 
53  EXTRACT_WORDS(hx,lx,x);
54  EXTRACT_WORDS(hp,lp,p);
55  sx = hx&0x80000000;
56  hp &= 0x7fffffff;
57  hx &= 0x7fffffff;
58 
59  /* purge off exception values */
60  if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
61  if((hx>=0x7ff00000)|| /* x not finite */
62  ((hp>=0x7ff00000)&& /* p is NaN */
63  (((hp-0x7ff00000)|lp)!=0)))
64  return (x*p)/(x*p);
65 
66 
67  if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */
68  if (((hx-hp)|(lx-lp))==0) return zero*x;
69  x = fabs(x);
70  p = fabs(p);
71  if (hp<0x00200000) {
72  if(x+x>p) {
73  x-=p;
74  if(x+x>=p) x -= p;
75  }
76  } else {
77  p_half = 0.5*p;
78  if(x>p_half) {
79  x-=p;
80  if(x>=p_half) x -= p;
81  }
82  }
83  GET_HIGH_WORD(hx,x);
84  SET_HIGH_WORD(x,hx^sx);
85  return x;
86 }
87 #endif
88