jcs's openbsd hax
openbsd
1/* $OpenBSD: bc.library,v 1.4 2012/03/14 07:35:53 otto Exp $ */
2
3/*
4 * Copyright (C) Caldera International Inc. 2001-2002.
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code and documentation must retain the above
11 * copyright notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * 3. All advertising materials mentioning features or use of this software
16 * must display the following acknowledgement:
17 * This product includes software developed or owned by Caldera
18 * International, Inc.
19 * 4. Neither the name of Caldera International, Inc. nor the names of other
20 * contributors may be used to endorse or promote products derived from
21 * this software without specific prior written permission.
22 *
23 * USE OF THE SOFTWARE PROVIDED FOR UNDER THIS LICENSE BY CALDERA
24 * INTERNATIONAL, INC. AND CONTRIBUTORS ``AS IS'' AND ANY EXPRESS OR
25 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
26 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
27 * IN NO EVENT SHALL CALDERA INTERNATIONAL, INC. BE LIABLE FOR ANY DIRECT,
28 * INDIRECT INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
29 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
30 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
31 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
32 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
33 * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34 * POSSIBILITY OF SUCH DAMAGE.
35 */
36
37/*
38 * @(#)bc.library 5.1 (Berkeley) 4/17/91
39 */
40
41scale = 20
42define e(x) {
43 auto a, b, c, d, e, g, t, w, y, r
44
45 r = ibase
46 ibase = A
47 t = scale
48 scale = 0
49 if (x > 0) scale = (0.435*x)/1
50 scale = scale + t + length(scale + t) + 1
51
52 w = 0
53 if (x < 0) {
54 x = -x
55 w = 1
56 }
57 y = 0
58 while (x > 2) {
59 x = x/2
60 y = y + 1
61 }
62
63 a = 1
64 b = 1
65 c = b
66 d = 1
67 e = 1
68 for (a = 1; 1 == 1; a++) {
69 b = b*x
70 c = c*a + b
71 d = d*a
72 g = c/d
73 if (g == e) {
74 g = g/1
75 while (y--) {
76 g = g*g
77 }
78 scale = t
79 ibase = r
80 if (w == 1) return (1/g)
81 return (g/1)
82 }
83 e = g
84 }
85}
86
87define l(x) {
88 auto a, b, c, d, e, f, g, u, s, t, r
89 r = ibase
90 ibase = A
91 if (x <= 0) {
92 a = (1 - 10^scale)
93 ibase = r
94 return (a)
95 }
96 t = scale
97
98 f = 1
99 if (x < 1) {
100 s = scale(x)
101 } else {
102 s = length(x)-scale(x)
103 }
104 scale = 0
105 a = (2.31*s)/1 /* estimated integer part of the answer */
106 s = t + length(a) + 2 /* estimated length of the answer */
107 while (x > 2) {
108 scale = 0
109 scale = (length(x) + scale(x))/2 + 1
110 if (scale < s) scale = s
111 x = sqrt(x)
112 f = f*2
113 }
114 while (x < .5) {
115 scale = 0
116 scale = scale(x)/2 + 1
117 if (scale < s) scale = s
118 x = sqrt(x)
119 f = f*2
120 }
121
122 scale = 0
123 scale = t + length(f) + length((1.05*(t+length(f))/1)) + 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}