usr.bin/bc/bc.library

   1 /*
   2  * $OpenBSD: bc.library,v 1.1 2003/09/25 19:34:22 otto Exp $
   3  * $DragonFly: src/usr.bin/bc/bc.library,v 1.1 2004/09/20 04:20:34 dillon 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
  47
  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                         if(w==1) return(1/g)
  79                         return(g/1)
  80                 }
  81                 e=g
  82         }
  83 }
  84
  85 define l(x){
  86         auto a, b, c, d, e, f, g, u, s, t
  87         if(x <=0) return(1-10^scale)
  88         t = scale
  89
  90         f=1
  91         scale = scale + scale(x) - length(x) + 1
  92         s=scale
  93         while(x > 2){
  94                 s = s + (length(x)-scale(x))/2 + 1
  95                 if(s>0) scale = s
  96                 x = sqrt(x)
  97                 f=f*2
  98         }
  99         while(x < .5){
 100                 s = s + (length(x)-scale(x))/2 + 1
 101                 if(s>0) scale = s
 102                 x = sqrt(x)
 103                 f=f*2
 104         }
 105
 106         scale = t + length(f) - scale(f) + 1
 107         u = (x-1)/(x+1)
 108
 109         scale = scale + 1.1*length(t) - 1.1*scale(t)
 110         s = u*u
 111         b = 2*f
 112         c = b
 113         d = 1
 114         e = 1
 115         for(a=3;1==1;a=a+2){
 116                 b=b*s
 117                 c=c*a+d*b
 118                 d=d*a
 119                 g=c/d
 120                 if(g==e){
 121                         scale = t
 122                         return(u*c/d)
 123                 }
 124                 e=g
 125         }
 126 }
 127
 128 define s(x){
 129         auto a, b, c, s, t, y, p, n, i
 130         t = scale
 131         y = x/.7853
 132         s = t + length(y) - scale(y)
 133         if(s<t) s=t
 134         scale = s
 135         p = a(1)
 136
 137         scale = 0
 138         if(x>=0) n = (x/(2*p)+1)/2
 139         if(x<0) n = (x/(2*p)-1)/2
 140         x = x - 4*n*p
 141         if(n%2!=0) x = -x
 142
 143         scale = t + length(1.2*t) - scale(1.2*t)
 144         y = -x*x
 145         a = x
 146         b = 1
 147         s = x
 148         for(i=3; 1==1; i=i+2){
 149                 a = a*y
 150                 b = b*i*(i-1)
 151                 c = a/b
 152                 if(c==0){scale=t; return(s/1)}
 153                 s = s+c
 154         }
 155 }
 156
 157 define c(x){
 158         auto t
 159         t = scale
 160         scale = scale+1
 161         x = s(x+2*a(1))
 162         scale = t
 163         return(x/1)
 164 }
 165
 166 define a(x){
 167         auto a, b, c, d, e, f, g, s, t
 168         if(x==0) return(0)
 169         if(x==1) {
 170                 if(scale<52) {
 171                         return(.7853981633974483096156608458198757210492923498437764/1)
 172                 }
 173         }
 174         t = scale
 175         f=1
 176         while(x > .5){
 177                 scale = scale + 1
 178                 x= -(1-sqrt(1.+x*x))/x
 179                 f=f*2
 180         }
 181         while(x < -.5){
 182                 scale = scale + 1
 183                 x = -(1-sqrt(1.+x*x))/x
 184                 f=f*2
 185         }
 186         s = -x*x
 187         b = f
 188         c = f
 189         d = 1
 190         e = 1
 191         for(a=3;1==1;a=a+2){
 192                 b=b*s
 193                 c=c*a+d*b
 194                 d=d*a
 195                 g=c/d
 196                 if(g==e){
 197                         scale = t
 198                         return(x*c/d)
 199                 }
 200                 e=g
 201         }
 202 }
 203
 204 define j(n,x){
 205         auto a,b,c,d,e,g,i,s,k,t
 206
 207         t = scale
 208         k = 1.36*x + 1.16*t - n
 209         k = length(k) - scale(k)
 210         if(k>0) scale = scale + k
 211
 212         s= -x*x/4
 213         if(n<0){
 214                 n= -n
 215                 x= -x
 216         }
 217         a=1
 218         c=1
 219         for(i=1;i<=n;i++){
 220                 a=a*x
 221                 c = c*2*i
 222         }
 223         b=a
 224         d=1
 225         e=1
 226         for(i=1;1;i++){
 227                 a=a*s
 228                 b=b*i*(n+i) + a
 229                 c=c*i*(n+i)
 230                 g=b/c
 231                 if(g==e){
 232                         scale = t
 233                         return(g/1)
 234                 }
 235                 e=g
 236         }
 237 }