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