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 /*
2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2425c28e83SPiotr Jasiukajtis */
2525c28e83SPiotr Jasiukajtis /*
2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
2725c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2825c28e83SPiotr Jasiukajtis */
2925c28e83SPiotr Jasiukajtis
30*ddc0e0b5SRichard Lowe #pragma weak __remainderl = remainderl
3125c28e83SPiotr Jasiukajtis
3225c28e83SPiotr Jasiukajtis #include "libm.h"
3325c28e83SPiotr Jasiukajtis #include "longdouble.h"
3425c28e83SPiotr Jasiukajtis
3525c28e83SPiotr Jasiukajtis /*
3625c28e83SPiotr Jasiukajtis * remainderl(x,p)
3725c28e83SPiotr Jasiukajtis * returns x REM p = x - [x/p]*p as if in infinite
3825c28e83SPiotr Jasiukajtis * precise arithmetic, where [x/p] is the (inifinite bit)
3925c28e83SPiotr Jasiukajtis * integer nearest x/p (in half way case choose the even one).
4025c28e83SPiotr Jasiukajtis * Method :
4125c28e83SPiotr Jasiukajtis * Based on fmodl() return x-[x/p]chopped*p exactly.
4225c28e83SPiotr Jasiukajtis */
4325c28e83SPiotr Jasiukajtis
4425c28e83SPiotr Jasiukajtis #define HFMAX 5.948657476786158825428796633140035080982e+4931L
4525c28e83SPiotr Jasiukajtis #define DBMIN 6.724206286224187012525355634643505205196e-4932L
4625c28e83SPiotr Jasiukajtis
4725c28e83SPiotr Jasiukajtis static const long double
4825c28e83SPiotr Jasiukajtis zero = 0.0L,
4925c28e83SPiotr Jasiukajtis half = 0.5L,
5025c28e83SPiotr Jasiukajtis hfmax = HFMAX, /* half of the maximum number */
5125c28e83SPiotr Jasiukajtis dbmin = DBMIN; /* double of the minimum (normal) number */
5225c28e83SPiotr Jasiukajtis
5325c28e83SPiotr Jasiukajtis long double
remainderl(long double x,long double p)5425c28e83SPiotr Jasiukajtis remainderl(long double x, long double p) {
5525c28e83SPiotr Jasiukajtis long double hp;
5625c28e83SPiotr Jasiukajtis int sx;
5725c28e83SPiotr Jasiukajtis
5825c28e83SPiotr Jasiukajtis if (isnanl(p))
5925c28e83SPiotr Jasiukajtis return (x + p);
6025c28e83SPiotr Jasiukajtis if (!finitel(x))
6125c28e83SPiotr Jasiukajtis return (x - x);
6225c28e83SPiotr Jasiukajtis p = fabsl(p);
6325c28e83SPiotr Jasiukajtis if (p <= hfmax)
6425c28e83SPiotr Jasiukajtis x = fmodl(x, p + p);
6525c28e83SPiotr Jasiukajtis sx = signbitl(x);
6625c28e83SPiotr Jasiukajtis x = fabsl(x);
6725c28e83SPiotr Jasiukajtis if (p < dbmin) {
6825c28e83SPiotr Jasiukajtis if (x + x > p) {
6925c28e83SPiotr Jasiukajtis if (x == p)
7025c28e83SPiotr Jasiukajtis x = zero;
7125c28e83SPiotr Jasiukajtis else
7225c28e83SPiotr Jasiukajtis x -= p; /* avoid x-x=-0 in RM mode */
7325c28e83SPiotr Jasiukajtis if (x + x >= p)
7425c28e83SPiotr Jasiukajtis x -= p;
7525c28e83SPiotr Jasiukajtis }
7625c28e83SPiotr Jasiukajtis } else {
7725c28e83SPiotr Jasiukajtis hp = half * p;
7825c28e83SPiotr Jasiukajtis if (x > hp) {
7925c28e83SPiotr Jasiukajtis if (x == p)
8025c28e83SPiotr Jasiukajtis x = zero;
8125c28e83SPiotr Jasiukajtis else
8225c28e83SPiotr Jasiukajtis x -= p; /* avoid x-x=-0 in RM mode */
8325c28e83SPiotr Jasiukajtis if (x >= hp)
8425c28e83SPiotr Jasiukajtis x -= p;
8525c28e83SPiotr Jasiukajtis }
8625c28e83SPiotr Jasiukajtis }
8725c28e83SPiotr Jasiukajtis return (sx == 0 ? x : -x);
8825c28e83SPiotr Jasiukajtis }
89