commit 894de25beb628716e1b7abe4fc35ca0255289160
Author: Sebastiano Tronto <sebastiano@tronto.net>
Date: Sun, 19 Mar 2023 22:34:16 +0100
Initial commit
Diffstat:
| A | README.md | | | 5 | +++++ |
| A | bc.library | | | 277 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
2 files changed, 282 insertions(+), 0 deletions(-)
diff --git a/README.md b/README.md
@@ -0,0 +1,5 @@
+This is a small collection of mathematical functions for bc(1),
+the standard UNIX basic calculator.
+
+To use them, copy `bc.library` in a convenient location and call
+bc with `bc /path/to/bc.library`.
diff --git a/bc.library b/bc.library
@@ -0,0 +1,277 @@
+/*
+ * The content of this file is free knowledge.
+ * No copyright applies. No warranty is provided.
+ *
+ * The following functions are available (aliases in parentheses):
+ * General utility: max, min, sgn, abs, floor, ceiling
+ * Arithmetic: factorial(fact), binomial(binom, bin), gcd, lcm, totient(phi)
+ * Calculus: exp, log(ln), pow, sin, cos, tan, cosh, sinh, tanh
+ * atan, atan2, asin, acos, atanh, asinh, acosh
+ * Constants (as functions): e(), pi()
+ *
+ * Approximation of inverse trigonometric functions and of pi is not good.
+ *
+ * For functions returning non-integer values, remember to set the
+ * desired scaled value before calling the function.
+ */
+
+/* General utility functions */
+
+define max(x, y) {
+ if (x >= y) return x
+ return y
+}
+
+define min(x, y) {
+ if (x <= y) return x
+ return y
+}
+
+define sgn(c) {
+ if (x > 0) return 1
+ if (x < 0) return -1
+ return 0
+}
+
+define abs(x) {
+ return x * sgn(x)
+}
+
+define floor(x) {
+ auto s, y
+ if (x < 0) return -ceiling(-x)
+ s = scale
+ scale = 0
+ y = x / 1
+ scale = s
+ return y
+}
+
+define ceiling(x) {
+ if (x < 0) return -floor(-x)
+ if (floor(x) == x) return x
+ return 1 + floor(x)
+}
+
+/* Arithmetic */
+
+define factorial(n) {
+ auto i, res, s
+ s = scale
+ scale = 0
+ res = 1
+ for (i = 1; i <= n; i++) res *= i
+ scale = s
+ return res
+}
+
+define fact(n) {
+ return factorial(n)
+}
+
+define binomial(n, k) {
+ auto i, j, told[], tnew[], s
+ if (k < 0) return -1
+ if (k > n) return 0
+ s = scale
+ scale = 0
+ tnew[0] = 1
+ for (i = 1; i <= n; i++) {
+ for (j = 0; j <= i; j++) told[j] = tnew[j]
+ for (j = 1; j <= i; j++) tnew[j] = told[j] + told[j-1]
+ }
+ scale = s
+ return tnew[k]
+}
+
+define binom(n, k) {
+ return binomial(n, k)
+}
+
+define bin(n, k) {
+ return binomial(n, k)
+}
+
+define gcd(a, b) {
+ auto s, aux
+ s = scale
+ scale = 0
+ a /= 1
+ b /= 1
+ while (b != 0) {
+ aux = a
+ a = b
+ b = aux % b
+ }
+ scale = s
+ return a
+}
+
+define lcm(a, b) {
+ if (a == 0 || b == 0) return 0
+ return a * b / gcd(a, b)
+}
+
+define phi(n) {
+ auto i, j, f[], r, s
+ s = scale
+ scale = 0
+ j = 0
+ r = n
+ for (i = 2; i*i <= r; i++) {
+ if (n % i == 0) {
+ f[j] = i
+ j += 1
+ while (n % i == 0) n /= i
+ }
+ }
+ for (i = 0; i < j; i++) r -= r / f[i]
+ scale = s
+ return r
+}
+
+define totient(n) {
+ return phi(n)
+}
+
+/* Calculus */
+
+define exp(x) {
+ auto i, n, d, series
+ scale += 10
+ i = 0
+ n = 1
+ d = 1
+ series = 0
+ while (n/d != 0) {
+ series += n/d
+ i += 1
+ n *= x
+ d *= i
+ }
+ scale -= 10
+ return series / 1
+}
+
+define e() {
+ return exp(1)
+}
+
+/* Log: first take some square roots to make x smaller, then
+ use Newton's method to solve e^y-x=0 */
+define log(x) {
+ auto m, s, t, n, d
+ scale += 10
+ while (x > 3) {
+ m += 1
+ x = sqrt(x)
+ }
+ t = x
+ d = exp(t)
+ n = d - x
+ while (n/d != 0) {
+ t -= n/d
+ d = exp(t)
+ n = d - x
+ }
+ scale -= 10
+ return (t * 2^m) / 1
+}
+
+define ln(x) {
+ return log(x)
+}
+
+define pow(a, b) {
+ if (b != floor(b)) return (2^floor(b)) * exp(log(a)*(b-floor(b)))
+ return a ^ b
+}
+
+define trig(x, ii, in, id) {
+ auto i, n, d, series
+ scale += 10
+ i = ii
+ n = in
+ d = id
+ while (n/d != 0) {
+ series += n/d
+ i += 2
+ n *= -x*x
+ d *= i*(i-1)
+ }
+ scale -= 10
+ return series / 1
+}
+
+define sin(x) {
+ return trig(x, 1, x, 1)
+}
+
+define cos(x) {
+ return trig(x, 0, 1, 1)
+}
+
+define tan(x) {
+ return sin(x) / cos(x)
+}
+
+define sinh(x) {
+ return 0.5*(exp(x) - exp(-x))
+}
+
+define cosh(x) {
+ return 0.5*(exp(x) + exp(-x))
+}
+
+define tanh(x) {
+ return sinh(x) / cosh(x)
+}
+
+/* Atan: approximate integral of 1/(1+x^2) */
+define atan(x) {
+ auto i, n, f1, f2, a, res, s
+ s = scale
+ scale = 10
+ n = 10^4
+ res = 0
+ for (i = 0; i < n; i++) {
+ f1 = 1 + (x*i/n)^2
+ f2 = 1 + (x*(i+1)/n)^2
+ a = (f1 + f2) / (2 * f1 * f2)
+ res += a * x / n
+ }
+ scale = s
+ return res/1
+}
+
+define pi() {
+ return 4 * atan(1)
+}
+
+define atan2(y, x) {
+ if (x > 0) return atan(y/x)
+ if (x < 0 && y >= 0) return atan(y/x) + pi()
+ if (x < 0 && y < 0) return atan(y/x) - pi()
+ if (x == 0 && y > 0) return pi() / 2
+ return -pi() / 2
+}
+
+define asin(x) {
+ return atan2(x, sqrt((1+x)*(1-x)))
+}
+
+define acos(x) {
+ return atan2(sqrt((1+x)*(1-x)), x)
+}
+
+define atanh(x) {
+ return 0.5 * log((1+x)/(1-x))
+}
+
+define asinh(x) {
+ return log(x + sqrt(x^2 + 1))
+}
+
+define acosh(x) {
+ return log(x + sqrt(x^2 - 1))
+}
