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 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak remainderl = __remainderl 31 32 #include "libm.h" 33 #include "longdouble.h" 34 35 /* 36 * remainderl(x,p) 37 * returns x REM p = x - [x/p]*p as if in infinite 38 * precise arithmetic, where [x/p] is the (inifinite bit) 39 * integer nearest x/p (in half way case choose the even one). 40 * Method : 41 * Based on fmodl() return x-[x/p]chopped*p exactly. 42 */ 43 44 #define HFMAX 5.948657476786158825428796633140035080982e+4931L 45 #define DBMIN 6.724206286224187012525355634643505205196e-4932L 46 47 static const long double 48 zero = 0.0L, 49 half = 0.5L, 50 hfmax = HFMAX, /* half of the maximum number */ 51 dbmin = DBMIN; /* double of the minimum (normal) number */ 52 53 long double 54 remainderl(long double x, long double p) { 55 long double hp; 56 int sx; 57 58 if (isnanl(p)) 59 return (x + p); 60 if (!finitel(x)) 61 return (x - x); 62 p = fabsl(p); 63 if (p <= hfmax) 64 x = fmodl(x, p + p); 65 sx = signbitl(x); 66 x = fabsl(x); 67 if (p < dbmin) { 68 if (x + x > p) { 69 if (x == p) 70 x = zero; 71 else 72 x -= p; /* avoid x-x=-0 in RM mode */ 73 if (x + x >= p) 74 x -= p; 75 } 76 } else { 77 hp = half * p; 78 if (x > hp) { 79 if (x == p) 80 x = zero; 81 else 82 x -= p; /* avoid x-x=-0 in RM mode */ 83 if (x >= hp) 84 x -= p; 85 } 86 } 87 return (sx == 0 ? x : -x); 88 } 89