bc.library revision 232994
1/*	$FreeBSD: head/usr.bin/bc/bc.library 232994 2012-03-15 01:43:44Z kevlo $							*/
2/*      $OpenBSD: bc.library,v 1.3 2007/02/03 21:15:06 otto Exp $	*/
3
4/*
5 * Copyright (C) Caldera International Inc.  2001-2002.
6 * All rights reserved.
7 *
8 * Redistribution and use in source and binary forms, with or without
9 * modification, are permitted provided that the following conditions
10 * are met:
11 * 1. Redistributions of source code and documentation must retain the above
12 *    copyright notice, this list of conditions and the following disclaimer.
13 * 2. Redistributions in binary form must reproduce the above copyright
14 *    notice, this list of conditions and the following disclaimer in the
15 *    documentation and/or other materials provided with the distribution.
16 * 3. All advertising materials mentioning features or use of this software
17 *    must display the following acknowledgement:
18 *      This product includes software developed or owned by Caldera
19 *      International, Inc.
20 * 4. Neither the name of Caldera International, Inc. nor the names of other
21 *    contributors may be used to endorse or promote products derived from
22 *    this software without specific prior written permission.
23 *
24 * USE OF THE SOFTWARE PROVIDED FOR UNDER THIS LICENSE BY CALDERA
25 * INTERNATIONAL, INC. AND CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR
26 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
27 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
28 * IN NO EVENT SHALL CALDERA INTERNATIONAL, INC. BE LIABLE FOR ANY DIRECT,
29 * INDIRECT INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
30 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
32 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
33 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
34 * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
35 * POSSIBILITY OF SUCH DAMAGE.
36 */
37
38/*
39 *	@(#)bc.library	5.1 (Berkeley) 4/17/91
40 */
41
42scale = 20
43define e(x) {
44	auto a, b, c, d, e, g, t, w, y, r
45
46	r = ibase
47	ibase = A
48	t = scale
49	scale = 0
50	if (x > 0) scale = (0.435*x)/1
51	scale = scale + t + 1
52
53	w = 0
54	if (x < 0) {
55		x = -x
56		w = 1
57	}
58	y = 0
59	while (x > 2) {
60		x = x/2
61		y = y + 1
62	}
63
64	a = 1
65	b = 1
66	c = b
67	d = 1
68	e = 1
69	for (a = 1; 1 == 1; a++) {
70		b = b*x
71		c = c*a + b
72		d = d*a
73		g = c/d
74		if (g == e) {
75			g = g/1
76			while (y--) {
77				g = g*g
78			}
79			scale = t
80			ibase = r
81			if (w == 1) return (1/g)
82			return (g/1)
83		}
84		e = g
85	}
86}
87
88define l(x) {
89	auto a, b, c, d, e, f, g, u, s, t, r
90	r = ibase
91	ibase = A
92	if (x <= 0) {
93		a = (1 - 10^scale)
94		ibase = r
95		return (a)
96	}
97	t = scale
98
99	f = 1
100	if (x < 1) {
101		s = scale(x)
102	} else {
103		s = length(x) - scale(x)
104	}
105	scale = 0
106	a = (2.31*s)/1 /* estimated integer part of the answer */
107	s = t + length(a) + 2 /* estimated length of the answer */
108	while (x > 2) {
109		scale=0
110		scale = (length(x) + scale(x))/2 + 1
111		if (scale < s) scale = s
112		x = sqrt(x)
113		f = f*2
114	}
115	while (x < .5) {
116		scale = 0
117		scale = scale(x)/2 + 1
118		if (scale < s) scale = s
119		x = sqrt(x)
120		f = f*2
121	}
122
123	scale = t + length(f) + length(t + length(f)) + 1
124	u = (x - 1)/(x + 1)
125	s = u*u
126	scale = t + 2
127	b = 2*f
128	c = b
129	d = 1
130	e = 1
131	for (a = 3; 1 == 1 ; a = a + 2) {
132		b = b*s
133		c = c*a + d*b
134		d = d*a
135		g = c/d
136		if (g == e) {
137			scale = t
138			ibase = r
139			return (u*c/d)
140		}
141		e = g
142	}
143}
144
145define s(x) {
146	auto a, b, c, s, t, y, p, n, i, r
147	r = ibase
148	ibase = A
149	t = scale
150	y = x/.7853
151	s = t + length(y) - scale(y)
152	if (s < t) s = t
153	scale = s
154	p = a(1)
155
156	scale = 0
157	if (x >= 0) n = (x/(2*p) + 1)/2
158	if (x < 0) n = (x/(2*p) - 1)/2
159	x = x - 4*n*p
160	if (n % 2 != 0) x = -x
161
162	scale = t + length(1.2*t) - scale(1.2*t)
163	y = -x*x
164	a = x
165	b = 1
166	s = x
167	for (i =3 ; 1 == 1; i = i + 2) {
168		a = a*y
169		b = b*i*(i - 1)
170		c = a/b
171		if (c == 0) {
172			scale = t
173			ibase = r
174			return (s/1)
175		}
176		s = s + c
177	}
178}
179
180define c(x) {
181	auto t, r
182	r = ibase
183	ibase = A
184	t = scale
185	scale = scale + 1
186	x = s(x + 2*a(1))
187	scale = t
188	ibase = r
189	return (x/1)
190}
191
192define a(x) {
193	auto a, b, c, d, e, f, g, s, t, r
194	if (x == 0) return(0)
195
196	r = ibase
197	ibase = A
198	if (x == 1) {
199		if (scale < 52) {
200			 a = .7853981633974483096156608458198757210492923498437764/1
201			 ibase = r
202			 return (a)
203		}
204	}
205	t = scale
206	f = 1
207	while (x > .5) {
208		scale = scale + 1
209		x = -(1 - sqrt(1. + x*x))/x
210		f = f*2
211	}
212	while (x < -.5) {
213		scale = scale + 1
214		x = -(1 - sqrt(1. + x*x))/x
215		f = f*2
216	}
217	s = -x*x
218	b = f
219	c = f
220	d = 1
221	e = 1
222	for (a = 3; 1 == 1; a = a + 2) {
223		b = b*s
224		c = c*a + d*b
225		d = d*a
226		g = c/d
227		if (g == e) {
228			ibase = r
229			scale = t
230			return (x*c/d)
231		}
232		e = g
233	}
234}
235
236define j(n,x) {
237	auto a, b, c, d, e, g, i, s, k, t, r
238
239	r = ibase
240	ibase = A
241	t = scale
242	k = 1.36*x + 1.16*t - n
243	k = length(k) - scale(k)
244	if (k > 0) scale = scale + k
245
246	s = -x*x/4
247	if (n < 0) {
248		n = -n
249		x = -x
250	}
251	a = 1
252	c = 1
253	for (i = 1; i <= n; i++) {
254		a = a*x
255		c = c*2*i
256	}
257	b = a
258	d = 1
259	e = 1
260	for (i = 1; 1; i++) {
261		a = a*s
262		b = b*i*(n + i) + a
263		c = c*i*(n + i)
264		g = b/c
265		if (g == e) {
266			ibase = r
267			scale = t
268			return (g/1)
269		}
270		e = g
271	}
272}
273/* vim: set filetype=bc shiftwidth=8 noexpandtab: */
274