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 __remquof = remquof
3125c28e8Piotr Jasiukajtis
3225c28e8Piotr Jasiukajtis/* INDENT OFF */
3325c28e8Piotr Jasiukajtis/*
3425c28e8Piotr Jasiukajtis * float remquof(float x, float y, int *quo) return remainderf(x,y) and an
3525c28e8Piotr Jasiukajtis * integer pointer quo such that *quo = N mod (2**31),  where N is the
3625c28e8Piotr Jasiukajtis * exact integeral part of x/y rounded to nearest even.
3725c28e8Piotr Jasiukajtis *
3825c28e8Piotr Jasiukajtis * remquof call internal fmodquof
3925c28e8Piotr Jasiukajtis */
4025c28e8Piotr Jasiukajtis
4125c28e8Piotr Jasiukajtis#include "libm.h"
4225c28e8Piotr Jasiukajtis#include "libm_protos.h"
4325c28e8Piotr Jasiukajtis#include <math.h>
4425c28e8Piotr Jasiukajtisextern float fabsf(float);
4525c28e8Piotr Jasiukajtis
4625c28e8Piotr Jasiukajtisstatic const int
4725c28e8Piotr Jasiukajtis	is = (int) 0x80000000,
4825c28e8Piotr Jasiukajtis	im = 0x007fffff,
4925c28e8Piotr Jasiukajtis	ii = 0x7f800000,
5025c28e8Piotr Jasiukajtis	iu = 0x00800000;
5125c28e8Piotr Jasiukajtis
5225c28e8Piotr Jasiukajtisstatic const float zero = 0.0F, half = 0.5F;
5325c28e8Piotr Jasiukajtis/* INDENT ON */
5425c28e8Piotr Jasiukajtis
5525c28e8Piotr Jasiukajtisstatic float
5625c28e8Piotr Jasiukajtisfmodquof(float x, float y, int *quo) {
5725c28e8Piotr Jasiukajtis	float w;
5825c28e8Piotr Jasiukajtis	int hx, ix, iy, iz, k, ny, nd, m, sq;
5925c28e8Piotr Jasiukajtis
6025c28e8Piotr Jasiukajtis	hx = *(int *) &x;
6125c28e8Piotr Jasiukajtis	ix = hx & 0x7fffffff;
6225c28e8Piotr Jasiukajtis	iy = *(int *) &y;
6325c28e8Piotr Jasiukajtis	sq = (iy ^ hx) & is;	/* sign of x/y */
6425c28e8Piotr Jasiukajtis	iy &= 0x7fffffff;
6525c28e8Piotr Jasiukajtis
6625c28e8Piotr Jasiukajtis	/* purge off exception values */
6725c28e8Piotr Jasiukajtis	*quo = 0;
6825c28e8Piotr Jasiukajtis	if (ix >= ii || iy > ii || iy == 0) {
6925c28e8Piotr Jasiukajtis		w = x * y;
7025c28e8Piotr Jasiukajtis		w = w / w;
7125c28e8Piotr Jasiukajtis	} else if (ix <= iy) {
7225c28e8Piotr Jasiukajtis		if (ix < iy)
7325c28e8Piotr Jasiukajtis			w = x;	/* return x if |x|<|y| */
7425c28e8Piotr Jasiukajtis		else {
7525c28e8Piotr Jasiukajtis			*quo = 1 + (sq >> 30);
7625c28e8Piotr Jasiukajtis			w = zero * x;	/* return sign(x)*0.0  */
7725c28e8Piotr Jasiukajtis		}
7825c28e8Piotr Jasiukajtis	} else {
7925c28e8Piotr Jasiukajtis		/* INDENT OFF */
8025c28e8Piotr Jasiukajtis		/*
8125c28e8Piotr Jasiukajtis		 * scale x,y to "normal" with
8225c28e8Piotr Jasiukajtis		 *	ny = exponent of y
8325c28e8Piotr Jasiukajtis		 *	nd = exponent of x minus exponent of y
8425c28e8Piotr Jasiukajtis		 */
8525c28e8Piotr Jasiukajtis		/* INDENT ON */
8625c28e8Piotr Jasiukajtis		ny = iy >> 23;
8725c28e8Piotr Jasiukajtis		k = ix >> 23;
8825c28e8Piotr Jasiukajtis
8925c28e8Piotr Jasiukajtis		/* special case for subnormal y or x */
9025c28e8Piotr Jasiukajtis		if (ny == 0) {
9125c28e8Piotr Jasiukajtis			ny = 1;
9225c28e8Piotr Jasiukajtis			while (iy < iu) {
9325c28e8Piotr Jasiukajtis				ny -= 1;
9425c28e8Piotr Jasiukajtis				iy += iy;
9525c28e8Piotr Jasiukajtis			}
9625c28e8Piotr Jasiukajtis			nd = k - ny;
9725c28e8Piotr Jasiukajtis			if (k == 0) {
9825c28e8Piotr Jasiukajtis				nd += 1;
9925c28e8Piotr Jasiukajtis				while (ix < iu) {
10025c28e8Piotr Jasiukajtis					nd -= 1;
10125c28e8Piotr Jasiukajtis					ix += ix;
10225c28e8Piotr Jasiukajtis				}
10325c28e8Piotr Jasiukajtis			} else
10425c28e8Piotr Jasiukajtis				ix = iu | (ix & im);
10525c28e8Piotr Jasiukajtis		} else {
10625c28e8Piotr Jasiukajtis			nd = k - ny;
10725c28e8Piotr Jasiukajtis			ix = iu | (ix & im);
10825c28e8Piotr Jasiukajtis			iy = iu | (iy & im);
10925c28e8Piotr Jasiukajtis		}
11025c28e8Piotr Jasiukajtis		/* INDENT OFF */
11125c28e8Piotr Jasiukajtis		/* fix point fmod for normalized ix and iy */
11225c28e8Piotr Jasiukajtis		/*
11325c28e8Piotr Jasiukajtis		 * while (nd--) {
11425c28e8Piotr Jasiukajtis		 *	iz = ix - iy;
11525c28e8Piotr Jasiukajtis		 *	if (iz < 0)
11625c28e8Piotr Jasiukajtis		 *		ix = ix + ix;
11725c28e8Piotr Jasiukajtis		 *	else if (iz == 0) {
11825c28e8Piotr Jasiukajtis		 *		*(int *) &w = is & hx;
11925c28e8Piotr Jasiukajtis		 *		return w;
12025c28e8Piotr Jasiukajtis		 *	} else
12125c28e8Piotr Jasiukajtis		 *		ix = iz + iz;
12225c28e8Piotr Jasiukajtis		 * }
12325c28e8Piotr Jasiukajtis		 */
12425c28e8Piotr Jasiukajtis		/* INDENT ON */
12525c28e8Piotr Jasiukajtis		/* unroll the above loop 4 times to gain performance */
12625c28e8Piotr Jasiukajtis		m = 0;
12725c28e8Piotr Jasiukajtis		k = nd >> 2;
12825c28e8Piotr Jasiukajtis		nd -= (k << 2);
12925c28e8Piotr Jasiukajtis		while (k--) {
13025c28e8Piotr Jasiukajtis			iz = ix - iy;
13125c28e8Piotr Jasiukajtis			if (iz >= 0) {
13225c28e8Piotr Jasiukajtis				m += 1;
13325c28e8Piotr Jasiukajtis				ix = iz + iz;
13425c28e8Piotr Jasiukajtis			} else
13525c28e8Piotr Jasiukajtis				ix += ix;
13625c28e8Piotr Jasiukajtis			m += m;
13725c28e8Piotr Jasiukajtis			iz = ix - iy;
13825c28e8Piotr Jasiukajtis			if (iz >= 0) {
13925c28e8Piotr Jasiukajtis				m += 1;
14025c28e8Piotr Jasiukajtis				ix = iz + iz;
14125c28e8Piotr Jasiukajtis			} else
14225c28e8Piotr Jasiukajtis				ix += ix;
14325c28e8Piotr Jasiukajtis			m += m;
14425c28e8Piotr Jasiukajtis			iz = ix - iy;
14525c28e8Piotr Jasiukajtis			if (iz >= 0) {
14625c28e8Piotr Jasiukajtis				m += 1;
14725c28e8Piotr Jasiukajtis				ix = iz + iz;
14825c28e8Piotr Jasiukajtis			} else
14925c28e8Piotr Jasiukajtis				ix += ix;
15025c28e8Piotr Jasiukajtis			m += m;
15125c28e8Piotr Jasiukajtis			iz = ix - iy;
15225c28e8Piotr Jasiukajtis			if (iz >= 0) {
15325c28e8Piotr Jasiukajtis				m += 1;
15425c28e8Piotr Jasiukajtis				ix = iz + iz;
15525c28e8Piotr Jasiukajtis			} else
15625c28e8Piotr Jasiukajtis				ix += ix;
15725c28e8Piotr Jasiukajtis			m += m;
15825c28e8Piotr Jasiukajtis			if (iz == 0) {
15925c28e8Piotr Jasiukajtis				iz = (k << 2) + nd;
16025c28e8Piotr Jasiukajtis				if (iz < 32)
16125c28e8Piotr Jasiukajtis					m <<= iz;
16225c28e8Piotr Jasiukajtis				else
16325c28e8Piotr Jasiukajtis					m = 0;
16425c28e8Piotr Jasiukajtis				m &= 0x7fffffff;
16525c28e8Piotr Jasiukajtis				*quo = sq >= 0 ? m : -m;
16625c28e8Piotr Jasiukajtis				*(int *) &w = is & hx;
16725c28e8Piotr Jasiukajtis				return (w);
16825c28e8Piotr Jasiukajtis			}
16925c28e8Piotr Jasiukajtis		}
17025c28e8Piotr Jasiukajtis		while (nd--) {
17125c28e8Piotr Jasiukajtis			iz = ix - iy;
17225c28e8Piotr Jasiukajtis			if (iz >= 0) {
17325c28e8Piotr Jasiukajtis				m += 1;
17425c28e8Piotr Jasiukajtis				ix = iz + iz;
17525c28e8Piotr Jasiukajtis			} else
17625c28e8Piotr Jasiukajtis				ix += ix;
17725c28e8Piotr Jasiukajtis			m += m;
17825c28e8Piotr Jasiukajtis		}
17925c28e8Piotr Jasiukajtis		/* end of unrolling */
18025c28e8Piotr Jasiukajtis
18125c28e8Piotr Jasiukajtis		iz = ix - iy;
18225c28e8Piotr Jasiukajtis		if (iz >= 0) {
18325c28e8Piotr Jasiukajtis			m += 1;
18425c28e8Piotr Jasiukajtis			ix = iz;
18525c28e8Piotr Jasiukajtis		}
18625c28e8Piotr Jasiukajtis		m &= 0x7fffffff;
18725c28e8Piotr Jasiukajtis		*quo = sq >= 0 ? m : -m;
18825c28e8Piotr Jasiukajtis
18925c28e8Piotr Jasiukajtis		/* convert back to floating value and restore the sign */
19025c28e8Piotr Jasiukajtis		if (ix == 0) {
19125c28e8Piotr Jasiukajtis			*(int *) &w = is & hx;
19225c28e8Piotr Jasiukajtis			return (w);
19325c28e8Piotr Jasiukajtis		}
19425c28e8Piotr Jasiukajtis		while (ix < iu) {
19525c28e8Piotr Jasiukajtis			ix += ix;
19625c28e8Piotr Jasiukajtis			ny -= 1;
19725c28e8Piotr Jasiukajtis		}
19825c28e8Piotr Jasiukajtis		while (ix > (iu + iu)) {
19925c28e8Piotr Jasiukajtis			ny += 1;
20025c28e8Piotr Jasiukajtis			ix >>= 1;
20125c28e8Piotr Jasiukajtis		}
20225c28e8Piotr Jasiukajtis		if (ny > 0)
20325c28e8Piotr Jasiukajtis			*(int *) &w = (is & hx) | (ix & im) | (ny << 23);
20425c28e8Piotr Jasiukajtis		else {		/* subnormal output */
20525c28e8Piotr Jasiukajtis			k = -ny + 1;
20625c28e8Piotr Jasiukajtis			ix >>= k;
20725c28e8Piotr Jasiukajtis			*(int *) &w = (is & hx) | ix;
20825c28e8Piotr Jasiukajtis		}
20925c28e8Piotr Jasiukajtis	}
21025c28e8Piotr Jasiukajtis	return (w);
21125c28e8Piotr Jasiukajtis}
21225c28e8Piotr Jasiukajtis
21325c28e8Piotr Jasiukajtisfloat
21425c28e8Piotr Jasiukajtisremquof(float x, float y, int *quo) {
21525c28e8Piotr Jasiukajtis	int hx, hy, sx, sq;
21625c28e8Piotr Jasiukajtis	float v;
21725c28e8Piotr Jasiukajtis
21825c28e8Piotr Jasiukajtis	hx = *(int *) &x;	/* high word of x */
21925c28e8Piotr Jasiukajtis	hy = *(int *) &y;	/* high word of y */
22025c28e8Piotr Jasiukajtis	sx = hx & is;		/* sign of x */
22125c28e8Piotr Jasiukajtis	sq = (hx ^ hy) & is;	/* sign of x/y */
22225c28e8Piotr Jasiukajtis	hx ^= sx;		/* |x| */
22325c28e8Piotr Jasiukajtis	hy &= 0x7fffffff;	/* |y| */
22425c28e8Piotr Jasiukajtis
22525c28e8Piotr Jasiukajtis	/* purge off exception values: y is 0 or NaN, x is Inf or NaN */
22625c28e8Piotr Jasiukajtis	*quo = 0;
22725c28e8Piotr Jasiukajtis	if (hx >= ii || hy > ii || hy == 0) {
22825c28e8Piotr Jasiukajtis		v = x * y;
22925c28e8Piotr Jasiukajtis		return (v / v);
23025c28e8Piotr Jasiukajtis	}
23125c28e8Piotr Jasiukajtis
23225c28e8Piotr Jasiukajtis	y = fabsf(y);
23325c28e8Piotr Jasiukajtis	x = fabsf(x);
23425c28e8Piotr Jasiukajtis	if (hy <= 0x7f7fffff) {
23525c28e8Piotr Jasiukajtis		x = fmodquof(x, y + y, quo);
23625c28e8Piotr Jasiukajtis		*quo = ((*quo) & 0x3fffffff) << 1;
23725c28e8Piotr Jasiukajtis	}
23825c28e8Piotr Jasiukajtis	if (hy < 0x01000000) {
23925c28e8Piotr Jasiukajtis		if (x + x > y) {
24025c28e8Piotr Jasiukajtis			*quo += 1;
24125c28e8Piotr Jasiukajtis			if (x == y)
24225c28e8Piotr Jasiukajtis				x = zero;
24325c28e8Piotr Jasiukajtis			else
24425c28e8Piotr Jasiukajtis				x -= y;
24525c28e8Piotr Jasiukajtis			if (x + x >= y) {
24625c28e8Piotr Jasiukajtis				x -= y;
24725c28e8Piotr Jasiukajtis				*quo += 1;
24825c28e8Piotr Jasiukajtis			}
24925c28e8Piotr Jasiukajtis		}
25025c28e8Piotr Jasiukajtis	} else {
25125c28e8Piotr Jasiukajtis		v = half * y;
25225c28e8Piotr Jasiukajtis		if (x > v) {
25325c28e8Piotr Jasiukajtis			*quo += 1;
25425c28e8Piotr Jasiukajtis			if (x == y)
25525c28e8Piotr Jasiukajtis				x = zero;
25625c28e8Piotr Jasiukajtis			else
25725c28e8Piotr Jasiukajtis				x -= y;
25825c28e8Piotr Jasiukajtis			if (x >= v) {
25925c28e8Piotr Jasiukajtis				x -= y;
26025c28e8Piotr Jasiukajtis				*quo += 1;
26125c28e8Piotr Jasiukajtis			}
26225c28e8Piotr Jasiukajtis		}
26325c28e8Piotr Jasiukajtis	}
26425c28e8Piotr Jasiukajtis	if (sq != 0)
26525c28e8Piotr Jasiukajtis		*quo = -(*quo);
26625c28e8Piotr Jasiukajtis	return (sx == 0 ? x : -x);
26725c28e8Piotr Jasiukajtis}
268