OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libquadmath/] [math/] [fmodq.c] - Blame information for rev 801

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 740 jeremybenn 
/* 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 
/*
16 
 * fmodq(x,y)
17 
 * Return x mod y in exact arithmetic
18 
 * Method: shift and subtract
19 
 */
20 
 
21 
#include "quadmath-imp.h"
22 
 
23 
static const __float128 one = 1.0, Zero[] = {0.0, -0.0,};
24 
 
25 
__float128
26 
fmodq (__float128 x, __float128 y)
27 
{
28 
  int64_t n,hx,hy,hz,ix,iy,sx,i;
29 
  uint64_t lx,ly,lz;
30 
 
31 
  GET_FLT128_WORDS64(hx,lx,x);
32 
  GET_FLT128_WORDS64(hy,ly,y);
33 
  sx = hx&0x8000000000000000ULL;        /* sign of x */
34 
  hx ^=sx;                              /* |x| */
35 
  hy &= 0x7fffffffffffffffLL;           /* |y| */
36 
 
37 
  /* purge off exception values */
38 
  if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
39 
    ((hy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* or y is NaN */
40 
      return (x*y)/(x*y);
41 
  if(hx<=hy) {
42 
      if((hx<hy)||(lx<ly)) return x;    /* |x|<|y| return x */
43 
      if(lx==ly)
44 
          return Zero[(uint64_t)sx>>63];        /* |x|=|y| return x*0*/
45 
  }
46 
 
47 
  /* determine ix = ilogb(x) */
48 
  if(hx<0x0001000000000000LL) { /* subnormal x */
49 
      if(hx==0) {
50 
          for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
51 
      } else {
52 
          for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
53 
      }
54 
  } else ix = (hx>>48)-0x3fff;
55 
 
56 
  /* determine iy = ilogb(y) */
57 
      if(hy<0x0001000000000000LL) {     /* subnormal y */
58 
          if(hy==0) {
59 
              for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
60 
          } else {
61 
              for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
62 
          }
63 
      } else iy = (hy>>48)-0x3fff;
64 
 
65 
  /* set up {hx,lx}, {hy,ly} and align y to x */
66 
      if(ix >= -16382)
67 
          hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
68 
      else {            /* subnormal x, shift x to normal */
69 
          n = -16382-ix;
70 
          if(n<=63) {
71 
              hx = (hx<<n)|(lx>>(64-n));
72 
              lx <<= n;
73 
          } else {
74 
              hx = lx<<(n-64);
75 
              lx = 0;
76 
          }
77 
      }
78 
      if(iy >= -16382)
79 
          hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
80 
      else {            /* subnormal y, shift y to normal */
81 
          n = -16382-iy;
82 
          if(n<=63) {
83 
              hy = (hy<<n)|(ly>>(64-n));
84 
              ly <<= n;
85 
          } else {
86 
              hy = ly<<(n-64);
87 
              ly = 0;
88 
          }
89 
      }
90 
 
91 
  /* fix point fmod */
92 
      n = ix - iy;
93 
      while(n--) {
94 
          hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
95 
          if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
96 
          else {
97 
              if((hz|lz)==0)             /* return sign(x)*0 */
98 
                  return Zero[(uint64_t)sx>>63];
99 
              hx = hz+hz+(lz>>63); lx = lz+lz;
100 
          }
101 
      }
102 
      hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
103 
      if(hz>=0) {hx=hz;lx=lz;}
104 
 
105 
  /* convert back to floating value and restore the sign */
106 
      if((hx|lx)==0)                     /* return sign(x)*0 */
107 
          return Zero[(uint64_t)sx>>63];
108 
      while(hx<0x0001000000000000LL) {  /* normalize x */
109 
          hx = hx+hx+(lx>>63); lx = lx+lx;
110 
          iy -= 1;
111 
      }
112 
      if(iy>= -16382) { /* normalize output */
113 
          hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
114 
          SET_FLT128_WORDS64(x,hx|sx,lx);
115 
      } else {          /* subnormal output */
116 
          n = -16382 - iy;
117 
          if(n<=48) {
118 
              lx = (lx>>n)|((uint64_t)hx<<(64-n));
119 
              hx >>= n;
120 
          } else if (n<=63) {
121 
              lx = (hx<<(64-n))|(lx>>n); hx = sx;
122 
          } else {
123 
              lx = hx>>(n-64); hx = sx;
124 
          }
125 
          SET_FLT128_WORDS64(x,hx|sx,lx);
126 
          x *= one;             /* create necessary signal */
127 
      }
128 
      return x;         /* exact output */
129 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.