#!/usr/local/bin/bc -l ### Funcs.BC - a large number of functions for use with GNU BC ## Not to be regarded as suitable for any purpose ## Not guaranteed to return correct answers scale=50; define pi() { auto s; if(scale==(s=scale(pi_)))return pi_ if(scale 1, 0.99 -> 0 define int(x) { auto os;os=scale;scale=0;x/=1;scale=os;return(x) } # Round down to integer below x define floor(x) { auto os,xx;os=scale;scale=0 xx=x/1;if(xx>x).=xx-- scale=os;return(xx) } # Round up to integer above x define ceil(x) { auto os,xx;x=-x;os=scale;scale=0 xx=x/1;if(xx>x).=xx-- scale=os;return(-xx) } # Fractional part of x: 12.345 -> 0.345 define frac(x) { auto os,xx;os=scale;scale=0 xx=x/1;if(xx>x).=xx-- scale=os;return(x-xx) } # Absolute value of x define abs(x) { if(x<0)return(-x)else return(x) } # Sign of x define sgn(x) { if(x<0)return(-1)else if(x>0)return(1);return(0) } # Round x up to next multiple of y define round_up( x,y) { return(y*ceil( x/y )) } # Round x down to previous multiple of y define round_down(x,y) { return(y*floor(x/y )) } # Round x to the nearest multiple of y define round( x,y) { auto os,oib; os=scale;oib=ibase .=scale++;ibase=A y*=floor(x/y+.5) ibase=oib;scale=os return y } # Find the remainder of x/y define int_remainder(x,y) { auto os; os=scale;scale=0 x/=1;y/=1;x%=y scale=os return(x) } define remainder(x,y) { os=scale;scale=0 if(x==x/1&&y==y/1){scale=os;return int_remainder(x,y)} scale=os return(x-round_down(x,y)) } # Greatest common divisor of x and y define int_gcd(x,y) { auto r,os; os=scale;scale=0 x/=1;y/=1 while(y>0){r=x%y;x=y;y=r} scale=os return(x) } define gcd(x,y) { auto r,os; os=scale;scale=0 if(x==x/1&&y==y/1){scale=os;return int_gcd(x,y)} scale=os while(y>0){r=remainder(x,y);x=y;y=r} return(x) } # Lowest common multiple of x and y define int_lcm(x,y) { auto r,m,os; os=scale;scale=0 x/=1;y/=1 m=x*y while(y>0){r=x%y;x=y;y=r} m/=x scale=os return(m) } define lcm(x,y) { return (x*y/gcd(x,y)) } # Remove largest possible power of 2 from x define oddpart(x){ auto os; os=scale;scale=0;x/=1 if(x==0){scale=os;return 1} while(!x%2)x/=2 scale=os;return x } # Largest power of 2 in x define evenpart(x) { auto os; os=scale;scale=0 x/=oddpart(x/1) scale=os;return x } ## Trig / Hyperbolic Trig # Sine define sin(x) { return s(x) } # alias for standard library # Cosine define c(x) { return s(x+pi()/2) } # as fast or faster than define cos(x) { return c(x) } # . standard library # Tangent define tan(x) { auto c;c=c(x);if(c==0)c=A^-scale;return(s(x)/c) } # Secant define sec(x) { auto c;c=c(x);if(c==0)c=A^-scale;return( 1/c) } # Cosecant define cosec(x) { auto s;s=s(x);if(s==0)s=A^-scale;return( 1/s) } # Cotangent define cotan(x) { auto s;s=s(x);if(s==0)s=A^-scale;return(c(x)/s) } # Arcsine define arcsin(x) { if(x==-1||x==1)return(pi()/2*x);return( a(x/sqrt(1-x*x)) ) } # Arccosine define arccos(x) { if(x==0)return(0);return pi()/2-arcsin(x) } # Arctangent (one argument) define arctan(x) { return a(x) } # alias for standard library # Arctangent (two arguments) define arctan2(x,y) { auto p; if(x==0&&y==0)return(0) p=(1-sgn(y))*pi()*(2*(x>=0)-1)/2 if(x==0||y==0)return(p) return(p+a(x/y)) } # Arcsecant define arcsec(x) { return( a(x/sqrt(x*x-1)) ) } # Arccosecant define arccosec(x) { return( a(x/sqrt(x*x-1))+pi()*(sgn(x)-1)/2 ) } # Arccotangent (one argument) define arccotan(x) { return( a(x)+pi()/2 ) } # Arccotangent (two arguments) define arccotan2(x,y) { return( arctan(x,y)+pi()/2 ) } # Hyperbolic Sine define sinh(x) { auto t;t=e(x);return((t-1/t)/2) } # Hyperbolic Cosine define cosh(x) { auto t;t=e(x);return((t+1/t)/2) } # Hyperbolic Tangent define tanh(x) { auto t;t=e(x+x)-1;return(t/(t+2)) } # Hyperbolic Secant define sech(x) { auto t;t=e(x);return(2/(t+1/t)) } # Hyperbolic Cosecant define cosech(x) { auto t;t=e(x);return(2/(t-1/t)) } # Hyperbolic Cotangent define coth(x) { auto t;t=e(x+x)-1;return((t+2)/t) } # Hyperbolic Arcsine define arcsinh(x) { return( l(x+sqrt(x*x+1)) ) } # Hyperbolic Arccosine define arccosh(x) { return( l(x+sqrt(x*x-1)) ) } # Hyperbolic Arctangent define arctanh(x) { return( l((1+x)/(1-x))/2 ) } # Hyperbolic Arcsecant define arcsech(x) { return( l((sqrt(1-x*x)+1)/x) ) } # Hyperbolic Arccosecant define arccosech(x) { return( l((sqrt(1+x*x)*sgn(x)+1)/x) ) } # Hyperbolic Arccotangent define arccoth(x) { return( l((x+1)/(x-1))/2 ) } # Length of the diagonal vector (0,0)-(x,y) [pythagoras] define pyth(x,y) { return(sqrt(x*x+y*y)) } define pyth3(x,y,z) { return(sqrt(x*x+y*y+z*z)) } # Gudermannian Function define gudermann(x) { return 2*(a(e(x))-a(1)) } # Inverse Gudermannian Function define arcgudermann(x) { return arctanh(s(x)) } # Bessel function define besselj(n,x) { return j(n,x) } # alias for standard library ## Exponential / Logs # Exponential e^x define exp(x) { return e(x) } # alias for standard library # Natural Logarithm (base e) define ln(x) { auto os,len,ln; if(x< 0){print "ln error: logarithm of a negative number\n";return 0} if(x==0)print "ln error: logarithm of zero; negative infinity\n" len=length(x)-scale(x)-1 if(lenlim||(maxmax2||3*length(t)>es+es)&&i>1){cf[i--]=0;break}#cheat: assume rational if(y==0)break;#completely rational } n=1;d=cf[i] if(i==0){print "id_frac2_: something is wrong; y=",y,", d=",d,"\n"} if(d!=-2&&i)while(--i){d=n+cf[i]*(t=d);n=t} if(d2){scale=os;return x^y} scale=os;x/=1;scale=0 x=fastintpow__(x,y); scale=os;return x; } # Raise x to a fractional power faster than e^(y*l(x)) define fastfracpow_(x,y) { auto f,yy,inv; inv=0;if(y<0){y=-y;inv=1} y-=int(y) if(y==0)return 1; if((yy=y*2^C)!=int(yy)){x=l(x);if(inv)x=-x;return e(y/1*x)} # faster using square roots for rational binary fractions # where denominator <= 8192 x=sqrt(x) for(f=1;y&&x!=1;x=sqrt(x))if(y+=y>=1){.=y--;f*=x} if(inv)f=1/f; return f; } # Find the yth root of x where y is integer define fastintroot_(x,y) { auto os,d,r,ys,eps; os=scale;scale=0;y/=1;scale=os if(y<0){x=1/x;y=-y} if(y==1){return x} if(y>=x-1){return fastfracpow_(x,1/y)} if(y*int((d=2^F)/y)==d){ r=1;while(r+=r<=y)x=sqrt(x) return x } scale=length(y)-scale(y);if(scale<5)scale=5;r=e(ln(x)/y) scale=os+5;if(scale<5)scale=5 d=1;eps=A^(3-scale) ys=y-1 while(d>eps){ d=r;r=(ys*r+x/fastintpow_(r,ys))/y d-=r;if(d<0)d=-d } scale=os return r/1 } # Raise x to the y-th power define pow(x,y) { auto os,p,ix,iy,fy,dn,s; if(y==0) return 1 if(x==0) return 0 if(00) print "even root" if(dn<0) print "irrational power" print " of a negative number\n" scale=os;return 0 } if(y==iy) { if(x==ix){p=fastintpow_(ix,iy);if(iy>0){scale=0;p/=1};scale=os;return p/1} scale+=scale;p=fastintpow_(x,iy);scale=os;return p/1 } if((dn=id_frac2_(y))!=-2){ #accurate rational roots (sometimes slower) if(dn<0)dn=-dn s=1;if(y<0){y=-y;s=-1} p=y*dn+1/2;scale=0;p/=1;scale=os if(p=A^3)x=fastintpow_(x,p) if(s<0)x=1/x return x } p=fastintpow_(ix,iy)*fastfracpow_(x,fy); scale=os+os if(ix)p*=fastintpow_(x/ix,iy) scale=os return p/1 #The above is usually faster and more accurate than # return( e(y*l(x)) ); } # y-th root of x [ x^(1/y) ] define root(x,y) { return pow(x,1/y) } # Specific cube root function # = stripped down version of fastintroot_(x,3) define cbrt(x) { auto os,d,r,eps; if(x<0)return -cbrt(-x) if(x==0)return 0 os=scale;scale=0;eps=A^(scale/3) if(xeps){ d=r;r=(r+r+x/(r*r))/3 d-=r;if(d<0)d=-d } scale=os return r/1 } # Logarithm of x in given base: log(2, 32) = 5 because 2^5 = 32 # tries to return a real answer where possible when given negative numbers # e.g. log(-2, 64) = 6 because (-2)^6 = 64 # likewise log(-2,-128) = 7 because (-2)^7 = -128 define log(base,x) { auto os,i,l,sx,dn,dnm2; if(base==x)return 1; if(x==0){print "log error: logarithm of zero; negative infinity\n"; return l(0)} if(x==1)return 0; if(base==0){print "log error: zero-based logarithm\n"; return 0 } if(base==1){print "log error: one-based logarithm; positive infinity\n";return -l(0)} scale+=6 if((-1A){c*=int_log(base,A);c-=2*(base<4)}else{c=0}}else{c/=length(base)+1} p=base^c;while(p<=x){.=c++;p*=base} scale=os;return(c-1) } # Lambert's W function 0 branch; Numerically solves w*e(w) = x for w # * is slow to converge near -1/e at high scales define lambertw0(x) { auto oib, a, b, w, ow, lx, ew, e1, eps; if(x==0) return 0; oib=ibase;ibase=A ew = -e(-1) if (x5000) { lx*=l(A);w=lx-(1-1/lx)*l(lx) } else { lx=l(x);w=l(x-4)-(1-1/lx)*l(lx) } # Iteration adapted from code found on Wikipedia # apparently by an anonymous user at 147.142.207.26 # and later another at 87.68.32.52 ow = 0 eps = A^-scale scale += 5 e1 = e(1) while(abs(ow-w)>eps&&w>-1){ ow = w if(x>0){ew=pow(e1,w)}else{ew=e(w)} a = w*ew b = a+ew a -= x; if(a==0)break b = b/a - 1 + 1/(w+1) w -= 1/b if(x<-0.367)w-=eps } scale -= 5 ibase=oib return w/1 } # Lambert's W function -1 branch; Numerically solves w*e(w) = x for w # * is slow to converge near -1/e at high scales define lambertw_1(x) { auto oib,os,oow,ow,w,ew,eps,d,iters; oib=ibase;ibase=A ew = -e(-1) if(ew>x||x>=0) { print "lambertw_1: expected argument in [-1/e,0)\n" ibase=oib if(x==0)return 1-A^scale if(x>0)return 0 return -1 } if(x==ew) return -1; os=scale eps=A^-os scale+=3 oow=ow=0 w=x w=l(-w) w-=l(-w) w+=sqrt(eps) iters=0 while(abs(ow-w)>eps){ oow=ow;ow=w if(w==-1)break w=(x*e(-w)+w*w)/(w+1) if(iters++==A+A||oow==w){iters=0;w-=A^-scale;scale+=2} } scale=os;ibase=oib return w/1 } # LambertW wrapper; takes most useful branch based on x # to pick a branch manually, use lambertw_1 or lambertw0 directly define w(x) { if(x<0)return lambertw_1(x) return lambertw0(x) } # Faster calculation of lambertw0(exp(x)) # . avoids large intermediate value and associated slowness # . numerically solves x = y+ln(y) for y define lambertw0_exp(x) { auto oy,y,eps; # Actual calculation is faster for x < 160 or thereabouts if(xeps)y=x-l(oy=y) return y } # Shorthand alias for the above define w_e(x){ return lambertw0_exp(x) } # Numerically solve pow(y,y) = x for y define powroot(x) { auto r; if(x==0) { print "powroot error: attempt to solve for zero\n" return 0 } if(x==1||x==-1) {return x} if(x<=r=e(-e(-1))){ print "powroot error: unimplemented for values\n <0";r return 0 } r = ln(x) r /= w(r) return r } ## Triangular numbers # xth triangular number define tri(x) { auto xx x=x*(x+1)/2;xx=int(x) if(x==xx)return(xx) return(x) } # 'triangular root' of x define trirt(x) { auto xx x=(sqrt(1+8*x)-1)/2;xx=int(x) if(x==xx)x=xx return(x) } # Workhorse for following 2 functions define tri_step_(t,s) { auto tt t=t+(1+s*sqrt(1+8*t))/2;tt=int(t) if(tt==t)return(tt) return(t) } # Turn tri(x) into tri(x+1) without knowing x define tri_succ(t) { return(tri_step_(t,0+1)) } # Turn tri(x) into tri(x-1) without knowing x define tri_pred(t) { return(tri_step_(t,0-1)) } ## Polygonal Numbers # the xth s-gonal number: # e.g. poly(3, 4) = tri(4) = 1+2+3+4 = 10; poly(4, x) = x*x, etc define poly(s, x) { auto xx x*=(s/2-1)*(x-1)+1;xx=int(x);if(x==xx)x=xx return x } # inverse of the above = polygonal root: # e.g. inverse_poly(3,x)=trirt(x); inverse_poly(4,x)=sqrt(x), etc define inverse_poly(s, r) { auto t,xx t=(s-=2)-2 r=(sqrt(8*s*r+t*t)+t)/s/2;xx=int(r);if(r==xx)r=xx return r } # converse of poly(); solves poly(s,x)=r for s # i.e. if the xth polygonal number is r, how many sides has the polygon? # e.g. if the 5th polygonal number is 15, converse_poly(5,15) = 3 # so the polygon must have 3 sides! (15 is the 5th triangular number) define converse_poly(x,r) { auto xx x=2*((r/x-1)/(x-1)+1);xx=int(x);if(x==xx)x=xx return x } ## Tetrahedral numbers # nth tetrahedral number define tet(n) { return n*(n+1)*(n+2)/6 } # tetrahedral root = inverse of the above define tetrt(t) { auto k,c3,w; if(t==0)return 0 if(t<0)return -2-tetrt(-t) k=3^5*t*t-1 if(k<0){print "tetrt: unimplemented for 0<|t|eps){ox=x;x+=x/n*(n-arigeomean(1,x))} } x+=5*eps scale-=6;return x*y/s }