Home | History | Annotate | Download | only in common
      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 }
    151