POK
e_fmod.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_fmod.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_fmod(x,y)
31  * Return x mod y in exact arithmetic
32  * Method: shift and subtract
33  */
34 
35 #ifdef POK_NEEDS_LIBMATH
36 
37 #include "math_private.h"
38 
39 static const double one = 1.0, Zero[] = {0.0, -0.0,};
40 
41 double
42 __ieee754_fmod(double x, double y)
43 {
44  int32_t n,hx,hy,hz,ix,iy,sx,i;
45  uint32_t lx,ly,lz;
46 
47  EXTRACT_WORDS(hx,lx,x);
48  EXTRACT_WORDS(hy,ly,y);
49  sx = hx&0x80000000; /* sign of x */
50  hx ^=sx; /* |x| */
51  hy &= 0x7fffffff; /* |y| */
52 
53  /* purge off exception values */
54  if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
55  ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
56  return (x*y)/(x*y);
57  if(hx<=hy) {
58  if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
59  if(lx==ly)
60  return Zero[(uint32_t)sx>>31]; /* |x|=|y| return x*0*/
61  }
62 
63  /* determine ix = ilogb(x) */
64  if(hx<0x00100000) { /* subnormal x */
65  if(hx==0) {
66  for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
67  } else {
68  for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
69  }
70  } else ix = (hx>>20)-1023;
71 
72  /* determine iy = ilogb(y) */
73  if(hy<0x00100000) { /* subnormal y */
74  if(hy==0) {
75  for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
76  } else {
77  for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
78  }
79  } else iy = (hy>>20)-1023;
80 
81  /* set up {hx,lx}, {hy,ly} and align y to x */
82  if(ix >= -1022)
83  hx = 0x00100000|(0x000fffff&hx);
84  else { /* subnormal x, shift x to normal */
85  n = -1022-ix;
86  if(n<=31) {
87  hx = (hx<<n)|(lx>>(32-n));
88  lx <<= n;
89  } else {
90  hx = lx<<(n-32);
91  lx = 0;
92  }
93  }
94  if(iy >= -1022)
95  hy = 0x00100000|(0x000fffff&hy);
96  else { /* subnormal y, shift y to normal */
97  n = -1022-iy;
98  if(n<=31) {
99  hy = (hy<<n)|(ly>>(32-n));
100  ly <<= n;
101  } else {
102  hy = ly<<(n-32);
103  ly = 0;
104  }
105  }
106 
107  /* fix point fmod */
108  n = ix - iy;
109  while(n--) {
110  hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
111  if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
112  else {
113  if((hz|lz)==0) /* return sign(x)*0 */
114  return Zero[(uint32_t)sx>>31];
115  hx = hz+hz+(lz>>31); lx = lz+lz;
116  }
117  }
118  hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
119  if(hz>=0) {hx=hz;lx=lz;}
120 
121  /* convert back to floating value and restore the sign */
122  if((hx|lx)==0) /* return sign(x)*0 */
123  return Zero[(uint32_t)sx>>31];
124  while(hx<0x00100000) { /* normalize x */
125  hx = hx+hx+(lx>>31); lx = lx+lx;
126  iy -= 1;
127  }
128  if(iy>= -1022) { /* normalize output */
129  hx = ((hx-0x00100000)|((iy+1023)<<20));
130  INSERT_WORDS(x,hx|sx,lx);
131  } else { /* subnormal output */
132  n = -1022 - iy;
133  if(n<=20) {
134  lx = (lx>>n)|((uint32_t)hx<<(32-n));
135  hx >>= n;
136  } else if (n<=31) {
137  lx = (hx<<(32-n))|(lx>>n); hx = sx;
138  } else {
139  lx = hx>>(n-32); hx = sx;
140  }
141  INSERT_WORDS(x,hx|sx,lx);
142  x *= one; /* create necessary signal */
143  }
144  return x; /* exact output */
145 }
146 #endif
147