1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 #pragma weak fmodl = __fmodl
  31 
  32 #include "libm.h"
  33 
  34 static const int
  35         is = -0x7fffffff - 1,
  36         im = 0x0000ffff,
  37         iu = 0x00010000;
  38 
  39 static const long double
  40         zero = 0.0L,
  41         one = 1.0L;
  42 
  43 #ifdef __LITTLE_ENDIAN
  44 #define __H0(x) *(3 + (int *) &x)
  45 #define __H1(x) *(2 + (int *) &x)
  46 #define __H2(x) *(1 + (int *) &x)
  47 #define __H3(x) *(0 + (int *) &x)
  48 #else
  49 #define __H0(x) *(0 + (int *) &x)
  50 #define __H1(x) *(1 + (int *) &x)
  51 #define __H2(x) *(2 + (int *) &x)
  52 #define __H3(x) *(3 + (int *) &x)
  53 #endif
  54 
  55 long double
  56 fmodl(long double x, long double y) {
  57         long double a, b;
  58         int n, ix, iy, k, sx;
  59         int hx;
  60         int x0, y0, z0, carry;
  61         unsigned x1, x2, x3, y1, y2, y3, z1, z2, z3;
  62 
  63         hx = __H0(x);
  64         x1 = __H1(x);
  65         x2 = __H2(x);
  66         x3 = __H3(x);
  67         y0 = __H0(y);
  68         y1 = __H1(y);
  69         y2 = __H2(y);
  70         y3 = __H3(y);
  71 
  72         sx = hx & 0x80000000;
  73         x0 = hx ^ sx;
  74         y0 &= 0x7fffffff;
  75 
  76         /* purge off exception values */
  77         if (x0 >= 0x7fff0000 ||      /* !finitel(x) */
  78             (y0 > 0x7fff0000) || (y0 == 0x7fff0000 && ((y1 | y2 | y3) != 0)) ||
  79             (y0 | y1 | y2 | y3) == 0)   /* isnanl(y) || y = 0 */
  80                 return ((x * y) / (x * y));
  81         a = fabsl(x);
  82         b = fabsl(y);
  83         if (a <= b) {
  84                 if (a < b)
  85                         return (x);
  86                 else
  87                         return (zero * x);
  88         }
  89         /* determine ix = ilogbl(x) */
  90         if (x0 < iu) {               /* subnormal x */
  91                 ix = -16382;
  92                 while (x0 == 0) {
  93                         ix -= 16;
  94                         x0 = x1 >> 16;
  95                         x1 = (x1 << 16) | (x2 >> 16);
  96                         x2 = (x2 << 16) | (x3 >> 16);
  97                         x3 = (x3 << 16);
  98                 }
  99                 while (x0 < iu) {
 100                         ix -= 1;
 101                         x0 = (x0 << 1) | (x1 >> 31);
 102                         x1 = (x1 << 1) | (x2 >> 31);
 103                         x2 = (x2 << 1) | (x3 >> 31);
 104                         x3 <<= 1;
 105                 }
 106         } else {
 107                 ix = (x0 >> 16) - 16383;
 108                 x0 = iu | (x0 & im);
 109         }
 110 
 111         /* determine iy = ilogbl(y) */
 112         if (y0 < iu) {               /* subnormal y */
 113                 iy = -16382;
 114                 while (y0 == 0) {
 115                         iy -= 16;
 116                         y0 = y1 >> 16;
 117                         y1 = (y1 << 16) | (y2 >> 16);
 118                         y2 = (y2 << 16) | (y3 >> 16);
 119                         y3 = (y3 << 16);
 120                 }
 121                 while (y0 < iu) {
 122                         iy -= 1;
 123                         y0 = (y0 << 1) | (y1 >> 31);
 124                         y1 = (y1 << 1) | (y2 >> 31);
 125                         y2 = (y2 << 1) | (y3 >> 31);
 126                         y3 <<= 1;
 127                 }
 128         } else {
 129                 iy = (y0 >> 16) - 16383;
 130                 y0 = iu | (y0 & im);
 131         }
 132 
 133         /* fix point fmod */
 134         n = ix - iy;
 135         while (n--) {
 136                 while (x0 == 0 && n >= 16) {
 137                         n -= 16;
 138                         x0 = x1 >> 16;
 139                         x1 = (x1 << 16) | (x2 >> 16);
 140                         x2 = (x2 << 16) | (x3 >> 16);
 141                         x3 = (x3 << 16);
 142                 }
 143                 while (x0 < iu && n >= 1) {
 144                         n -= 1;
 145                         x0 = (x0 << 1) | (x1 >> 31);
 146                         x1 = (x1 << 1) | (x2 >> 31);
 147                         x2 = (x2 << 1) | (x3 >> 31);
 148                         x3 = (x3 << 1);
 149                 }
 150                 carry = 0;
 151                 z3 = x3 - y3;
 152                 carry = (z3 > x3);
 153                 if (carry == 0) {
 154                         z2 = x2 - y2;
 155                         carry = (z2 > x2);
 156                 } else {
 157                         z2 = x2 - y2 - 1;
 158                         carry = (z2 >= x2);
 159                 }
 160                 if (carry == 0) {
 161                         z1 = x1 - y1;
 162                         carry = (z1 > x1);
 163                 } else {
 164                         z1 = x1 - y1 - 1;
 165                         carry = (z1 >= x1);
 166                 }
 167                 z0 = x0 - y0 - carry;
 168                 if (z0 < 0) {        /* double x */
 169                         x0 = x0 + x0 + ((x1 & is) != 0);
 170                         x1 = x1 + x1 + ((x2 & is) != 0);
 171                         x2 = x2 + x2 + ((x3 & is) != 0);
 172                         x3 = x3 + x3;
 173                 } else {
 174                         if (z0 == 0) {
 175                                 if ((z1 | z2 | z3) == 0) {      /* 0: done */
 176                                         __H0(a) = hx & is;
 177                                         __H1(a) = __H2(a) = __H3(a) = 0;
 178                                         return (a);
 179                                 }
 180                         }
 181                         /* x = z << 1 */
 182                         z0 = z0 + z0 + ((z1 & is) != 0);
 183                         z1 = z1 + z1 + ((z2 & is) != 0);
 184                         z2 = z2 + z2 + ((z3 & is) != 0);
 185                         z3 = z3 + z3;
 186                         x0 = z0;
 187                         x1 = z1;
 188                         x2 = z2;
 189                         x3 = z3;
 190                 }
 191         }
 192 
 193         carry = 0;
 194         z3 = x3 - y3;
 195         carry = (z3 > x3);
 196         if (carry == 0) {
 197                 z2 = x2 - y2;
 198                 carry = (z2 > x2);
 199         } else {
 200                 z2 = x2 - y2 - 1;
 201                 carry = (z2 >= x2);
 202         }
 203         if (carry == 0) {
 204                 z1 = x1 - y1;
 205                 carry = (z1 > x1);
 206         } else {
 207                 z1 = x1 - y1 - 1;
 208                 carry = (z1 >= x1);
 209         }
 210         z0 = x0 - y0 - carry;
 211         if (z0 >= 0) {
 212                 x0 = z0;
 213                 x1 = z1;
 214                 x2 = z2;
 215                 x3 = z3;
 216         }
 217         /* convert back to floating value and restore the sign */
 218         if ((x0 | x1 | x2 | x3) == 0) {
 219                 __H0(a) = hx & is;
 220                 __H1(a) = __H2(a) = __H3(a) = 0;
 221                 return (a);
 222         }
 223         while (x0 < iu) {
 224                 if (x0 == 0) {
 225                         iy -= 16;
 226                         x0 = x1 >> 16;
 227                         x1 = (x1 << 16) | (x2 >> 16);
 228                         x2 = (x2 << 16) | (x3 >> 16);
 229                         x3 = (x3 << 16);
 230                 } else {
 231                         x0 = x0 + x0 + ((x1 & is) != 0);
 232                         x1 = x1 + x1 + ((x2 & is) != 0);
 233                         x2 = x2 + x2 + ((x3 & is) != 0);
 234                         x3 = x3 + x3;
 235                         iy -= 1;
 236                 }
 237         }
 238 
 239         /* normalize output */
 240         if (iy >= -16382) {
 241                 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
 242                 __H1(a) = x1;
 243                 __H2(a) = x2;
 244                 __H3(a) = x3;
 245         } else {                /* subnormal output */
 246                 n = -16382 - iy;
 247                 k = n & 31;
 248                 if (k != 0) {
 249                         if (k <= 16) {
 250                                 x3 = (x2 << (32 - k)) | (x3 >> k);
 251                                 x2 = (x1 << (32 - k)) | (x2 >> k);
 252                                 x1 = (x0 << (32 - k)) | (x1 >> k);
 253                                 x0 >>= k;
 254                         } else {
 255                                 x3 = (x2 << (32 - k)) | (x3 >> k);
 256                                 x2 = (x1 << (32 - k)) | (x2 >> k);
 257                                 x1 = (x0 << (32 - k)) | (x1 >> k);
 258                                 x0 = 0;
 259                         }
 260                 }
 261                 while (n >= 32) {
 262                         n -= 32;
 263                         x3 = x2;
 264                         x2 = x1;
 265                         x1 = x0;
 266                         x0 = 0;
 267                 }
 268                 __H0(a) = x0 | sx;
 269                 __H1(a) = x1;
 270                 __H2(a) = x2;
 271                 __H3(a) = x3;
 272                 a *= one;
 273         }
 274         return (a);
 275 }