functions.bc

### functions.bc

################################################################################
### Helpful Constants and Functions

pi=4*a(1)
ex=e(1)
define sgn(x) { if(x>0) return 1; if(x<0) return -1; }
define abs(x) { return sgn(x)*x }
define heavyside(x) { return (x>0) }

# Return the integer-part of a number (not the floor)
define int(x) { auto s; s=scale; scale=0; x/=1; scale=s; return x; }

################################################################################
### Logarithms & Exponentials

define ln(x) { return l(x) }
define log(x) { return l(x)/l(10) }
define logb(x,b) { return l(x)/l(b) }
define pow(x,n) { return e(n*l(x)) }

################################################################################
### Trigonometry

define rad2deg(x) { return x*(45/a(1)) }
define deg2rad(x) { return x*(a(1)/45) }

define cos(x) { return c(x) }
define sin(x) { return s(x) }
define tan(x) { return s(x)/c(x) }
define sec(x) { return 1/c(x) }
define csc(x) { return 1/s(x) }
define cot(x) { return c(x)/s(x) }
define arccos(x) {
    pi = 4*a(1)
    if(x == 1) return 0
    if(x == -1) return pi() 
    return pi/2-a(x/sqrt(1-(x^2)))
}
define arcsin(x) { 
    pi = 4*a(1)
    if(x == 1) return pi/2
    if(x == -1) return -pi/2
    return sgn(x)*a(sqrt((1/(1-(x^2)))-1))
}
define arctan(x) { return a(x) }
define atan2(y,x) { 
    pi = 4*a(1)
    if(x == 0) return sgn(y)*pi/2
    if(x < 0) {
        if(y < 0)
            return a(y/x)-pi
        return a(y/x)+pi
    }
    return a(y/x)
}
define arcsec(x) { return arccos(1/x) }
define arccsc(x) { return arcsin(1/x) }
define arccot(x) { return (4*a(1))/2-a(x) }

define cosh(x) { return (e(x)+e(-x))/2 }
define sinh(x) { return (e(x)-e(-x))/2 }
define tanh(x) { return sinh(x)/cosh(x) }
define sech(x) { return 1/cosh(x) }
define csch(x) { return 1/sinh(x) }
define coth(x) { return 1/tanh(x) }

define arcosh(x) { return ln(x+sqrt(x^2-1)) }
define arsinh(x) { return ln(x+sqrt(x^2+1)) }
define artanh(x) { return (1/2)*ln((1+x)/(1-x))}
define arsech(x) { return arcosh(1/x) }
define arcsch(x) { return arsinh(1/x) }
define arcoth(x) { return artanh(1/x) }

################################################################################
### Combinatorics

define factorial(n) {
    auto i
    if ( (n < 0) || (int(n) != n) ) {
        print "Error: factorials defined for positive integers only\n"
        return
    }
    factorial[0] = factorial[1] = 1
    for (i=n; factorial[i] == 0; --i) {}
    for (++i; i<=n; ++i)
        factorial[i] = i*factorial[i-1]
    return factorial[n]
}

# Return nPk, the number of ways to permute k of n objects
define pick(n,k) {
    auto i, r
    if ( (n < 0) || (k < 0) || (int(n) != n) || (int(k) != k) ) {
        print "Error: permutations defined for positive integers only\n"
        return
    }
    r = n-k
    for (i=1; n > r; --n)
        i*=n
    return i
}

# Return nCk, the number of ways to choose k of n objects
define choose(n,k) {
    if ( (n < 0) || (k < 0) || (int(n) != n) || (int(k) != k) ) {
        print "Error: combinations defined for positive integers only\n"
        return
    }
    if ( n*(n+1)/2 + k >= 2^24 ) {
        return int(factorial(n)/(factorial(n-k)*factorial(k)))
    }
    return choose_(n,k)
}

# Dynamic, recursive helper function for choose(n,k)
define choose_(n,k) {
    auto i
    if ( (k == 0) || (k == n) )
        return 1
    i = n*(n+1)/2 + k
    if ( choose[i] != 0)
        return choose[i]
    choose[i] = choose_(n-1,k) + choose_(n-1,k-1)
    return choose[i]
}

# Return the nth Fibonacci number
define fibonacci(n) {
    auto i
    if ( (n < 1) || (int(n) != n) ) {
        print "Error: Fibonacci numbers indexed by positive integers\n"
        return
    }
    fibonacci[1] = fibonacci[2] = 1
    for (i=n; fibonacci[i] == 0; --i) {}
    for (++i; i<=n; ++i)
        fibonacci[i] = fibonacci[i-1] + fibonacci[i-2]
    return fibonacci[n]
}

################################################################################
### Calculus

# Recursively applies Newton's Method with initial parameter x
# and (global) functions f and ff=f' until the process converges up to scale
define newton(x) { 
    print x, "\n"
    if (f(x)/9)
        return newton(x - f(x)/ff(x))
    return x
}

# Numerically calculate the definite integral 
# of f (global) from a to b 
define integrate(a,b) {
    return boole_(a,b,(b-a)*(e((scale/5)*l(10))))
}

# Composite Boole's Rule
# The Definite Integral of f from a to b using n subdivisions
define boole_(a,b,n) {
    auto h, i, left, mid, right, sum
    h = (b-a)/(4*n)
    right=f(a)
    for (i=a+4*h; i<=b; i+=4*h) {
        left = right
        mid[1] = f(i-3*h)
        mid[2] = f(i-2*h)
        mid[3] = f(i-h)
        right = f(i)
        sum += 7*left+32*mid[1]+12*mid[2]+32*mid[3]+7*right
    }
    return 2*sum*h/45
}

### Trapezoid Rule 
# The Definite Integral of f from a to b using n subdivisions
define trapezoid_(a,b,n) {
    auto h, i, sum
    h = (b-a)/n
    for (i=1; i<n; i++) {
        sum += f(a+i*h)
    }
    return (h/2)*( f(a) + 2*sum + f(b) )
}

################################################################################
### Number Theory

# Returns the nth prime integer
define prime(n) {
    auto i, s 
    prime[1]=2
    prime[2]=3
    for(i=n; prime[i]==0; --i) {}
    return prime_(i,n)
}

# Fill in the array of prime numbers prime[]
# beyond index j, up through k
define prime_(j,k) {
    auto i, n, flag, p, s
    s = scale
    scale = 0
    for(i=j+1; i<=k; ++i) {
        flag=0;
        for(n=prime[i-1]+2; !flag; n=n+2) {
            flag=1
            for(p=1; p<i; ++p) {
                if (n % prime[p] == 0) {flag=0; break;}
            }
            if(flag) {prime[i]=n; break;}
        }
    }
    scale = s
    return prime[k]
}