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 #pragma ident "%Z%%M% %I% %E% SMI" 23 24 /* 25 * Copyright (c) 1988 by Sun Microsystems, Inc. 26 */ 27 28 /* Special version adapted from libm for use in libc. */ 29 30 #ifdef i386 31 static n0 = 1, n1 = 0; 32 #else 33 static n0 = 0, n1 = 1; 34 #endif 35 36 static double two52 = 4.503599627370496000E+15; 37 static double twom52 = 2.220446049250313081E-16; 38 39 static double 40 setexception(n, x) 41 int n; 42 double x; 43 { 44 } 45 46 double 47 copysign(x, y) 48 double x, y; 49 { 50 long *px = (long *) &x; 51 long *py = (long *) &y; 52 px[n0] = (px[n0] & 0x7fffffff) | (py[n0] & 0x80000000); 53 return x; 54 } 55 56 static double 57 fabs(x) 58 double x; 59 { 60 long *px = (long *) &x; 61 #ifdef i386 62 px[1] &= 0x7fffffff; 63 #else 64 px[0] &= 0x7fffffff; 65 #endif 66 return x; 67 } 68 69 static int 70 finite(x) 71 double x; 72 { 73 long *px = (long *) &x; 74 return ((px[n0] & 0x7ff00000) != 0x7ff00000); 75 } 76 77 static int 78 ilogb(x) 79 double x; 80 { 81 long *px = (long *) &x, k; 82 k = px[n0] & 0x7ff00000; 83 if (k == 0) { 84 if ((px[n1] | (px[n0] & 0x7fffffff)) == 0) 85 return 0x80000001; 86 else { 87 x *= two52; 88 return ((px[n0] & 0x7ff00000) >> 20) - 1075; 89 } 90 } else if (k != 0x7ff00000) 91 return (k >> 20) - 1023; 92 else 93 return 0x7fffffff; 94 } 95 96 static double 97 scalbn(x, n) 98 double x; 99 int n; 100 { 101 long *px = (long *) &x, k; 102 double twom54 = twom52 * 0.25; 103 k = (px[n0] & 0x7ff00000) >> 20; 104 if (k == 0x7ff) 105 return x + x; 106 if ((px[n1] | (px[n0] & 0x7fffffff)) == 0) 107 return x; 108 if (k == 0) { 109 x *= two52; 110 k = ((px[n0] & 0x7ff00000) >> 20) - 52; 111 } 112 k = k + n; 113 if (n > 5000) 114 return setexception(2, x); 115 if (n < -5000) 116 return setexception(1, x); 117 if (k > 0x7fe) 118 return setexception(2, x); 119 if (k <= -54) 120 return setexception(1, x); 121 if (k > 0) { 122 px[n0] = (px[n0] & 0x800fffff) | (k << 20); 123 return x; 124 } 125 k += 54; 126 px[n0] = (px[n0] & 0x800fffff) | (k << 20); 127 return x * twom54; 128 } 129 130 double 131 fmod(x, y) 132 double x, y; 133 { 134 int ny, nr; 135 double r, z, w; 136 int finite(), ilogb(); 137 double fabs(), scalbn(), copysign(); 138 139 /* purge off exception values */ 140 if (!finite(x) || y != y || y == 0.0) { 141 return (x * y) / (x * y); 142 } 143 /* scale and subtract to get the remainder */ 144 r = fabs(x); 145 y = fabs(y); 146 ny = ilogb(y); 147 while (r >= y) { 148 nr = ilogb(r); 149 if (nr == ny) 150 w = y; 151 else { 152 z = scalbn(y, nr - ny - 1); 153 w = z + z; 154 } 155 if (r >= w) 156 r -= w; 157 else 158 r -= z; 159 } 160 161 /* restore sign */ 162 return copysign(r, x); 163 } 164