sys/lib/bclib

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