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