sysdeps/ieee754/ldbl-128ibm/e_fmodl.c

   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 /*
  16  * __ieee754_fmodl(x,y)
  17  * Return x mod y in exact arithmetic
  18  * Method: shift and subtract
  19  */
  20
  21 #include "math.h"
  22 #include "math_private.h"
  23 #include <ieee754.h>
  24
  25 #ifdef __STDC__
  26 static const long double one = 1.0, Zero[] = {0.0, -0.0,};
  27 #else
  28 static long double one = 1.0, Zero[] = {0.0, -0.0,};
  29 #endif
  30
  31 #ifdef __STDC__
  32         long double __ieee754_fmodl(long double x, long double y)
  33 #else
  34         long double __ieee754_fmodl(x,y)
  35         long double x,y;
  36 #endif
  37 {
  38         int64_t n,hx,hy,hz,ix,iy,sx,i;
  39         u_int64_t lx,ly,lz;
  40         int temp;
  41
  42         GET_LDOUBLE_WORDS64(hx,lx,x);
  43         GET_LDOUBLE_WORDS64(hy,ly,y);
  44         sx = hx&0x8000000000000000ULL;          /* sign of x */
  45         hx ^=sx;                                /* |x| */
  46         hy &= 0x7fffffffffffffffLL;             /* |y| */
  47
  48     /* purge off exception values */
  49         if((hy|(ly&0x7fffffffffffffff))==0||(hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
  50           (hy>0x7ff0000000000000LL))    /* or y is NaN */
  51             return (x*y)/(x*y);
  52         if(hx<=hy) {
  53             if((hx<hy)||(lx<ly)) return x;      /* |x|<|y| return x */
  54             if(lx==ly)
  55                 return Zero[(u_int64_t)sx>>63]; /* |x|=|y| return x*0*/
  56         }
  57
  58     /* determine ix = ilogb(x) */
  59         if(hx<0x0010000000000000LL) {   /* subnormal x */
  60             if(hx==0) {
  61                 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
  62             } else {
  63                 for (ix = -1022, i=hx<<19; i>0; i<<=1) ix -=1;
  64             }
  65         } else ix = (hx>>52)-0x3ff;
  66
  67     /* determine iy = ilogb(y) */
  68         if(hy<0x0010000000000000LL) {   /* subnormal y */
  69             if(hy==0) {
  70                 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
  71             } else {
  72                 for (iy = -1022, i=hy<<19; i>0; i<<=1) iy -=1;
  73             }
  74         } else iy = (hy>>52)-0x3ff;
  75
  76     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
  77        bit mantissa so the following operatations will give the correct
  78        result.  */
  79         ldbl_extract_mantissa(&hx, &lx, &temp, x);
  80         ldbl_extract_mantissa(&hy, &ly, &temp, y);
  81
  82     /* set up {hx,lx}, {hy,ly} and align y to x */
  83         if(ix >= -1022)
  84             hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
  85         else {          /* subnormal x, shift x to normal */
  86             n = -1022-ix;
  87             if(n<=63) {
  88                 hx = (hx<<n)|(lx>>(64-n));
  89                 lx <<= n;
  90             } else {
  91                 hx = lx<<(n-64);
  92                 lx = 0;
  93             }
  94         }
  95         if(iy >= -1022)
  96             hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
  97         else {          /* subnormal y, shift y to normal */
  98             n = -1022-iy;
  99             if(n<=63) {
 100                 hy = (hy<<n)|(ly>>(64-n));
 101                 ly <<= n;
 102             } else {
 103                 hy = ly<<(n-64);
 104                 ly = 0;
 105             }
 106         }
 107
 108     /* fix point fmod */
 109         n = ix - iy;
 110         while(n--) {
 111             hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 112             if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
 113             else {
 114                 if((hz|(lz&0x7fffffffffffffff))==0)             /* return sign(x)*0 */
 115                     return Zero[(u_int64_t)sx>>63];
 116                 hx = hz+hz+(lz>>63); lx = lz+lz;
 117             }
 118         }
 119         hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 120         if(hz>=0) {hx=hz;lx=lz;}
 121
 122     /* convert back to floating value and restore the sign */
 123         if((hx|(lx&0x7fffffffffffffff))==0)                     /* return sign(x)*0 */
 124             return Zero[(u_int64_t)sx>>63];
 125         while(hx<0x0001000000000000LL) {        /* normalize x */
 126             hx = hx+hx+(lx>>63); lx = lx+lx;
 127             iy -= 1;
 128         }
 129         if(iy>= -1022) {        /* normalize output */
 130             x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
 131         } else {                /* subnormal output */
 132             n = -1022 - iy;
 133             if(n<=48) {
 134                 lx = (lx>>n)|((u_int64_t)hx<<(64-n));
 135                 hx >>= n;
 136             } else if (n<=63) {
 137                 lx = (hx<<(64-n))|(lx>>n); hx = sx;
 138             } else {
 139                 lx = hx>>(n-64); hx = sx;
 140             }
 141             x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
 142             x *= one;           /* create necessary signal */
 143         }
 144         return x;               /* exact output */
 145 }