usr.bin/bc/bc.library

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