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