125c28e8Piotr Jasiukajtis/*
225c28e8Piotr Jasiukajtis * CDDL HEADER START
325c28e8Piotr Jasiukajtis *
425c28e8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e8Piotr Jasiukajtis *
825c28e8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e8Piotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e8Piotr Jasiukajtis * and limitations under the License.
1225c28e8Piotr Jasiukajtis *
1325c28e8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e8Piotr Jasiukajtis *
1925c28e8Piotr Jasiukajtis * CDDL HEADER END
2025c28e8Piotr Jasiukajtis */
2125c28e8Piotr Jasiukajtis
2225c28e8Piotr Jasiukajtis/*
2325c28e8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
2425c28e8Piotr Jasiukajtis */
2525c28e8Piotr Jasiukajtis/*
2625c28e8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
2725c28e8Piotr Jasiukajtis * Use is subject to license terms.
2825c28e8Piotr Jasiukajtis */
2925c28e8Piotr Jasiukajtis
30ddc0e0bRichard Lowe#pragma weak __remquo = remquo
3125c28e8Piotr Jasiukajtis
3225c28e8Piotr Jasiukajtis/* INDENT OFF */
3325c28e8Piotr Jasiukajtis/*
3425c28e8Piotr Jasiukajtis * double remquo(double x, double y, int *quo) return remainder(x,y) and an
3525c28e8Piotr Jasiukajtis * integer pointer quo such that *quo = N mod {2**31}, where N is the
3625c28e8Piotr Jasiukajtis * exact integral part of x/y rounded to nearest even.
3725c28e8Piotr Jasiukajtis *
3825c28e8Piotr Jasiukajtis * remquo call internal fmodquo
3925c28e8Piotr Jasiukajtis */
4025c28e8Piotr Jasiukajtis/* INDENT ON */
4125c28e8Piotr Jasiukajtis
4225c28e8Piotr Jasiukajtis#include "libm.h"
4325c28e8Piotr Jasiukajtis#include "libm_protos.h"
4425c28e8Piotr Jasiukajtis#include <math.h>		/* fabs() */
4525c28e8Piotr Jasiukajtis#include <sys/isa_defs.h>
4625c28e8Piotr Jasiukajtis
4725c28e8Piotr Jasiukajtis#if defined(_BIG_ENDIAN)
4825c28e8Piotr Jasiukajtis#define	HIWORD	0
4925c28e8Piotr Jasiukajtis#define	LOWORD	1
5025c28e8Piotr Jasiukajtis#else
5125c28e8Piotr Jasiukajtis#define	HIWORD	1
5225c28e8Piotr Jasiukajtis#define	LOWORD	0
5325c28e8Piotr Jasiukajtis#endif
5425c28e8Piotr Jasiukajtis#define	__HI(x)	((int *) &x)[HIWORD]
5525c28e8Piotr Jasiukajtis#define	__LO(x)	((int *) &x)[LOWORD]
5625c28e8Piotr Jasiukajtis
5725c28e8Piotr Jasiukajtisstatic const double one = 1.0, Zero[] = {0.0, -0.0};
5825c28e8Piotr Jasiukajtis
5925c28e8Piotr Jasiukajtisstatic double
6025c28e8Piotr Jasiukajtisfmodquo(double x, double y, int *quo) {
6125c28e8Piotr Jasiukajtis	int n, hx, hy, hz, ix, iy, sx, sq, i, m;
6225c28e8Piotr Jasiukajtis	unsigned lx, ly, lz;
6325c28e8Piotr Jasiukajtis
6425c28e8Piotr Jasiukajtis	hx = __HI(x);		/* high word of x */
6525c28e8Piotr Jasiukajtis	lx = __LO(x);		/* low  word of x */
6625c28e8Piotr Jasiukajtis	hy = __HI(y);		/* high word of y */
6725c28e8Piotr Jasiukajtis	ly = __LO(y);		/* low  word of y */
6825c28e8Piotr Jasiukajtis	sx = hx & 0x80000000;	/* sign of x */
6925c28e8Piotr Jasiukajtis	sq = (hx ^ hy) & 0x80000000;	/* sign of x/y */
7025c28e8Piotr Jasiukajtis	hx ^= sx;		/* |x| */
7125c28e8Piotr Jasiukajtis	hy &= 0x7fffffff;	/* |y| */
7225c28e8Piotr Jasiukajtis
7325c28e8Piotr Jasiukajtis	/* purge off exception values */
7425c28e8Piotr Jasiukajtis	*quo = 0;
7525c28e8Piotr Jasiukajtis	if ((hy | ly) == 0 || hx >= 0x7ff00000 ||	/* y=0, or x !finite */
7625c28e8Piotr Jasiukajtis	    (hy | ((ly | -ly) >> 31)) > 0x7ff00000)	/* or y is NaN */
7725c28e8Piotr Jasiukajtis		return ((x * y) / (x * y));
7825c28e8Piotr Jasiukajtis	if (hx <= hy) {
7925c28e8Piotr Jasiukajtis		if (hx < hy || lx < ly)
8025c28e8Piotr Jasiukajtis			return (x);	/* |x|<|y| return x */
8125c28e8Piotr Jasiukajtis		if (lx == ly) {
8225c28e8Piotr Jasiukajtis			*quo = 1 + (sq >> 30);
8325c28e8Piotr Jasiukajtis			/* |x|=|y| return x*0 */
8425c28e8Piotr Jasiukajtis			return (Zero[(unsigned) sx >> 31]);
8525c28e8Piotr Jasiukajtis		}
8625c28e8Piotr Jasiukajtis	}
8725c28e8Piotr Jasiukajtis
8825c28e8Piotr Jasiukajtis	/* determine ix = ilogb(x) */
8925c28e8Piotr Jasiukajtis	if (hx < 0x00100000) {	/* subnormal x */
9025c28e8Piotr Jasiukajtis		if (hx == 0) {
9125c28e8Piotr Jasiukajtis			for (ix = -1043, i = lx; i > 0; i <<= 1)
9225c28e8Piotr Jasiukajtis				ix -= 1;
9325c28e8Piotr Jasiukajtis		} else {
9425c28e8Piotr Jasiukajtis			for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
9525c28e8Piotr Jasiukajtis				ix -= 1;
9625c28e8Piotr Jasiukajtis		}
9725c28e8Piotr Jasiukajtis	} else
9825c28e8Piotr Jasiukajtis		ix = (hx >> 20) - 1023;
9925c28e8Piotr Jasiukajtis
10025c28e8Piotr Jasiukajtis	/* determine iy = ilogb(y) */
10125c28e8Piotr Jasiukajtis	if (hy < 0x00100000) {	/* subnormal y */
10225c28e8Piotr Jasiukajtis		if (hy == 0) {
10325c28e8Piotr Jasiukajtis			for (iy = -1043, i = ly; i > 0; i <<= 1)
10425c28e8Piotr Jasiukajtis				iy -= 1;
10525c28e8Piotr Jasiukajtis		} else {
10625c28e8Piotr Jasiukajtis			for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
10725c28e8Piotr Jasiukajtis				iy -= 1;
10825c28e8Piotr Jasiukajtis		}
10925c28e8Piotr Jasiukajtis	} else
11025c28e8Piotr Jasiukajtis		iy = (hy >> 20) - 1023;
11125c28e8Piotr Jasiukajtis
11225c28e8Piotr Jasiukajtis	/* set up {hx,lx}, {hy,ly} and align y to x */
11325c28e8Piotr Jasiukajtis	if (ix >= -1022)
11425c28e8Piotr Jasiukajtis		hx = 0x00100000 | (0x000fffff & hx);
11525c28e8Piotr Jasiukajtis	else {			/* subnormal x, shift x to normal */
11625c28e8Piotr Jasiukajtis		n = -1022 - ix;
11725c28e8Piotr Jasiukajtis		if (n <= 31) {
11825c28e8Piotr Jasiukajtis			hx = (hx << n) | (lx >> (32 - n));
11925c28e8Piotr Jasiukajtis			lx <<= n;
12025c28e8Piotr Jasiukajtis		} else {
12125c28e8Piotr Jasiukajtis			hx = lx << (n - 32);
12225c28e8Piotr Jasiukajtis			lx = 0;
12325c28e8Piotr Jasiukajtis		}
12425c28e8Piotr Jasiukajtis	}
12525c28e8Piotr Jasiukajtis	if (iy >= -1022)
12625c28e8Piotr Jasiukajtis		hy = 0x00100000 | (0x000fffff & hy);
12725c28e8Piotr Jasiukajtis	else {			/* subnormal y, shift y to normal */
12825c28e8Piotr Jasiukajtis		n = -1022 - iy;
12925c28e8Piotr Jasiukajtis		if (n <= 31) {
13025c28e8Piotr Jasiukajtis			hy = (hy << n) | (ly >> (32 - n));
13125c28e8Piotr Jasiukajtis			ly <<= n;
13225c28e8Piotr Jasiukajtis		} else {
13325c28e8Piotr Jasiukajtis			hy = ly << (n - 32);
13425c28e8Piotr Jasiukajtis			ly = 0;
13525c28e8Piotr Jasiukajtis		}
13625c28e8Piotr Jasiukajtis	}
13725c28e8Piotr Jasiukajtis
13825c28e8Piotr Jasiukajtis	/* fix point fmod */
13925c28e8Piotr Jasiukajtis	n = ix - iy;
14025c28e8Piotr Jasiukajtis	m = 0;
14125c28e8Piotr Jasiukajtis	while (n--) {
14225c28e8Piotr Jasiukajtis		hz = hx - hy;
14325c28e8Piotr Jasiukajtis		lz = lx - ly;
14425c28e8Piotr Jasiukajtis		if (lx < ly)
14525c28e8Piotr Jasiukajtis			hz -= 1;
14625c28e8Piotr Jasiukajtis		if (hz < 0) {
14725c28e8Piotr Jasiukajtis			hx = hx + hx + (lx >> 31);
14825c28e8Piotr Jasiukajtis			lx = lx + lx;
14925c28e8Piotr Jasiukajtis		} else {
15025c28e8Piotr Jasiukajtis			m += 1;
15125c28e8Piotr Jasiukajtis			if ((hz | lz) == 0) {	/* return sign(x)*0 */
15225c28e8Piotr Jasiukajtis				if (n < 31)
15325c28e8Piotr Jasiukajtis					m <<= 1 + n;
15425c28e8Piotr Jasiukajtis				else
15525c28e8Piotr Jasiukajtis					m = 0;
15625c28e8Piotr Jasiukajtis				m &= 0x7fffffff;
15725c28e8Piotr Jasiukajtis				*quo = sq >= 0 ? m : -m;
15825c28e8Piotr Jasiukajtis				return (Zero[(unsigned) sx >> 31]);
15925c28e8Piotr Jasiukajtis			}
16025c28e8Piotr Jasiukajtis			hx = hz + hz + (lz >> 31);
16125c28e8Piotr Jasiukajtis			lx = lz + lz;
16225c28e8Piotr Jasiukajtis		}
16325c28e8Piotr Jasiukajtis		m += m;
16425c28e8Piotr Jasiukajtis	}
16525c28e8Piotr Jasiukajtis	hz = hx - hy;
16625c28e8Piotr Jasiukajtis	lz = lx - ly;
16725c28e8Piotr Jasiukajtis	if (lx < ly)
16825c28e8Piotr Jasiukajtis		hz -= 1;
16925c28e8Piotr Jasiukajtis	if (hz >= 0) {
17025c28e8Piotr Jasiukajtis		hx = hz;
17125c28e8Piotr Jasiukajtis		lx = lz;
17225c28e8Piotr Jasiukajtis		m += 1;
17325c28e8Piotr Jasiukajtis	}
17425c28e8Piotr Jasiukajtis	m &= 0x7fffffff;
17525c28e8Piotr Jasiukajtis	*quo = sq >= 0 ? m : -m;
17625c28e8Piotr Jasiukajtis
17725c28e8Piotr Jasiukajtis	/* convert back to floating value and restore the sign */
17825c28e8Piotr Jasiukajtis	if ((hx | lx) == 0) {	/* return sign(x)*0 */
17925c28e8Piotr Jasiukajtis		return (Zero[(unsigned) sx >> 31]);
18025c28e8Piotr Jasiukajtis	}
18125c28e8Piotr Jasiukajtis	while (hx < 0x00100000) {	/* normalize x */
18225c28e8Piotr Jasiukajtis		hx = hx + hx + (lx >> 31);
18325c28e8Piotr Jasiukajtis		lx = lx + lx;
18425c28e8Piotr Jasiukajtis		iy -= 1;
18525c28e8Piotr Jasiukajtis	}
18625c28e8Piotr Jasiukajtis	if (iy >= -1022) {	/* normalize output */
18725c28e8Piotr Jasiukajtis		hx = (hx - 0x00100000) | ((iy + 1023) << 20);
18825c28e8Piotr Jasiukajtis		__HI(x) = hx | sx;
18925c28e8Piotr Jasiukajtis		__LO(x) = lx;
19025c28e8Piotr Jasiukajtis	} else {			/* subnormal output */
19125c28e8Piotr Jasiukajtis		n = -1022 - iy;
19225c28e8Piotr Jasiukajtis		if (n <= 20) {
19325c28e8Piotr Jasiukajtis			lx = (lx >> n) | ((unsigned) hx << (32 - n));
19425c28e8Piotr Jasiukajtis			hx >>= n;
19525c28e8Piotr Jasiukajtis		} else if (n <= 31) {
19625c28e8Piotr Jasiukajtis			lx = (hx << (32 - n)) | (lx >> n);
19725c28e8Piotr Jasiukajtis			hx = sx;
19825c28e8Piotr Jasiukajtis		} else {
19925c28e8Piotr Jasiukajtis			lx = hx >> (n - 32);
20025c28e8Piotr Jasiukajtis			hx = sx;
20125c28e8Piotr Jasiukajtis		}
20225c28e8Piotr Jasiukajtis		__HI(x) = hx | sx;
20325c28e8Piotr Jasiukajtis		__LO(x) = lx;
20425c28e8Piotr Jasiukajtis		x *= one;	/* create necessary signal */
20525c28e8Piotr Jasiukajtis	}
20625c28e8Piotr Jasiukajtis	return (x);		/* exact output */
20725c28e8Piotr Jasiukajtis}
20825c28e8Piotr Jasiukajtis
20925c28e8Piotr Jasiukajtis#define	zero	Zero[0]
21025c28e8Piotr Jasiukajtis
21125c28e8Piotr Jasiukajtisdouble
21225c28e8Piotr Jasiukajtisremquo(double x, double y, int *quo) {
21325c28e8Piotr Jasiukajtis	int hx, hy, sx, sq;
21425c28e8Piotr Jasiukajtis	double v;
21525c28e8Piotr Jasiukajtis	unsigned ly;
21625c28e8Piotr Jasiukajtis
21725c28e8Piotr Jasiukajtis	hx = __HI(x);		/* high word of x */
21825c28e8Piotr Jasiukajtis	hy = __HI(y);		/* high word of y */
21925c28e8Piotr Jasiukajtis	ly = __LO(y);		/* low  word of y */
22025c28e8Piotr Jasiukajtis	sx = hx & 0x80000000;	/* sign of x */
22125c28e8Piotr Jasiukajtis	sq = (hx ^ hy) & 0x80000000;	/* sign of x/y */
22225c28e8Piotr Jasiukajtis	hx ^= sx;		/* |x| */
22325c28e8Piotr Jasiukajtis	hy &= 0x7fffffff;	/* |y| */
22425c28e8Piotr Jasiukajtis
22525c28e8Piotr Jasiukajtis	/* purge off exception values */
22625c28e8Piotr Jasiukajtis	*quo = 0;
22725c28e8Piotr Jasiukajtis	if ((hy | ly) == 0 || hx >= 0x7ff00000 ||	/* y=0, or x !finite */
22825c28e8Piotr Jasiukajtis	    (hy | ((ly | -ly) >> 31)) > 0x7ff00000)	/* or y is NaN */
22925c28e8Piotr Jasiukajtis		return ((x * y) / (x * y));
23025c28e8Piotr Jasiukajtis
23125c28e8Piotr Jasiukajtis	y = fabs(y);
23225c28e8Piotr Jasiukajtis	x = fabs(x);
23325c28e8Piotr Jasiukajtis	if (hy <= 0x7fdfffff) {
23425c28e8Piotr Jasiukajtis		x = fmodquo(x, y + y, quo);
23525c28e8Piotr Jasiukajtis		*quo = ((*quo) & 0x3fffffff) << 1;
23625c28e8Piotr Jasiukajtis	}
23725c28e8Piotr Jasiukajtis	if (hy < 0x00200000) {
23825c28e8Piotr Jasiukajtis		if (x + x > y) {
23925c28e8Piotr Jasiukajtis			*quo += 1;
24025c28e8Piotr Jasiukajtis			if (x == y)
24125c28e8Piotr Jasiukajtis				x = zero;
24225c28e8Piotr Jasiukajtis			else
24325c28e8Piotr Jasiukajtis				x -= y;
24425c28e8Piotr Jasiukajtis			if (x + x >= y) {
24525c28e8Piotr Jasiukajtis				x -= y;
24625c28e8Piotr Jasiukajtis				*quo += 1;
24725c28e8Piotr Jasiukajtis			}
24825c28e8Piotr Jasiukajtis		}
24925c28e8Piotr Jasiukajtis	} else {
25025c28e8Piotr Jasiukajtis		v = 0.5 * y;
25125c28e8Piotr Jasiukajtis		if (x > v) {
25225c28e8Piotr Jasiukajtis			*quo += 1;
25325c28e8Piotr Jasiukajtis			if (x == y)
25425c28e8Piotr Jasiukajtis				x = zero;
25525c28e8Piotr Jasiukajtis			else
25625c28e8Piotr Jasiukajtis				x -= y;
25725c28e8Piotr Jasiukajtis			if (x >= v) {
25825c28e8Piotr Jasiukajtis				x -= y;
25925c28e8Piotr Jasiukajtis				*quo += 1;
26025c28e8Piotr Jasiukajtis			}
26125c28e8Piotr Jasiukajtis		}
26225c28e8Piotr Jasiukajtis	}
26325c28e8Piotr Jasiukajtis	if (sq != 0)
26425c28e8Piotr Jasiukajtis		*quo = -(*quo);
26525c28e8Piotr Jasiukajtis	return (sx == 0 ? x : -x);
26625c28e8Piotr Jasiukajtis}
267