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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak remainder = __remainder 30 31 /* 32 * remainder(x,p) 33 * Code originated from 4.3bsd. 34 * Modified by K.C. Ng for SUN 4.0 libm. 35 * Return : 36 * returns x REM p = x - [x/p]*p as if in infinite precise arithmetic, 37 * where [x/p] is the (inifinite bit) integer nearest x/p (in half way 38 * case choose the even one). 39 * Method : 40 * Based on fmod() return x-[x/p]chopped*p exactly. 41 */ 42 43 #include "libm.h" 44 45 static const double zero = 0.0, half = 0.5; 46 47 double 48 remainder(double x, double p) { 49 double halfp; 50 int ix, hx, hp; 51 52 ix = ((int *)&x)[HIWORD]; 53 hx = ix & ~0x80000000; 54 hp = ((int *)&p)[HIWORD] & ~0x80000000; 55 56 if (hp > 0x7ff00000 || (hp == 0x7ff00000 && ((int *)&p)[LOWORD] != 0)) 57 return (x * p); 58 if (hx > 0x7ff00000 || (hx == 0x7ff00000 && ((int *)&x)[LOWORD] != 0)) 59 return (x * p); 60 61 if ((hp | ((int *)&p)[LOWORD]) == 0 || hx == 0x7ff00000) 62 return (_SVID_libm_err(x, p, 28)); 63 64 p = fabs(p); 65 if (hp < 0x7fe00000) 66 x = fmod(x, p + p); 67 x = fabs(x); 68 if (hp < 0x00200000) { 69 if (x + x > p) { 70 if (x == p) /* avoid x-x=-0 in RM mode */ 71 return ((ix < 0)? -zero : zero); 72 x -= p; 73 if (x + x >= p) 74 x -= p; 75 } 76 } else { 77 halfp = half * p; 78 if (x > halfp) { 79 if (x == p) /* avoid x-x=-0 in RM mode */ 80 return ((ix < 0)? -zero : zero); 81 x -= p; 82 if (x >= halfp) 83 x -= p; 84 } 85 } 86 return ((ix < 0)? -x : x); 87 } 88