# 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