125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis /*
2225c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2325c28e83SPiotr Jasiukajtis */
2425c28e83SPiotr Jasiukajtis /*
2525c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
2625c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2725c28e83SPiotr Jasiukajtis */
2825c28e83SPiotr Jasiukajtis
29*ddc0e0b5SRichard Lowe #pragma weak __remainder = remainder
3025c28e83SPiotr Jasiukajtis
3125c28e83SPiotr Jasiukajtis /*
3225c28e83SPiotr Jasiukajtis * remainder(x,p)
3325c28e83SPiotr Jasiukajtis * Code originated from 4.3bsd.
3425c28e83SPiotr Jasiukajtis * Modified by K.C. Ng for SUN 4.0 libm.
3525c28e83SPiotr Jasiukajtis * Return :
3625c28e83SPiotr Jasiukajtis * returns x REM p = x - [x/p]*p as if in infinite precise arithmetic,
3725c28e83SPiotr Jasiukajtis * where [x/p] is the (inifinite bit) integer nearest x/p (in half way
3825c28e83SPiotr Jasiukajtis * case choose the even one).
3925c28e83SPiotr Jasiukajtis * Method :
4025c28e83SPiotr Jasiukajtis * Based on fmod() return x-[x/p]chopped*p exactly.
4125c28e83SPiotr Jasiukajtis */
4225c28e83SPiotr Jasiukajtis
4325c28e83SPiotr Jasiukajtis #include "libm.h"
4425c28e83SPiotr Jasiukajtis
4525c28e83SPiotr Jasiukajtis static const double zero = 0.0, half = 0.5;
4625c28e83SPiotr Jasiukajtis
4725c28e83SPiotr Jasiukajtis double
remainder(double x,double p)4825c28e83SPiotr Jasiukajtis remainder(double x, double p) {
4925c28e83SPiotr Jasiukajtis double halfp;
5025c28e83SPiotr Jasiukajtis int ix, hx, hp;
5125c28e83SPiotr Jasiukajtis
5225c28e83SPiotr Jasiukajtis ix = ((int *)&x)[HIWORD];
5325c28e83SPiotr Jasiukajtis hx = ix & ~0x80000000;
5425c28e83SPiotr Jasiukajtis hp = ((int *)&p)[HIWORD] & ~0x80000000;
5525c28e83SPiotr Jasiukajtis
5625c28e83SPiotr Jasiukajtis if (hp > 0x7ff00000 || (hp == 0x7ff00000 && ((int *)&p)[LOWORD] != 0))
5725c28e83SPiotr Jasiukajtis return (x * p);
5825c28e83SPiotr Jasiukajtis if (hx > 0x7ff00000 || (hx == 0x7ff00000 && ((int *)&x)[LOWORD] != 0))
5925c28e83SPiotr Jasiukajtis return (x * p);
6025c28e83SPiotr Jasiukajtis
6125c28e83SPiotr Jasiukajtis if ((hp | ((int *)&p)[LOWORD]) == 0 || hx == 0x7ff00000)
6225c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, p, 28));
6325c28e83SPiotr Jasiukajtis
6425c28e83SPiotr Jasiukajtis p = fabs(p);
6525c28e83SPiotr Jasiukajtis if (hp < 0x7fe00000)
6625c28e83SPiotr Jasiukajtis x = fmod(x, p + p);
6725c28e83SPiotr Jasiukajtis x = fabs(x);
6825c28e83SPiotr Jasiukajtis if (hp < 0x00200000) {
6925c28e83SPiotr Jasiukajtis if (x + x > p) {
7025c28e83SPiotr Jasiukajtis if (x == p) /* avoid x-x=-0 in RM mode */
7125c28e83SPiotr Jasiukajtis return ((ix < 0)? -zero : zero);
7225c28e83SPiotr Jasiukajtis x -= p;
7325c28e83SPiotr Jasiukajtis if (x + x >= p)
7425c28e83SPiotr Jasiukajtis x -= p;
7525c28e83SPiotr Jasiukajtis }
7625c28e83SPiotr Jasiukajtis } else {
7725c28e83SPiotr Jasiukajtis halfp = half * p;
7825c28e83SPiotr Jasiukajtis if (x > halfp) {
7925c28e83SPiotr Jasiukajtis if (x == p) /* avoid x-x=-0 in RM mode */
8025c28e83SPiotr Jasiukajtis return ((ix < 0)? -zero : zero);
8125c28e83SPiotr Jasiukajtis x -= p;
8225c28e83SPiotr Jasiukajtis if (x >= halfp)
8325c28e83SPiotr Jasiukajtis x -= p;
8425c28e83SPiotr Jasiukajtis }
8525c28e83SPiotr Jasiukajtis }
8625c28e83SPiotr Jasiukajtis return ((ix < 0)? -x : x);
8725c28e83SPiotr Jasiukajtis }
88