usr.bin/bc/bc.library

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