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