1179882Sdas/* @(#)e_fmod.c 1.3 95/01/18 */ 2179882Sdas/*- 3179882Sdas * ==================================================== 4179882Sdas * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5179882Sdas * 6179882Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 7179882Sdas * Permission to use, copy, modify, and distribute this 8179882Sdas * software is freely granted, provided that this notice 9179882Sdas * is preserved. 10179882Sdas * ==================================================== 11179882Sdas */ 12179882Sdas 13179882Sdas#include <sys/cdefs.h> 14179882Sdas__FBSDID("$FreeBSD$"); 15179882Sdas 16179882Sdas#include <float.h> 17179882Sdas#include <stdint.h> 18179882Sdas 19179882Sdas#include "fpmath.h" 20179882Sdas#include "math.h" 21179882Sdas#include "math_private.h" 22179882Sdas 23179882Sdas#define BIAS (LDBL_MAX_EXP - 1) 24179882Sdas 25179882Sdas#if LDBL_MANL_SIZE > 32 26179882Sdastypedef uint64_t manl_t; 27179882Sdas#else 28179882Sdastypedef uint32_t manl_t; 29179882Sdas#endif 30179882Sdas 31179882Sdas#if LDBL_MANH_SIZE > 32 32179882Sdastypedef uint64_t manh_t; 33179882Sdas#else 34179882Sdastypedef uint32_t manh_t; 35179882Sdas#endif 36179882Sdas 37179882Sdas/* 38179882Sdas * These macros add and remove an explicit integer bit in front of the 39179882Sdas * fractional mantissa, if the architecture doesn't have such a bit by 40179882Sdas * default already. 41179882Sdas */ 42179882Sdas#ifdef LDBL_IMPLICIT_NBIT 43179882Sdas#define SET_NBIT(hx) ((hx) | (1ULL << LDBL_MANH_SIZE)) 44179882Sdas#define HFRAC_BITS LDBL_MANH_SIZE 45179882Sdas#else 46179882Sdas#define SET_NBIT(hx) (hx) 47179882Sdas#define HFRAC_BITS (LDBL_MANH_SIZE - 1) 48179882Sdas#endif 49179882Sdas 50179882Sdas#define MANL_SHIFT (LDBL_MANL_SIZE - 1) 51179882Sdas 52179882Sdasstatic const long double one = 1.0, Zero[] = {0.0, -0.0,}; 53179882Sdas 54179882Sdas/* 55179882Sdas * fmodl(x,y) 56179882Sdas * Return x mod y in exact arithmetic 57179882Sdas * Method: shift and subtract 58179882Sdas * 59179882Sdas * Assumptions: 60179882Sdas * - The low part of the mantissa fits in a manl_t exactly. 61179882Sdas * - The high part of the mantissa fits in an int64_t with enough room 62179882Sdas * for an explicit integer bit in front of the fractional bits. 63179882Sdas */ 64179882Sdaslong double 65179882Sdasfmodl(long double x, long double y) 66179882Sdas{ 67179882Sdas union IEEEl2bits ux, uy; 68179882Sdas int64_t hx,hz; /* We need a carry bit even if LDBL_MANH_SIZE is 32. */ 69179882Sdas manh_t hy; 70179882Sdas manl_t lx,ly,lz; 71179882Sdas int ix,iy,n,sx; 72179882Sdas 73179882Sdas ux.e = x; 74179882Sdas uy.e = y; 75179882Sdas sx = ux.bits.sign; 76179882Sdas 77179882Sdas /* purge off exception values */ 78179882Sdas if((uy.bits.exp|uy.bits.manh|uy.bits.manl)==0 || /* y=0 */ 79179882Sdas (ux.bits.exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */ 80179882Sdas (uy.bits.exp == BIAS + LDBL_MAX_EXP && 81179882Sdas ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)!=0)) /* or y is NaN */ 82179882Sdas return (x*y)/(x*y); 83179882Sdas if(ux.bits.exp<=uy.bits.exp) { 84179882Sdas if((ux.bits.exp<uy.bits.exp) || 85179882Sdas (ux.bits.manh<=uy.bits.manh && 86179882Sdas (ux.bits.manh<uy.bits.manh || 87179882Sdas ux.bits.manl<uy.bits.manl))) { 88179882Sdas return x; /* |x|<|y| return x or x-y */ 89179882Sdas } 90179882Sdas if(ux.bits.manh==uy.bits.manh && ux.bits.manl==uy.bits.manl) { 91179882Sdas return Zero[sx]; /* |x|=|y| return x*0*/ 92179882Sdas } 93179882Sdas } 94179882Sdas 95179882Sdas /* determine ix = ilogb(x) */ 96179882Sdas if(ux.bits.exp == 0) { /* subnormal x */ 97179882Sdas ux.e *= 0x1.0p512; 98179882Sdas ix = ux.bits.exp - (BIAS + 512); 99179882Sdas } else { 100179882Sdas ix = ux.bits.exp - BIAS; 101179882Sdas } 102179882Sdas 103179882Sdas /* determine iy = ilogb(y) */ 104179882Sdas if(uy.bits.exp == 0) { /* subnormal y */ 105179882Sdas uy.e *= 0x1.0p512; 106179882Sdas iy = uy.bits.exp - (BIAS + 512); 107179882Sdas } else { 108179882Sdas iy = uy.bits.exp - BIAS; 109179882Sdas } 110179882Sdas 111179882Sdas /* set up {hx,lx}, {hy,ly} and align y to x */ 112179882Sdas hx = SET_NBIT(ux.bits.manh); 113179882Sdas hy = SET_NBIT(uy.bits.manh); 114179882Sdas lx = ux.bits.manl; 115179882Sdas ly = uy.bits.manl; 116179882Sdas 117179882Sdas /* fix point fmod */ 118179882Sdas n = ix - iy; 119179882Sdas 120179882Sdas while(n--) { 121179882Sdas hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 122179882Sdas if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;} 123179882Sdas else { 124179882Sdas if ((hz|lz)==0) /* return sign(x)*0 */ 125179882Sdas return Zero[sx]; 126179882Sdas hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz; 127179882Sdas } 128179882Sdas } 129179882Sdas hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 130179882Sdas if(hz>=0) {hx=hz;lx=lz;} 131179882Sdas 132179882Sdas /* convert back to floating value and restore the sign */ 133179882Sdas if((hx|lx)==0) /* return sign(x)*0 */ 134179882Sdas return Zero[sx]; 135181063Sdas while(hx<(1ULL<<HFRAC_BITS)) { /* normalize x */ 136179882Sdas hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx; 137179882Sdas iy -= 1; 138179882Sdas } 139179882Sdas ux.bits.manh = hx; /* The mantissa is truncated here if needed. */ 140179882Sdas ux.bits.manl = lx; 141179882Sdas if (iy < LDBL_MIN_EXP) { 142179882Sdas ux.bits.exp = iy + (BIAS + 512); 143179882Sdas ux.e *= 0x1p-512; 144179882Sdas } else { 145179882Sdas ux.bits.exp = iy + BIAS; 146179882Sdas } 147179882Sdas x = ux.e * one; /* create necessary signal */ 148179882Sdas return x; /* exact output */ 149179882Sdas} 150