routines.bc

### routines.bc

# Print the Pythagorean triple generated 
# by parameters m and n, and return the hypotenuse.
define pythagtriple(m,n) {
    print abs(n*n-m*m), "\n"
    print abs(2*m*n), "\n"
    return m*m+n*n
}

# Print the Pythagorean quadruple generated 
# by parameters m, n, p, and q, and return the hypotenuse.
define pythagquadruple(m,n,p,q) {
    print abs(n*n+m*m-p*p-q*q), "\n"
    print abs(2*(m*q+n*p)), "\n"
    print abs(2*(n*q-m*p)), "\n"
    return m*m+n*n+p*p+q*q
}

# Prints the angle x, presumed as a decimal number of degrees,
# in DMS (degrees, minutes, seconds) notation.
define void dd2dms(x) { 
    d = int(x)
    m = abs(int(frac(x)*60))
    s = abs(frac(frac(x)*60)*60)
    print d,"°",m,"′",s,"″\n"
}

# Iteratively applies Newton's Method, printing each result,
# with initial guess x and (global) function f.
define void newtoniter(x) { 
    print x/1, "\n"
    if (f(x)/9)
        newtoniter(x - f(x)/derivative(x))
}

# Print the roots and vertex coordinates of ax²+bx+c.
define void quadratic(a,b,c) {
    print "For vertex (h,k) and roots r[1] and r[2],\n"

    h = -b/(2*a)
    k = a*h^2 + b*h + c
    print "h = ", h, "\n"
    print "k = ", k, "\n"

    if(b^2-4*a*c < 0) {
        print "The roots are non-real: \n"
        scale /= 2
        print h/1, " ± i", sqrt(k/a), "\n"
        scale *= 2
    } else {
        r[1] = h - sqrt(-k/a)
        r[2] = h + sqrt(-k/a)
        print "r[1] = ", r[1], "\n"
        print "r[2] = ", r[2], "\n"
    }
}

# # Print the roots and local extrema coordinates of ax³+bx²+cx+d
# define void cubic(a,b,c,d) {
#     # 3a*x^2 + b*x + c
#     # 6a*x + b 
#     print "For inflection point at i,\n"
#     print "local extrema at h[1] and h[2],\n"
#     print "and roots r[1] and r[2] and r[3],\n"
# }
# 
# # Print the roots and local extrema coordinates of ax⁴+bx³+cx²+dx+e
# define void quartic(a,b,c,d,e) {
#     #print "For TK jerk points at ,\n"
#     print "For inflection points at i[1] and i[2],\n"
#     print "local extrema at h[1] and h[2] and h[3],\n"
#     print "and roots r[1] and r[2] and r[3] and r[4],\n"
# }

# Print the coefficients of the simple continued fraction representation of x
# as well as each convergent obtained by truncating at that coefficient,
# which are successively better rational approximations to the number x
define void contfrac(x) {
    auto i, j, k, l, a[], h[], k[], r

    if(x == int(x)) { print x, " = ", x, "/", 1, "\n"; return }

    r = x
    a[0] = int(r)
    h[0] = a[0]
    k[0] = 1
    l = length(a[0])

    r = 1/(r-a[0])
    a[1] = int(r)
    h[1] = a[1]*h[0]+1
    k[1] = a[1]
    l = max(l,length(a[1]))

    for(i=2; x!=h[i-1]/k[i-1]; ++i) {
        r = 1/(r-a[i-1])
        a[i] = int(r)
        h[i] = a[i]*h[i-1]+h[i-2]
        k[i] = a[i]*k[i-1]+k[i-2]
        l = max(l,length(a[i]))
    }
    
    for(j=0; j<i; ++j) {
        for(k=l-length(a[j]); k; k--) { print " " }
        print a[j], " | ", h[j]/k[j], " = ", h[j], "/", k[j], "\n"
    }
}

# Print n expressed in bases 2, 3, ..., 36.
define void bases(n) {
    auto i, base
    base = obase
    for(i=2; i<=36; i++) {
        obase = A
        if (i<A) print " "
        print "  ", i, " |"
        if (i<H) print " "
        obase = i
        print n, "\n"
    }
    obase = base
}

# Print the prime integer factorization of n.
define void factor(n) {
    auto i, s, p
    s = scale
    scale = 0
    for(i=1; n!=1; i++) {
        p=prime(i)
        while(n % p == 0) {
            print p, " "
            n/=p 
        }
    } 
    print "✓\n"
    scale = s
    return
}

# Convert two-dimensional rectangular coordinates 
# to polar coordinates and vice-versa respectively.
define void rect2polar(x,y) {
    r = sqrt(x^2+y^2)
    theta = atan2(y,x)
    print "r = ", r, "\n"
    print "theta = ", theta, "\n"
    return
}
define void polar2rect(r,theta) {
    x = r*c(theta)
    y = r*s(theta)
    print "x = ", x, "\n"
    print "y = ", y, "\n"
    return
}

# Convert three-dimensional rectangular coordinates 
# to cylindrical coordinates and vice-versa respectively.
define void rect2cyl(x,y,zed) {
    rect2polar(x,y)
    z = zed
    print "z = ", z, "\n"
    return
}
define void cyl2rect(r,theta,zed) {
    polar2rect(r,theta)
    z = zed
    print "z = ", z, "\n"
    return
}

# Convert three-dimensional rectangular coordinates 
# to spherical coordinates and vice-versa respectively,
# following the mathematician’s convention 
# with the zenith angle θ in the range [0,2π) with θ = 0 along the positive x-axis,
# and the azimuth angle φ in the range [0,π] with φ = 0 along the positive z-axis.
define void rect2sphere(x,y,z) {
    rho = sqrt(x^2+y^2+z^2)
    theta = atan2(y,x)
    phi = arccos(z/rho)
    print "rho = ", rho, "\n"
    print "theta = ", theta, "\n"
    print "phi = ", phi, "\n"
    return
}
define void sphere2rect(rho,theta,phi) {
    x = rho*s(phi)*c(theta)
    y = rho*s(phi)*s(theta)
    z = rho*c(phi)
    print "x = ", x, "\n"
    print "y = ", y, "\n"
    print "z = ", z, "\n"
    return
}

# Print the sequence of positive integers that results from 
# iteratively applying the function prescribed by the Collatz Conjecture.
define void collatz(n) { 
    auto s
    s = scale
    scale = 0
    collatz_(n)
    scale = s
}
define collatz_(n) {
    if(n==1)
        return 1
    print n, " → " 
    if(n%2) {
        return collatz_(3*n+1)
    } else {
        return collatz_(n/2)
    }
}

# Print every way that m can be written as a sum 
# of powers of consecutive positive integers.
define void sumofpowers(m) {
    auto s, n, i, j, cap, powers[], any
    s = scale
    for(n=1; 2^n<m; n++) {
        cap = pow(m,1/n)+1
        scale = 0
        for(i=1; i<cap; i++)
            powers[i] = i^n
        i = j = sum = any = 1
        while(j < cap) {
            if(sum >= m) {
                if(sum == m) {
                    if(any) {print "For n = ", n, "\n"; any = 0; }
                    print "  ", m, " = ", i, "ⁿ + ⋯  + ", j, "ⁿ\n"
                } 
                sum -= powers[i++]
            } else {
                sum += powers[++j]
            }
        }
        scale = s
    }
    print "✓\n"
}

# Prints the digits of the integer part of x expressed in base obase.
define void intdigits(x) {
    auto i
    i = intdigits_(x);
    for(i=1; obase-intdigits[i]; i++)
        print "intdigits[", i, "] = ", intdigits[i], "\n";
    print "fracdigits[", i, "] = ", fracdigits[i], " (obase) \n";
}

# Creates an array intdigits[] containing the digits 
# of the integer part of x expressed in base obase
# for which intdigits[i] is the digit in the ith position left of the radix point
# and which terminates with a value of obase.
define void intdigits_(x) {
    auto s, i
    s = scale;
    scale = 0;
    x/=1;
    for(i=1; x; x=(x-intdigits[i++])/obase)
        intdigits[i] = x % obase;
    intdigits[i] = obase;
    scale = s;
}

# Prints the digits of the fractional part of x expressed in base obase.
define void fracdigits(x) {
    auto i
    i = fracdigits_(x);
    for(i=1; obase-fracdigits[i]; i++)
        print "fracdigits[", i, "] = ", fracdigits[i], "\n";
    print "fracdigits[", i, "] = ", fracdigits[i], " (obase) \n";
}

# Creates an array fracdigits[] containing the digits 
# of the fractional part of x expressed in base obase
# for which fracdigits[i] is the digit in the ith position right of the radix point
# and which terminates with a value of obase.
define void fracdigits_(x) {
    auto i, d
    d = int(1+scale*l(A)/l(obase))
    for(i=1; i<=d; x -= fracdigits[i++]) {
        x = (frac(x))*obase
        fracdigits[i] = int(x);
    }
    fracdigits[i] = obase
}

# Returns 1 if the digits of x expressed in obase 
# are only those digits allowed in *digits[],
# indicated with a value at that digits's index in *digits[] of 1.
define onlydigits(x, *digits[]) {
    auto i, j, flag
    intdigits_(x)
    fracdigits_(x)
    for(j=1; obase-intdigits[j]; j++) {
        flag = 0
        for(i=0; i<obase; i++) {
            if(digits[i] && intdigits[j]==i) {
                flag = 1
                break
            }
        }
        if(!flag) { return 0 }
    }
    for(j=1; obase-fracdigits[j]; j++) {
        flag = 0
        for(i=0; i<obase; i++) {
            if(digits[i] && fracdigits[j]==i) {
                flag = 1
                break
            }
        }
        if(!flag) { return 0 }
    }
    return 1
}

# Prints the unique sum of non-adjacent Fibonacci numbers equal to n.
define void zeckendorf(n) {
    auto f
    f = fibonacci(int(l(n*sqrt(5)+1/2)/l((1+sqrt(5))/2)))
    print f
    if(n-f) {
        print " + "
        zeckendorf(n-f)
    } else {
        print " ✓\n"
    }
}