2 files changed, 282 insertions, 0 deletions

diff --git a/README.md b/README.md
new file mode 100644
index 0000000..86f57d8
--- /dev/null
+++ 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
new file mode 100644
index 0000000..6b94a01
--- /dev/null
+++ 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))
+}