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, 2011 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  }