h8rt3rmin8r

funcs.bc

Jun 17th, 2018
251
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #!/usr/local/bin/bc -l
  2.  
  3. ### Funcs.BC - a large number of functions for use with GNU BC
  4.  
  5.   ## Not to be regarded as suitable for any purpose
  6.   ## Not guaranteed to return correct answers
  7.  
  8. scale=50;
  9. define pi() {
  10.   auto s;
  11.   if(scale==(s=scale(pi_)))return pi_
  12.   if(scale<s)return pi_/1
  13.   scale+=5;pi_=a(1)*4;scale-=5
  14.   return pi_/1
  15. }
  16. e = e(1);
  17. define phi(){return((1+sqrt(5))/2)} ; phi = phi()
  18. define psi(){return((1-sqrt(5))/2)} ; psi = psi()
  19.  
  20. # Reset base to ten
  21. obase=ibase=A;
  22.  
  23. ## Integer and Rounding
  24.  
  25. # Round to next integer nearest 0:  -1.99 -> 1, 0.99 -> 0
  26. define int(x)   { auto os;os=scale;scale=0;x/=1;scale=os;return(x) }
  27.  
  28. # Round down to integer below x
  29. define floor(x) {
  30.   auto os,xx;os=scale;scale=0
  31.   xx=x/1;if(xx>x).=xx--
  32.   scale=os;return(xx)
  33. }
  34.  
  35. # Round up to integer above x
  36. define ceil(x) {
  37.   auto os,xx;x=-x;os=scale;scale=0
  38.   xx=x/1;if(xx>x).=xx--
  39.   scale=os;return(-xx)
  40. }
  41.  
  42. # Fractional part of x:  12.345 -> 0.345
  43. define frac(x) {
  44.   auto os,xx;os=scale;scale=0
  45.   xx=x/1;if(xx>x).=xx--
  46.   scale=os;return(x-xx)
  47. }
  48.  
  49. # Absolute value of x
  50. define abs(x) { if(x<0)return(-x)else return(x) }
  51.  
  52. # Sign of x
  53. define sgn(x) { if(x<0)return(-1)else if(x>0)return(1);return(0) }
  54.  
  55. # Round x up to next multiple of y
  56. define round_up(  x,y) { return(y*ceil( x/y )) }
  57.  
  58. # Round x down to previous multiple of y
  59. define round_down(x,y) { return(y*floor(x/y )) }
  60.  
  61. # Round x to the nearest multiple of y
  62. define round(     x,y) {
  63.   auto os,oib;
  64.   os=scale;oib=ibase
  65.   .=scale++;ibase=A
  66.     y*=floor(x/y+.5)
  67.   ibase=oib;scale=os
  68.   return y
  69. }
  70.  
  71. # Find the remainder of x/y
  72. define int_remainder(x,y) {
  73.   auto os;
  74.   os=scale;scale=0
  75.    x/=1;y/=1;x%=y
  76.   scale=os
  77.   return(x)
  78. }
  79. define remainder(x,y) {
  80.   os=scale;scale=0
  81.    if(x==x/1&&y==y/1){scale=os;return int_remainder(x,y)}
  82.   scale=os
  83.   return(x-round_down(x,y))
  84. }
  85.  
  86. # Greatest common divisor of x and y
  87. define int_gcd(x,y) {
  88.   auto r,os;
  89.   os=scale;scale=0
  90.   x/=1;y/=1
  91.   while(y>0){r=x%y;x=y;y=r}
  92.   scale=os
  93.   return(x)
  94. }
  95. define gcd(x,y) {
  96.   auto r,os;
  97.   os=scale;scale=0
  98.    if(x==x/1&&y==y/1){scale=os;return int_gcd(x,y)}
  99.   scale=os
  100.   while(y>0){r=remainder(x,y);x=y;y=r}
  101.   return(x)
  102. }
  103.  
  104. # Lowest common multiple of x and y
  105. define int_lcm(x,y) {
  106.   auto r,m,os;
  107.   os=scale;scale=0
  108.   x/=1;y/=1
  109.   m=x*y
  110.   while(y>0){r=x%y;x=y;y=r}
  111.   m/=x
  112.   scale=os
  113.   return(m)
  114. }
  115. define lcm(x,y) { return (x*y/gcd(x,y)) }
  116.  
  117. # Remove largest possible power of 2 from x
  118. define oddpart(x){
  119.   auto os;
  120.   os=scale;scale=0;x/=1
  121.   if(x==0){scale=os;return 1}
  122.   while(!x%2)x/=2
  123.   scale=os;return x
  124. }
  125.  
  126. # Largest power of 2 in x
  127. define evenpart(x) {
  128.   auto os;
  129.   os=scale;scale=0
  130.   x/=oddpart(x/1)
  131.   scale=os;return x
  132. }
  133.  
  134. ## Trig / Hyperbolic Trig
  135.  
  136. # Sine
  137. define sin(x) { return s(x) } # alias for standard library
  138. # Cosine
  139. define c(x)   { return s(x+pi()/2) } # as fast or faster than
  140. define cos(x) { return c(x)        } # . standard library
  141. # Tangent
  142. define tan(x)   { auto c;c=c(x);if(c==0)c=A^-scale;return(s(x)/c) }
  143.  
  144. # Secant
  145. define sec(x)   { auto c;c=c(x);if(c==0)c=A^-scale;return(   1/c) }
  146. # Cosecant
  147. define cosec(x) { auto s;s=s(x);if(s==0)s=A^-scale;return(   1/s) }
  148. # Cotangent
  149. define cotan(x) { auto s;s=s(x);if(s==0)s=A^-scale;return(c(x)/s) }
  150.  
  151. # Arcsine
  152. define arcsin(x) { if(x==-1||x==1)return(pi()/2*x);return( a(x/sqrt(1-x*x)) ) }
  153. # Arccosine
  154. define arccos(x) { if(x==0)return(0);return pi()/2-arcsin(x) }
  155.  
  156. # Arctangent (one argument)
  157. define arctan(x)  { return a(x) } # alias for standard library
  158.  
  159. # Arctangent (two arguments)
  160. define arctan2(x,y) {
  161.   auto p;
  162.   if(x==0&&y==0)return(0)
  163.   p=(1-sgn(y))*pi()*(2*(x>=0)-1)/2
  164.   if(x==0||y==0)return(p)
  165.   return(p+a(x/y))
  166. }
  167.  
  168. # Arcsecant
  169. define arcsec(x)      { return( a(x/sqrt(x*x-1)) ) }
  170. # Arccosecant
  171. define arccosec(x)    { return( a(x/sqrt(x*x-1))+pi()*(sgn(x)-1)/2 ) }
  172. # Arccotangent (one argument)
  173. define arccotan(x)    { return( a(x)+pi()/2 ) }
  174. # Arccotangent (two arguments)
  175. define arccotan2(x,y) { return( arctan(x,y)+pi()/2 ) }
  176.  
  177. # Hyperbolic Sine
  178. define sinh(x) { auto t;t=e(x);return((t-1/t)/2) }
  179. # Hyperbolic Cosine
  180. define cosh(x) { auto t;t=e(x);return((t+1/t)/2) }
  181. # Hyperbolic Tangent
  182. define tanh(x) { auto t;t=e(x+x)-1;return(t/(t+2)) }
  183.  
  184. # Hyperbolic Secant
  185. define sech(x)   { auto t;t=e(x);return(2/(t+1/t)) }
  186. # Hyperbolic Cosecant
  187. define cosech(x) { auto t;t=e(x);return(2/(t-1/t)) }
  188. # Hyperbolic Cotangent
  189. define coth(x)   { auto t;t=e(x+x)-1;return((t+2)/t) }
  190.  
  191. # Hyperbolic Arcsine
  192. define arcsinh(x) { return( l(x+sqrt(x*x+1)) ) }
  193. # Hyperbolic Arccosine
  194. define arccosh(x) { return( l(x+sqrt(x*x-1)) ) }
  195. # Hyperbolic Arctangent
  196. define arctanh(x) { return( l((1+x)/(1-x))/2 ) }
  197.  
  198. # Hyperbolic Arcsecant
  199. define arcsech(x)   { return( l((sqrt(1-x*x)+1)/x) ) }
  200. # Hyperbolic Arccosecant
  201. define arccosech(x) { return( l((sqrt(1+x*x)*sgn(x)+1)/x) ) }
  202. # Hyperbolic Arccotangent
  203. define arccoth(x)   { return( l((x+1)/(x-1))/2 ) }
  204.  
  205. # Length of the diagonal vector (0,0)-(x,y) [pythagoras]
  206. define pyth(x,y) { return(sqrt(x*x+y*y)) }
  207. define pyth3(x,y,z) { return(sqrt(x*x+y*y+z*z)) }
  208.  
  209. # Gudermannian Function
  210. define gudermann(x)    { return 2*(a(e(x))-a(1)) }
  211. # Inverse Gudermannian Function
  212. define arcgudermann(x) {
  213.   return arctanh(s(x))
  214. }
  215.  
  216. # Bessel function
  217. define besselj(n,x) { return j(n,x) } # alias for standard library
  218.  
  219. ## Exponential / Logs
  220.  
  221. # Exponential e^x
  222. define exp(x) { return e(x) } # alias for standard library
  223.  
  224. # Natural Logarithm (base e)
  225. define ln(x) {
  226.   auto os,len,ln;
  227.   if(x< 0){print "ln error: logarithm of a negative number\n";return 0}
  228.   if(x==0)print "ln error: logarithm of zero; negative infinity\n"
  229.   len=length(x)-scale(x)-1
  230.   if(len<A)return l(x);
  231.   os=scale;scale+=length(len)+1
  232.   ln=l(x/A^len)+len*l(A)
  233.   scale=os
  234.   return ln/1
  235. } # speed improvement on standard library
  236.  
  237. # workhorse function for pow and log - new, less clever version
  238. # Helps determine whether a fractional power is legitimate for a negative number
  239. # . expects to be fed a positive value
  240. # . returns -odd for even/odd; odd2 for odd1/odd2;
  241. #           even for odd/even;   -2 for irrational
  242. # . note that the return value is the denominator of the fraction if the
  243. #   fraction is rational, and the sign of the return value states whether
  244. #   the numerator is odd (positive) or even (negative)
  245. # . since even/even is not possible, -2 is used to signify irrational
  246. define id_frac2_(y){
  247.   auto os,oib,es,eps,lim,max,p,max2,i,cf[],f[],n,d,t;
  248.   os=scale
  249.   if(cf_max){
  250.     # cf.bc is present!
  251.     .=cf_new(cf[],y);if(scale(cf[0]))return -2;
  252.     .=frac_from_cf(f[],cf[],1)
  253.     d=f[0];scale=0;if(f[1]%2==0)d=-d;scale=os
  254.    return d
  255.   }
  256.   oib=ibase;ibase=A
  257.   scale=0
  258.    es=3*os/4
  259.   scale=os
  260.    eps=A^-es
  261.    y+=eps/A
  262.   scale=es
  263.    y/=1
  264.   scale=0
  265.   if(y<0)y=-y
  266.   d=y-(n=y/1)
  267.   if(d<eps){t=2*(n%2)-1;scale=os;ibase=oib;return t}#integers are x/1
  268.   t=y/2;t=y-t-t
  269.   # Find numerator and denominator of fraction, if any
  270.   lim=A*A;max2=A^5*(max=A^int(os/2));p=1
  271.   i=0;y=t
  272.   while(1) {
  273.     scale=es;y=1/y;scale=0
  274.     y-=(t=cf[++i]=y/1);p*=1+t
  275.     if(i>lim||(max<p&&p<max2)){cf[i=1]=-2;break}#escape if number seems irrational    
  276.     if((p>max2||3*length(t)>es+es)&&i>1){cf[i--]=0;break}#cheat: assume rational
  277.     if(y==0)break;#completely rational
  278.   }
  279.   n=1;d=cf[i]
  280.   if(i==0){print "id_frac2_: something is wrong; y=",y,", d=",d,"\n"}
  281.   if(d!=-2&&i)while(--i){d=n+cf[i]*(t=d);n=t}
  282.   if(d<A^os){d*=2*(n%2)-1}else{d=-2}
  283.   scale=os;ibase=oib
  284.   return d;
  285. }
  286.  
  287. # raise x to integer power y faster than bc's x^y
  288. # . it seems bc (at time of writing) uses
  289. # . an O(n) repeated multiplication algorithm
  290. # . for the ^ operator, which is inefficient given
  291. # . that there is a simple O(log n) alternative:
  292. define fastintpow__(x,y) {
  293.   auto r,hy;
  294.   if(y==0)return(1)
  295.   if(y==1)return(x)
  296.   r=fastintpow__(x,hy=y/2)
  297.   r*=r;if(hy+hy<y)r*=x
  298.   return( r )
  299. }
  300. define fastintpow_(x,y) {
  301.   auto ix,os;
  302.   if(y<0)return fastintpow_(1/x,-y)
  303.   if(y==0)return(1)
  304.   if(y==1)return(x)
  305.   if(x==1)return(1)
  306.   os=scale;scale=0
  307.   if(x==-1){y%=2;y+=y;scale=os;return 1-y}
  308.   # bc is still faster for integers
  309.   if(x==(ix=x/1)){scale=os;return ix^y}
  310.   # ...and small no. of d.p.s, but not for values <= 2
  311.   if(scale(x)<3&&x>2){scale=os;return x^y}
  312.   scale=os;x/=1;scale=0
  313.   x=fastintpow__(x,y);
  314.   scale=os;return x;
  315. }
  316.  
  317. # Raise x to a fractional power faster than e^(y*l(x))
  318. define fastfracpow_(x,y) {
  319.   auto f,yy,inv;
  320.   inv=0;if(y<0){y=-y;inv=1}
  321.   y-=int(y)
  322.   if(y==0)return 1;
  323.   if((yy=y*2^C)!=int(yy)){x=l(x);if(inv)x=-x;return e(y/1*x)}
  324.   # faster using square roots for rational binary fractions
  325.   # where denominator <= 8192
  326.   x=sqrt(x)
  327.   for(f=1;y&&x!=1;x=sqrt(x))if(y+=y>=1){.=y--;f*=x}
  328.   if(inv)f=1/f;
  329.   return f;
  330. }
  331.  
  332. # Find the yth root of x where y is integer
  333. define fastintroot_(x,y) {
  334.   auto os,d,r,ys,eps;
  335.   os=scale;scale=0;y/=1;scale=os
  336.   if(y<0){x=1/x;y=-y}
  337.   if(y==1){return x}
  338.   if(y>=x-1){return fastfracpow_(x,1/y)}
  339.   if(y*int((d=2^F)/y)==d){
  340.     r=1;while(r+=r<=y)x=sqrt(x)
  341.     return x
  342.   }
  343.   scale=length(y)-scale(y);if(scale<5)scale=5;r=e(ln(x)/y)
  344.   scale=os+5;if(scale<5)scale=5
  345.   d=1;eps=A^(3-scale)
  346.   ys=y-1
  347.   while(d>eps){
  348.     d=r;r=(ys*r+x/fastintpow_(r,ys))/y
  349.     d-=r;if(d<0)d=-d
  350.   }
  351.   scale=os
  352.   return r/1
  353. }
  354.  
  355. # Raise x to the y-th power
  356. define pow(x,y) {
  357.  auto os,p,ix,iy,fy,dn,s;
  358.  if(y==0) return 1
  359.  if(x==0) return 0
  360.  if(0<x&&x<1){x=1/x;y=-y}
  361.  os=scale;scale=0
  362.   ix=x/1;iy=y/1;fy=y-iy;dn=0
  363.  scale=os;#scale=length(x/1)
  364.  if(y!=iy&&x<0){
  365.    dn=id_frac2_(y)# -ve implies even numerator
  366.    scale=0;if(dn%2){# odd denominator
  367.      scale=os
  368.      if(dn<0)return  pow(-x,y) # even/odd
  369.      /*else*/return -pow(-x,y) #  odd/odd
  370.    }
  371.    print "pow error: "
  372.    if(dn>0) print "even root"
  373.    if(dn<0) print "irrational power"
  374.    print " of a negative number\n"
  375.    scale=os;return 0
  376.  }
  377.  if(y==iy) {
  378.    if(x==ix){p=fastintpow_(ix,iy);if(iy>0){scale=0;p/=1};scale=os;return p/1}
  379.    scale+=scale;p=fastintpow_(x,iy);scale=os;return p/1
  380.  }
  381.  if((dn=id_frac2_(y))!=-2){ #accurate rational roots (sometimes slower)
  382.    if(dn<0)dn=-dn
  383.    s=1;if(y<0){y=-y;s=-1}
  384.    p=y*dn+1/2;scale=0;p/=1;scale=os
  385.    if(p<A^3)x=fastintpow_(x,p)
  386.    x=fastintroot_(x,dn)
  387.    if(p>=A^3)x=fastintpow_(x,p)
  388.    if(s<0)x=1/x
  389.    return x
  390.  }
  391.  p=fastintpow_(ix,iy)*fastfracpow_(x,fy);
  392.  scale=os+os
  393.  if(ix)p*=fastintpow_(x/ix,iy)
  394.  scale=os
  395.  return p/1
  396.  #The above is usually faster and more accurate than
  397.  # return( e(y*l(x)) );
  398. }
  399.  
  400. # y-th root of x [ x^(1/y) ]
  401. define root(x,y) {
  402.   return pow(x,1/y)
  403. }
  404.  
  405. # Specific cube root function
  406. # = stripped down version of fastintroot_(x,3)
  407. define cbrt(x) {
  408.   auto os,d,r,eps;
  409.   if(x<0)return -cbrt(-x)
  410.   if(x==0)return 0
  411.   os=scale;scale=0;eps=A^(scale/3)
  412.   if(x<eps){scale=os;return 1/cbrt(1/x)}
  413.   scale=5;r=e(ln(x)/3)
  414.   scale=os+5;if(scale<5)scale=5
  415.   d=1;eps=A^(3-scale)
  416.   while(d>eps){
  417.     d=r;r=(r+r+x/(r*r))/3
  418.     d-=r;if(d<0)d=-d
  419.   }
  420.   scale=os
  421.   return r/1
  422. }
  423.  
  424. # Logarithm of x in given base:  log(2, 32) = 5 because 2^5 = 32
  425. #  tries to return a real answer where possible when given negative numbers
  426. #  e.g.     log(-2,  64) = 6 because (-2)^6 =   64
  427. #  likewise log(-2,-128) = 7 because (-2)^7 = -128
  428. define log(base,x) {
  429.   auto os,i,l,sx,dn,dnm2;
  430.   if(base==x)return 1;
  431.   if(x==0){print "log error: logarithm of zero; negative infinity\n";     return  l(0)}
  432.   if(x==1)return 0;
  433.   if(base==0){print "log error: zero-based logarithm\n";                  return    0 }
  434.   if(base==1){print "log error: one-based logarithm; positive infinity\n";return -l(0)}
  435.   scale+=6
  436.   if((-1<base&&base<0)||(0<base&&base<1)){x=-log(1/base,x);scale-=6;return x/1}
  437.   if((-1<x   &&   x<0)||(0<x   &&   x<1)){x=-log(base,1/x);scale-=6;return x/1}
  438.   if(base<0){
  439.     sx=1;if(x<0){x=-x;sx=-1}
  440.     l=log(-base,x)
  441.     dn=id_frac2_(l)
  442.     os=scale;scale=0;dnm2=dn%2;scale=os
  443.     if(dnm2&&dn*sx<0){scale-=6;return l/1}
  444.     print "log error: -ve base: "
  445.     if(dnm2)print "wrong sign for "
  446.     print "implied "
  447.     if(dnm2)print "odd root/integer power\n"
  448.     if(!dnm2){
  449.       if(dn!=-2)print "even root\n"
  450.       if(dn==-2)print "irrational power\n"
  451.     }
  452.     scale-=6;return 0;
  453.   }
  454.   if(x<0){
  455.     print "log error: +ve base: logarithm of a negative number\n"
  456.     scale-=6;return 0;
  457.   }
  458.   x=ln(x)/ln(base);scale-=6;return x/1
  459. }
  460.  
  461. # Integer-only logarithm of x in given base
  462. # (compare digits function in digits.bc)
  463. define int_log(base,x) {
  464.  auto os,p,c;
  465.  if(0<x&&x<1) {return -int_log(base,1/x)}
  466.  os=scale;scale=0;base/=1;x/=1
  467.   if(base<2)base=ibase;
  468.   if(x==0)    {scale=os;return  1-base*A^os}
  469.   if(x<base)  {scale=os;return  0    }
  470.   c=length(x) # cheat and use what bc knows about decimal length
  471.   if(base==A){scale=os;return c-1}
  472.   if(base<A){if(x>A){c*=int_log(base,A);c-=2*(base<4)}else{c=0}}else{c/=length(base)+1}
  473.   p=base^c;while(p<=x){.=c++;p*=base}
  474.   scale=os;return(c-1)
  475. }
  476.  
  477. # Lambert's W function 0 branch; Numerically solves w*e(w) = x for w
  478. # * is slow to converge near -1/e at high scales
  479. define lambertw0(x) {
  480.   auto oib, a, b, w, ow, lx, ew, e1, eps;
  481.   if(x==0) return 0;
  482.   oib=ibase;ibase=A
  483.   ew = -e(-1)
  484.   if (x<ew) {
  485.     print "lambertw0: expected argument in range [-1/e,oo)\n"
  486.     ibase=oib
  487.     return -1
  488.   }
  489.   if (x==ew) {ibase=oib;return -1}
  490.   # First approximation from :
  491.   #   http://www.desy.de/~t00fri/qcdins/texhtml/lambertw/
  492.   #   (A. Ringwald and F. Schrempp)
  493.   # via Wikipedia
  494.   if(x < 0){
  495.     w = x/ew
  496.   } else if(x < 500){
  497.     lx=l(x+1);w=0.665*(1+0.0195*lx)*lx+0.04
  498.   } else if((lx=length(x)-scale(x))>5000) {
  499.     lx*=l(A);w=lx-(1-1/lx)*l(lx)
  500.   } else {
  501.     lx=l(x);w=l(x-4)-(1-1/lx)*l(lx)
  502.   }
  503.   # Iteration adapted from code found on Wikipedia
  504.   #   apparently by an anonymous user at 147.142.207.26
  505.   #   and later another at 87.68.32.52
  506.   ow = 0
  507.   eps = A^-scale
  508.   scale += 5
  509.   e1 = e(1)
  510.   while(abs(ow-w)>eps&&w>-1){
  511.     ow = w
  512.     if(x>0){ew=pow(e1,w)}else{ew=e(w)}
  513.     a = w*ew
  514.     b = a+ew
  515.     a -= x;
  516.     if(a==0)break
  517.     b = b/a - 1 + 1/(w+1)
  518.     w -= 1/b
  519.     if(x<-0.367)w-=eps
  520.   }
  521.   scale -= 5
  522.   ibase=oib
  523.   return w/1
  524. }
  525.  
  526. # Lambert's W function -1 branch; Numerically solves w*e(w) = x for w
  527. # * is slow to converge near -1/e at high scales
  528. define lambertw_1(x) {
  529.   auto oib,os,oow,ow,w,ew,eps,d,iters;
  530.   oib=ibase;ibase=A
  531.   ew = -e(-1)
  532.   if(ew>x||x>=0) {
  533.     print "lambertw_1: expected argument in [-1/e,0)\n"
  534.     ibase=oib
  535.     if(x==0)return 1-A^scale
  536.     if(x>0)return 0
  537.     return -1
  538.   }
  539.   if(x==ew) return -1;
  540.   os=scale
  541.   eps=A^-os
  542.   scale+=3
  543.   oow=ow=0
  544.   w=x
  545.   w=l(-w)
  546.   w-=l(-w)
  547.   w+=sqrt(eps)
  548.   iters=0
  549.   while(abs(ow-w)>eps){
  550.     oow=ow;ow=w
  551.     if(w==-1)break
  552.     w=(x*e(-w)+w*w)/(w+1)
  553.     if(iters++==A+A||oow==w){iters=0;w-=A^-scale;scale+=2}
  554.   }
  555.   scale=os;ibase=oib
  556.   return w/1
  557. }
  558.  
  559. # LambertW wrapper; takes most useful branch based on x
  560. # to pick a branch manually, use lambertw_1 or lambertw0 directly
  561. define w(x) {
  562.   if(x<0)return lambertw_1(x)
  563.   return lambertw0(x)
  564. }
  565.  
  566. # Faster calculation of lambertw0(exp(x))
  567. # . avoids large intermediate value and associated slowness
  568. # . numerically solves x = y+ln(y) for y
  569. define lambertw0_exp(x) {
  570.   auto oy,y,eps;
  571.   # Actual calculation is faster for x < 160 or thereabouts
  572.   if(x<C*D)return lambertw0(e(x));
  573.   oy=0;y=l(x);y=x-y+y/x;eps=A^-scale
  574.   while(abs(oy-y)>eps)y=x-l(oy=y)
  575.   return y
  576. }
  577.  
  578. # Shorthand alias for the above
  579. define w_e(x){ return lambertw0_exp(x) }
  580.  
  581. # Numerically solve pow(y,y) = x for y
  582. define powroot(x) {
  583.   auto r;
  584.   if(x==0) {
  585.     print "powroot error: attempt to solve for zero\n"
  586.     return 0
  587.   }
  588.   if(x==1||x==-1) {return x}
  589.   if(x<=r=e(-e(-1))){
  590.     print "powroot error: unimplemented for values\n  <0";r
  591.     return 0
  592.   }
  593.   r = ln(x)
  594.   r /= w(r)
  595.   return r
  596. }
  597.  
  598. ## Triangular numbers
  599.  
  600. # xth triangular number
  601. define tri(x) {
  602.   auto xx
  603.   x=x*(x+1)/2;xx=int(x)
  604.   if(x==xx)return(xx)
  605.   return(x)
  606. }
  607.  
  608. # 'triangular root' of x
  609. define trirt(x) {
  610.   auto xx
  611.   x=(sqrt(1+8*x)-1)/2;xx=int(x)
  612.   if(x==xx)x=xx
  613.   return(x)
  614. }
  615.  
  616. # Workhorse for following 2 functions
  617. define tri_step_(t,s) {
  618.   auto tt
  619.   t=t+(1+s*sqrt(1+8*t))/2;tt=int(t)
  620.   if(tt==t)return(tt)
  621.   return(t)
  622. }
  623.  
  624. # Turn tri(x) into tri(x+1) without knowing x
  625. define tri_succ(t) {
  626.   return(tri_step_(t,0+1))
  627. }
  628.  
  629. # Turn tri(x) into tri(x-1) without knowing x
  630. define tri_pred(t) {
  631.   return(tri_step_(t,0-1))
  632. }
  633.  
  634. ## Polygonal Numbers
  635.  
  636. # the xth s-gonal number:
  637. #   e.g. poly(3, 4) = tri(4) = 1+2+3+4 = 10; poly(4, x) = x*x, etc
  638. define poly(s, x) {
  639.   auto xx
  640.   x*=(s/2-1)*(x-1)+1;xx=int(x);if(x==xx)x=xx
  641.   return x
  642. }
  643.  
  644. # inverse of the above = polygonal root:
  645. #   e.g. inverse_poly(3,x)=trirt(x); inverse_poly(4,x)=sqrt(x), etc
  646. define inverse_poly(s, r) {
  647.   auto t,xx
  648.   t=(s-=2)-2
  649.   r=(sqrt(8*s*r+t*t)+t)/s/2;xx=int(r);if(r==xx)r=xx
  650.   return r
  651. }
  652.  
  653. # converse of poly(); solves poly(s,x)=r for s
  654. #   i.e. if the xth polygonal number is r, how many sides has the polygon?
  655. #   e.g. if the 5th polygonal number is 15, converse_poly(5,15) = 3
  656. #     so the polygon must have 3 sides! (15 is the 5th triangular number)
  657. define converse_poly(x,r) {
  658.   auto xx
  659.   x=2*((r/x-1)/(x-1)+1);xx=int(x);if(x==xx)x=xx
  660.   return x
  661. }
  662.  
  663. ## Tetrahedral numbers
  664.  
  665. # nth tetrahedral number
  666. define tet(n) { return n*(n+1)*(n+2)/6 }
  667.  
  668. # tetrahedral root = inverse of the above
  669. define tetrt(t) {
  670.   auto k,c3,w;
  671.   if(t==0)return 0
  672.   if(t<0)return -2-tetrt(-t)
  673.   k=3^5*t*t-1
  674.   if(k<0){print "tetrt: unimplemented for 0<|t|<sqrt(3^-5)\n"; return 0}
  675.   c3=cbrt(3)
  676.   k=cbrt(sqrt(3*k)+3^3*t)
  677.   return k/c3^2+1/(c3*k)-1
  678. }
  679.  
  680. ## Arithmetic-Geometric mean
  681.  
  682. define arigeomean(a,b) {
  683.   auto c,s;
  684.   if(a==b)return a;
  685.   s=1;if(a<0&&b<0){s=-1;a=-a;b=-b}
  686.   if(a<0||b<0){print "arigeomean: mismatched signs\n";return 0}
  687.   while(a!=b){c=(a+b)/2;a=sqrt(a*b);b=c}
  688.   return s*a
  689. }
  690.  
  691. # solve n = arigeomean(x,y)
  692. define inv_arigeomean(n, y){
  693.   auto ns,ox,x,b,c,d,i,s,eps;
  694.   if(n==y)return n;
  695.   s=1;if(n<0&&y<0){s=-1;n=-n;y=-y}
  696.   if(n<0||y<0){print "inv_arigeomean: mismatched signs\n";return 0}  
  697.   if(n<y){x=y;y=n;n=x}
  698.   n/=y
  699.   scale+=2;eps=A^-scale;scale+=4
  700.   ns=scale
  701.   x=n*(1+ln(n));ox=-1
  702.   for(i=0;i<A;i++){
  703.     # try to force quadratic convergence
  704.     if(abs(x-ox)<eps){i=-1;break}
  705.     ox=x;scale+=scale
  706.     b=x+x/n*(n-arigeomean(1,x));
  707.     c=b+b/n*(n-arigeomean(1,b));
  708.     d=b+b-c-x
  709.     if(d){x=(b*b-c*x)/d}else{x=b;i=-1;break}
  710.     scale=ns
  711.   }
  712.   if(i!=-1){
  713.     # give up and converge linearly
  714.     x=(x+ox)/2
  715.     while(abs(x-ox)>eps){ox=x;x+=x/n*(n-arigeomean(1,x))}
  716.   }
  717.   x+=5*eps
  718.   scale-=6;return x*y/s
  719. }
RAW Paste Data