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 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 #pragma weak __remquol = remquol
  32 
  33 #include "libm.h"
  34 #if defined(__SUNPRO_C)
  35 #include <sunmath.h>                      /* fabsl */
  36 #endif
  37 
  38 static const int is = -0x7fffffff - 1, im = 0x0000ffff, iu = 0x00010000;
  39 static const long double zero = 0.0L, one = 1.0L;
  40 
  41 
  42 #if defined(__sparc)
  43 #define __H0(x)         ((int *)&x)[0]
  44 #define __H1(x)         ((int *)&x)[1]
  45 #define __H2(x)         ((int *)&x)[2]
  46 #define __H3(x)         ((int *)&x)[3]
  47 #else
  48 #error Unsupported architecture
  49 #endif
  50 
  51 /*
  52  * On entrance: *quo is initialized to 0, x finite and y non-zero & ordered
  53  */
  54 static long double
  55 fmodquol(long double x, long double y, int *quo)
  56 {
  57         long double a, b;
  58         int n, ix, iy, k, sx, sq, m;
  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 & is;
  73         sq = (hx ^ y0) & is;
  74         x0 = hx ^ sx;
  75         y0 &= ~0x80000000;
  76 
  77         a = fabsl(x);
  78         b = fabsl(y);
  79 
  80         if (a <= b) {
  81                 if (a < b) {
  82                         return (x);
  83                 } else {
  84                         *quo = 1 + (sq >> 30);
  85                         return (zero * x);
  86                 }
  87         }
  88 
  89         /* determine ix = ilogbl(x) */
  90         if (x0 < iu) {                       /* subnormal x */
  91                 ix = 0;
  92                 ix = -16382;
  93 
  94                 while (x0 == 0) {
  95                         ix -= 16;
  96                         x0 = x1 >> 16;
  97                         x1 = (x1 << 16) | (x2 >> 16);
  98                         x2 = (x2 << 16) | (x3 >> 16);
  99                         x3 = (x3 << 16);
 100                 }
 101 
 102                 while (x0 < iu) {
 103                         ix -= 1;
 104                         x0 = (x0 << 1) | (x1 >> 31);
 105                         x1 = (x1 << 1) | (x2 >> 31);
 106                         x2 = (x2 << 1) | (x3 >> 31);
 107                         x3 <<= 1;
 108                 }
 109         } else {
 110                 ix = (x0 >> 16) - 16383;
 111                 x0 = iu | (x0 & im);
 112         }
 113 
 114         /* determine iy = ilogbl(y) */
 115         if (y0 < iu) {                       /* subnormal y */
 116                 iy = -16382;
 117 
 118                 while (y0 == 0) {
 119                         iy -= 16;
 120                         y0 = y1 >> 16;
 121                         y1 = (y1 << 16) | (y2 >> 16);
 122                         y2 = (y2 << 16) | (y3 >> 16);
 123                         y3 = (y3 << 16);
 124                 }
 125 
 126                 while (y0 < iu) {
 127                         iy -= 1;
 128                         y0 = (y0 << 1) | (y1 >> 31);
 129                         y1 = (y1 << 1) | (y2 >> 31);
 130                         y2 = (y2 << 1) | (y3 >> 31);
 131                         y3 <<= 1;
 132                 }
 133         } else {
 134                 iy = (y0 >> 16) - 16383;
 135                 y0 = iu | (y0 & im);
 136         }
 137 
 138         /* fix point fmod */
 139         n = ix - iy;
 140         m = 0;
 141 
 142         while (n--) {
 143                 while (x0 == 0 && n >= 16) {
 144                         m <<= 16;
 145                         n -= 16;
 146                         x0 = x1 >> 16;
 147                         x1 = (x1 << 16) | (x2 >> 16);
 148                         x2 = (x2 << 16) | (x3 >> 16);
 149                         x3 = (x3 << 16);
 150                 }
 151 
 152                 while (x0 < iu && n >= 1) {
 153                         m += m;
 154                         n -= 1;
 155                         x0 = (x0 << 1) | (x1 >> 31);
 156                         x1 = (x1 << 1) | (x2 >> 31);
 157                         x2 = (x2 << 1) | (x3 >> 31);
 158                         x3 = (x3 << 1);
 159                 }
 160 
 161                 carry = 0;
 162                 z3 = x3 - y3;
 163                 carry = z3 > x3;
 164 
 165                 if (carry == 0) {
 166                         z2 = x2 - y2;
 167                         carry = z2 > x2;
 168                 } else {
 169                         z2 = x2 - y2 - 1;
 170                         carry = z2 >= x2;
 171                 }
 172 
 173                 if (carry == 0) {
 174                         z1 = x1 - y1;
 175                         carry = z1 > x1;
 176                 } else {
 177                         z1 = x1 - y1 - 1;
 178                         carry = z1 >= x1;
 179                 }
 180 
 181                 z0 = x0 - y0 - carry;
 182 
 183                 if (z0 < 0) {                /* double x */
 184                         x0 = x0 + x0 + ((x1 & is) != 0);
 185                         x1 = x1 + x1 + ((x2 & is) != 0);
 186                         x2 = x2 + x2 + ((x3 & is) != 0);
 187                         x3 = x3 + x3;
 188                         m += m;
 189                 } else {
 190                         m += 1;
 191 
 192                         if (z0 == 0) {
 193                                 if ((z1 | z2 | z3) == 0) {
 194                                         /* 0: we are done */
 195                                         if (n < 31)
 196                                                 m <<= (1 + n);
 197                                         else
 198                                                 m = 0;
 199 
 200                                         m &= ~0x80000000;
 201                                         *quo = sq >= 0 ? m : -m;
 202                                         __H0(a) = hx & is;
 203                                         __H1(a) = __H2(a) = __H3(a) = 0;
 204                                         return (a);
 205                                 }
 206                         }
 207 
 208                         /* x = z << 1 */
 209                         z0 = z0 + z0 + ((z1 & is) != 0);
 210                         z1 = z1 + z1 + ((z2 & is) != 0);
 211                         z2 = z2 + z2 + ((z3 & is) != 0);
 212                         z3 = z3 + z3;
 213                         x0 = z0;
 214                         x1 = z1;
 215                         x2 = z2;
 216                         x3 = z3;
 217                         m += m;
 218                 }
 219         }
 220 
 221         carry = 0;
 222         z3 = x3 - y3;
 223         carry = z3 > x3;
 224 
 225         if (carry == 0) {
 226                 z2 = x2 - y2;
 227                 carry = z2 > x2;
 228         } else {
 229                 z2 = x2 - y2 - 1;
 230                 carry = z2 >= x2;
 231         }
 232 
 233         if (carry == 0) {
 234                 z1 = x1 - y1;
 235                 carry = z1 > x1;
 236         } else {
 237                 z1 = x1 - y1 - 1;
 238                 carry = z1 >= x1;
 239         }
 240 
 241         z0 = x0 - y0 - carry;
 242 
 243         if (z0 >= 0) {
 244                 x0 = z0;
 245                 x1 = z1;
 246                 x2 = z2;
 247                 x3 = z3;
 248                 m += 1;
 249         }
 250 
 251         m &= ~0x80000000;
 252         *quo = sq >= 0 ? m : -m;
 253 
 254         /* convert back to floating value and restore the sign */
 255         if ((x0 | x1 | x2 | x3) == 0) {
 256                 __H0(a) = hx & is;
 257                 __H1(a) = __H2(a) = __H3(a) = 0;
 258                 return (a);
 259         }
 260 
 261         while (x0 < iu) {
 262                 if (x0 == 0) {
 263                         iy -= 16;
 264                         x0 = x1 >> 16;
 265                         x1 = (x1 << 16) | (x2 >> 16);
 266                         x2 = (x2 << 16) | (x3 >> 16);
 267                         x3 = (x3 << 16);
 268                 } else {
 269                         x0 = x0 + x0 + ((x1 & is) != 0);
 270                         x1 = x1 + x1 + ((x2 & is) != 0);
 271                         x2 = x2 + x2 + ((x3 & is) != 0);
 272                         x3 = x3 + x3;
 273                         iy -= 1;
 274                 }
 275         }
 276 
 277         /* normalize output */
 278         if (iy >= -16382) {
 279                 __H0(a) = sx | (x0 - iu) | ((iy + 16383) << 16);
 280                 __H1(a) = x1;
 281                 __H2(a) = x2;
 282                 __H3(a) = x3;
 283         } else {                        /* subnormal output */
 284                 n = -16382 - iy;
 285                 k = n & 31;
 286 
 287                 if (k <= 16) {
 288                         x3 = (x2 << (32 - k)) | (x3 >> k);
 289                         x2 = (x1 << (32 - k)) | (x2 >> k);
 290                         x1 = (x0 << (32 - k)) | (x1 >> k);
 291                         x0 >>= k;
 292                 } else {
 293                         x3 = (x2 << (32 - k)) | (x3 >> k);
 294                         x2 = (x1 << (32 - k)) | (x2 >> k);
 295                         x1 = (x0 << (32 - k)) | (x1 >> k);
 296                         x0 = 0;
 297                 }
 298 
 299                 while (n >= 32) {
 300                         n -= 32;
 301                         x3 = x2;
 302                         x2 = x1;
 303                         x1 = x0;
 304                         x0 = 0;
 305                 }
 306 
 307                 __H0(a) = x0 | sx;
 308                 __H1(a) = x1;
 309                 __H2(a) = x2;
 310                 __H3(a) = x3;
 311                 a *= one;
 312         }
 313 
 314         return (a);
 315 }
 316 
 317 long double
 318 remquol(long double x, long double y, int *quo)
 319 {
 320         int hx, hy, sx, sq;
 321         long double v;
 322 
 323         hx = __H0(x);                   /* high word of x */
 324         hy = __H0(y);                   /* high word of y */
 325         sx = hx & is;                       /* sign of x */
 326         sq = (hx ^ hy) & is;                /* sign of x/y */
 327         hx ^= sx;                       /* |x| */
 328         hy &= ~0x80000000;
 329 
 330         /* purge off exception values */
 331         *quo = 0;
 332 
 333         /* y=0, y is NaN, x is NaN or inf */
 334         if (y == 0.0L || y != y || hx >= 0x7fff0000)
 335                 return ((x * y) / (x * y));
 336 
 337         y = fabsl(y);
 338         x = fabsl(x);
 339 
 340         if (hy <= 0x7ffdffff) {
 341                 x = fmodquol(x, y + y, quo);
 342                 *quo = ((*quo) & 0x3fffffff) << 1;
 343         }
 344 
 345         if (hy < 0x00020000) {
 346                 if (x + x > y) {
 347                         *quo += 1;
 348 
 349                         if (x == y)
 350                                 x = zero;
 351                         else
 352                                 x -= y;
 353 
 354                         if (x + x >= y) {
 355                                 x -= y;
 356                                 *quo += 1;
 357                         }
 358                 }
 359         } else {
 360                 v = 0.5L * y;
 361 
 362                 if (x > v) {
 363                         *quo += 1;
 364 
 365                         if (x == y)
 366                                 x = zero;
 367                         else
 368                                 x -= y;
 369 
 370                         if (x >= v) {
 371                                 x -= y;
 372                                 *quo += 1;
 373                         }
 374                 }
 375         }
 376 
 377         if (sq != 0)
 378                 *quo = -(*quo);
 379 
 380         return (sx == 0 ? x : -x);
 381 }