Annotation of src/usr.bin/bc/bc.library, Revision 1.4
1.4 ! otto 1: /* $OpenBSD: bc.library,v 1.3 2007/02/03 21:15:06 otto Exp $ */
1.1 otto 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:
41: scale = 20
1.2 otto 42: define e(x) {
1.3 otto 43: auto a, b, c, d, e, g, t, w, y, r
1.1 otto 44:
1.3 otto 45: r = ibase
46: ibase = A
1.1 otto 47: t = scale
1.4 ! otto 48: scale = 0
! 49: if (x > 0) scale = (0.435*x)/1
! 50: scale = scale + t + length(scale + t) + 1
1.1 otto 51:
52: w = 0
1.2 otto 53: if (x < 0) {
1.1 otto 54: x = -x
55: w = 1
56: }
57: y = 0
1.2 otto 58: while (x > 2) {
1.1 otto 59: x = x/2
60: y = y + 1
61: }
62:
1.2 otto 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
1.1 otto 72: g = c/d
1.2 otto 73: if (g == e) {
1.1 otto 74: g = g/1
1.2 otto 75: while (y--) {
1.1 otto 76: g = g*g
77: }
78: scale = t
1.3 otto 79: ibase = r
80: if (w == 1) return (1/g)
81: return (g/1)
1.1 otto 82: }
1.2 otto 83: e = g
1.1 otto 84: }
85: }
86:
1.2 otto 87: define l(x) {
1.3 otto 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: }
1.1 otto 96: t = scale
97:
1.2 otto 98: f = 1
1.4 ! otto 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 */
1.2 otto 107: while (x > 2) {
1.4 ! otto 108: scale = 0
! 109: scale = (length(x) + scale(x))/2 + 1
! 110: if (scale < s) scale = s
1.1 otto 111: x = sqrt(x)
1.2 otto 112: f = f*2
1.1 otto 113: }
1.2 otto 114: while (x < .5) {
1.4 ! otto 115: scale = 0
! 116: scale = scale(x)/2 + 1
! 117: if (scale < s) scale = s
1.1 otto 118: x = sqrt(x)
1.2 otto 119: f = f*2
1.1 otto 120: }
121:
1.4 ! otto 122: scale = 0
! 123: scale = t + length(f) + length((1.05*(t+length(f))/1)) + 1
1.2 otto 124: u = (x - 1)/(x + 1)
1.1 otto 125: s = u*u
1.4 ! otto 126: scale = t + 2
1.1 otto 127: b = 2*f
128: c = b
129: d = 1
130: e = 1
1.2 otto 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) {
1.1 otto 137: scale = t
1.3 otto 138: ibase = r
1.2 otto 139: return (u*c/d)
1.1 otto 140: }
1.2 otto 141: e = g
1.1 otto 142: }
143: }
144:
1.2 otto 145: define s(x) {
1.3 otto 146: auto a, b, c, s, t, y, p, n, i, r
147: r = ibase
148: ibase = A
1.1 otto 149: t = scale
150: y = x/.7853
151: s = t + length(y) - scale(y)
1.2 otto 152: if (s < t) s = t
1.1 otto 153: scale = s
154: p = a(1)
155:
156: scale = 0
1.2 otto 157: if (x >= 0) n = (x/(2*p) + 1)/2
158: if (x < 0) n = (x/(2*p) - 1)/2
1.1 otto 159: x = x - 4*n*p
1.2 otto 160: if (n % 2 != 0) x = -x
1.1 otto 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
1.2 otto 167: for (i =3 ; 1 == 1; i = i + 2) {
1.1 otto 168: a = a*y
1.2 otto 169: b = b*i*(i - 1)
1.1 otto 170: c = a/b
1.2 otto 171: if (c == 0) {
172: scale = t
1.3 otto 173: ibase = r
1.2 otto 174: return (s/1)
175: }
176: s = s + c
1.1 otto 177: }
178: }
179:
1.2 otto 180: define c(x) {
1.3 otto 181: auto t, r
182: r = ibase
183: ibase = A
1.1 otto 184: t = scale
1.2 otto 185: scale = scale + 1
186: x = s(x + 2*a(1))
1.1 otto 187: scale = t
1.3 otto 188: ibase = r
1.2 otto 189: return (x/1)
1.1 otto 190: }
191:
1.2 otto 192: define a(x) {
1.3 otto 193: auto a, b, c, d, e, f, g, s, t, r
1.2 otto 194: if (x == 0) return(0)
1.3 otto 195:
196: r = ibase
197: ibase = A
1.2 otto 198: if (x == 1) {
199: if (scale < 52) {
1.3 otto 200: a = .7853981633974483096156608458198757210492923498437764/1
201: ibase = r
202: return (a)
1.1 otto 203: }
204: }
205: t = scale
1.2 otto 206: f = 1
207: while (x > .5) {
1.1 otto 208: scale = scale + 1
1.2 otto 209: x = -(1 - sqrt(1. + x*x))/x
210: f = f*2
1.1 otto 211: }
1.2 otto 212: while (x < -.5) {
1.1 otto 213: scale = scale + 1
1.2 otto 214: x = -(1 - sqrt(1. + x*x))/x
215: f = f*2
1.1 otto 216: }
217: s = -x*x
218: b = f
219: c = f
220: d = 1
221: e = 1
1.2 otto 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) {
1.3 otto 228: ibase = r
1.1 otto 229: scale = t
1.2 otto 230: return (x*c/d)
1.1 otto 231: }
1.2 otto 232: e = g
1.1 otto 233: }
234: }
235:
1.2 otto 236: define j(n,x) {
1.3 otto 237: auto a, b, c, d, e, g, i, s, k, t, r
1.1 otto 238:
1.3 otto 239: r = ibase
240: ibase = A
1.1 otto 241: t = scale
242: k = 1.36*x + 1.16*t - n
243: k = length(k) - scale(k)
1.2 otto 244: if (k > 0) scale = scale + k
1.1 otto 245:
1.2 otto 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
1.1 otto 255: c = c*2*i
256: }
1.2 otto 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) {
1.3 otto 266: ibase = r
1.1 otto 267: scale = t
1.2 otto 268: return (g/1)
1.1 otto 269: }
1.2 otto 270: e = g
1.1 otto 271: }
272: }