xref: /illumos-gate/usr/src/lib/libm/common/C/jn.c (revision 685c1a21)
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 
30ddc0e0b5SRichard Lowe #pragma weak __jn = jn
31ddc0e0b5SRichard Lowe #pragma weak __yn = yn
3225c28e83SPiotr Jasiukajtis 
3325c28e83SPiotr Jasiukajtis /*
3425c28e83SPiotr Jasiukajtis  * floating point Bessel's function of the 1st and 2nd kind
3525c28e83SPiotr Jasiukajtis  * of order n: jn(n,x),yn(n,x);
3625c28e83SPiotr Jasiukajtis  *
3725c28e83SPiotr Jasiukajtis  * Special cases:
3825c28e83SPiotr Jasiukajtis  *	y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
3925c28e83SPiotr Jasiukajtis  *	y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
4025c28e83SPiotr Jasiukajtis  * Note 2. About jn(n,x), yn(n,x)
4125c28e83SPiotr Jasiukajtis  *	For n=0, j0(x) is called,
4225c28e83SPiotr Jasiukajtis  *	for n=1, j1(x) is called,
4325c28e83SPiotr Jasiukajtis  *	for n<x, forward recursion us used starting
4425c28e83SPiotr Jasiukajtis  *	from values of j0(x) and j1(x).
4525c28e83SPiotr Jasiukajtis  *	for n>x, a continued fraction approximation to
4625c28e83SPiotr Jasiukajtis  *	j(n,x)/j(n-1,x) is evaluated and then backward
4725c28e83SPiotr Jasiukajtis  *	recursion is used starting from a supposed value
4825c28e83SPiotr Jasiukajtis  *	for j(n,x). The resulting value of j(0,x) is
4925c28e83SPiotr Jasiukajtis  *	compared with the actual value to correct the
5025c28e83SPiotr Jasiukajtis  *	supposed value of j(n,x).
5125c28e83SPiotr Jasiukajtis  *
5225c28e83SPiotr Jasiukajtis  *	yn(n,x) is similar in all respects, except
5325c28e83SPiotr Jasiukajtis  *	that forward recursion is used for all
5425c28e83SPiotr Jasiukajtis  *	values of n>1.
5525c28e83SPiotr Jasiukajtis  *
5625c28e83SPiotr Jasiukajtis  */
5725c28e83SPiotr Jasiukajtis 
5825c28e83SPiotr Jasiukajtis #include "libm.h"
5925c28e83SPiotr Jasiukajtis #include <float.h>	/* DBL_MIN */
6025c28e83SPiotr Jasiukajtis #include <values.h>	/* X_TLOSS */
6125c28e83SPiotr Jasiukajtis #include "xpg6.h"	/* __xpg6 */
6225c28e83SPiotr Jasiukajtis 
6325c28e83SPiotr Jasiukajtis #define	GENERIC double
6425c28e83SPiotr Jasiukajtis 
6525c28e83SPiotr Jasiukajtis static const GENERIC
6625c28e83SPiotr Jasiukajtis 	invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
6725c28e83SPiotr Jasiukajtis 	two	= 2.0,
6825c28e83SPiotr Jasiukajtis 	zero	= 0.0,
6925c28e83SPiotr Jasiukajtis 	one	= 1.0;
7025c28e83SPiotr Jasiukajtis 
7125c28e83SPiotr Jasiukajtis GENERIC
jn(int n,GENERIC x)72*685c1a21SRichard Lowe jn(int n, GENERIC x)
73*685c1a21SRichard Lowe {
7425c28e83SPiotr Jasiukajtis 	int i, sgn;
7525c28e83SPiotr Jasiukajtis 	GENERIC a, b, temp = 0;
7625c28e83SPiotr Jasiukajtis 	GENERIC z, w, ox, on;
7725c28e83SPiotr Jasiukajtis 
7825c28e83SPiotr Jasiukajtis 	/*
7925c28e83SPiotr Jasiukajtis 	 * J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
8025c28e83SPiotr Jasiukajtis 	 * Thus, J(-n,x) = J(n,-x)
8125c28e83SPiotr Jasiukajtis 	 */
82*685c1a21SRichard Lowe 	ox = x;
83*685c1a21SRichard Lowe 	on = (GENERIC)n;
84*685c1a21SRichard Lowe 
8525c28e83SPiotr Jasiukajtis 	if (n < 0) {
8625c28e83SPiotr Jasiukajtis 		n = -n;
8725c28e83SPiotr Jasiukajtis 		x = -x;
8825c28e83SPiotr Jasiukajtis 	}
8925c28e83SPiotr Jasiukajtis 	if (isnan(x))
9025c28e83SPiotr Jasiukajtis 		return (x*x);	/* + -> * for Cheetah */
91*685c1a21SRichard Lowe 	if (!((int)_lib_version == libm_ieee ||
92*685c1a21SRichard Lowe 	    (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
93*685c1a21SRichard Lowe 		if (fabs(x) > X_TLOSS)
9425c28e83SPiotr Jasiukajtis 			return (_SVID_libm_err(on, ox, 38));
9525c28e83SPiotr Jasiukajtis 	}
9625c28e83SPiotr Jasiukajtis 	if (n == 0)
9725c28e83SPiotr Jasiukajtis 		return (j0(x));
9825c28e83SPiotr Jasiukajtis 	if (n == 1)
9925c28e83SPiotr Jasiukajtis 		return (j1(x));
10025c28e83SPiotr Jasiukajtis 	if ((n&1) == 0)
101*685c1a21SRichard Lowe 		sgn = 0;			/* even n */
10225c28e83SPiotr Jasiukajtis 	else
10325c28e83SPiotr Jasiukajtis 		sgn = signbit(x);	/* old n  */
10425c28e83SPiotr Jasiukajtis 	x = fabs(x);
10525c28e83SPiotr Jasiukajtis 	if (x == zero||!finite(x)) b = zero;
10625c28e83SPiotr Jasiukajtis 	else if ((GENERIC)n <= x) {
10725c28e83SPiotr Jasiukajtis 					/*
10825c28e83SPiotr Jasiukajtis 					 * Safe to use
10925c28e83SPiotr Jasiukajtis 					 *  J(n+1,x)=2n/x *J(n,x)-J(n-1,x)
11025c28e83SPiotr Jasiukajtis 					 */
111*685c1a21SRichard Lowe 		if (x > 1.0e91) {
11225c28e83SPiotr Jasiukajtis 				/*
11325c28e83SPiotr Jasiukajtis 				 * x >> n**2
11425c28e83SPiotr Jasiukajtis 				 *    Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
11525c28e83SPiotr Jasiukajtis 				 *   Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
11625c28e83SPiotr Jasiukajtis 				 *   Let s=sin(x), c=cos(x),
11725c28e83SPiotr Jasiukajtis 				 *	xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
11825c28e83SPiotr Jasiukajtis 				 *
11925c28e83SPiotr Jasiukajtis 				 *	   n	sin(xn)*sqt2	cos(xn)*sqt2
12025c28e83SPiotr Jasiukajtis 				 *	----------------------------------
12125c28e83SPiotr Jasiukajtis 				 *	   0	 s-c		 c+s
122*685c1a21SRichard Lowe 				 *	   1	-s-c		-c+s
12325c28e83SPiotr Jasiukajtis 				 *	   2	-s+c		-c-s
12425c28e83SPiotr Jasiukajtis 				 *	   3	 s+c		 c-s
12525c28e83SPiotr Jasiukajtis 				 */
126*685c1a21SRichard Lowe 			switch (n&3) {
127*685c1a21SRichard Lowe 			case 0:
128*685c1a21SRichard Lowe 				temp =  cos(x)+sin(x);
129*685c1a21SRichard Lowe 				break;
130*685c1a21SRichard Lowe 			case 1:
131*685c1a21SRichard Lowe 				temp = -cos(x)+sin(x);
132*685c1a21SRichard Lowe 				break;
133*685c1a21SRichard Lowe 			case 2:
134*685c1a21SRichard Lowe 				temp = -cos(x)-sin(x);
135*685c1a21SRichard Lowe 				break;
136*685c1a21SRichard Lowe 			case 3:
137*685c1a21SRichard Lowe 				temp =  cos(x)-sin(x);
138*685c1a21SRichard Lowe 				break;
139*685c1a21SRichard Lowe 			}
140*685c1a21SRichard Lowe 			b = invsqrtpi*temp/sqrt(x);
141*685c1a21SRichard Lowe 		} else {
14225c28e83SPiotr Jasiukajtis 			a = j0(x);
14325c28e83SPiotr Jasiukajtis 			b = j1(x);
14425c28e83SPiotr Jasiukajtis 			for (i = 1; i < n; i++) {
145*685c1a21SRichard Lowe 				temp = b;
146*685c1a21SRichard Lowe 				/* avoid underflow */
147*685c1a21SRichard Lowe 				b = b*((GENERIC)(i+i)/x) - a;
148*685c1a21SRichard Lowe 				a = temp;
14925c28e83SPiotr Jasiukajtis 			}
15025c28e83SPiotr Jasiukajtis 		}
151*685c1a21SRichard Lowe 	} else {
152*685c1a21SRichard Lowe 		if (x < 1e-9) {	/* use J(n,x) = 1/n!*(x/2)^n */
153*685c1a21SRichard Lowe 			b = pow(0.5*x, (GENERIC) n);
154*685c1a21SRichard Lowe 			if (b != zero) {
155*685c1a21SRichard Lowe 				for (a = one, i = 1; i <= n; i++)
156*685c1a21SRichard Lowe 					a *= (GENERIC)i;
157*685c1a21SRichard Lowe 				b = b/a;
158*685c1a21SRichard Lowe 			}
15925c28e83SPiotr Jasiukajtis 		} else {
160*685c1a21SRichard Lowe 			/*
161*685c1a21SRichard Lowe 			 * use backward recurrence
162*685c1a21SRichard Lowe 			 *			x	  x^2	  x^2
163*685c1a21SRichard Lowe 			 *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
164*685c1a21SRichard Lowe 			 *			2n  - 2(n+1) - 2(n+2)
165*685c1a21SRichard Lowe 			 *
166*685c1a21SRichard Lowe 			 *			1	  1	    1
167*685c1a21SRichard Lowe 			 *  (for large x)   =  ----  ------   ------   .....
168*685c1a21SRichard Lowe 			 *			2n   2(n+1)   2(n+2)
169*685c1a21SRichard Lowe 			 *			-- - ------ - ------ -
170*685c1a21SRichard Lowe 			 *			 x	 x		 x
171*685c1a21SRichard Lowe 			 *
172*685c1a21SRichard Lowe 			 * Let w = 2n/x and h = 2/x, then the above quotient
173*685c1a21SRichard Lowe 			 * is equal to the continued fraction:
174*685c1a21SRichard Lowe 			 *		    1
175*685c1a21SRichard Lowe 			 *	= -----------------------
176*685c1a21SRichard Lowe 			 *			   1
177*685c1a21SRichard Lowe 			 *	   w - -----------------
178*685c1a21SRichard Lowe 			 *			  1
179*685c1a21SRichard Lowe 			 *			w+h - ---------
180*685c1a21SRichard Lowe 			 *			   w+2h - ...
181*685c1a21SRichard Lowe 			 *
182*685c1a21SRichard Lowe 			 * To determine how many terms needed, let
183*685c1a21SRichard Lowe 			 * Q(0) = w, Q(1) = w(w+h) - 1,
184*685c1a21SRichard Lowe 			 * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
185*685c1a21SRichard Lowe 			 * When Q(k) > 1e4	good for single
186*685c1a21SRichard Lowe 			 * When Q(k) > 1e9	good for double
187*685c1a21SRichard Lowe 			 * When Q(k) > 1e17	good for quaduple
188*685c1a21SRichard Lowe 			 */
189*685c1a21SRichard Lowe 			/* determine k */
190*685c1a21SRichard Lowe 			GENERIC t, v;
191*685c1a21SRichard Lowe 			double q0, q1, h, tmp;
192*685c1a21SRichard Lowe 			int k, m;
193*685c1a21SRichard Lowe 			w  = (n+n)/(double)x;
194*685c1a21SRichard Lowe 			h = 2.0/(double)x;
195*685c1a21SRichard Lowe 			q0 = w;
196*685c1a21SRichard Lowe 			z = w + h;
197*685c1a21SRichard Lowe 			q1 = w*z - 1.0;
198*685c1a21SRichard Lowe 			k = 1;
199*685c1a21SRichard Lowe 
200*685c1a21SRichard Lowe 			while (q1 < 1.0e9) {
201*685c1a21SRichard Lowe 				k += 1;
202*685c1a21SRichard Lowe 				z += h;
203*685c1a21SRichard Lowe 				tmp = z*q1 - q0;
204*685c1a21SRichard Lowe 				q0 = q1;
205*685c1a21SRichard Lowe 				q1 = tmp;
206*685c1a21SRichard Lowe 			}
207*685c1a21SRichard Lowe 			m = n+n;
208*685c1a21SRichard Lowe 			for (t = zero, i = 2*(n+k); i >= m; i -= 2)
209*685c1a21SRichard Lowe 				t = one/(i/x-t);
210*685c1a21SRichard Lowe 			a = t;
211*685c1a21SRichard Lowe 			b = one;
212*685c1a21SRichard Lowe 			/*
213*685c1a21SRichard Lowe 			 * estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
214*685c1a21SRichard Lowe 			 * hence, if n*(log(2n/x)) > ...
215*685c1a21SRichard Lowe 			 *  single:
216*685c1a21SRichard Lowe 			 *    8.8722839355e+01
217*685c1a21SRichard Lowe 			 *  double:
218*685c1a21SRichard Lowe 			 *    7.09782712893383973096e+02
219*685c1a21SRichard Lowe 			 *  long double:
220*685c1a21SRichard Lowe 			 *    1.1356523406294143949491931077970765006170e+04
221*685c1a21SRichard Lowe 			 * then recurrent value may overflow and the result is
222*685c1a21SRichard Lowe 			 * likely underflow to zero
223*685c1a21SRichard Lowe 			 */
224*685c1a21SRichard Lowe 			tmp = n;
225*685c1a21SRichard Lowe 			v = two/x;
226*685c1a21SRichard Lowe 			tmp = tmp*log(fabs(v*tmp));
227*685c1a21SRichard Lowe 			if (tmp < 7.09782712893383973096e+02) {
228*685c1a21SRichard Lowe 				for (i = n-1; i > 0; i--) {
229*685c1a21SRichard Lowe 					temp = b;
230*685c1a21SRichard Lowe 					b = ((i+i)/x)*b - a;
231*685c1a21SRichard Lowe 					a = temp;
232*685c1a21SRichard Lowe 				}
233*685c1a21SRichard Lowe 			} else {
23425c28e83SPiotr Jasiukajtis 				for (i = n-1; i > 0; i--) {
235*685c1a21SRichard Lowe 					temp = b;
236*685c1a21SRichard Lowe 					b = ((i+i)/x)*b - a;
237*685c1a21SRichard Lowe 					a = temp;
23825c28e83SPiotr Jasiukajtis 					if (b > 1e100) {
23925c28e83SPiotr Jasiukajtis 						a /= b;
24025c28e83SPiotr Jasiukajtis 						t /= b;
24125c28e83SPiotr Jasiukajtis 						b  = 1.0;
24225c28e83SPiotr Jasiukajtis 					}
24325c28e83SPiotr Jasiukajtis 				}
244*685c1a21SRichard Lowe 			}
24525c28e83SPiotr Jasiukajtis 			b = (t*j0(x)/b);
246*685c1a21SRichard Lowe 		}
24725c28e83SPiotr Jasiukajtis 	}
248*685c1a21SRichard Lowe 	if (sgn != 0)
24925c28e83SPiotr Jasiukajtis 		return (-b);
25025c28e83SPiotr Jasiukajtis 	else
25125c28e83SPiotr Jasiukajtis 		return (b);
25225c28e83SPiotr Jasiukajtis }
25325c28e83SPiotr Jasiukajtis 
25425c28e83SPiotr Jasiukajtis GENERIC
yn(int n,GENERIC x)255*685c1a21SRichard Lowe yn(int n, GENERIC x)
256*685c1a21SRichard Lowe {
25725c28e83SPiotr Jasiukajtis 	int i;
25825c28e83SPiotr Jasiukajtis 	int sign;
25925c28e83SPiotr Jasiukajtis 	GENERIC a, b, temp = 0, ox, on;
26025c28e83SPiotr Jasiukajtis 
261*685c1a21SRichard Lowe 	ox = x;
262*685c1a21SRichard Lowe 	on = (GENERIC)n;
26325c28e83SPiotr Jasiukajtis 	if (isnan(x))
26425c28e83SPiotr Jasiukajtis 		return (x*x);	/* + -> * for Cheetah */
26525c28e83SPiotr Jasiukajtis 	if (x <= zero) {
26625c28e83SPiotr Jasiukajtis 		if (x == zero) {
26725c28e83SPiotr Jasiukajtis 			/* return -one/zero; */
26825c28e83SPiotr Jasiukajtis 			return (_SVID_libm_err((GENERIC)n, x, 12));
26925c28e83SPiotr Jasiukajtis 		} else {
27025c28e83SPiotr Jasiukajtis 			/* return zero/zero; */
27125c28e83SPiotr Jasiukajtis 			return (_SVID_libm_err((GENERIC)n, x, 13));
27225c28e83SPiotr Jasiukajtis 		}
27325c28e83SPiotr Jasiukajtis 	}
274*685c1a21SRichard Lowe 	if (!((int)_lib_version == libm_ieee ||
275*685c1a21SRichard Lowe 	    (__xpg6 & _C99SUSv3_math_errexcept) != 0)) {
276*685c1a21SRichard Lowe 		if (x > X_TLOSS)
27725c28e83SPiotr Jasiukajtis 			return (_SVID_libm_err(on, ox, 39));
27825c28e83SPiotr Jasiukajtis 	}
27925c28e83SPiotr Jasiukajtis 	sign = 1;
28025c28e83SPiotr Jasiukajtis 	if (n < 0) {
28125c28e83SPiotr Jasiukajtis 		n = -n;
28225c28e83SPiotr Jasiukajtis 		if ((n&1) == 1) sign = -1;
28325c28e83SPiotr Jasiukajtis 	}
28425c28e83SPiotr Jasiukajtis 	if (n == 0)
28525c28e83SPiotr Jasiukajtis 		return (y0(x));
28625c28e83SPiotr Jasiukajtis 	if (n == 1)
28725c28e83SPiotr Jasiukajtis 		return (sign*y1(x));
28825c28e83SPiotr Jasiukajtis 	if (!finite(x))
28925c28e83SPiotr Jasiukajtis 		return (zero);
29025c28e83SPiotr Jasiukajtis 
29125c28e83SPiotr Jasiukajtis 	if (x > 1.0e91) {
29225c28e83SPiotr Jasiukajtis 				/*
29325c28e83SPiotr Jasiukajtis 				 * x >> n**2
29425c28e83SPiotr Jasiukajtis 				 *  Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
29525c28e83SPiotr Jasiukajtis 				 *  Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
29625c28e83SPiotr Jasiukajtis 				 *  Let s = sin(x), c = cos(x),
29725c28e83SPiotr Jasiukajtis 				 *  xn = x-(2n+1)*pi/4, sqt2 = sqrt(2), then
29825c28e83SPiotr Jasiukajtis 				 *
29925c28e83SPiotr Jasiukajtis 				 *    n	sin(xn)*sqt2	cos(xn)*sqt2
30025c28e83SPiotr Jasiukajtis 				 *	----------------------------------
30125c28e83SPiotr Jasiukajtis 				 *	 0	 s-c		 c+s
302*685c1a21SRichard Lowe 				 *	 1	-s-c		-c+s
30325c28e83SPiotr Jasiukajtis 				 *	 2	-s+c		-c-s
30425c28e83SPiotr Jasiukajtis 				 *	 3	 s+c		 c-s
30525c28e83SPiotr Jasiukajtis 				 */
30625c28e83SPiotr Jasiukajtis 		switch (n&3) {
307*685c1a21SRichard Lowe 		case 0:
308*685c1a21SRichard Lowe 			temp =  sin(x)-cos(x);
309*685c1a21SRichard Lowe 			break;
310*685c1a21SRichard Lowe 		case 1:
311*685c1a21SRichard Lowe 			temp = -sin(x)-cos(x);
312*685c1a21SRichard Lowe 			break;
313*685c1a21SRichard Lowe 		case 2:
314*685c1a21SRichard Lowe 			temp = -sin(x)+cos(x);
315*685c1a21SRichard Lowe 			break;
316*685c1a21SRichard Lowe 		case 3:
317*685c1a21SRichard Lowe 			temp =  sin(x)+cos(x);
318*685c1a21SRichard Lowe 			break;
31925c28e83SPiotr Jasiukajtis 		}
32025c28e83SPiotr Jasiukajtis 		b = invsqrtpi*temp/sqrt(x);
32125c28e83SPiotr Jasiukajtis 	} else {
32225c28e83SPiotr Jasiukajtis 		a = y0(x);
32325c28e83SPiotr Jasiukajtis 		b = y1(x);
32425c28e83SPiotr Jasiukajtis 		/*
32525c28e83SPiotr Jasiukajtis 		 * fix 1262058 and take care of non-default rounding
32625c28e83SPiotr Jasiukajtis 		 */
32725c28e83SPiotr Jasiukajtis 		for (i = 1; i < n; i++) {
32825c28e83SPiotr Jasiukajtis 			temp = b;
32925c28e83SPiotr Jasiukajtis 			b *= (GENERIC) (i + i) / x;
33025c28e83SPiotr Jasiukajtis 			if (b <= -DBL_MAX)
33125c28e83SPiotr Jasiukajtis 				break;
33225c28e83SPiotr Jasiukajtis 			b -= a;
33325c28e83SPiotr Jasiukajtis 			a = temp;
33425c28e83SPiotr Jasiukajtis 		}
33525c28e83SPiotr Jasiukajtis 	}
33625c28e83SPiotr Jasiukajtis 	if (sign > 0)
33725c28e83SPiotr Jasiukajtis 		return (b);
33825c28e83SPiotr Jasiukajtis 	else
33925c28e83SPiotr Jasiukajtis 		return (-b);
34025c28e83SPiotr Jasiukajtis }
341