#!/bin/sh

# The MIT License (MIT)
#
# Copyright © 2024 pacman64
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the “Software”), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.


# ca [expression]
#
#
# CAlculator is an easier-to-use version of bc (basic calculator) where
#
#   - you give the expression as an argument, while bc reads it from stdin
#   - you can use either ** or ^ to raise powers
#   - you can use [ and ] or ( and ) interchangeably
#   - the number of decimals is 20 by default
#   - automatically includes the extended bc math library via option -l
#   - there are several extra predefined values and functions (see below)
#
#
# Predefined functions and values, besides the ones `bc` already comes with
#
#
# nano   a billionth
# micro  a millionth
# milli  a thousandth
#
# kilo   1000
# mega   a million
# giga   a billion
# tera   a thousand billion / trillion
# peta   a million billion / quadrillion
# exa    a quintillion
# zetta  a sixtillion
#
# binkilo   1024
# binmega   1024^2
# bingiga   1024^3
# bintera   1024^4
# binpeta   1024^5
# binexa    1024^6
# binzetta  1024^7
#
# kb  1024
# mb  1024^2
# gb  1024^3
# tb  1024^4
# pb  1024^5
#
# kib  1024
# mib  1024^2
# gib  1024^3
# tib  1024^4
# pib  1024^5
#
# mol      1 mole, the exact avogadro number
# mole     the exact avogadro number
#
# cup2l    liters in 1 cup
# floz2ml  milliliters in 1 fluid ounce
# gal2l    liters in 1 gallon
# in2cm    centimeters in 1 inch
# lb2kg    kilograms in 1 pound
# mi2km    kilometers in 1 mile
# mpg2kpl  kilometers per liter in 1 mile per gallon
# nmi2km   kilometers in 1 nautical mile
# oz2g     grams in 1 weight ounce
# psi2pa   pascals in 1 pound per square-inch
# ton2kg   kilograms in 1 american ton
# yd2m     meters in 1 yard
#
# ftin(feet, inches)  turn a feet/inches pair into meters
# lboz(feet, inches)  turn a pound/ounces pair into kilograms
#
# eu()          approximate Euler's number
# euler()       approximate Euler's number
# pi()          approximate any circle's circumference/diameter ratio
# tau()         2*pi; approximate any circle's circumference/radius ratio
# abs(x)        the absolute value function
# exp(x)        the exponential function; alias for bc's e
# ln(x)         natural logarithm; an alias for bc's l
# log(x)        natural logarithm; an alias for bc's l
# log2(x)       base-2 logarithm
# log10(x)      base-10 logarithm
# sin(x)        the sine function
# cos(x)        the cosine function
# atan(x)       the arc-tangent function
# hypot(x, y)   length of the hypothenuse
# mod(x, y)     modulus function, which is different from the % operator
# r(x, d)       round to the given number of decimal digits; can be negative
# round(x, d)   round to the given number of decimal digits; can be negative
# r0(x)         round to integer
# r1(x)         round to 1 decimal place
# ...
# r20(x)        round to 20 decimal places
# round0(x)     round to integer
# round1(x)     round to 1 decimal place
# ...
# round20(x)    round to 20 decimal places
# f(x)          the factorial function
# fac(x)        the factorial function
# fact(x)       the factorial function
# factorial(x)  the factorial function
# com(n, k)     number of combinations
# comb(n, k)    number of combinations
# p(n, k)       number of permutations
# per(n, k)     number of permutations
# perm(n, k)    number of permutations
# lcm(x, y)     least-common multiple
# gcd(x, y)     greatest-common divisor
# min(x, y)     the minimum of 2 numbers
# max(x, y)     the maximum of 2 numbers
# sgn(x)        the sign (-1, 0, or +1) of a number
# sign(x)       the sign (-1, 0, or +1) of a number
# isprime(x)    check if number is prime
# lgamma(x)     lanczos approximation of the log-gamma function
# lbeta(x, y)   lanczos approximation of the log-beta function
# gamma(x)      approximation of the gamma function
# beta(x)       approximation of the beta function
# logistic(x)   the logistic function
# mix(x, y, k)  interpolate 2 numbers; extrapolates when k is outside [0, 1]
# clamp(x, min, max)   trap a number in the range given
# wrap(x, min, max)    normalize a number in the range given
# unwrap(x, min, max)  expand a normalized number into the range given


# handle help options
case "$1" in
    -h|--h|-help|--help)
        awk '/^# +ca/, /^$/ { gsub(/^# ?/, ""); print }' "$0"
        exit 0
    ;;
esac

# handle optional scale-override parameter
scale=20
if [ $# -eq 2 ]; then
    scale="$1"
    [ $# -gt 0 ] && shift
fi

# ensure output is all on 1 line, define several funcs and values, then
# inject the expression given as this script's arguments, transformed
# according to the rules described above
{
BC_LINE_LENGTH=0 bc -l << ENDOFSCRIPT
scale = ${scale};

femto = 0.000000000000001;
pico = 0.000000000001;
nano = 0.000000001;
micro = 0.000001;
milli = 0.001;

kilo = 1000;
mega = 1000 * kilo;
giga = 1000 * mega;
tera = 1000 * giga;
peta = 1000 * tera;
exa =  1000 * peta;
zetta =  1000 * exa;

binkilo = 1024;
binmega = 1024 * binkilo;
bingiga = 1024 * binmega;
bintera = 1024 * bingiga;
binpeta = 1024 * bintera;
binexa = 1024 * binpeta;
binzetta = 1024 * binexa;

kb = 1024;
mb = 1024 * kb;
gb = 1024 * mb;
tb = 1024 * gb;
pb = 1024 * tb;
eb = 1024 * pb;
zb = 1024 * eb;

kib = 1024;
mib = 1024 * kib;
gib = 1024 * mib;
tib = 1024 * gib;
pib = 1024 * tib;
zib = 1024 * pib;

mol = 602214076000000000000000;
mole = 602214076000000000000000;

cup2l = 0.23658824;
floz2l = 0.0295735295625;
floz2ml = 29.5735295625;
ft2m = 0.3048;
gal2l = 3.785411784;
in2cm = 2.54;
lb2kg = 0.45359237;
mi2km = 1.609344;
mpg2kpl = 0.425143707;
nmi2km = 1.852;
oz2g = 28.349523125
psi2pa = 6894.757293168;
ton2kg = 907.18474;
yd2m = 0.9144;

ga2l = gal2l;
nm2km = nmi2km;
tn2kg = ton2kg;

hour = 3600;
day = 24 * hour;
week = 7 * day;

hr = hour;
wk = week;

define ftin(f, i) {
    return (0.3048 * f + 0.0254 * i);
}

define lboz(l, o) {
    return (0.45359237 * l + 0.028349523 * o);
}

define eu() {
    return (e(1));
}

define euler() {
    return (e(1));
}

define pi() {
    return (4*a(1));
}

define tau() {
    return (8*a(1));
}

define deg(x) {
    return (180 * x / pi());
}

define rad(x) {
    return (pi() * x / 180);
}

define abs(x) {
    if (x >= 0) return (x);
    return (-x);
}

define exp(x) {
    return (e(x));
}

define ln(x) {
    return (l(x));
}

define log(x) {
    return (l(x));
}

define log2(x) {
    return (l(x) / l(2));
}

define log10(x) {
    return (l(x) / l(10));
}

define sin(x) {
    return (s(x));
}

define cos(x) {
    return (c(x));
}

define tan(x) {
    return (s(x) / c(x));
}

define cot(x) {
    return (c(x) / s(x));
}

define atan(x) {
    return (a(x));
}

define sinh(x) {
    return ((e(x) - e(-x)) / 2);
}

define cosh(x) {
    return ((e(x) + e(-x)) / 2);
}

define tanh(x) {
    return ((e(x) - e(-x)) / (e(x) + e(-x)));
}

define coth(x) {
    return ((e(x) + e(-x)) / (e(x) - e(-x)));
}

define hypot(x, y) {
    return (sqrt(x*x + y*y));
}

define sinc(x) {
    if (x == 0) return (1);
    return (s(x) / x);
}

define min(x, y) {
    if (x <= y) return (x);
    return (y);
}

define max(x, y) {
    if (x >= y) return (x);
    return (y);
}

define mod(x, y) {
    auto s, m;
    s = scale;
    scale = 0;
    m = x % y;
    scale = s;
    return (m);
}

define mod1(x) {
    return (mod(x, 1));
}

define modf(x) {
    return (mod(x, 1));
}

define round0(x) {
    auto i;
    i = x - mod(x, 1);
    if (x - i >= 0.5) {
        return (i + 1);
    }
    return (i);
}

define round(x, d) {
    auto k;
    k = 10^d;
    return (round0(x * k) / k);
}

define r(x, d) {
    return (round(x, d));
}

define r0(x) {
    return (round0(x));
}

define r1(x) {
    return (round(x, 1));
}

define r2(x) {
    return (round(x, 2));
}

define r3(x) {
    return (round(x, 3));
}

define r4(x) {
    return (round(x, 4));
}

define r5(x) {
    return (round(x, 5));
}

define r6(x) {
    return (round(x, 6));
}

define r7(x) {
    return (round(x, 7));
}

define r8(x) {
    return (round(x, 8));
}

define r9(x) {
    return (round(x, 9));
}

define r10(x) {
    return (round(x, 10));
}

define r11(x) {
    return (round(x, 11));
}

define r12(x) {
    return (round(x, 12));
}

define r13(x) {
    return (round(x, 13));
}

define r14(x) {
    return (round(x, 14));
}

define r15(x) {
    return (round(x, 15));
}

define r16(x) {
    return (round(x, 16));
}

define r17(x) {
    return (round(x, 17));
}

define r18(x) {
    return (round(x, 18));
}

define r19(x) {
    return (round(x, 19));
}

define r20(x) {
    return (round(x, 20));
}

define round1(x) {
    return (round(x, 1));
}

define round2(x) {
    return (round(x, 2));
}

define round3(x) {
    return (round(x, 3));
}

define round4(x) {
    return (round(x, 4));
}

define round5(x) {
    return (round(x, 5));
}

define round6(x) {
    return (round(x, 6));
}

define round7(x) {
    return (round(x, 7));
}

define round8(x) {
    return (round(x, 8));
}

define round9(x) {
    return (round(x, 9));
}

define round10(x) {
    return (round(x, 10));
}

define round11(x) {
    return (round(x, 11));
}

define round12(x) {
    return (round(x, 12));
}

define round13(x) {
    return (round(x, 13));
}

define round14(x) {
    return (round(x, 14));
}

define round15(x) {
    return (round(x, 15));
}

define round16(x) {
    return (round(x, 16));
}

define round17(x) {
    return (round(x, 17));
}

define round18(x) {
    return (round(x, 18));
}

define round19(x) {
    return (round(x, 19));
}

define round20(x) {
    return (round(x, 20));
}

define fac(x) {
    auto f, i;
    if (x < 0) return (0);
    f = 1;
    for (i = x; i >= 2; i--) { f *= i; }
    return (f);
}

define f(x) {
    return (fac(x));
}

define fact(x) {
    return (fac(x));
}

define factorial(x) {
    return (fac(x));
}

define per(n, k) {
    auto p, i;
    if (n < k) return (0);
    p = 1;
    for (i = n; i >= n - k + 1; i--) { p *= i; }
    return (p);
}

define p(n, k) {
    return (per(n, k));
}

define perm(n, k) {
    return (per(n, k));
}

define com(n, k) {
    if (n < k) return (0);
    return (per(n, k) / fac(k));
}

define comb(n, k) {
    return (com(n, k));
}

define gcd(x, y) {
    return (x * y / lcm(x, y));
}

define lcm(x, y) {
    auto a, b, z;

    /* the LCM is defined only for positive integers */
    # if (mod(x, 1) != 0 || x < 1 || mod(y, 1) != 0 || y < 1) { return 0; }
    if (mod(x, 1) != 0) return (0);
    if (x < 1) return (0);
    if (mod(y, 1) != 0) return (0);
    if (y < 1) return (0);

    a = min(x, y);
    b = max(x, y);

    z = b;
    while (mod(z, a) != 0) { z += b; }
    return (z);
}

define sgn(x) {
    if (x > 0) return (1);
    if (x < 0) return (-1);
    return (0);
}

define sign(x) {
    if (x > 0) return (1);
    if (x < 0) return (-1);
    return (0);
}

define isprime(n) {
    auto div;

    # if (mod(n, 1) != 0 || n < 2) { return 0; }
    if (mod(n, 1) != 0) return (0);
    if (n < 2) return (0);

    /* 2 is the only even prime number */
    if (mod(n, 2) == 0) {
        if (n == 2) return (1);
        return (0);
    }

    /*
        even numbers have already been handled, so only odd numbers
        make it here: only check up to the square-root, so the loop
        has O(n**1.5) time-complexity
    */
    for (div = 3; div <= sqrt(n); div += 2) {
        if (mod(n, div) == 0) return (0);
    }
    return (1);
}

/*
I adapted an implementation of Lanczos' log-gamma approximation I found at
http://introcs.cs.princeton.edu/java/91float/Gamma.java.html
*/
define lgamma(x) {
    auto tmp, a, b, c, d;
    tmp = (x - 0.5) * l(x + 4.5) - (x + 4.5);
    a = 76.18009173 / (x + 0) - 86.50532033 / (x + 1);
    b = 24.01409822 / (x + 2) - 1.231739516 / (x + 3);
    c = 0.00120858003 / (x + 4) - 0.00000536382 / (x + 5);
    d = sqrt(2 * 3.14159265358979323844);
    return (tmp + l((1.0 + a + b + c) * d));
}

define lbeta(x, y) {
    return (lgamma(x) + lgamma(y) - lgamma(x + y));
}

define gamma(x) {
    return (e(lgamma(x)));
}

define beta(x, y) {
    return (e(lbeta(x, y)));
}

define logistic(x) {
    return (1 / (1 + e(-x)));
}

define mix(x, y, k) {
    # return (x + (1 - k) * (y - x));
    return (x + k * y);
}

define clamp(x, min, max) {
    return (max(min(max, x), min));
}

define wrap(x, min, max) {
    return ((x - min) / (max - min));
}

define unwrap(x, min, max) {
    return ((max - min) * x + min);
}

$(echo "$@" | sed 's-^+--g; s-_--g; s-\*\*-^-g; s-\[-(-g; s-\]-)-g')
ENDOFSCRIPT
} |
# ensure the result shows at least a zero before the decimal dot, then
# rid the result of trailing zero decimals and/or trailing decimal dots
sed -E 's-^\.-0.-; s/^-\./-0./; s-(\.[0-9]+[1-9]+)0+$-\1-; s-\.0*$--'