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 = 0
  52         if (x > 0) scale = (0.435*x)/1
  53         scale = scale + t + length(scale + t) + 1
  54
  55         w = 0
  56         if (x < 0) {
  57                 x = -x
  58                 w = 1
  59         }
  60         y = 0
  61         while (x > 2) {
  62                 x = x/2
  63                 y = y + 1
  64         }
  65
  66         a = 1
  67         b = 1
  68         c = b
  69         d = 1
  70         e = 1
  71         for (a = 1; 1 == 1; a++) {
  72                 b = b*x
  73                 c = c*a + b
  74                 d = d*a
  75                 g = c/d
  76                 if (g == e) {
  77                         g = g/1
  78                         while (y--) {
  79                                 g = g*g
  80                         }
  81                         scale = t
  82                         ibase = r
  83                         if (w == 1) return (1/g)
  84                         return (g/1)
  85                 }
  86                 e = g
  87         }
  88 }
  89
  90 define l(x) {
  91         auto a, b, c, d, e, f, g, u, s, t, r
  92         r = ibase
  93         ibase = A
  94         if (x <= 0) {
  95                 a = (1 - 10^scale)
  96                 ibase = r
  97                 return (a)
  98         }
  99         t = scale
 100
 101         f = 1
 102         if (x < 1) {
 103                 s = scale(x)
 104         } else {
 105                 s = length(x)-scale(x)
 106         }
 107         scale = 0
 108         a = (2.31*s)/1 /* estimated integer part of the answer */
 109         s = t + length(a) + 2 /* estimated length of the answer */
 110         while (x > 2) {
 111                 scale = 0
 112                 scale = (length(x) + scale(x))/2 + 1
 113                 if (scale < s) scale = s
 114                 x = sqrt(x)
 115                 f = f*2
 116         }
 117         while (x < .5) {
 118                 scale = 0
 119                 scale = scale(x)/2 + 1
 120                 if (scale < s) scale = s
 121                 x = sqrt(x)
 122                 f = f*2
 123         }
 124
 125         scale = 0
 126         scale = t + length(f) + length((1.05*(t+length(f))/1)) + 1
 127         u = (x - 1)/(x + 1)
 128         s = u*u
 129         scale = t + 2
 130         b = 2*f
 131         c = b
 132         d = 1
 133         e = 1
 134         for (a = 3; 1 == 1 ; a = a + 2) {
 135                 b = b*s
 136                 c = c*a + d*b
 137                 d = d*a
 138                 g = c/d
 139                 if (g == e) {
 140                         scale = t
 141                         ibase = r
 142                         return (u*c/d)
 143                 }
 144                 e = g
 145         }
 146 }
 147
 148 define s(x) {
 149         auto a, b, c, s, t, y, p, n, i, r
 150         r = ibase
 151         ibase = A
 152         t = scale
 153         y = x/.7853
 154         s = t + length(y) - scale(y)
 155         if (s < t) s = t
 156         scale = s
 157         p = a(1)
 158
 159         scale = 0
 160         if (x >= 0) n = (x/(2*p) + 1)/2
 161         if (x < 0) n = (x/(2*p) - 1)/2
 162         x = x - 4*n*p
 163         if (n % 2 != 0) x = -x
 164
 165         scale = t + length(1.2*t) - scale(1.2*t)
 166         y = -x*x
 167         a = x
 168         b = 1
 169         s = x
 170         for (i =3 ; 1 == 1; i = i + 2) {
 171                 a = a*y
 172                 b = b*i*(i - 1)
 173                 c = a/b
 174                 if (c == 0) {
 175                         scale = t
 176                         ibase = r
 177                         return (s/1)
 178                 }
 179                 s = s + c
 180         }
 181 }
 182
 183 define c(x) {
 184         auto t, r
 185         r = ibase
 186         ibase = A
 187         t = scale
 188         scale = scale + 1
 189         x = s(x + 2*a(1))
 190         scale = t
 191         ibase = r
 192         return (x/1)
 193 }
 194
 195 define a(x) {
 196         auto a, b, c, d, e, f, g, s, t, r
 197         if (x == 0) return(0)
 198
 199         r = ibase
 200         ibase = A
 201         if (x == 1) {
 202                 if (scale < 52) {
 203                          a = .7853981633974483096156608458198757210492923498437764/1
 204                          ibase = r
 205                          return (a)
 206                 }
 207         }
 208         t = scale
 209         f = 1
 210         while (x > .5) {
 211                 scale = scale + 1
 212                 x = -(1 - sqrt(1. + x*x))/x
 213                 f = f*2
 214         }
 215         while (x < -.5) {
 216                 scale = scale + 1
 217                 x = -(1 - sqrt(1. + x*x))/x
 218                 f = f*2
 219         }
 220         s = -x*x
 221         b = f
 222         c = f
 223         d = 1
 224         e = 1
 225         for (a = 3; 1 == 1; a = a + 2) {
 226                 b = b*s
 227                 c = c*a + d*b
 228                 d = d*a
 229                 g = c/d
 230                 if (g == e) {
 231                         ibase = r
 232                         scale = t
 233                         return (x*c/d)
 234                 }
 235                 e = g
 236         }
 237 }
 238
 239 define j(n,x) {
 240         auto a, b, c, d, e, g, i, s, k, t, r
 241
 242         r = ibase
 243         ibase = A
 244         t = scale
 245         k = 1.36*x + 1.16*t - n
 246         k = length(k) - scale(k)
 247         if (k > 0) scale = scale + k
 248
 249         s = -x*x/4
 250         if (n < 0) {
 251                 n = -n
 252                 x = -x
 253         }
 254         a = 1
 255         c = 1
 256         for (i = 1; i <= n; i++) {
 257                 a = a*x
 258                 c = c*2*i
 259         }
 260         b = a
 261         d = 1
 262         e = 1
 263         for (i = 1; 1; i++) {
 264                 a = a*s
 265                 b = b*i*(n + i) + a
 266                 c = c*i*(n + i)
 267                 g = b/c
 268                 if (g == e) {
 269                         ibase = r
 270                         scale = t
 271                         return (g/1)
 272                 }
 273                 e = g
 274         }
 275 }