bc.library (4292B)
1 /* 2 * The content of this file is free knowledge. 3 * No copyright applies. No warranty is provided. 4 * 5 * The following functions are available (aliases in parentheses): 6 * General utility: max, min, sgn, abs, floor, ceiling 7 * Arithmetic: factorial(fact), binomial(binom, bin), gcd, lcm, totient(phi) 8 * Calculus: exp, log(ln), pow, root, sin, cos, tan, cosh, sinh, tanh 9 * atan, atan2, asin, acos, atanh, asinh, acosh 10 * Constants (as functions): e(), pi() 11 * 12 * Approximation of inverse trigonometric functions and of pi is not good. 13 * 14 * For functions returning non-integer values, remember to set the 15 * desired scaled value before calling the function. 16 */ 17 18 /* General utility functions */ 19 20 define max(x, y) { 21 if (x >= y) return x 22 return y 23 } 24 25 define min(x, y) { 26 if (x <= y) return x 27 return y 28 } 29 30 define sgn(c) { 31 if (x > 0) return 1 32 if (x < 0) return -1 33 return 0 34 } 35 36 define abs(x) { 37 return x * sgn(x) 38 } 39 40 define floor(x) { 41 auto s, y 42 if (x < 0) return -ceiling(-x) 43 s = scale 44 scale = 0 45 y = x / 1 46 scale = s 47 return y 48 } 49 50 define ceiling(x) { 51 if (x < 0) return -floor(-x) 52 if (floor(x) == x) return x 53 return 1 + floor(x) 54 } 55 56 /* Arithmetic */ 57 58 define factorial(n) { 59 auto i, res, s 60 s = scale 61 scale = 0 62 res = 1 63 for (i = 1; i <= n; i++) res *= i 64 scale = s 65 return res 66 } 67 68 define fact(n) { 69 return factorial(n) 70 } 71 72 define binomial(n, k) { 73 auto i, j, told[], tnew[], s 74 if (k < 0) return -1 75 if (k > n) return 0 76 s = scale 77 scale = 0 78 tnew[0] = 1 79 for (i = 1; i <= n; i++) { 80 for (j = 0; j <= i; j++) told[j] = tnew[j] 81 for (j = 1; j <= i; j++) tnew[j] = told[j] + told[j-1] 82 } 83 scale = s 84 return tnew[k] 85 } 86 87 define binom(n, k) { 88 return binomial(n, k) 89 } 90 91 define bin(n, k) { 92 return binomial(n, k) 93 } 94 95 define gcd(a, b) { 96 auto s, aux 97 s = scale 98 scale = 0 99 a /= 1 100 b /= 1 101 while (b != 0) { 102 aux = a 103 a = b 104 b = aux % b 105 } 106 scale = s 107 return a 108 } 109 110 define lcm(a, b) { 111 if (a == 0 || b == 0) return 0 112 return a * b / gcd(a, b) 113 } 114 115 define phi(n) { 116 auto i, j, f[], r, s 117 s = scale 118 scale = 0 119 j = 0 120 r = n 121 for (i = 2; i*i <= r; i++) { 122 if (n % i == 0) { 123 f[j] = i 124 j += 1 125 while (n % i == 0) n /= i 126 } 127 } 128 for (i = 0; i < j; i++) r -= r / f[i] 129 scale = s 130 return r 131 } 132 133 define totient(n) { 134 return phi(n) 135 } 136 137 /* Calculus */ 138 139 define exp(x) { 140 auto i, n, d, series 141 scale += 10 142 i = 0 143 n = 1 144 d = 1 145 series = 0 146 while (n/d != 0) { 147 series += n/d 148 i += 1 149 n *= x 150 d *= i 151 } 152 scale -= 10 153 return series / 1 154 } 155 156 define e() { 157 return exp(1) 158 } 159 160 /* Log: first take some square roots to make x smaller, then 161 use Newton's method to solve e^y-x=0 */ 162 define log(x) { 163 auto m, s, t, n, d 164 scale += 10 165 while (x > 3) { 166 m += 1 167 x = sqrt(x) 168 } 169 t = x 170 d = exp(t) 171 n = d - x 172 while (n/d != 0) { 173 t -= n/d 174 d = exp(t) 175 n = d - x 176 } 177 scale -= 10 178 return (t * 2^m) / 1 179 } 180 181 define ln(x) { 182 return log(x) 183 } 184 185 define pow(a, b) { 186 if (b != floor(b)) return (2^floor(b)) * exp(log(a)*(b-floor(b))) 187 return a ^ b 188 } 189 190 define root(a, b) { 191 auto res 192 scale += 5 193 res = pow(b, 1/a) 194 scale -= 5 195 return res / 1 196 } 197 198 define trig(x, ii, in, id) { 199 auto i, n, d, series 200 scale += 10 201 i = ii 202 n = in 203 d = id 204 while (n/d != 0) { 205 series += n/d 206 i += 2 207 n *= -x*x 208 d *= i*(i-1) 209 } 210 scale -= 10 211 return series / 1 212 } 213 214 define sin(x) { 215 return trig(x, 1, x, 1) 216 } 217 218 define cos(x) { 219 return trig(x, 0, 1, 1) 220 } 221 222 define tan(x) { 223 return sin(x) / cos(x) 224 } 225 226 define sinh(x) { 227 return 0.5*(exp(x) - exp(-x)) 228 } 229 230 define cosh(x) { 231 return 0.5*(exp(x) + exp(-x)) 232 } 233 234 define tanh(x) { 235 return sinh(x) / cosh(x) 236 } 237 238 /* Atan: approximate integral of 1/(1+x^2) */ 239 define atan(x) { 240 auto i, n, f1, f2, a, res, s 241 s = scale 242 scale = 10 243 n = 10^4 244 res = 0 245 for (i = 0; i < n; i++) { 246 f1 = 1 + (x*i/n)^2 247 f2 = 1 + (x*(i+1)/n)^2 248 a = (f1 + f2) / (2 * f1 * f2) 249 res += a * x / n 250 } 251 scale = s 252 return res/1 253 } 254 255 define pi() { 256 return 4 * atan(1) 257 } 258 259 define atan2(y, x) { 260 if (x > 0) return atan(y/x) 261 if (x < 0 && y >= 0) return atan(y/x) + pi() 262 if (x < 0 && y < 0) return atan(y/x) - pi() 263 if (x == 0 && y > 0) return pi() / 2 264 return -pi() / 2 265 } 266 267 define asin(x) { 268 return atan2(x, sqrt((1+x)*(1-x))) 269 } 270 271 define acos(x) { 272 return atan2(sqrt((1+x)*(1-x)), x) 273 } 274 275 define atanh(x) { 276 return 0.5 * log((1+x)/(1-x)) 277 } 278 279 define asinh(x) { 280 return log(x + sqrt(x^2 + 1)) 281 } 282 283 define acosh(x) { 284 return log(x + sqrt(x^2 - 1)) 285 }