Lib/software/bc/funcs.bc

   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 }