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