usr/src/lib/libbc/libc/gen/common/fmod.c

   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, Version 1.0 only
   6  * (the "License").  You may not use this file except in compliance
   7  * with the License.
   8  *
   9  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
  10  * or http://www.opensolaris.org/os/licensing.
  11  * See the License for the specific language governing permissions
  12  * and limitations under the License.
  13  *
  14  * When distributing Covered Code, include this CDDL HEADER in each
  15  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  16  * If applicable, add the following below this CDDL HEADER, with the
  17  * fields enclosed by brackets "[]" replaced with your own identifying
  18  * information: Portions Copyright [yyyy] [name of copyright owner]
  19  *
  20  * CDDL HEADER END
  21  */
  22 /*
  23  * Copyright 1988 Sun Microsystems, Inc.  All rights reserved.
  24  * Use is subject to license terms.
  25  */
  26
  27 #pragma ident   "%Z%%M% %I%     %E% SMI"
  28
  29 /* Special version adapted from libm for use in libc. */
  30
  31 static int      n0 = 0, n1 = 1;
  32
  33 static double   two52 = 4.503599627370496000E+15;
  34 static double   twom52 = 2.220446049250313081E-16;
  35
  36 static double
  37 setexception(int n, double x)
  38 {
  39         return (0.0);
  40 }
  41
  42 double
  43 copysign(double x, double y)
  44 {
  45         long           *px = (long *) &x;
  46         long           *py = (long *) &y;
  47         px[n0] = (px[n0] & 0x7fffffff) | (py[n0] & 0x80000000);
  48         return (x);
  49 }
  50
  51 static double
  52 fabs(double x)
  53 {
  54         long           *px = (long *) &x;
  55         px[0] &= 0x7fffffff;
  56
  57         return (x);
  58 }
  59
  60 static int
  61 finite(double x)
  62 {
  63         long           *px = (long *) &x;
  64         return ((px[n0] & 0x7ff00000) != 0x7ff00000);
  65 }
  66
  67 static int
  68 ilogb(double x)
  69 {
  70         long           *px = (long *) &x, k;
  71         k = px[n0] & 0x7ff00000;
  72         if (k == 0) {
  73                 if ((px[n1] | (px[n0] & 0x7fffffff)) == 0)
  74                         return (0x80000001);
  75                 else {
  76                         x *= two52;
  77                         return ((px[n0] & 0x7ff00000) >> 20) - 1075;
  78                 }
  79         } else if (k != 0x7ff00000)
  80                 return (k >> 20) - 1023;
  81         else
  82                 return (0x7fffffff);
  83 }
  84
  85 static double
  86 scalbn(double x, int n)
  87 {
  88         long           *px = (long *) &x, k;
  89         double          twom54 = twom52 * 0.25;
  90         k = (px[n0] & 0x7ff00000) >> 20;
  91         if (k == 0x7ff)
  92                 return (x + x);
  93         if ((px[n1] | (px[n0] & 0x7fffffff)) == 0)
  94                 return (x);
  95         if (k == 0) {
  96                 x *= two52;
  97                 k = ((px[n0] & 0x7ff00000) >> 20) - 52;
  98         }
  99         k = k + n;
 100         if (n > 5000)
 101                 return (setexception(2, x));
 102         if (n < -5000)
 103                 return (setexception(1, x));
 104         if (k > 0x7fe)
 105                 return (setexception(2, x));
 106         if (k <= -54)
 107                 return (setexception(1, x));
 108         if (k > 0) {
 109                 px[n0] = (px[n0] & 0x800fffff) | (k << 20);
 110                 return (x);
 111         }
 112         k += 54;
 113         px[n0] = (px[n0] & 0x800fffff) | (k << 20);
 114         return (x * twom54);
 115 }
 116
 117 double
 118 fmod(double x, double y)
 119 {
 120         int             ny, nr;
 121         double          r, z, w;
 122
 123         int finite(), ilogb();
 124         double fabs(), scalbn(), copysign();
 125
 126         /* purge off exception values */
 127         if (!finite(x) || y != y || y == 0.0) {
 128                 return ((x * y) / (x * y));
 129         }
 130         /* scale and subtract to get the remainder */
 131         r = fabs(x);
 132         y = fabs(y);
 133         ny = ilogb(y);
 134         while (r >= y) {
 135                 nr = ilogb(r);
 136                 if (nr == ny)
 137                         w = y;
 138                 else {
 139                         z = scalbn(y, nr - ny - 1);
 140                         w = z + z;
 141                 }
 142                 if (r >= w)
 143                         r -= w;
 144                 else
 145                         r -= z;
 146         }
 147
 148         /* restore sign */
 149         return (copysign(r, x));
 150 }