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