Mathematics.rex
Introduction
This year I started posting on RosettaCode and found the copy-and-paste way to get a working program annoying. Therefore I'll gather all procedures I have until now and present them here in a set libraries and stand-alone files. Other entries will link to this page.
See documentation for help and tips.
Complex
/* Complex numbers library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2024 */
/* Numbers are coded as string 're im' (rectangular) */
/* Example '123 456' = '123+456i' */
Cabs:
/* Absolute value (modulus, magnitude) */
procedure expose glob.
arg xx
/* Modulus */
return Cmod(xx)
Cadd:
/* Add */
procedure
if Arg() < 2 then say argument
do n = 1 to Arg()
xx = Arg(n) 0 0; parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Sum */
if n = 1 then do
zr = xr; zi = xi
end
else do
zr = zr+xr; zi = zi+xi
end
end
return (zr) (zi)
Carg:
/* Argument (phase) */
procedure expose glob.
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Argument between -pi and pi */
if xr = 0 & xi = 0 then
return 0
else
return Atan2(xi,xr)
Cconj:
/* Conjugate */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Conjugate */
return (xr) (-xi)
Ccueq:
/* Roots of cubic ax^3+bx^2+cx+d = 0 */
procedure expose cueq. glob.
arg xx; xx = xx 0 0 0 0
parse var xx a b c d .
if \ Number(a) then say argument
if \ Number(b) then say argument
if \ Number(c) then say argument
if \ Number(d) then say argument
if a = 0 then say argument
/* Method Cardano et al */
numeric digits Digits()*2
/* Reduce to x^3+ax^2+bx+c */
parse value b/a c/a d/a with a1 b1 c1
/* Reduce to x^3+ax+b */
h1 = a1*a1; h2 = h1*a1
a2 = b1-h1/3; b2 = -h2/27+h2/9-a1*b1/3+c1
/* Discriminant */
s = Round(b2*b2/4+a2*a2*a2/27,Digits()/2)
/* Roots */
cueq. = 0
select
when s < 0 then do
h1 = Abs(a2)/3
f = Acos(-b2/(2*Sqrt(h1*h1*h1)))
h1 = 2*Sqrt(h1); h2 = f/3; h3 = Pi()/3; h4 = a1/3
cueq.root.1 = h1*Cos(h2)-h4
cueq.root.2 = -h1*Cos(h2-h3)-h4
cueq.root.3 = -h1*Cos(h2+h3)-h4
end
when s = 0 then do
h1 = Curt(-b2/2); h2 = a1/3
cueq.root.1 = 2*h1-h2
cueq.root.2 = -h1-h2
cueq.root.3 = cueq.root.2
end
when s > 0 then do
h1 = b2/2
f = Sqrt((h1)**2+(a2*a2*a2/27)); g = Curt(-h1+f); h = Curt(-h1-f)
h1 = g+h; h2 = -h1/2; h3 = a1/3; h4 = Sqrt(3)*(g-h)/2
cueq.root.1 = h1-h3
cueq.root.2 = (h2-h3) (h4)
cueq.root.3 = (h2-h3) (-h4)
end
end
/* Discard roundoff digits */
do i = 1 to 3
cueq.root.i = Cround(cueq.root.i,Digits()-4)
end
numeric digits Digits()/2
/* Normalize */
do i = 1 to 3
cueq.root.i = Cnorm(cueq.root.i)
end
/* Number of roots */
cueq.0 = 3
return 3
Cdiv:
/* Divide */
procedure
if Arg() < 2 then say argument
do n = 1 to Arg()
xx = Arg(n) 0 0; parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Quotient */
if n = 1 then do
zr = xr; zi = xi
end
else do
if xr = 0 & xi = 0 then say argument
d = xr*xr+xi*xi
r = (zr*xr+zi*xi)/d; i = (zi*xr-zr*xi)/d
zr = r; zi = i
end
end
return (zr) (zi)
Cequ:
/* Equality */
procedure
arg xx,yy; xx = xx 0 0; yy = yy 0 0
parse var xx xr xi .
parse var yy yr yi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
if \ Number(yr) then say argument
if \ Number(yi) then say argument
/* Equal? */
if xr = yr & xi = yi then
return 1
else
return 0
Ceval:
/* Evaluate polynomial with real coefficients */
procedure
arg xx,yy; yy = yy 0 0
parse var yy yr yi .
if xx = '' then say argument
if \ Number(yr) then say argument
if \ Number(yi) then say argument
/* Evaluate */
zz = 0 0
numeric digits Digits()+2
do n = 1 to Words(xx)
zz = Cadd(Cmul(zz,yy),Word(xx,n))
end
numeric digits Digits()-2
/* Normalize */
return Cnorm(zz)
Cim:
/* Imaginary part */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Imagery part */
return (xi+0)
Cinv:
/* Invert */
procedure
arg xx
/* Inversion */
return Cdiv(1,xx)
Cmod:
/* Modulus (absolute value, magnitude) */
procedure expose glob.
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Modulus */
return Sqrt(xr*xr+xi*xi)
Cmul:
/* Multiply */
procedure
if Arg() < 2 then say argument
do n = 1 to Arg()
xx = Arg(n) 0 0; parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Product */
if n = 1 then do
zr = xr; zi = xi
end
else do
r = zr*xr-zi*xi; i = zr*xi+zi*xr
zr = r; zi = i
end
end
return (zr) (zi)
Cneg:
/* Negate */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Negation */
return (-xr) (-xi)
Cnorm:
/* Normalize */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Round to numeric digits and trim zeroes */
return (xr/1) (xi/1)
Complex:
/* Complex? */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
if xi = 0 then
return 0
else
return 1
Cpol2form:
/* Polar -> formula */
procedure
arg xx; xx = xx 0 0
parse var xx xm xa .
if \ Number(xm) then say argument
if \ Number(xa) then say argument
/* Formula */
if xm = 0 then
return 0
za = xm'r'
select
when xa < 0 then
zb = '-'Abs(xa)'th'
when xa = 0 then
zb = ''
when xa > 0 then
zb = ','xa'th'
end
return za||zb
Crect:
/* Polar -> rectangular */
procedure expose glob.
arg xx; xx = xx 0 0
parse var xx xm xa .
if \ Number(xm) then say argument
if \ Number(xa) then say argument
/* Rectangular */
return (xm*Cos(xa)) (xm*Sin(xa))
Cpower:
/* Integer power */
procedure
arg xx,yy; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
/* Repeated divide or multiply */
numeric digits Digits()+2
select
when yy < 0 then do
zz = Cinv(xx)
do Abs(yy)-1
zz = Cdiv(zz,xx)
end
end
when yy = 0 then
return 1
when yy > 0 then do
zz = xx
do yy-1
zz = Cmul(zz,xx)
end
end
end
/* Normalize */
numeric digits Digits()-2
return Cnorm(zz)
Cqueq:
/* Roots of quartic ax^4+bx^3+cx^2+dx+e = 0 */
procedure expose queq. glob.
arg xx; xx = xx 0 0 0 0 0
parse var xx a b c d e .
if \ Number(a) then say argument
if \ Number(b) then say argument
if \ Number(c) then say argument
if \ Number(d) then say argument
if \ Number(e) then say argument
if a = 0 then say argument
/* Method Abramowitz and Stegun, courtesy PlanetCalc */
numeric digits Digits()*2
/* Reduce to x^4+ax^3+bx^2+cx+d */
parse value b/a c/a d/a e/a with a3 a2 a1 a0
/* Solve resolvent cubic x^3+ax^2+bx+c for real roots */
n = Cueq((1) (-a2) (a1*a3-4*a0) (4*a0*a2-a1**2-a0*a3**2))
/* Use max root */
u1 = cueq.root.1
do i = 2 to n
u1 = Max(cueq.root.i,u1)
end
/* Compute coefficients square equations */
h = Sqrt(Round(a3**2/4+u1-a2,Digits()-2)); p1 = a3/2+h; p2 = a3/2-h
h = Sqrt(Round(u1**2/4-a0,Digits()-2)); q1 = u1/2+h; q2 = u1/2-h
if Round(p1*q2+p2*q1-a1,Digits()/2) <> 0 then
parse value q2 q1 with q1 q2
/* Find roots by solving ax^2+bx+c for complex roots */
queq. = 0
call Csqeq 1 p1 q1
queq.root.1 = sqeq.root.1
queq.root.2 = sqeq.root.2
call Csqeq 1 p2 q2
queq.root.3 = sqeq.root.1
queq.root.4 = sqeq.root.2
/* Discard roundoff digits */
do i = 1 to 4
queq.root.i = Cround(queq.root.i,Digits()-4)
end
/* Normalize */
numeric digits Digits()/2
do i = 1 to 4
queq.root.i = Cnorm(queq.root.i)
end
/* Number of roots */
queq.0 = 4
return 4
Cre:
/* Real part */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Real part */
return xr+0
Crect2form:
/* Rectangular -> formula */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Formula */
if xr = 0 then
za = ''
else
za = xr
select
when xi < 0 then
if xi = -1 then
zb = '-i'
else
zb = xi'i'
when xi = 0 then
if xr = 0 then
zb = 0
else
zb = ''
when xi > 0 then
if xi = 1 then
zb = '+i'
else
zb = '+'xi'i'
end
return za||zb
Cpol:
/* Rectangular -> polar */
procedure expose glob.
arg xx
/* Polar */
return (Cmod(xx)) (Carg(xx))
Cround:
/* Round to specified decimals */
procedure
arg xx,yy; xx = xx 0 0
parse var xx xr xi .
if yy = '' then
yy = Digits()
if \ Number(xr) then say argument
if \ Number(xi) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
/* Round */
return (Round(xr,yy)) (Round(xi,yy))
Cscale:
/* Scale */
procedure
arg xx,yy; yy = yy 0 0
parse var yy yr yi .
if \ Number(xx) then say argument
if \ Number(yr) then say argument
if \ Number(yi) then say argument
/* Multiply by constant */
return (xx*yr) (xx*yi)
Csci:
/* Convert to scientific format */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Force format */
return (Format(xr,,,,0)) (Format(xi,,,,0))
Csqeq:
/* Roots of square ax^2+bx+c = 0 */
procedure expose sqeq. glob.
arg xx; xx = xx 0 0 0
parse var xx a b c .
if \ Number(a) then say argument
if \ Number(b) then say argument
if \ Number(c) then say argument
if a = 0 then say argument
/* Abc formula */
numeric digits Digits()+5
/* Discriminant */
s = Round(b*b-4*a*c,Digits()-3)
/* Roots */
sqeq. = 0; h1 = -1/(2*a)
select
when s < 0 then do
h2 = Sqrt(-s)
sqeq.root.1 = (h1*b) (h1*h2)
sqeq.root.2 = (h1*b) (-h1*h2)
end
when s = 0 then do
sqeq.root.1 = h1*b
sqeq.root.2 = sqeq.root.1
end
when s > 0 then do
h2 = Sqrt(s)
sqeq.root.1 = h1*(b+h2)
sqeq.root.2 = h1*(b-h2)
end
end
/* Discard roundoff digits */
do i = 1 to 2
sqeq.root.i = Cround(sqeq.root.i,Digits()-3)
end
numeric digits Digits()-5
/* Normalize */
do i = 1 to 2
sqeq.root.i = Cnorm(sqeq.root.i)
end
/* Number of roots */
sqeq.0 = 2
return 2
Csqrt:
/* Square root */
procedure expose glob.
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Principal value */
numeric digits Digits()+2
m = Cmod(xx)
if xi < 0 then
s = -1
else
s = 1
a = Round(m+xr,Digits()-2); b = Round(m-xr,Digits()-2)
r = Sqrt(a/2); i = Sqrt(b/2)*s
numeric digits Digits()-2
return (r/1) (i/1)
Cstd:
/* Convert to standard format */
procedure
arg xx; xx = xx 0 0
parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Force format */
return (Format(xr,,,0)) (Format(xi,,,0))
Csub:
/* Subtract */
procedure
if Arg() < 2 then say argument
do n = 1 to Arg()
xx = Arg(n) 0 0; parse var xx xr xi .
if \ Number(xr) then say argument
if \ Number(xi) then say argument
/* Difference */
if n = 1 then do
zr = xr; zi = xi
end
else do
zr = zr-xr; zi = zi-xi
end
end
return (zr) (zi)
Constants
/* Constants library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */
Apery:
/* Apery constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.20205690315959428539973816151144999076498629234049888179227155534183820578631309018645587360933526
else do
/* Wedeniwski recurring */
numeric digits Digits()+2
numeric fuzz 2
/* Term k = 0 */
f10 = 1; f20 = 1; f21 = 1; f32 = 2; f43 = 6; p = 12463; s = 1
zz = 0; v = 0
do k = 1
/* Add term */
zz = zz+s*((f10*f20*f21)**3*p)/(f32*f43**3)
if zz = v then
leave
v = zz
/* Prepare next term */
k2 = k+k; k3 = k2+k; k4 = k3+k
/* Factorials */
f10 = f10*k; f20 = f20*k2*(k2-1); f21 = f21*k2*(k2+1); f32 = f32*k3*(k3+1)*(k3+2); f43 = f43*k4*(k4+1)*(k4+2)*(k4+3)
/* Polynomial */
p = ((((126392*k+412708)*k+531578)*k+336367)*k+104000)*k+12463
/* Sign */
s = -s
end
zz = zz/24
numeric digits Digits()-2
end
return zz/1
Artin:
/* Artin constant */
procedure
/* Fast value */
zz = 0.373955813619202288054728054346416415111629248606150042094742802417350182040028082344304317087250568981603
return zz/1
Backhouse:
/* Backhouse constant */
procedure
/* Fast value */
zz = 1.45607494858268967139959535111654355765317837484713154027070243741400150626538989559964531940186030910992
return zz/1
Bernstein:
/* Bernstein constant */
procedure
/* Fast value */
zz = 0.28016949902386913303643649123067200004248213981236
return zz/1
Bloch:
/* Bloch constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.471861653452681784874468793611316149077012621739443243631460266888006805241859656902448655624594840993
else do
/* Formula */
numeric digits Digits()+1
zz = Sqrt((Sqrt(3)-1)*0.5)*Gamma(1/3)*Gamma(11/12)/Gamma(0.25)
numeric digits Digits()-1
end
return zz/1
Bronze:
/* Bronze ratio constant = positive root of x^2-3x-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 3.302775637731994646559610633735247973125648286922623106355226528113583474146505222602309541009245359
else
/* Formula */
zz = 0.5*(3+Sqrt(13))
return zz/1
Brun:
/* Brun constant */
procedure
/* Fast value */
zz = 1.902160583104
return zz/1
Cahen:
/* Cahen constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 0.6434105462883380261822543077575647632865878602682395059870309203074927764618326108484408955504634320
else do
numeric digits Digits()+2
/* Sylvester */
zz = 0; v = 0; s = 2
do n = 0
zz = zz+(-1)**n/(s-1)
if zz = v
then leave
v = zz; s = s*s-s+1
end
numeric digits Digits()-2
end
return zz/1
Capacity:
/* Logarithmic capacity of the unit disk constant */
procedure
/* Fast value */
zz = 0.5901702995080481130226689702792442936168583174407236497579321997021520903603578974892293080979039771
return zz/1
Catalan1:
/* Catalan constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062
else do
numeric digits Digits()+2
numeric fuzz 2
/* Broadhurst */
s1 = 0; v = 0
do n = 0 until s1 = v
v = s1
a = 1/2**(4*n); b = 8*n
c1 = (b+1)*(b+1); c2 = (b+2)*(b+2); c3 = (b+3)*(b+3); c5 = (b+5)*(b+5); c6 = (b+6)*(b+6); c7 = (b+7)*(b+7)
s1 = s1+a*(1/(2*c1)-1/(2*c2)+1/(4*c3)-1/(8*c5)+1/(8*c6)-1/(16*c7))
end
s2 = 0; v = 0
do n = 0 until s2 = v
v = s2
a = 1/2**(12*n); b = 8*n
c1 = (b+1)*(b+1); c2 = (b+2)*(b+2); c3 = (b+3)*(b+3); c5 = (b+5)*(b+5); c6 = (b+6)*(b+6); c7 = (b+7)*(b+7)
s2 = s2+a*(1/(8*c1)+1/(16*c2)+1/(64*c3)-1/(512*c5)-1/(1024*c6)-1/(4096*c7))
end
zz = 3*s1-2*s2
numeric digits Digits()-2
end
return zz/1
Chaitin:
/* Chaitin constant */
procedure
/* Fast value */
zz = 0.0078749969978123844
return zz/1
Champernowne:
/* Champernowne constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 0.12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565
else do
p = Digits()+4
/* Concatenate integers */
zz = '0.'
do n = 1 while Length(zz) <= p
zz = zz||n
end
end
return zz/1
Conway2:
/* Conway constant = real root of x^71-x^69-2x^68...-6x^2+3x-6 */
procedure
/* Fast value or first guess */
zz = 1.30357726903429639125709911215255189073070250465940487575486139062855088785246155712681576686442522555
if Digits() > 100 then do
/* Newton */
numeric digits Digits()+3
numeric fuzz 3
do until zz = v
v = zz
/* Powers */
y2 = zz*zz; y3 = zz*y2; y4 = zz*y3; y5 = zz*y4; y6 = zz*y5; y7 = zz*y6; y8 = zz*y7; y9 = zz*y8
y10 = zz*y9; y11 = zz*y10; y12 = zz*y11; y13 = zz*y12; y14 = zz*y13; y15 = zz*y14; y16 = zz*y15; y17 = zz*y16; y18 = zz*y17; y19 = zz*y18
y20 = zz*y19; y21 = zz*y20; y22 = zz*y21; y23 = zz*y22; y24 = zz*y23; y25 = zz*y24; y26 = zz*y25; y27 = zz*y26; y28 = zz*y27; y29 = zz*y28
y30 = zz*y29; y31 = zz*y30; y32 = zz*y31; y33 = zz*y32; y34 = zz*y33; y35 = zz*y34; y36 = zz*y35; y37 = zz*y36; y38 = zz*y37; y39 = zz*y38
y40 = zz*y39; y41 = zz*y40; y42 = zz*y41; y43 = zz*y42; y44 = zz*y43; y45 = zz*y44; y46 = zz*y45; y47 = zz*y46; y48 = zz*y47; y49 = zz*y48
y50 = zz*y49; y51 = zz*y50; y52 = zz*y51; y53 = zz*y52; y54 = zz*y53; y55 = zz*y54; y56 = zz*y55; y57 = zz*y56; y58 = zz*y57; y59 = zz*y58
y60 = zz*y59; y61 = zz*y60; y62 = zz*y61; y63 = zz*y62; y64 = zz*y63; y65 = zz*y64; y66 = zz*y65; y67 = zz*y66; y68 = zz*y67; y69 = zz*y68
y70 = zz*y69; y71 = zz*y70
/* Polynomial and derivative */
f = y71-y69-2*y68-y67+2*y66+2*y65+y64-y63-y62-y61-y60,
-y59+2*y58+5*y57+3*y56-2*y55-10*y54-3*y53-2*y52+6*y51+6*y50,
+y49+9*y48-3*y47-7*y46-8*y45-8*y44+10*y43+6*y42+8*y41-5*y40,
-12*y39+7*y38-7*y37+7*y36+y35-3*y34+10*y33+y32-6*y31-2*y30,
-10*y29-3*y28+2*y27+9*y26-3*y25+14*y24-8*y23-7*y21+9*y20,
+3*y19-4*y18-10*y17-7*y16+12*y15+7*y14+2*y13-12*y12-4*y11-2*y10,
+5*y9+y7-7*y6+7*y5-4*y4+12*y3-6*y2+3*zz-6
f1 = 71*y70-69*y68-136*y67-67*y66+132*y65+130*y64+64*y63-63*y62-62*y61-61*y60-60*y59,
-59*y58+116*y57+285*y56+168*y55-110*y54-540*y53-159*y52-104*y51+306*y50+300*y49,
+49*y48+432*y47-141*y46-322*y45-360*y44-352*y43+430*y42+252*y41+328*y40-200*y39,
-468*y38+266*y37-259*y36+252*y35+35*y34-102*y33+330*y32+ 32*y31-186*y30-60*y29,
-290*y28-84*y27+54*y26+234*y25-75*y24+336*y23-184*y22-147*y20+180*y19,
+57*y18-72*y17-170*y16-112*y15+170*y14+98*y13+ 26*y12-144*y11-44*y10-20*y9,
+45*y8+7*y6-42*y5+35*y4-16*y3+36*y2-12*zz+3
/* Next value */
zz = zz-f/f1
end
numeric digits Digits()-3
end
return zz/1
Copeland:
/* Copeland-Erdos constant */
procedure
/* Fast value */
zz = 0.23571113171923293137414347535961677173798389971011031071091131271311371391491511571631671731791811911
return zz/1
Devicci:
/* DeVicci constant = positive root of 4x^8-28x^6-7x^4+16x^2+16 */
procedure
/* Fast value or first guess */
zz = 1.007434756884279376098253595231099141925690114113669770234963798571152313280286777962520551474635923942
if Digits() > 100 then do
/* Newton */
numeric digits Digits()+2
numeric fuzz 2
zz = 1.007434756884279376098253595231099141925690114113669770234963798571152313280286777962520551474635923942
do until zz = v
v = zz
/* Powers */
y2 = zz*zz; y3 = y2*zz; y4 = y3*zz; y5 = y4*zz; y6 = y5*zz; y7 = y6*zz; y8 = y7*zz
/* Polynomial and derivatives */
f = 4*y8-28*y6-7*y4+16*y2+16
f1 = 32*y7-168*y5-28*y3+32*zz
/* Next value */
zz = zz-f/f1
end
numeric digits Digits()-2
end
return zz/1
Dottie:
/* Dottie constant = root of Cos(x) - x */
procedure expose glob.
/* Fast value or first guess */
zz = 0.7390851332151606416553120876738734040134117589007574649656806357732846548835475945993761069317665318
if Digits() > 100 then do
numeric digits Digits()+2
/* Halley */
do until zz = v
numeric fuzz
v = zz
c = Cos(zz); s = Sin(zz); f = c-zz; f1 = -s-1; f2 = -c
zz = zz - (2*f*f1) / (2*f1*f1-f*f2)
numeric fuzz 2
end
numeric fuzz
numeric digits Digits()-2
end
return zz/1
Dubois:
/* Second Dubois-Reymond constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.194528049465325113615213730287503906590157785275923662043563911261286898039528881692156242539560897386876
else
/* Formula */
zz = (E()**2-7)*0.5
return zz/1
E:
/* Euler constant */
procedure expose glob.
p = Digits()
/* In memory? */
if glob.e.p <> '' then
return glob.e.p
if p < 101 then
/* Fast value */
glob.e.p = 2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516643+0
else do
numeric digits Digits()+2
numeric fuzz 2
/* Taylor */
zz = 2; t = 1; v = zz
do n = 2 until zz = v
v = zz
t = t/n; zz = zz+t
end
numeric digits Digits()-2
glob.e.p = zz+0
end
return glob.e.p
Embree:
/* Embree-Trefethen constant */
procedure
/* Fast value */
zz = 0.70258
return zz/1
Erdos1:
/* Erdos-Borwein constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.606695152415291763783301523190924580480579671505756435778079553691418420743486690565711801670155575897045429
else do
numeric digits Digits()+2
numeric fuzz 2
/* Series */
zz = 0; v = 0
do n = 1 until zz = v
v = zz
a = 2**n; zz = zz+(1/a**n)*(a+1)/(a-1)
end
numeric digits Digits()-2
end
return zz/1
Erdos2:
/* Erdos-Tenenbaum-Ford constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.08607133205593420688757309877692267776059110953033317349202023666542263581462287979938053460252876807163
else
/* Formula */
zz = 1-(1+Ln(Ln(2)))/Ln(2)
return zz/1
Euler:
/* Euler-Mascheroni constant */
procedure expose glob.
p = Digits()
/* In memory? */
if glob.euler.p <> '' then
return glob.euler.p
if p < 101 then
/* Fast value */
glob.euler.p = 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495+0
else do
numeric digits Digits()+2
/* Brent McMillan */
n = Ceil((Digits()*Ln(10)+Ln(Pi()))*0.25); m = Ceil(2.07*Digits())
n2 = n*n; ak = -Ln(n); bk = 1; s = ak; v = 1
do k = 1 to m
bk = bk*n2/(k*k); ak = (ak*n2/k+bk)/k
s = s+ak; v = v+bk
end
numeric digits Digits()-2
glob.euler.p = s/v
end
return glob.euler.p
Feigenbaum1:
/* Feigenbaum first constant */
procedure
/* Fast value */
zz = 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716
return zz/1
Feigenbaum2:
/* Feigenbaum second constant */
procedure
/* Fast value */
zz = 2.50290787509589282228390287321821578638127137672714997733619205677923546317959020670329964974643383412959
return zz/1
Feller:
/* Feller-Tornier constant */
procedure
/* Fast value */
zz = 0.66131704946962233528976584627411853328547528983291635498090562622662503174312230494226174078428187
return zz/1
Foias:
/* Foias constant */
procedure
/* Fast value */
zz = 1.18745235112650105459548015839651935121569268158586035301010412619878041872352540738702465760608657943378
return zz/1
Fraction:
/* First continued fraction constant */
procedure
/* Fast value */
zz = 0.697774657964007982006790592551752599486658262998021232368630082816530852764641112996965654182676568723982
return zz/1
Fransen:
/* Fransen-Robertson constant */
procedure
/* Fast value */
zz = 2.80777024202851936522150118655777293230808592093019829122005480959710088912190166551018530816819663814187
return zz/1
Gauss:
/* Gauss constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 0.8346268416740731862814297327990468089939930134903470024498273701036819927095264118696911603512753241
else
/* Formula */
zz = 1/Agmean(1,Sqrt(2))
return zz/1
Gelfond1:
/* Gelfond constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 23.1406926327792690057290863679485473802661062426002119934450464095243423506904527835169719970675492
else
/* Formula */
zz = Pow(E(),Pi())
return zz/1
Gelfond2:
/* Gelfond-Schneider constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 2.66514414269022518865029724987313984827421131371465949283597959336492044617870595486760918000519642
else
/* Formula */
zz = Pow(2,Sqrt(2))
return zz/1
Gieseking:
/* Gieseking constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 1.01494160640965362502120255427452028594168930753029979201748910677659747625824402213647035422825669
else
/* Formula */
zz = (Trigamma(1/3)-Trigamma(2/3))/(4*Sqrt(3))
return zz/1
Glaisher:
/* Glaisher-Kinkelin constant */
procedure
/* Fast value */
zz = 1.28242712910062263687534256886979172776768892732500119206374002174040630885882646112973649195820237
return zz/1
Golden:
/* Golden ratio constant = root of x^2-x-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939114
else
/* Formula */
zz = (Sqrt(5)+1)*0.5
return zz/1
Goldenangle:
/* Golden angle constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 2.39996322972865332223155550663361385312499901105811504293511275073130733823943879077996206066058396
else
/* Formula */
zz = Pi()*(3-Sqrt(5))
return zz/1
Golomb1:
/* Golomb constant */
procedure
/* Fast value */
zz = 0.624329988543550870992936383100837244179642620180529286973551902495638088855113254462460276195539868869
return zz/1
Gompertz:
/* Gompertz constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.596347362323194074341078499369279376074177860152548781573484910482327219114874417470430497093612760344237
else do
numeric digits Digits()+2
/* Taylor */
zz = 0; s = 1; f = 1
do n = 1 until zz = v
numeric fuzz
v = zz
s = -s; f = f*n; zz = zz+s/(n*f)
numeric fuzz 2
end
numeric fuzz
zz = -E()*(Euler()+zz)
numeric digits Digits()-2
end
return zz/1
Hafner:
/* Hafner-Sarnak-McCurley constant */
procedure
/* Fast value */
zz = 0.353236371854995984543516550432682011280164778566690446416085942814238325002669003483672078334335498967
return zz/1
Heath:
/* Heath-Brown-Moroz constant */
procedure
/* Fast value */
zz = 0.00131764115485317810981735232251358595125073432329525167879254742178602344409610895090834
return zz/1
Hermite:
/* Second Hermite constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.15470053837925152901829756100391491129520350254025375203720465296795534460586669138743079117149905
else
/* Formula */
zz = 2/Sqrt(3)
return zz/1
Kepler:
/* Kepler-Bouwkamp constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.1149420448532962007010401574695987428307953372008635168440233965189660128253530511779407724849858370
else do
/* Many extra digits */
numeric digits Digits()*3
/* Asymptotic */
zz = 0; v = 0
do k = 1 to 100 until zz = v
v = zz
a = 2*k; b = Zeta(a); c = 2**a
zz = zz+((c-1)/a)*b*(b-1-1/c)
end
zz = Exp(-2*zz)
numeric digits Digits()/3
end
return zz/1
Khinchin:
/* Khinchin constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 2.685452001065306445309714835481795693820382293994462953051152345557218859537152002801141174931847698
else do
numeric digits Digits()+2
/* Asymptotic */
y1 = 0; v = 0
do n = 1 until y1 = v
numeric fuzz
v = y1
y2 = 0
do k = 1 to 2*n-1
y2 = y2+(-1)**(k+1)/k
end
y1 = y1+((Zeta(2*n)-1)/n)*y2
numeric fuzz 2
end
numeric fuzz
zz = Exp(y1/Ln(2))
numeric digits Digits()-2
end
return zz/1
Komornik:
/* Komornik-Loreti constant */
procedure
/* Fast value */
zz = 1.78723165018296593301327489033700838533793140296181099778478147050555749175060568699131001863407533302
return zz/1
Landau1:
/* Landau-Ramanujan constant */
procedure
/* Fast value */
zz = 0.7642236535892206629906987312500923281167905413934095147216866737496146416587328588384015050131312337
return zz/1
Landau2:
/* Landau constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.54325896534297670695272829530061323113886329375835698895573256910043468701530783731285212710225971085753+0
else
/* Formula */
zz = Pow(2,5/3)*Pi()**2/(3*Gamma(1/3)**3)
return zz/1
Laplace:
/* Laplace limit constant */
procedure
/* Fast value */
zz = 0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168
return zz/1
Lebesgue2:
/* Lebesgue constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.9894312738311469517416488090188667149242014060911110482841224326443725315845786546346329831401895579
else do
numeric digits Digits()+2
/* Series */
zz = 0; v = 0
do n = 0 until zz = v
numeric fuzz
v = zz
zz = zz+(Lambda1(2*n+2)-1)/(2*n+1)
numeric fuzz 2
end
numeric fuzz
a = Pi()**2
zz = 8*zz/a+4*(2*Ln(2)+Euler())/a
numeric digits Digits()-2
end
return zz/1
Lemniscate:
/* Lemniscate constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 2.6220575542921198104648395898911194136827549514316231628168217038007905870704142502302955329614291
else
/* Formula */
zz = Pi()/Agmean(1,Sqrt(2))
return zz/1
Levy1:
/* Levy-Khinchin constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 1.18656911041562545282172297594723712056835653647205433595425429865280963205625444330034830110848687594663
else
/* Formula */
zz = Pi()**2/(12*Ln(2))
return zz/1
Levy2:
/* Second Levy constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 3.27582291872181115978768188245384386360847552598237414940519892419072321564496035518127754047917452949269
else
/* Formula */
zz = Exp(Levy1())
return zz/1
Lieb:
/* Lieb constant */
procedure
if Digits() < 101 then
zz = 1.53960071783900203869106341467188654839360467005367166938293953729060712614115558851657438822866540060055
/* Fast value */
else
/* Formula */
zz = 8/(3*Sqrt(3))
return zz/1
Ln2:
/* Natural log of 2 constant */
procedure
/* Fast value */
zz = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875
return zz/1
Ln4:
/* Natural log of 4 constant */
procedure
/* Fast value */
zz = 1.386294361119890618834464242916353136151000268720510508241360018986787243939389431211726653992837375
return zz/1
Ln8:
/* Natural log of 8 constant */
procedure
/* Fast value */
zz = 2.079441541679835928251696364374529704226500403080765762362040028480180865909084146817589980989256063
return zz/1
Ln10:
/* Natural log of 10 constant */
procedure
/* Fast value */
zz = 2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959830
return zz/1
Loch:
/* Loch constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.97027011439203392574025601921001083378128470478516128661035052993125419989173704803621267490802902646924
else
/* Formula */
zz = 6*Ln(2)*Ln(10)/Pi()**2
return zz/1
Magicangle:
/* Magic angle constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.9553166181245092781638571025157577542434146950100054909596981293219120459039764553873916025856280734
else
/* Formula */
zz = Acos(1/Sqrt(3))
return zz/1
Meissel:
/* Meissel-Mertens constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.261497212847642783755426838608695859051566648261199206192064213924924510897368209714142631434246651051617
else do
numeric digits Digits()+2
zz = 0; v = 0
do n = 2 to 1000
t = Moebius(n)*Ln(Zeta(n))/n
if t <> 0 then do
zz = zz+t
if zz = v then
leave
v = zz
end
end
zz = Euler()+zz
numeric digits Digits()-2
end
return zz/1
Mills:
/* Mills constant */
procedure
/* Fast value */
zz = 1.30637788386308069046861449260260571291678458515671364436805375996643405376682659882150140370119739570729
return zz/1
Mrb:
/* MRB constant */
procedure
/* Fast value */
zz = 0.187859642462067120248517934054273230055903094900138786172004684089477231564660213703296654433107496903842
return zz/1
Nielsen:
/* Nielsen-Ramanujan constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.82246703342411321823620758332301259460947495060339921886777911468500373520160043691681445030987935
else
/* Formula */
zz = Pi()**2/12
return zz/1
Niven:
/* Niven constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 1.70521114010536776428855145343450816076202765165346909999428490654731319216812249193424413210087100179
else do
numeric digits Digits()+2
/* Series */
zz = 0; v = 0
do k = 2 to 100 until zz = v
numeric fuzz
v = zz
zz = zz+1-1/Zeta(k)
numeric fuzz 2
end
numeric fuzz
zz = zz+1
numeric digits Digits()-2
end
return zz/1
Omega:
/* Omega constant = LambertW0(1) = root of x*e^x-1 */
procedure expose glob.
/* Fast value or first guess */
zz = 0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194469617522946
if Digits() > 100 then do
numeric digits Digits()+3
/* Halley */
do until zz = v
numeric fuzz
v = zz
a = Exp(zz); b = a*zz-1
zz = zz-b/(a*(zz+1)-((zz+2)*b/(2*zz+2)))
numeric fuzz 3
end
numeric fuzz
numeric digits Digits()-3
end
return zz/1
Paperfolding:
/* Paper folding constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 0.850736188201867260367797760532066604411399493082710647728168261031830158459319445348545982642193923996091
else do
numeric digits Digits()+2
numeric fuzz 2
/* Definition */
zz = 0; v = 0
do n = 0 until zz = v
v = zz
zz = zz+8**(2**n)/(2**(2**(n+2))-1)
end
numeric fuzz
numeric digits Digits()-2
end
return zz/1
Parabolic:
/* Universal parabolic constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 2.29558714939263807403429804918949038759783220363858348392997534664410966268413312668409442623789762
else do
/* Formula */
a = Sqrt(2); zz = Ln(a+1)+a
end
return zz/1
Pi:
/* Pi constant */
procedure expose glob.
p = Digits()
/* In memory? */
if glob.pi.p <> '' then
return glob.pi.p
if p < 101 then
/* Fast value */
glob.pi.p = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211707+0
else do
numeric digits Digits()+2
if p < 201 then do
/* Chudnovsky */
zz = 0
do n = 0 until zz = v
numeric fuzz
v = zz; zz = zz + Fact(6*n)*(13591409+545140134*n)/(Fact(3*n)*Fact(n)**3*-640320**(3*n))
numeric fuzz 2
end
numeric fuzz
zz = 4270934400/(Sqrt(10005)*zz)
end
else do
/* Agmean */
zz = 0.25; a = 1; g = Sqrt(0.5); n = 1
do until a = v
numeric fuzz
v = a
x = (a+g)*0.5; g = Sqrt(a*g)
zz = zz-n*(x-a)**2; n = n+n; a = x
numeric fuzz 2
end
numeric fuzz
zz = a*a/zz
end
numeric digits Digits()-2
glob.pi.p = zz+0
end
return glob.pi.p
Plastic:
/* Plastic ratio constant = real root of x^3-x-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.32471795724474602596090885447809734073440405690173336453401505030282785124554759405469934798178728
else do
/* Formula */
a = Sqrt(69)/18; zz = Curt(0.5+a)+Curt(0.5-a)
end
return zz/1
Porter:
/* Porter constant */
procedure
/* Fast value */
zz = 1.46707807943397547289779848470722995344990332241488777739968581761660674432904480843036932751117401521266
return zz/1
Prouhet:
/* Prouhet-Thue-Morse constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 0.4124540336401075977833613682584552830894783744557695575733794153487935923657825889638045404862121334
else do
numeric digits Digits()+2
/* Infinite product */
zz = 1; v = 0
do n = 0
zz = zz*(1-1/2**(2**n))
if zz = v
then leave
v = zz
end
zz = (2-zz)*0.25
numeric digits Digits()-2
end
return zz/1
Ramanujan:
/* Ramanujan constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 262537412640768743.999999999999250072597198185688879353856337336990862707537410378210647910118607313
else
/* Formula */
zz = Exp(Pi()*Sqrt(163))
return zz/1
Recifibo:
/* Reciprocal Fibonacci constant */
procedure
if Digits() < 101 then
/* Fast value */
zz = 3.35988566624317755317201130291892717968890513373196848649555381532513031899668338361541621645679008729704
else do
numeric digits Digits()+2
numeric fuzz 2
/* Recurring */
zz = 1; v = 0; a = 0; b = 1
do n = 2 until zz = v
v = zz
c = a+b; zz = zz+1/c
a = b; b = c
end
numeric digits Digits()-2
end
return zz/1
Robbins:
/* Robbins constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.661707182267176235155831133248413581746400135790953604808944229479584646138597631306652480768107120151709
else
/* Formula */
zz = (4+17*Sqrt(2)-6*Sqrt(3)-7*Pi())/105+(Ln(1+Sqrt(2)))*0.2+(Ln(2+Sqrt(3)))*0.4
return zz/1
Salem:
/* Salem constant = largest positive root of x^10+x^9-x^7-x^6-x^5-x^4-x^3+x+1 */
procedure
/* Fast value or first guess */
zz = 1.176280818259917506544070338474035050693415806564695259830106347029688376548549962096830115581815395
if Digits() > 100 then do
/* Newton */
numeric digits Digits()+2; numeric fuzz 2
do until zz = v
v = zz
/* Powers */
y2 = zz*zz; y3 = y2*zz; y4 = y3*zz; y5 = y4*zz; y6 = y5*zz; y7 = y6*zz; y8 = y7*zz; y9 = y8*zz; y10 = y9*zz
/* Polynomial and derivative */
f = y10+y9-y7-y6-y5-y4-y3+zz+1
f1 = 10*y9+9*y8-7*y6-6*y5-5*y4-4*y3-3*y2+1
/* Next value */
zz = zz-f/f1
end
numeric digits Digits()-2
end
return zz/1
Sierpinski:
/* Sierpinski constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 2.584981759579253217065893587383171160088051651852630917321544987971932044001157120211117724527064
else
/* Formula */
zz = Pi()*Ln(Exp(2*Euler())/(2*Gauss()**2))
return zz/1
Silver:
/* Silver ratio constant = positive root of x^2-2x-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 2.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157
else
/* Formula */
zz = Sqrt(2)+1
return zz/1
Soldner:
/* Soldner constant */
procedure
/* Fast value */
zz = 1.45136923488338105028396848589202744949303228364801586309300455766242559575451783565953135771108682884
return zz/1
Somos:
/* Somos quadratic recurrence constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 1.66168794963359412129581892274995074996441863502506820818971116802560902982638372790836917641146116715528
else do
numeric digits Digits()+2
/* Definition */
zz = 1; v = 0; a = 2
do n = 1
zz = zz*Pow(n,1/a)
if zz = v
then leave
v = zz; a = 2*a
end
numeric digits Digits()-2
end
return zz/1
Sqrt2:
/* Square root of 2 constant */
procedure
/* Fast value */
zz = 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573
return zz/1
Sqrt3:
/* Square root of 3 constant */
procedure
/* Fast value */
zz = 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576
return zz/1
Sqrt5:
/* Square root of 5 constant */
procedure
/* Fast value */
zz = 2.236067977499789696409173668731276235440618359611525724270897245410520925637804899414414408378782275
return zz/1
Stephens:
/* Stephens constant */
procedure
/* Fast value */
zz = 0.57595996889294543964316337549249669250651396717649236360064079866537255169886852843640987209172618
return zz/1
Supergolden:
/* Super golden ratio constant = real solution of x^3-x^2-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.46557123187676802665673122521993910802557756847228570164318311124926299668501784047812580119490927
else do
/* Formula */
a = 3*Sqrt(93); zz = (1+Curt((29+a)*0.5)+Curt((29-a)*0.5))/3
end
return zz/1
Taniguchi:
/* Taniguchi constant */
procedure
/* Fast value */
zz = 0.678234491917391978035538279482894814096332239189440103036460415964983370740123233213762122933484616888328
return zz/1
Tau1:
/* Tau constant = 2*pi */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 6.28318530717958647692528676655900576839433879875021164194988918461563281257241799725606965068423414
else
/* Formula */
zz = 2*Pi()
return zz/1
Tetranacci2:
/* Tetranacci constant = root of x^4-x^3-x^2-x-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.92756197548292530426190586173662216869855425516338472714664703800966606229781555914981825346189065325
else do
/* Formula */
a = Curt((Sqrt(177304464)+7020)); zz = 1/4+Sqrt(11/48-a/72+7/a)+Sqrt(11/24+a/72-7/a+1/Sqrt(704/507-128*a/1521+7168/(169*a)))
end
return zz/1
Tribonacci2:
/* Tribonacci constant = root of x^3-x^2-x-1 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 1.839286755214161132551852564653286600424178746097592246778758639404203222081966425738435419428307014
else do
/* Formula */
a = 3*Sqrt(33); zz = (1+Curt(19+a)+Curt(19-a))/3
end
return zz/1
Twinprimes:
/* Twin primes constant */
procedure
/* Fast value */
zz = 0.660161815846869573927812110014555778432623360284733413319448423335405642304495277143760031413839867911779
return zz/1
Vanderpauw:
/* Van der Pauw constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 4.532360141827193809627682945716666810171861467723795584186016547940600953721305102259083879604016089653
else
/* Formula */
zz = Pi()/Ln(2)
return zz/1
Viswanath:
/* Viswanath constant */
procedure
/* Fast value */
zz = 1.1319882487943
return zz/1
Wallis:
/* Wallis constant = real root of x^3-2x-5 */
procedure
if Digits() < 101 then
/* Fast value */
zz = 2.09455148154232659148238654057930296385730610562823918030412852904531218998348366714626728177715776
else do
/* Formula */
a = Sqrt(1929); zz = Curt((45-a)/18)+Curt((45+a)/18)
end
return zz/1
Weierstrass:
/* Weierstrass constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
zz = 0.474949379987920650332504636327982968559549373217202982283331024864557929174883860274275641250502144418903
else
/* Formula */
zz = Pow(2,1.25)*Sqrt(Pi())*Exp(Pi()*0.125)/Gamma(0.25)**2
return zz/1
Zscore:
/* Z-score 97.5 percentile constant */
procedure
/* Fast value */
zz = 1.95996398454005423552459443052055152795555007786954839847695264636163527414488266779825470949281420601799
return zz/1
Functions
/* Functions library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */
/* Only real arguments and results supported */
Ackermann:
/* Ackermann's function */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
/* Fast values */
if xx = 0 then
return yy+1
if xx = 1 then
return yy+2
if xx = 2 then
return yy+yy+3
if xx = 3 then
return 2**(yy+3)-3
if xx = 4 then do
if yy = 0 then
return 2**(2**2)-3
if yy = 1 then
return 2**(2**(2**2))-3
if yy = 2 then
return 2**(2**(2**(2**2)))-3
end
/* Cf definition */
if yy = 0 then
return Ackermann(xx-1,1)
else
return Ackermann(xx-1,Ackermann(xx,yy-1))
Acos:
/* Arc cosine */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if Abs(xx) > 1 then say argument
/* Formula */
return Pi()/2-Asin(xx)
Acosh:
/* Arc cosine hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < 1 then say argument
/* Formula */
return Ln(xx+Sqrt(xx*xx-1))
Acot:
/* Arc cotangent */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return Acos(xx/Sqrt(xx*xx+1))
Acoth:
/* Arc cotangent hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < 1 then say argument
/* Formula */
return Ln((xx+1)/(xx-1))*0.5
Acsc:
/* Arc cosecant */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < 1 then say argument
/* Formulas */
if xx > 0 then
return Asin(1/xx)
else
return -(Pi()+Asin(1/xx))
Acsch:
/* Arc cosecant hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx = 0 then say argument
/* Formula */
zz = 1/xx
return Ln(zz+Sqrt(zz*zz+1))
Agmean:
/* Arithmetic geometric mean */
procedure
arg xx,yy
if \ Number(xx) then say argument
if xx <= 0 then say argument
if \ Number(yy) then say argument
if yy < 0 then say argument
/* Fast values */
if yy = 0 then
return 0
if xx = yy then
return xx
/* Series */
numeric digits Digits()+2
vx = xx+1
do until xx = vx
numeric fuzz 0
vx = xx; vy = yy
xx = (vx+vy)*0.5; yy = Sqrt(vx*vy)
numeric fuzz 2
end
numeric fuzz 0
numeric digits Digits()-2
return xx/1
Aliquot:
/* Aliquot function = Sum of all divisors of xx including 1 */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx = 1 then
return 1
/* Euclid's method */
return Sigma(xx)-xx
Alog:
/* Antilog = 10^x */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return Pow(10,xx)
Amean:
/* Arithmetic mean */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
/* Formula */
return (xx+yy)/2
Asec:
/* Arc secant */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < 1 then say argument
/* Formula */
return Sign(xx)*Acos(1/xx)
Asech:
/* Arc secant hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx <= 0 then say argument
if xx > 1 then say argument
/* Formula */
zz = 1/xx
return Ln(zz+Sqrt(zz*zz-1))
Asin:
/* Arc sine */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if Abs(xx) > 1 then say argument
/* Acos is faster for xx in [3/4,1) */
if Abs(xx) >= 0.75 then
return Sign(xx)*Acos(Sqrt(1-xx*xx))
/* Taylor series */
numeric digits Digits()+2
numeric fuzz 2
t = xx; zz = xx; xx = xx*xx
do n = 2 by 2 until zz = v
v = zz
t = t*xx*(n-1)/n; zz = zz+t/(n+1)
end
numeric digits Digits()-2
return zz/1
Asinh:
/* Arc sine hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return Ln(xx+Sqrt(xx*xx+1))
Atan:
/* Arc tangent */
procedure
arg xx
if \ Number(xx) then say argument
/* Twice argument reduction */
numeric digits Digits()+2
do 2
xx = xx/(1+Sqrt(xx*xx+1))
end
/* Taylor series */
numeric fuzz 2
t = xx; zz = xx; xx = xx*xx
do n = 3 by 2 until zz = v
v = zz
t = -t*xx; zz = zz+t/n
end
zz = 4*zz
numeric digits Digits()-2
return zz/1
Atanh:
/* Arc tangent hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if Abs(xx) >= 1 then say argument
/* Formula */
return Ln((1+xx)/(1-xx))*0.5
Atan2:
/* Arc tangent correct quadrant */
procedure expose glob.
arg yy,xx
if \ Number(xx) then say argument
if \ Number(yy) then say argument
if xx = 0 & yy = 0 then say argument
/* Formulas */
if xx > 0 then
return Atan(yy/xx)
if xx < 0 then
if yy < 0 then
return Atan(yy/xx)-Pi()
else
return Atan(yy/xx)+Pi()
if yy < 0 then
return -Pi()*0.5
else
return Pi()*0.5
Beta:
/* Beta */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
/* Formula */
return (Gamma(xx)*Gamma(yy))/Gamma(xx+yy)
Ceil:
/* Ceiling */
procedure
arg xx
if \ Number(xx) then say argument
/* Formulas */
if Whole(xx) then
return xx
else
return Trunc(xx)+(xx>=0)
Comb:
/* Combination = binomial coefficient */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
/* Fast values */
if yy = 0 then
return 1
if xx = yy then
return 1
if xx < yy then
return 0
numeric digits Digits()+2
/* Reduce yy */
if xx-yy < yy then
yy = xx-yy
/* Recurring loop */
zz = xx
do i = 1 to yy-1
zz = zz*(xx-i)/(i+1)
end
numeric digits Digits()-2
return zz/1
Cos:
/* Cosine */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Digits */
d1 = Digits(); d2 = Abs(Xpon(xx))+2; d3 = d1+d2
numeric digits d3
/* Reduce xx */
xx = xx//(2*Pi())
if Abs(xx) > Pi() then
xx = xx-Sign(xx)*2*Pi()
/* Taylor */
numeric fuzz d2
t = 1; zz = 1; xx = xx*xx
do n = 2 by 2 until zz = v
v = zz
t = -t*xx/(n*(n-1)); zz = zz+t
end
numeric digits d1
return zz/1
Cosh:
/* Cosine hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
zz = Exp(xx)
return (zz+1/zz)*0.5
Cot:
/* Cotangent */
procedure expose glob.
arg xx
zz = Sin(xx)
if zz = 0 then say argument
/* Formula */
return Cos(xx)/zz
Coth:
/* Cotangent hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx = 0 then say argument
/* Formula */
return 1/Tanh(xx)
Csc:
/* Cosecant */
procedure expose glob.
arg xx
zz = Sin(xx)
if zz = 0 then say argument
/* Formula */
return 1/zz
Csch:
/* Cosecant hyperbolic */
procedure expose glob.
arg xx
if xx = 0 then say argument
/* Formula */
return 1/Sinh(xx)
Curt:
/* Cubic root = x^(1/3) */
procedure
arg xx
if \ Number(xx) then say argument
/* Fast values */
if xx = 0 then
return 0
if xx = 1 then
return 1
/* Negative argument allowed */
if xx < 0 then do
s = -1; xx = Abs(xx)
end
else
s = 1
p = Digits()
numeric digits p+2
/* Argument reduction to [0...1000) */
i = Xpon(xx); i = (i-(2*(i<0)))%3; xx = xx/1000**i
/* First guess 1 decimal accurate */
t = '4.5 17.5 45.5 94.5 170.5 279.5 427.6 620.5 864.5 1000'
do zz = 1 until Word(t,zz) > xx
end
/* Dynamic precision */
d = Digits()
do n = 1 while d > 2
d.n = d; d = d%2+1
end
d.n = 2
/* Halley */
do k = n to 1 by -1
numeric digits d.k
zz = (xx+2*zz*zz*zz)/(3*zz*zz)
end
/* Inverse reduction */
zz = zz*10**i
numeric digits p
return s*zz
Deg:
/* Radians->degrees */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return xx*180/Pi()
Dfact:
/* Double factorial = n!! */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* In memory? */
if glob.dfact.xx <> '' then
return glob.dfact.xx
/* Fast values */
if xx < 2 then
return 1
/* Loops cf definition */
if xx//2 = 1 then do
zz = 1
do n = 3 to xx by 2
zz = zz*n
end
end
else do
zz = 2
do n = 4 to xx by 2
zz = zz*n
end
end
glob.dfact.xx = zz
return zz
Digamma:
/* Digamma */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx <= 0 & Whole(xx) then say argument
/* Fast formulas */
if xx = 0.125 then
return -0.5*Pi()-4*Ln(2)-( Pi()+Ln(Sqrt(2)+1)-Ln(Sqrt(2)-1))/Sqrt(2)-Euler()
if xx = 0.25 then
return -0.5*Pi()-3*Ln(2)-Euler()
if xx = 0.5 then
return -2*Ln(2)-Euler()
if xx = 1 then
return -Euler()
if xx > 0 & xx < 500 & Whole(xx) then
return Harmonic(xx-1)-Euler()
/* Fast series */
if xx > 0 & Half(xx) then do
xx = Trunc(xx)
zz = 0
do n = 1 to xx
zz = zz+2/(2*n-1)
end
zz = zz-2*Ln(2)-Euler()
return zz
end
numeric digits Digits()+2
/* Argument reduction */
p = Digits()
a = 0
do while xx < p
a = a+1/xx; xx = xx+1
end
/* Series using Zeta function */
zz = -0.5/xx; v = 0
do n = 2 by 2 until zz = v
numeric fuzz
v = zz
zz = zz+Zeta(1-n)/(xx**n)
numeric fuzz 2
end
/* Inverse reduction */
numeric fuzz
zz = Ln(xx)+zz-a
numeric digits Digits()-2
return zz/1
Digitprod:
/* Digital product = product(digits) */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Multiply digits */
zz = 1
do n = 1 to Length(xx) while zz <> 0
zz = zz*Substr(xx,n,1)
end
return zz
Digitroot:
/* Digital root */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
return 1+(xx-1)//9
Digitsum:
/* digitsum = sum(digits) */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Sum digits */
zz = 0
do n = 1 to Length(xx)
zz = zz+Substr(xx,n,1)
end
return zz
Divisor:
/* Divisor = sum(Divisors(x)^y) */
procedure expose glob. divi.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
if yy = '' then
yy = 1
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
/* Get divisors */
numeric digits Digits()+2
d = Divisors(xx)
/* Loop cf definition */
select
when yy = 0 then
zz = d
when yy = 1 then do
zz = 0
do n = 1 to d
zz = zz+divi.divisor.n
end
end
when Whole(yy) then do
zz = 0
do n = 1 to d
zz = zz+divi.divisor.n**yy
end
end
otherwise do
zz = 0
do n = 1 to d
zz = zz+Pow(divi.divisor.n,yy)
end
end
end
numeric digits Digits()-2
return zz/1
Dxgamma:
/* First derivative of Gamma */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx <= 0 & Whole(xx) then say argument
/* From derivative definition */
p = Digits()
numeric digits p*2+2
d = 10**(-p)
zz = (Gamma(xx+d)-Gamma(xx))/d
numeric digits p
return zz/1
Dxzeta:
/* First derivative of Zeta */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx = 1 then say argument
/* Fast values */
if xx < 0 then
if Whole(xx) then
if Even(xx) then
return (-1)**(-0.5*xx)*Fact(-xx)*Zeta(1-xx)/(2*(2*Pi())**-xx)
if xx = -2 then
return -Zeta(3)/(4*Pi()**2)
if xx = -1 then
return 1/12-Ln(Glaisher())
if xx = 0 then
return -0.5*Ln(2*Pi())
if xx = 0.5 then
return 0.25*(Pi()+2*Euler()+6*Ln(2)+2*Ln(Pi()))*Zeta(0.5)
if xx = 2 then
return Pi()*Pi()*(Euler()+Ln(2)-12*Ln(Glaisher())+Ln(Pi()))/6
if xx > 5 then do
/* Taylor series for the Zeta derivative */
numeric digits Digits()+2
zz = 0
a = Pow(2,xx); b = Ln(4); c = a-2
do n = 1 to 10000 until zz = v
numeric fuzz
v = zz
if n//2 = 0 then
t1 = 1
else
t1 = -1
zz = zz+(t1*a*Pow(n,-xx))*(b+c*Ln(n))/(c*c)
numeric fuzz 2
end
numeric fuzz
numeric digits Digits()-2
end
else do
/* From derivative definition */
p = Digits()
numeric digits p*4
d = 10**(-p+p%5); zz = (Zeta(xx+d)-Zeta(xx))/d
numeric digits p
end
return zz/1
Erf:
/* Error */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Asymptotic series */
numeric digits Digits()+2
a = (2*xx*Exp(-xx*xx))/Sqrt(Pi()); b = 2*xx*xx
numeric fuzz 2
zz = 0; f = 1; g = 1
do n = 0 until zz = v
v = zz
zz = zz+f/g
f = f*b; g = g*(2*n+3)
end
zz = a*zz
numeric digits Digits()-2
return zz/1
Erfc:
/* Error complementary */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return 1-Erf(xx)
Even:
/* Is a number even? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
/* Formula */
return \ Abs(xx//2)
Exp:
/* Exponential = e^xx */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Fast values */
numeric digits Digits()+2
if Whole(xx) then
return E()**xx
/* Argument reduction */
i = xx%1
if Abs(xx-i) > 0.5 then
i = i+Sign(xx)
/* Taylor series */
numeric fuzz 2
xx = xx-i; zz = 1; t = 1
do n = 1 until zz = v
v = zz
t = (t*xx)/n; zz = zz+t
end
/* Inverse reduction */
numeric fuzz
zz = zz*E()**i
numeric digits Digits()-2
return zz/1
Fact:
/* Integer factorial = n! */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Current in memory? */
if glob.fact.xx <> '' then
return glob.fact.xx
w = xx-1
/* Previous in memory? */
if glob.fact.w = '' then do
/* Loop cf definition */
zz = 1
do n = 2 to xx
zz = zz*n
end
glob.fact.xx = zz
end
else
/* Multiply */
glob.fact.xx = glob.fact.w*xx
return glob.fact.xx
Factorial:
/* Real factorial = x! */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < 0 & Whole(xx) then say argument
/* Fast values */
if Whole(xx) then
return Fact(xx)
/* Formula */
return Gamma(xx+1)
Ffact:
/* Falling factorial */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
/* In memory? */
if glob.ffact.xx.yy <> '' then
return glob.ffact.xx.yy
/* Fast values */
if xx < 2 then
return 1
if yy = 0 then
return 1
/* Calculate cf definition */
zz = 1
do k = 0 to yy-1
zz = zz*(xx-k)
end
glob.ffact.xx.yy = zz
return zz
Floor:
/* Floor */
procedure
arg xx
if \ Number(xx) then say argument
/* Formula */
if Whole(xx) then
return xx
else
return Trunc(xx)-(xx<0)
Frac:
/* Fractional part */
procedure
arg xx
if \ Number(xx) then say argument
/* Formula */
return xx-xx%1
Gamma:
/* Gamma */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < 1 & Whole(xx) then say argument
/* Fast values */
if xx < 0 then do
if Half(xx) then do
i = Abs(Floor(xx)); zz = (-1)**i*2**(2*i)*Fact(i)*Sqrt(Pi())/Fact(2*i)
return zz/1
end
end
if xx > 0 then do
if Whole(xx) then
return Fact(xx-1)
if Half(xx) then do
i = Floor(xx); zz = Fact(2*i)*Sqrt(Pi())/(2**(2*i)*Fact(i))
return zz/1
end
end
p = Digits()
if p < 61 then do
/* Lanczos with predefined coefficients */
/* Map negative xx to positive xx */
if xx < 0 then
return Pi()/(Gamma(1-xx)*Sin(Pi()*xx))
/* Argument reduction to (0.5...1.5) */
numeric digits p+2
c = Trunc(xx); xx = xx-c
if xx < 0.5 then do
xx = xx+1; c = c-1
end
/* Series coefficients 1/Gamma(xx) in 80 digits Fransen & Wrigge */
c.1 = 1.00000000000000000000000000000000000000000000000000000000000000000000000000000000
c.2 = 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467
c.3 = -0.65587807152025388107701951514539048127976638047858434729236244568387083835372210
c.4 = -0.04200263503409523552900393487542981871139450040110609352206581297618009687597599
c.5 = 0.16653861138229148950170079510210523571778150224717434057046890317899386605647425
c.6 = -0.04219773455554433674820830128918739130165268418982248637691887327545901118558900
c.7 = -0.00962197152787697356211492167234819897536294225211300210513886262731167351446074
c.8 = 0.00721894324666309954239501034044657270990480088023831800109478117362259497415854
c.9 = -0.00116516759185906511211397108401838866680933379538405744340750527562002584816653
c.10 = -0.00021524167411495097281572996305364780647824192337833875035026748908563946371678
c.11 = 0.00012805028238811618615319862632816432339489209969367721490054583804120355204347
c.12 = -0.00002013485478078823865568939142102181838229483329797911526116267090822918618897
c.13 = -0.00000125049348214267065734535947383309224232265562115395981534992315749121245561
c.14 = 0.00000113302723198169588237412962033074494332400483862107565429550539546040842730
c.15 = -0.00000020563384169776071034501541300205728365125790262933794534683172533245680371
c.16 = 0.00000000611609510448141581786249868285534286727586571971232086732402927723507435
c.17 = 0.00000000500200764446922293005566504805999130304461274249448171895337887737472132
c.18 = -0.00000000118127457048702014458812656543650557773875950493258759096189263169643391
c.19 = 0.00000000010434267116911005104915403323122501914007098231258121210871073927347588
c.20 = 0.00000000000778226343990507125404993731136077722606808618139293881943550732692987
c.21 = -0.00000000000369680561864220570818781587808576623657096345136099513648454655443000
c.22 = 0.00000000000051003702874544759790154813228632318027268860697076321173501048565735
c.23 = -0.00000000000002058326053566506783222429544855237419746091080810147188058196444349
c.24 = -0.00000000000000534812253942301798237001731872793994898971547812068211168095493211
c.25 = 0.00000000000000122677862823826079015889384662242242816545575045632136601135999606
c.26 = -0.00000000000000011812593016974587695137645868422978312115572918048478798375081233
c.27 = 0.00000000000000000118669225475160033257977724292867407108849407966482711074006109
c.28 = 0.00000000000000000141238065531803178155580394756670903708635075033452562564122263
c.29 = -0.00000000000000000022987456844353702065924785806336992602845059314190367014889830
c.30 = 0.00000000000000000001714406321927337433383963370267257066812656062517433174649858
c.31 = 0.00000000000000000000013373517304936931148647813951222680228750594717618947898583
c.32 = -0.00000000000000000000020542335517666727893250253513557337960820379352387364127301
c.33 = 0.00000000000000000000002736030048607999844831509904330982014865311695836363370165
c.34 = -0.00000000000000000000000173235644591051663905742845156477979906974910879499841377
c.35 = -0.00000000000000000000000002360619024499287287343450735427531007926413552145370486
c.36 = 0.00000000000000000000000001864982941717294430718413161878666898945868429073668232
c.37 = -0.00000000000000000000000000221809562420719720439971691362686037973177950067567580
c.38 = 0.00000000000000000000000000012977819749479936688244144863305941656194998646391332
c.39 = 0.00000000000000000000000000000118069747496652840622274541550997151855968463784158
c.40 = -0.00000000000000000000000000000112458434927708809029365467426143951211941179558301
c.41 = 0.00000000000000000000000000000012770851751408662039902066777511246477487720656005
c.42 = -0.00000000000000000000000000000000739145116961514082346128933010855282371056899245
c.43 = 0.00000000000000000000000000000000001134750257554215760954165259469306393008612196
c.44 = 0.00000000000000000000000000000000004639134641058722029944804907952228463057968680
c.45 = -0.00000000000000000000000000000000000534733681843919887507741819670989332090488591
c.46 = 0.00000000000000000000000000000000000032079959236133526228612372790827943910901464
c.47 = -0.00000000000000000000000000000000000000444582973655075688210159035212464363740144
c.48 = -0.00000000000000000000000000000000000000131117451888198871290105849438992219023663
c.49 = 0.00000000000000000000000000000000000000016470333525438138868182593279063941453996
c.50 = -0.00000000000000000000000000000000000000001056233178503581218600561071538285049997
c.51 = 0.00000000000000000000000000000000000000000026784429826430494783549630718908519485
c.52 = 0.00000000000000000000000000000000000000000002424715494851782689673032938370921241
/* Series expansion */
xx = xx-1; s = 0
do k = 52 by -1 to 1
s = s*xx+c.k
end
zz = 1/s
/* Undo reduction */
if c = -1 then
zz = zz/xx
else do
do i = 1 to c
zz = (xx+i)*zz
end
end
end
else do
xx = xx-1
/* Spouge */
/* Estimate digits and iterations */
q = Floor(p*1.5); a = Floor(p*1.1)
numeric digits q
/* Series */
s = 0
do k = 1 to a-1
s = s+((-1)**(k-1)*Pow(a-k,k-0.5)*Exp(a-k))/(Fact(k-1)*(xx+k))
end
s = s+Sqrt(2*Pi()); zz = Pow(xx+a,xx+0.5)*Exp(-a-xx)*s
end
/* Normalize */
numeric digits p
return zz/1
Gcd:
/* Greatest common divisor */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
/* Fast values */
xx = Abs(xx); yy = Abs(yy)
if yy = 0 then
return xx
/* Euclidian division */
do until yy = 0
parse value xx//yy yy with yy xx
end
return xx
Gmean:
/* Geometric mean */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
/* Formula */
return Sqrt(xx*yy)
Half:
/* Is a number half integer? */
procedure
arg xx
if \ Number(xx) then say argument
/* Formula */
return (Frac(Abs(xx))=0.5)
Isqrt:
/* Integer square root = Floor(xx^(1/2)) */
procedure
arg xx
if \ Number(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 1 then
return 0
if xx < 4 then
return 1
if xx < 9 then
return 2
/* Uses only integers */
xx = xx%1; q = 1
do until q > xx
q = q*4
end
zz = 0
do while q > 1
q = q%4; t = xx-zz-q; zz = zz%2
if t >= 0 then do
xx = t; zz = zz+q
end
end
return zz
Lambda1:
/* Dirichlet lambda */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return (1-Pow(2,-xx))*Zeta(xx)
Lambertw0:
/* Lambert W0 */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Fast values */
if xx = 0 then
return 0
numeric digits Digits()+2
/* Argument reduction */
e = E()
select
when xx > E() then do
a = Ln(xx); zz = a-Ln(a)
end
when xx > 0 then
zz = xx/E()
when xx >= -1/E() then do
a = E()*xx; b = a+1; c = Sqrt(b); zz = (a*Ln(1+c))/(b+c)
end
otherwise
say xx argument
end
/* Halley */
numeric fuzz 2
do until zz = v
v = zz
a = Exp(zz); b = a*zz; zz = zz-(b-xx)/(a+b-((zz+2)*(b-xx))/(2*zz+2))
end
numeric fuzz
numeric digits Digits()-2
return zz/1
Lambertwm1:
/* Lambert W-1 */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx < -1/E() then say argument
if xx >= 0 then say argument
/* Argument reduction */
numeric digits Digits()+2
select
when xx > -0.25 then do
a = Ln(-xx); zz = a-Ln(-a)
end
when xx > -1/E() then do
zz = -1-Sqrt(2)*Sqrt(E()*xx+1)
end
end
/* Halley */
numeric fuzz 2
do until zz = v
v = zz
a = Exp(zz); b = a*zz; zz = zz-(b-xx)/(a+b-((zz+2)*(b-xx))/(2*zz+2))
end
numeric fuzz
numeric digits Digits()-2
return zz/1
Lcm:
/* Least common multiple function */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
/* Fast values */
xx = Abs(xx); yy = Abs(yy)
if yy = 0 then
return 0
/* Euclid */
a = xx*yy
do until yy = 0
parse value xx//yy yy with yy xx
end
return a%xx
Ln:
/* Natural logarithm */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx <= 0 then say argument
/* Fast values */
if xx = 1 then
return 0
p = Digits()
/* In memory? */
if glob.ln.xx.p <> '' then
return glob.ln.xx.p
/* Precalculated values */
if xx = 2 & p < 101 then do
glob.ln.xx.p = Ln2()
return glob.ln.xx.p
end
if xx = 4 & p < 101 then do
glob.ln.xx.p = Ln4()
return glob.ln.xx.p
end
if xx = 8 & p < 101 then do
glob.ln.xx.p = Ln8()
return glob.ln.xx.p
end
if xx = 10 & p < 101 then do
glob.ln.xx.p = Ln10()
return glob.ln.xx.p
end
numeric digits p+2
/* Argument reduction */
yy = xx; i = 0; e = 1/E()
if yy < 0.5 then do
zz = 1/yy
do while zz > 1.5
zz = zz*e; i = i-1
end
yy = 1/zz
end
if yy > 1.5 then do
do while yy > 1.5
yy = yy*e; i = i+1
end
end
/* Taylor series */
q = (yy-1)/(yy+1); f = q; zz = q; v = q; q = q*q
do n = 3 by 2
f = f*q; zz = zz+f/n
if zz = v then
leave
v = zz
end
numeric digits p
/* Inverse reduction */
glob.ln.xx.p = 2*zz+i
return glob.ln.xx.p
Log:
/* Briggs logarithm */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return Ln(xx)/Ln(10)
Logxy:
/* Logarithm xx base yy */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
if yy = 1 then say argument
/* Formula */
return Ln(xx)/Ln(yy)
Mant:
/* Mantissa */
procedure
arg xx
if \ Number(xx) then say argument
/* Fast values */
if xx = 0 then
return 0
/* Formula */
zz = xx*1E+99999
return Left(zz,Pos('E',zz)-1)
Nroot:
/* Nth root = x^(1/y) */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 1 then say argument
if xx < 0 & Even(yy) then say argument
/* Fast values */
if xx = 0 then
return 0
if xx = 1 then
return 1
/* Formulas using faster methods */
if yy = 2 then
return Sqrt(xx)
if yy = 3 then
return Curt(xx)
if yy = 4 then
return Qurt(xx)
/* Calculate */
sx = Sign(xx); xx = Abs(xx)
p1 = Digits(); p2 = p1+2
numeric digits 3
/* First guess low accuracy */
zz = 1/Exp(Ln(xx)/yy)
numeric digits p2
/* Dynamic precision */
d = p2
do k = 1 while d > 4
d.k = d; d = d%2+1
end
d.k = 4
/* Halley */
a = 1/yy; b = yy+1
do j = k to 1 by -1
numeric digits d.j
zz = zz*a*(b-xx*zz**yy)
end
zz = 1/zz
if sx < 0 then
zz = -zz
numeric digits p1
return zz/1
Number:
/* Is xx a number? */
procedure
arg xx
return Datatype(xx,'n')
Odd:
/* Is a number odd? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
/* Formula */
return Abs(xx//2)
Perm:
/* Permutations */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if xx < 0 then say argument
if yy < 0 then say argument
if xx < yy then say argument
/* Cf definition */
zz = 1
do n = (xx-yy)+1 to xx
zz = zz*n
end
return zz
Pow:
/* Power = x^y */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Number(yy) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 0
if yy = 0 then
return 1
if Whole(yy) then
return xx**yy
/* Formulas */
if Abs(yy//1) = 0.5 then
return Sqrt(xx)**Sign(yy)*xx**(yy%1)
else
return Exp(yy*Ln(xx))
Powmod:
/* Power modulus = x^(y mod z) */
procedure
arg xx,yy,zz
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
if \ Number(zz) then say argument
if \ Whole(zz) then say argument
if zz < 1 then say argument
/* Fast values */
if zz = 1 then
return 0
/* Binary method */
numeric digits Max(Length(Format(xx,,,0)),Length(Format(yy,,,0)),2*Length(Format(zz,,,0)))
b = xx//zz; rr = 1
do while yy > 0
if yy//2 then
rr = rr*xx//zz
xx = xx*xx//zz; yy = yy%2
end
return rr
Prim:
/* Primorial = p# */
procedure expose glob. prim.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
/* Get primes up to the xth one */
p = Primes(Primecount(xx))
/* Calculate product */
zz = 1
numeric digits Digits()+2
do n = 1 to p
zz = zz*prim.prime.n
end
numeric digits Digits()-2
return zz/1
Primecount:
/* Prime count = Pi() function */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx = 1 then
return 0
/* Estimate number of primes up to x */
return Round(xx/(Ln(xx)-1))
Primeest:
/* Nth prime number estimate */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Estimate */
select
when xx <= 5 then do
zz = xx+xx+1
end
when xx <= 12 then do
a = Ln(xx); b = Ln(a); zz = xx*(a+b)
end
when xx <= 200 then do
a = Ln(xx); b = Ln(a); zz = xx*(a+b-1+1.8*b/a)
end
otherwise do
a = Ln(xx); b = Ln(a); zz = xx*(a+b-0.9385)+100
end
end
return Trunc(zz)
Primorial:
/* Primorial = p# */
procedure expose prim.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 2 then
return 1
/* Get primes up to xx */
p = Primes(xx)
/* Calculate product */
zz = 1
numeric digits Digits()+2
do n = 1 to p
zz = zz*prim.prime.n
end
numeric digits Digits()-2
return zz/1
Qurt:
/* Quartic root = x^(1/4) */
procedure
arg xx
if \ Number(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 0
if xx = 1 then
return 1
p = Digits()
numeric digits p+2
/* Argument reduction to [0...10000) */
i = Xpon(xx); i = (i+(i<0))%4-(i<0); xx = xx/10000**i
/* First guess 1 decimal accurate */
t = '8.5 48.5 168.5 440.5 960.5 1848.5 3248.5 5328.5 8280.5 10000'
do zz = 1 until Word(t,zz) > xx
end
/* Dynamic precision */
d = Digits()
do n = 1 while d > 2
d.n = d; d = Trunc(0.75*d)
end
d.n = 2
/* Halley */
f1 = 0.75; f2 = 0.25*xx
do k = n to 1 by -1
numeric digits d.k
zz = f1*zz+f2/(zz*zz*zz)
end
numeric digits p
/* Inverse reduction */
return zz*10**i
Rad:
/* Degrees->radians */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return xx*Pi()/180
Radical:
/* Radical = product of unique prime factors */
procedure expose ufac.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get unique factors */
n = Ufactors(xx)
/* Calculate product */
zz = 1
do i = 1 to n
zz = zz*ufac.factor.i
end
return zz
Rand:
/* Random number */
procedure expose glob.
arg xx,yy
if xx <> '' then do
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if xx >= yy then say argument
end
/* Fixed precision */
p = Digits(); numeric digits 30
/* Get and save first seed from system Date and Time */
if glob.rand = '' then do
a = Date('s'); b = Time('l')
glob.rand = Substr(b,10,3)||Substr(b,7,2)||Substr(b,4,2)||Substr(b,1,2)||Substr(a,7,2)||Substr(a,5,2)||Substr(a,3,2)
end
/* Calculate next random number cf HP calculators */
glob.rand = Right((glob.rand*2851130928467),15)
/* Uniform deviate [0,1) */
zz = '0.'Left(glob.rand,Length(glob.rand)-3)
if xx <> '' then
zz = Floor(zz*(yy+1-xx)+xx)
numeric digits 12
return zz/1
Repunit:
/* Repeated 1's */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Formula */
return Copies('1',xx)
Rfact:
/* Rising factorial */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx = 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
/* In memory? */
if glob.Rfact.xx.yy <> '' then
return glob.Rfact.xx.yy
/* Fast values */
if xx = 1 then
return 1
if yy = 0 then
return 1
/* Cf definition */
zz = 1
do k = 0 to yy-1
zz = zz*(xx+k)
end
glob.Rfact.xx.yy = zz
return zz
Round:
/* Round to specified decimals */
procedure
arg xx,yy
if yy = '' then
yy = Digits()
if \ Number(xx) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
/* Format */
if yy < 0 then do
yy = Abs(yy); a = 10**yy
zz = Format(xx/a,,0,0)*a
end
else
zz = Format(xx,,yy,0)
return zz/1
Sci:
/* Convert to scientific format */
procedure
arg xx
if \ Number(xx) then say argument
/* Force format */
return Format(xx,,,,0)
Sec:
/* Secant */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
zz = Cos(xx)
if zz = 0 then say argument
/* Formula */
return 1/zz
Sech:
/* Secant hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
return 1/Cosh(xx)
Sigma:
/* Sigma = Sum of all divisors of x including 1 and x */
procedure expose divi.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx = 1 then do
divi.0 = 1
return 1
end
/* Euclid's method */
m = xx//2; zz = 1+xx; n = 2
do j = 2+m by 1+m while j*j < xx
if xx//j = 0 then do
zz = zz+j+xx%j; n = n+2
end
end
if j*j = xx then do
zz = zz+j; n = n+1
end
/* Store number of divisors */
divi.0 = n
/* Return sum */
return zz
Sin:
/* Sine */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Digits */
d1 = Digits(); d2 = Abs(Xpon(xx))+2; d3 = d1+d2
numeric digits d3
/* Argument reduction */
p = Pi(); xx = xx//(2*p)
if Abs(xx) > p then
xx = xx-Sign(xx)*2*p
/* Taylor series */
t = xx; zz = xx; xx = xx*xx
numeric fuzz d2
do n = 2 by 2 until zz = v
v = zz
t = -t*xx/(n*(n+1)); zz = zz+t
end
numeric digits d1
return zz/1
Sinh:
/* Sine hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
zz = Exp(xx)
return (zz-1/zz)*0.5
Sqrt:
/* Square root = x^(1/2) */
procedure
arg xx
if \ Number(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 0
if xx = 1 then
return 1
p = Digits()
/* Predefined values */
if xx = 2 & p < 101 then
return Sqrt2()
if xx = 3 & p < 101 then
return Sqrt3()
if xx = 5 & p < 101 then
return Sqrt5()
numeric digits p+2
/* Argument reduction to [0...100) */
i = Xpon(xx); i = (i-(i<0))%2; xx = xx/100**i
/* First guess 1 decimal accurate */
t = '2.5 6.5 12.5 20.5 30.5 42.5 56.5 72.5 90.5 100'
do zz = 1 until Word(t,zz) > xx
end
/* Dynamic precision */
d = Digits()
do n = 1 while d > 2
d.n = d; d = d%2+1
end
d.n = 2
/* Heron */
do k = n to 1 by -1
numeric digits d.k
zz = (zz+xx/zz)*0.5
end
numeric digits p
return zz*10**i
Std:
/* Convert to standard format */
procedure
arg xx
if \ Number(xx) then say argument
/* Force format */
return Format(xx,,,0)
Tan:
/* Tangent */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
zz = Cos(xx)
if zz = 0 then say argument
/* Formula */
return Sin(xx)/zz
Tanh:
/* Tangent hyperbolic */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
/* Formula */
zz = Exp(2*xx)
return (zz-1)/(zz+1)
Tau2:
/* Ramanujan Tau */
procedure expose divi.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Niebur series */
numeric digits Digits()+2
s = 0; a = 18*xx*xx
do n = 1 to xx-1
b = n*n
s = s+b*(35*b-52*n*xx+a)*Sigma(n)*Sigma(xx-n)
end
zz = xx**4*Sigma(xx)-24*s
numeric digits Digits()-2
return zz/1
Totient:
/* Euler's totient function = Phi */
procedure expose fact.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 3 then
return 1
if xx < 5 then
return 2
/* Multiplicative property using Factors */
f = Factors(xx)+1; fact.factor.f = 0
zz = 1; v = fact.factor.1; n = 1
do i = 2 to f
a = fact.factor.i
if a = v then
n = n+1
else do
zz = zz*v**(n-1)*(v-1)
v = a; n = 1
end
end
return zz
Trigamma:
/* Trigamma */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx <= 0 & Whole(xx) then say argument
/* Fast values */
if xx = 1 then
return Zeta(2)
numeric digits Digits()+2
p = Digits()
/* Argument reduction */
a = 0
do while xx < p
a = a+1/(xx*xx); xx = xx+1
end
/* Series using Bernoulli numbers */
zz = 0; v = 0
do n = 2 by 2
zz = zz+Bernoulli(n)/xx**(n+1)
if zz = v then
leave
v = zz
end
/* Inverse reduction and result */
zz = zz+1/xx+1/(2*xx*xx)+a
numeric digits Digits()-2
return zz/1
Whole:
/* Is a number integer? */
procedure
arg xx
/* Formula */
return Datatype(xx,'w')
Xpon:
/* Exponent */
procedure
arg xx
if \ Number(xx) then say argument
/* Formula */
numeric digits 9
if xx = 0 then
return 0
else
return Right(xx*1E+99999999,9)-99999999
Zeta:
/* Zeta */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if xx = 1 then say argument
/* Fast values */
if xx < 0 & Whole(xx) then
if Even(xx) then
return 0
if xx = 0 then
return -0.5
if xx = 3 then
return Apery()
if xx < 0 & Whole(xx) then
if Odd(xx) then
return ((-1)**-xx*Bernoulli(1-xx))/(1-xx)
if xx > 0 & Whole(xx) then
if Even(xx) then
return ((-1)**(xx/2+1)*Bernoulli(xx)*(2*Pi())**xx)/(2*Fact(xx))
if xx < 0 then
return -2*Factorial(-xx)*Sin(Pi()*-xx/2)*Zeta(1-xx)/Pow(Pi()*2,1-xx)
/* Selected odd integers fast cf Plouffe */
p = Digits()
numeric digits p+2
if xx > 2 & xx < 22 & Whole(xx) then do
if Odd(xx) then do
select
when xx = 3 then do
a = 180; b = 7; c = 360; d = 0
end
when xx = 5 then do
a = 1470; b = 5; c = 3024; d = 84
end
when xx = 7 then do
a = 56700; b = 19; c = 113400; d = 0
end
when xx = 9 then do
a = 18523890; b = 625; c = 37122624; d = 74844
end
when xx = 11 then do
a = 425675250; b = 1453; c = 851350500; d = 0
end
when xx = 13 then do
a = 257432175; b = 89; c = 514926720; d = 62370
end
when xx = 15 then do
a = 390769879500; b = 13687; c = 781539759000; d = 0
end
when xx = 17 then do
a = 1904417007743250; b = 6758333; c = 3808863131673600; d = 29116187100
end
when xx = 19 then do
a = 21438612514068750; b = 7708537; c = 42877225028137500; d = 0
end
when xx = 21 then do
a = 1881063815762259253125; b = 68529640373; c = 3762129424572110592000; d = 1793047592085750
end
end
s1 = 0; s2 = 0; v1 = 0; v2 = 0
do n = 1 to 1000
m = n**xx; e = Exp(2*Pi()*n); q = m*e
s1 = s1+1/(q-m); s2 = s2+1/(q+m)
if s1 = v1 & s2 = v2 then
leave
v1 = s1; v2 = s2
end
zz = (b*Pi()**xx-c*s1-d*s2)/a
numeric digits p
return zz/1
end
end
/* Standard for all other xx */
/* Max iterations estimate */
select
when p < 15 then
n = p*2
when p < 89 then
n = p*3
otherwise
n = p*4
end
/* Argument reduction */
a = (-1)/(2**n*(1-Pow(2,1-xx)))
/* Series */
s = 0; zz = 0
do j = 0 to 2*n-1
b = (-1)**j; c = Pow(j+1,xx)
if j < n then
s = 0
else
s = s+Comb(n,j-n)
zz = zz+(s-2**n)*b/c
end
/* Inverse reduction */
zz = a*zz
numeric digits Digits()-2
return zz/1
Numbers
/* Numbers library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */
Abundant:
/* Is a number abundant? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 101 then do
a = '12 18 20 24 30 36 40 42 48 54 56 60 66 70 72 78 80 84 88 90 96 100'
if Wordpos(xx,a) = 0 then
return 0
else
return 1
end
/* Cf definition */
if Sigma(xx) > 2*xx then
return 1
else
return 0
Alcuin:
/* Alcuin numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 4 then
return 0
/* Formula */
if Even(xx) then
return Round(xx*xx/48)
else
return Round((xx+3)**2/48)
Almostprime:
/* Is a number almost k-prime? */
procedure expose glob. fact.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 1 then say argument
/* Fast values */
if yy = 1 then
return Prime(xx)
if yy = 2 then
return Semiprime(xx)
/* Cf definition */
if Factors(xx) = yy then
return 1
else
return 0
Arithmetic:
/* Is a number arithmetic? */
procedure expose divi.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 101 then do
a = '1 3 5 6 7 11 13 14 15 17 19 20 21 22 23 27 29 30 31 33 ',
|| '35 37 38 39 41 42 43 44 45 46 47 49 51 53 54 55 56 57 59 60 61 62 65 66 ',
|| '67 68 69 70 71 73 77 78 79 83 85 86 87 89 91 92 93 94 95 96 97 99 '
if Wordpos(xx,a) = 0 then
return 0
else
return 1
end
/* Cf definition */
s = Sigma(xx)
if Whole(s/divi.0) then
return 1
else
return 0
Bell:
/* Bell numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
if xx = 1 then
return 1
/* Recurring */
numeric digits Digits()+2
if xx < 36 then do
zz = 0
do k = 0 to xx
zz = zz+Stirling2(xx,k)
end
end
else do
zz. = 0; zz.0 = 1; zz.1 = 1
do n = 2 to xx
s = 0
do k = 0 to n-1
s = s+Comb(n-1,k)*zz.k
end
zz.n = s
end
zz = zz.xx
end
numeric digits Digits()-2
return zz/1
Beraha:
/* Beraha numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 1 then
return 4
if xx = 2 then
return 0
if xx = 3 then
return 1
if xx = 4 then
return 2
if xx = 6 then
return 3
if xx = 8 then
return 2+Sqrt(2)
if xx = 10 then
return 0.5*(5+Sqrt(5))
/* Formula */
return 2+2*Cos(2*Pi()/xx)
Bernoulli:
/* Bernoulli numbers */
procedure expose glob. work.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
if xx = 1 then
return 0.5
if xx//2 = 1 then
return 0
p = Digits()
/* In memory? */
if glob.bernoulli.xx.p <> '' then
return glob.bernoulli.xx.p
/* Take extra digits */
numeric digits p+xx
/* Recurring method Akiyama-Tanigawa */
do m = 0 by 1 to xx
work.m = 1/(m+1)
do j = m by -1 to 1
k = j-1; work.k = j*(work.k-work.j)
end
end
numeric digits p
/* Result */
zz = work.0+0
glob.bernoulli.xx.p = zz
return zz
Cake:
/* Cake numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
return (xx*xx*xx+5*xx+6)/6
Catalan2:
/* Catalan numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Recurring */
numeric digits Digits()+2
zz = 1
do n = 1 to xx
zz = zz*(4*n-2)/(n+1)
end
numeric digits Digits()-2
return zz/1
Caterer:
/* Lazy caterer's numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
return (xx*xx+xx+2)/2
Central:
/* Central binomial coefficient numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
return Comb(2*xx,xx)
Chowla:
/* Chowla numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Formula */
return Aliquot(xx)-1
Composite:
/* Is a number composite? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 2 then say argument
/* Fast values */
if xx < 101 then do
c = '4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 33 34 35 36 38 39 ',
|| '40 42 44 45 46 48 49 50 51 52 54 55 56 57 58 60 62 63 64 65 66 68 69 70 72 ',
|| '74 75 76 77 78 80 81 82 84 85 86 87 88 90 91 92 93 94 95 96 98 99 100 '
if Wordpos(xx,c) = 0 then
return 0
else
return 1
end
if xx//2 = 0 then
return 1
if xx//3 = 0 then
return 1
if Right(xx,1) = 5 then
return 1
/* Method trial division */
do n = 5 by 6 while n*n <= xx
if xx//n = 0 then
return 1
if xx//(n+2) = 0 then
return 1
end
return 0
Conway1:
/* Conway look-and-say numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Series cf definition */
zz = '1'
do n = 1 to xx
ns = ''; c = 1; m = Length(zz)
do j = 1 to m
if j = m | Substr(zz,j,1) <> Substr(zz,j+1,1) then do
ns = ns||c||Substr(zz,j,1); c = 1
end
else
c = c+1
end
zz = ns
end
return zz
Coprime:
/* Are 2 numbers coprime? */
procedure
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
/* Formula */
return (Gcd(xx,yy) = 1)
Cullen:
/* Cullen numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
return xx*2**xx+1
Deficient:
/* Is a number deficient? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 101 then do
d = '1 2 3 4 5 7 8 9 10 11 13 14 15 16 17 19 21 22 23 25 26 27 29 31 32 33 34 35 ',
|| '37 38 39 41 43 44 45 46 47 49 50 51 52 53 55 57 58 59 61 62 63 64 65 67 68 ',
|| '69 71 73 74 75 76 77 79 81 82 83 85 86 87 89 91 92 93 94 95 97 98 99 '
if Wordpos(xx,d) = 0 then
return 0
else
return 1
end
/* Cf definition */
if Sigma(xx) < 2*xx then
return 1
else
return 0
Derangement:
/* Derangement numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Recurring cf definition */
numeric digits Digits()+2
a = 1; b = 0; zz = 0
do n = 2 to xx
zz = (a+b)*(n-1)
a = b; b = zz
end
numeric digits Digits()-2
return zz/1
Favard:
/* Favard numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
if xx = 1 then
return Pi()*0.5
if xx = 2 then
return Pi()**2*0.125
if xx = 3 then
return Pi()**3/24
if xx = 4 then
return 5*Pi()**4/384
if xx = 5 then
return Pi()**5/240
if xx = 6 then
return 61*Pi()**6/46080
if xx = 7 then
return 17*Pi()**7/40320
if xx = 8 then
return 277*Pi()**8/2064384
if xx = 9 then
return 31*Pi()**9/725760
if xx = 10 then
return 50521*Pi()**10/3715891200
if xx = 11 then
return 691*Pi()**11/159667200
if xx = 12 then
return 540553*Pi()**12/392398110720
numeric digits Digits()+2
/* Recurring */
zz.0 = 1; zz.1 = Pi()*0.5
do j = 2 to xx
s = 0
do k = 1 to Floor(j*0.5)
a = 2*k; b = a-1; c = j-a
s = s+zz.b*zz.c
end
zz.j = s*Pi()/(2*j)
end
numeric digits Digits()-2
return zz.xx/1
Fermat:
/* Fermat numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
numeric digits Digits()+2
/* Formula */
zz = 2**(2**xx)
numeric digits Digits()-2
return zz/1
Ffactor:
/* First prime factor */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 2 then say argument
/* Fast values */
if xx < 4 then
return xx
/* Check low factors */
n = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i)
if xx//p = 0 then
return p
end
/* Check higher factors */
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then
if xx//j = 0 then
return j
if p = 3 then
iterate
k = j+2
if xx//k = 0 then
return k
end
/* Last factor */
if xx > 1 then
return xx
return 0
Fibonacci:
/* Fibonacci numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 2 then
return xx
/* Recurring */
numeric digits Digits()+2
a = 0; b = 1
do n = 2 to xx
zz = a+b; a = b; b = zz
end
numeric digits Digits()-2
return zz/1
Goldbach:
/* Goldbach numbers */
procedure expose prim.
arg xx
if \ Number(xx) then say argument
if \ Even(xx) then say argument
if xx < 4 then say argument
/* Scan primes */
p = Primes(xx); zz = 0
do i = 2 to xx%2
if prim.flag.i then do
j = xx-i
if prim.flag.j then do
zz = zz+1
end
end
end
return zz/1
Golomb2:
/* Golomb-Silverman numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 1 then
return 1
/* Recurring */
numeric digits Digits()+2
zz. = 0; zz.1 = 1
do n = 2 to xx
a = n-1; b = zz.a; c = zz.b; d = n-c
zz.n = 1+zz.d
end
numeric digits Digits()-2
return zz.xx/1
Gould:
/* Gould numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
/* Recurring */
numeric digits Digits()+2
zz. = 0; zz.0 = 1
do n = 1 to xx
if Even(n) then do
a = n%2; zz = zz.a
end
else do
a = n%2; zz = 2*zz.a
end
zz.n = zz
end
numeric digits Digits()-2
return zz.xx/1
Gregory1:
/* Gregory numbers 1 */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx = 1 then
return 0.5
/* Recurring */
numeric digits Digits()+2
zz.1 = 0.5
if xx = 1 then
return zz.1
do n = 2 to xx
s = 0
do k = 1 to n-1
s = s+zz.k/(n+1-k)
end
zz.n = -s+1/(n+1)
end
numeric digits Digits()-2
return (-1)**(xx-1)*zz.xx
Gregory2:
/* Gregory numbers 2 */
procedure
arg xx
if \ Number(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx = 1 then
return 0.25*Pi()
/* Formula */
return Atan(1/xx)
Harmonic:
/* Harmonic numbers */
procedure
arg xx,yy
if yy = '' then
yy = 1
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
/* Fast values */
if xx = 0 then
return 0
if xx < 2 then
return 1
numeric digits Digits()+2
/* From definition */
zz = 0
if yy = 1 then do
do n = 1 to xx
zz = zz+1/n
end
end
else do
do n = 1 to xx
zz = zz+1/n**yy
end
end
numeric digits Digits()-2
return zz/1
Highcomposite:
/* Is a number highly composite? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* ... to do */
return 0
Home:
/* Home prime numbers */
procedure expose fact. glob. work.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Run thru chain */
numeric digits 25
do while Factors(xx) > 1
xx = ''
do i = 1 to fact.0
xx = xx||fact.factor.i
end
end
return xx
Jacobsthal:
/* Jacobsthal numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 3 then
return 1
numeric digits Digits()+2
/* Formula */
zz = (2**xx-(-1)**xx)/3
numeric digits Digits()-2
return zz/1
Kaprekar:
/* Is a number Kaprekar? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 101 then do
k = '1 9 45 55 99'
if Wordpos(xx,k) = 0 then
return 0
else
return 1
end
/* Method casting out nines */
p = Max(Digits(),2*Length(xx))
numeric digits p
a = xx//9
if a > 2 then
return 0
s = xx*xx
if a = s//9 then do
do i = 1 to Length(s)%2
parse var s l +(i) r
if xx = l+r then
return 1
end
end
return 0
Lah:
/* Lah numbers */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 1 then say argument
if xx < yy then say argument
/* Fast values */
if xx = yy then
return 1
/* Formulas */
s = (-1)**xx
if yy = 1 then
return s*Fact(xx)
else
return s*Comb(xx-1,yy-1)*Fact(xx)/Fact(yy)
Lambda2:
/* Lambda numbers */
procedure expose efac.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 3 then
return 1
if xx < 5 then
return 2
/* Using Efactors and Lcm */
f = Efactors(xx)
zz = 1
do i = 1 to f
p = efac.factor.i
e = efac.exponent.i
if p = 2 & e > 2 then
zz = Lcm(zz,2**(e-2))
else
zz = Lcm(zz,(p-1)*p**(e-1))
end
return zz
Lebesgue1:
/* Lebesgue numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
if xx = 1 then
return 1/3+2*Sqrt(3)/Pi()
if xx = 2 then
return 0.2+Sqrt(25-2*Sqrt(5))/Pi()
if xx = 3 then
return 1/7+(22*Sin(Pi()/7)-2*Cos(Pi()/14)+10*Cos(3*Pi()/14))/(3*Pi())
if xx = 4 then
return 13/(2*Sqrt(3)*Pi())+1/9+(7*Sin(2*Pi()/9)-5*Sin(Pi()/9)-Cos(Pi()/18))/Pi()
numeric digits Digits()+2
/* Series expansion */
zz = 0
do n = 1 to xx
zz = zz+Tan(n*Pi()/(2*xx+1))/n
end
zz = 1/(2*xx+1)+2*zz/Pi()
numeric digits Digits()-2
return zz/1
Leonardo:
/* Leonardo numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 3 then
return 1
/* From definition */
a = 1; b = 1
do n = 3 to xx
zz = a+b+1
a = b; b = zz
end
return zz/1
Lucas:
/* Lucas numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 1 then
return 2
if xx = 2 then
return 1
/* Recurring */
numeric digits Digits()+2
a = 2; b = 1; zz = 0
do n = 2 to xx
zz = a+b
a = b; b = zz
end
numeric digits Digits()-2
return zz/1
Magic:
/* Magic numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
return (xx*xx*xx+xx)/2
Metallic:
/* Metallic ratio numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
if xx = 1 then
return Golden()
if xx = 2 then
return Silver()
if xx = 3 then
return Bronze()
/* Formula */
return 0.5*(xx+Sqrt(xx*xx+4))
Moebius:
/* Moebius numbers */
procedure expose glob. fact. ufac.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Using number of (unique) prime factors */
call Factors(xx)
call Ufactors(xx)
if fact.0 = ufac.0 then
if Even(fact.0) then
return 1
else
return -1
else
return 0
Motzkin:
/* Motzkin numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 3 then
return 1
numeric digits Digits()+2
/* Recurring */
a = 1; b = 1; zz = 0
do n = 2 to xx
zz = ((3*n-3)*a+(2*n+1)*b)/(n+2)
a = b; b = zz
end
numeric digits Digits()-2
return zz/1
Narayana:
/* Narayana cows numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 3 then
return 1
/* Recurring */
numeric digits Digits()+2
a = 1; b = 1; c = 1; zz = 0
do n = 3 to xx
zz = a+c
a = b; b = c; c = zz
end
numeric digits Digits()-2
return zz/1
Padovan:
/* Padovan numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 3 then
return 1
numeric digits Digits()+2
/* Recurring */
a = 1; b = 1; c = 1; zz = 0
do n = 3 to xx
zz = a+b
a = b; b = c; c = zz
end
numeric digits Digits()-2
return zz/1
Partition:
/* Partition numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
numeric digits Digits()+2
/* Recurring */
yy. = 0; yy.0 = 1
do n = 1 to xx
s = 0
do k = 1 to n
a = n-k*(3*k-1)/2; b = n-k*(3*k+1)/2
s = s+(-1)**(k+1)*(yy.a+yy.b)
end
yy.n = s
end
numeric digits Digits()-2
return yy.xx/1
Pell1:
/* Pell numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 2 then
return xx
numeric digits Digits()+2
/* Recurring */
a = 0; b = 1; zz = 0
do n = 2 to xx
zz = a+2*b
a = b; b = zz
end
numeric digits Digits()-2
return zz/1
Pell2:
/* Pell-Lucas companion numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 2 then
return 2
/* Recurring */
numeric digits Digits()+2
a = 2; b = 2; zz = 0
do n = 2 to xx
zz = a+2*b
a = b; b = zz
end
numeric digits Digits()-2
return zz/1
Perfect:
/* Is a number perfect? */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
p = '6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128 ',
|| '2658455991569831744654692615953842176 ',
|| '191561942608236107294793378084303638130997321548169216 ',
|| '13164036458569648337239753460458722910223472318386943117783728128 ',
|| '14474011154664524427946373126085988481573677491474835889066354349131199152128 '
if Wordpos(xx,p) = 0 then
return 0
else
return 1
Perrin:
/* Perrin numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 3
if xx = 1 then
return 0
if xx = 2 then
return 2
numeric digits Digits()+2
/* Recurring */
a = 3; b = 0; c = 2; zz = 0
do n = 3 to xx
zz = a+b
a = b; b = c; c = zz
end
numeric digits Digits()-2
return zz/1
Persistence:
/* Additive persistence */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast value */
if xx = 0 then
return 0
/* Cf definition */
do zz = 1 until xx < 10
xx = Digitsum(xx)
end
return zz
Pollardrho:
/* Find one prime factor */
procedure
arg xx
/* Init */
numeric digits Length(xx)+2
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 2 then say argument
/* Method Pollard Rho */
a = 2; b = 2; yy = 1
do while yy = 1
a = (a*a+1)//xx; b = (b*b+1)//xx; b = (b*b+1)//xx
yy = Gcd(Abs(a-b),xx)
end
if xx = yy then
return yy
else
return Min(yy,Trunc(xx/yy))
Practical:
/* Is a number practical? */
procedure expose fact.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 101 then do
p = '1 2 4 6 8 12 16 18 20 24 28 30 32 36 40 42 ',
|| '48 54 56 60 64 66 72 78 80 84 88 90 96 100 '
if Wordpos(xx,p) = 0 then
return 0
else
return 1
end
if Odd(xx) then
return 0
/* Method Srinivasan-Stewart-Sierpinski */
n = Factors(xx)
f = 2; s = 2
do i = 2 to n
if fact.factor.i > s+1
then return 0
f = f*fact.factor.i; s = Sigma(f)
end
return 1
Prime:
/* Is a number prime? */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Low primes also used as deterministic witnesses */
lp = '2 3 5 7 11 13 17 19 23 29 31 37 41'
/* Fast values */
if xx < 2 then
return 0
if xx < 42 then do
if Wordpos(xx,lp) = 0 then
return 0
else
return 1
end
if xx//2 = 0 then
return 0
if xx//3 = 0 then
return 0
if Right(xx,1) = 5 then
return 0
/* Miller-Rabin primality test */
numeric digits 2*Length(xx)
d = xx-1; e = d
/* Reduce xx-1 by factors of 2 */
do s = -1 while d//2 = 0
d = d%2
end
/* Thresholds deterministic witnesses */
th = '2047 1373653 25326001 3215031751 2152302898747 3474749660383 341550071728321 0 ',
|| '3825123056546413051 0 0 318665857834031151167461 3317044064679887385961981 '
w = Words(th)
/* Up to 13 deterministic trials */
if xx < Word(th,w) then do
do k = 1 to w
if xx < Word(th,k) then
leave
end
end
/* or 3 probabilistic trials */
else do
lp = ''
do k = 1 to 3
r = Rand(2,e)/1; lp = lp||r||' '
end
k = k-1
end
/* Algorithm using generated witnesses */
do i = 1 to k
a = Word(lp,i); zz = Powmod(a,d,xx)
if zz = 1 then
iterate
if zz = e then
iterate
do s
zz = (zz*zz)//xx
if zz = 1 then
return 0
if zz = e then
leave
end
if zz <> e then
return 0
end
return 1
Pronic:
/* Pronic numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Formula */
numeric digits Digits()+2
zz = xx*xx+xx
numeric digits Digits()-2
return zz/1
Recaman:
/* Recaman numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 0
/* Recurring cf definition */
zz. = 0; zz.0 = 0; s. = 0; s.0 = 1
do n = 1 to xx
m = n-1; a = zz.m-n
if a > 0 & s.a = 0 then do
zz.n = a; s.a = 1
end
else do
b = zz.m+n
zz.n = b; s.b = 1
end
end
return zz.xx
Semiprime:
/* Is a number semi prime? */
procedure expose prim.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Low semiprimes */
s = '4 6 9 10 14 15 21 22 25 26 33 34 35 38 39 46 49 ',
|| '51 55 57 58 62 65 69 74 77 82 85 86 87 91 93 94 95 '
/* Fast values */
if xx < 4 then
return 0
if xx < 101 then do
if Wordpos(xx,s) = 0 then
return 0
else
return 1
end
/* Wheeled scan */
do i = 2 for 2
if xx//i = 0 then
if Prime(xx%i) then
return 1
else
return 0
end
do i = 5 by 6 until i*i > xx
do j = i by 2 for 2
if xx//j = 0 then
if Prime(xx%j) then
return 1
else
return 0
end
end
return 0
Sorting:
/* Sorting numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 3
if xx = 1 then
return 0
if xx = 2 then
return 2
/* Recurring */
zz = 0; i = xx-1; z = 1
do while i >= 0
zz = zz+i; i = i-z; z = z+z
end
return zz/1
Stirling1:
/* Stirling numbers 1 */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
/* Fast values */
if yy = xx then
return 1
if yy > xx then
return 0
if yy = 0 then
return 0
/* Formulas */
s = (-1)**(xx-yy)
if yy = 1 then
return s*Fact(xx-1)
if yy = 2 then
return s*Round(Fact(xx-1)*Harmonic(xx-1))
if yy = 3 then
return s*Round(Fact(xx-1)*(Harmonic(xx-1)**2-Harmonic(xx-1,2))/2)
if yy = xx-3 then
return s*Comb(xx,2)*Comb(xx,4)
if yy = xx-2 then
return s*Round((3*xx-1)*Comb(xx,3)/4)
if yy = xx-1 then
return s*Comb(xx,2)
numeric digits Digits()+2
/* Definition */
zz. = 0; zz.1 = 1
do i = 1 to xx-1
a = i
do j = i+1 to 1 by -1
b = j-1; zz.j = zz.b+a*zz.j
end
end
numeric digits Digits()-2
return s*zz.yy
Stirling2:
/* Stirling numbers 2 */
procedure expose glob.
arg xx,yy
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
if \ Number(yy) then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
if xx < yy then say argument
/* Fast values */
if xx = 0 then
return 0
if yy = 0 then
return 0
if xx = yy then
return 1
if yy = 1 then
return 1
if yy = 2 then
return 2**(xx-1)-1
if yy = 3 then
return (3**xx-3*2**xx+3)/6
if yy = xx-1 then
return Comb(xx,2)
numeric digits Digits()+2
/* Definition */
zz. = 1; zz = xx-yy
do i = 2 to yy
do j = 1 to zz
a = j-1; zz.j = zz.j+i*zz.a
end
end
numeric digits Digits()-2
return zz.zz/1
Sylvester:
/* Sylvester numbers */
procedure
arg xx
if \ Whole(xx) then say argument
if \ Number(xx) then say argument
if xx < 1 then say argument
if xx > 33 then say argument
numeric digits Digits()+2
/* Definition */
zz = 2
do n = 2 to xx
zz = zz*zz-zz+1
end
numeric digits Digits()-2
return zz/1
Tetranacci1:
/* Tetranacci numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 4 then
return 0
/* Recurring */
a = 0; b = 0; c = 0; d = 1
do n = 3 to xx
zz = a+b+c; a = b; b = c; c = d; d = zz
end
return zz/1
Thue:
/* Thue-Morse numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Concatenate 2s complements */
zz = 0
do i = 1 to xx
zz = zz||Translate(zz,10,01)
end
return zz
Tribonacci1:
/* Tribonacci numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 3 then
return 0
/* Recurring */
a = 0; b = 0; c = 1
do n = 3 to xx
zz = a+b+c; a = b; b = c; c = zz
end
return zz/1
Ulam:
/* Ulam numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx < 3 then
return xx
/* Series */
numeric digits Digits()+2
zz.1 = 1; zz.2 = 2; d = 2; e.= 0; e.1 = 1; e.2 = 1; z = 3
do until d = xx
n = 0
do j = 1 to d
b = z-zz.j
if e.b then do
if zz.j <> b then do
n = n+1
if n > 2 then
leave j
end
end
end
if n = 2 then do
d = d+1; zz.d = z; e.z = 1
end
z = z+1
end
zz = zz.d
numeric digits Digits()-2
return zz/1
Wedderburn:
/* Wedderburn numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx < 2 then
return xx
/* Recurring */
numeric digits Digits()+2
zz. = 0; zz.1 = 1
do n = 2 to xx
if Even(n) then do
a = n%2; s = zz.a*(zz.a+1)/2
end
else do
a = n%2+1; s = 0
end
do i = 1 to a-1
j = n-i; s = s+zz.i*zz.j
end
zz.n = s
end
numeric digits Digits()-2
return zz.xx/1
Woodall:
/* Woodall numbers */
procedure
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Formula */
numeric digits Digits()+2
zz = xx*2**xx-1
numeric digits Digits()-2
return zz/1
ZigZag:
/* Zigzag alternating permutations numbers */
procedure expose glob.
arg xx
if \ Number(xx) then say argument
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Fast values */
if xx = 0 then
return 1
if xx = 1 then
return 1
/* Recurring */
numeric digits Digits()+2
zz. = 0; zz.0 = 1; zz.1 = 1
do n = 2 to xx
s = 0
do k = 0 to n-1
a = Comb(n-1,k); b = n-k-1
s = s+a*zz.k*zz.b
end
zz.n = s/2
end
zz = zz.xx
numeric digits Digits()-2
return zz/1
Polynomial
/* Polynomial library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2024 */
/* Coefficients are coded as string, highest first */
/* Example ax^3+cx = 'a 0 c 0' */
Padd:
/* Addition */
procedure expose poly.
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
/* Add */
wx = Words(xx); wy = Words(yy); zz = Max(wx,wy)
poly. = 0; poly.0 = zz
a = zz-wx; b = zz-wy
do i = 1 to wx
j = a+i; poly.coef.j = Word(xx,i)
end
do i = 1 to wy
j = b+i; poly.coef.j = poly.coef.j+Word(yy,i)
end
return zz
Parr2form:
/* Array -> Formula */
procedure expose poly.
/* Generate formula */
zz = ''; wm = poly.0
do i = 1 to wm
a = poly.coef.i
if a <> 0 then do
select
when a < 0 then
s = '-'
when i > 1 then
s = '+'
otherwise
s = ''
end
a = Abs(a); e = wm-i
if a = 1 & e > 0 then
a = ''
select
when e > 1 then
b = 'x^'e
when e = 1 then
b = 'x'
otherwise
b = ''
end
zz = zz||s||a||b
end
end
if zz = '' then
zz = 0
return Strip(zz)
Parr2lst:
/* Array -> List */
procedure expose poly.
zz = ''
do i = 1 to poly.0
zz = zz poly.coef.i
end
return Strip(zz)
Pdif:
/* Differentiation */
procedure expose poly.
arg xx
if xx = '' then say argument
/* dif */
wx = Words(xx); zz = wx-1
poly. = 0; poly.0 = zz
do i = 1 to zz
poly.coef.i = (zz-i+1)*Word(xx,i)
end
return zz
Pdiv:
/* Long division (in progress) */
procedure
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
/* Divide... to do */
return 0
Peval:
/* Evaluation */
procedure
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
if \ Number(yy) then say argument
/* Evaluate */
zz = 0
numeric digits Digits()+2
do n = 1 to Words(xx)
zz = zz*yy+Word(xx,n)
end
numeric digits Digits()-2
/* Normalize */
return zz/1
Plst2arr:
/* List -> Array */
procedure
arg xx
if xx = '' then say argument
/* Generate array */
zz = Words(xx)
poly. = 0; poly.0 = zz
do i = 1 to zz
poly.coef.i = Word(xx,i)
end
return zz
Plst2form:
/* List -> Formula */
procedure
arg xx
if xx = '' then say argument
/* Generate formula */
zz = ''; wm = Words(xx)
do i = 1 to wm
a = Word(xx,i)
if a <> 0 then do
select
when a < 0 then
s = '-'
when i > 1 then
s = '+'
otherwise
s = ''
end
a = Abs(a); e = wm-i
if a = 1 & e > 0 then
a = ''
select
when e > 1 then
b = 'x^'e
when e = 1 then
b = 'x'
otherwise
b = ''
end
zz = zz||s||a||b
end
end
if zz = '' then
zz = 0
return zz
Pmul:
/* Multiplication */
procedure expose poly. work.
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
/* Multiply */
wx = Words(xx); wy = Words(yy)
do i = 1 to wx
work.coef.1.i = Word(xx,i)
end
do i = 1 to wy
work.coef.2.i = Word(yy,i)
end
numeric digits Digits()+2
zz = wx+wy-1
poly. = 0; poly.0 = zz
do i = 1 to wx
do j = 1 to wy
k = i+j-1; poly.coef.k = poly.coef.k+work.coef.1.i*work.coef.2.j
end
end
numeric digits Digits()-2
/* Normalize coefs */
call Pnorm
/* Convert to list */
return zz
Pnorm:
/* Normalize coefs */
procedure expose poly.
do i = 1 to poly.0
poly.coef.i = poly.coef.i/1
end
return poly.0
Ppow:
/* Exponentiation */
procedure expose poly. comb. work.
arg xx,yy
if xx = '' then say argument
if \ Whole(yy) then say argument
if yy < 0 then say argument
/* Exponentiate */
numeric digits Digits()+2
wx = Words(xx); zz = wx*yy-yy+1
poly. = 0; poly.0 = zz
select
when wx = 1 then
/* Power of a number */
poly.coef.1 = xx**yy
when wx = 2 then do
/* Newton's binomial */
call Combinations -yy
a = Word(xx,1); b = Word(xx,2)
do i = 1 to zz
j = yy-i+1; k = i-1
poly.coef.i = comb.yy.k*a**j*b**k
end
end
otherwise do
/* Repeated multiplication */
do i = 1 to wx
poly.coef.1.i = Word(xx,i)
poly.coef.2.i = poly.coef.1.i
end
wy = wx
do i = 2 to yy
work. = 0
do j = 1 to wx
do k = 1 to wy
l = j+k-1; work.coef.l = work.coef.l+poly.coef.1.j*poly.coef.2.k
end
end
if i = yy then
leave i
wx = wx+wy-1
do j = 1 to wx
poly.coef.1.j = work.coef.j
end
end
do i = 1 to zz
poly.coef.i = work.coef.i
end
end
end
numeric digits Digits()-2
/* Normalize coefs */
call Pnorm
return zz
Psub:
/* Polynomial subtraction */
procedure expose poly.
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
/* Subtract */
wx = Words(xx); wy = Words(yy); zz = Max(wx,wy)
poly. = 0; poly.0 = zz
a = zz-wx; b = zz-wy
do i = 1 to wx
j = a+i; poly.coef.j = Word(xx,i)
end
do i = 1 to wy
j = b+i; poly.coef.j = poly.coef.j-Word(yy,i)
end
return zz
Rational
/* Rational library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2024 */
/* Rationals are coded as string 'num den' */
/* Example 123/456 = '123 456' */
Rabs:
/* Absolute value */
procedure
arg xx
parse var xx nx dx
if dx = '' then
dx = 1
nx = Abs(nx); dx = Abs(dx)
return nx dx
Radd:
/* Add */
procedure
arg xx,yy
parse var xx nx dx
if dx = '' then
dx = 1
parse var yy ny dy
if dy = '' then
dy = 1
return nx*dy+ny*dx dx*dy
Rdiv:
/* Divide */
procedure
arg xx,yy
parse var xx nx dx
if dx = '' then
dx = 1
parse var yy ny dy
if dy = '' then
dy = 1
return nx*dy dx*ny
Rinv:
/* Invert */
procedure
arg xx
parse var xx nx dx
if dx = '' then
dx = 1
return dx nx
Rmul:
/* Multiply */
procedure
arg xx,yy
parse var xx nx dx
if dx = '' then
dx = 1
parse var yy ny dy
if dy = '' then
dy = 1
return nx*ny dx*dy
Rneg:
/* Negate */
procedure
arg xx
parse var xx nx dx
if dx = '' then
dx = 1
return -nx dx
Rpow:
/* Power */
procedure
arg xx,yy
parse var xx nx dx
if dx = '' then
dx = 1
return nx**yy dx**yy
Rrat:
/* Float -> Rational */
procedure
arg xx
in = xx; h = 10**(Digits()-1)
/* Starting values */
num = 0; den = 1; tnum = 1; tden = 0
/* Limit precision */
do while tnum <= h & tden <= h
/* Integer part */
n = Trunc(xx)
/* Calculate numerator and denominator */
z = tnum; tnum = n*tnum+num; num = z
z = tden; tden = n*tden+den; den = z
/* Check precision */
if n = xx | tnum/tden = in then do
if tnum > h | tden > h then
leave
/* Possible final values */
num = tnum; den = tden
leave
end
/* Next round */
xx = 1/(xx-n)
end
return num den
Rsimp:
/* Simplify */
procedure
arg xx
parse var xx nx dx
if dx = '' then
dx = 1
g = Gcd(nx,dx)
return nx/g dx/g
Rsub:
/* Subtract */
procedure
arg xx,yy
parse var xx nx dx
if dx = '' then
dx = 1
parse var yy ny dy
if dy = '' then
dy = 1
return nx*dy-ny*dx dx*dy
Roots
/* Roots library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2024 */
/* Real coefficients and real roots */
/* Coefficients are coded as string, highest first */
/* Example ax^3+cx = 'a 0 c' */
Cueq:
/* Roots of cubic ax^3+bx^2+cx+d = 0 */
procedure expose cueq. glob.
arg xx; xx = xx 0 0 0 0
parse var xx a b c d .
if \ Number(a) then say argument
if \ Number(b) then say argument
if \ Number(c) then say argument
if \ Number(d) then say argument
if a = 0 then say argument
/* Method Cardano et al */
numeric digits Digits()*2+2
/* Reduce to x^3+ax^2+bx+c */
parse value b/a c/a d/a with a1 b1 c1
/* Reduce to x^3+ax+b */
h1 = a1*a1; h2 = h1*a1
a2 = b1-h1/3; b2 = -h2/27+h2/9-a1*b1/3+c1
/* Discriminant */
s = Round(b2*b2/4+a2*a2*a2/27,digits()-2)
/* Roots */
cueq. = 0
select
when s < 0 then do
h1 = Abs(a2)/3
f = Acos(-b2/(2*Sqrt(h1*h1*h1)))
h1 = 2*Sqrt(h1); h2 = f/3; h3 = Pi()/3; h4 = a1/3
x1 = h1*Cos(h2)-h4; x2 = -h1*Cos(h2-h3)-h4; x3 = -h1*Cos(h2+h3)-h4
zz = 1; cueq.root.zz = x1
if x2 <> x1 then do
zz = zz+1; cueq.root.zz = x2
end
if x3 <> x1 & x3 <> x2 then do
zz = zz+1; cueq.root.zz = x3
end
end
when s = 0 then do
h1 = Curt(-0.5*b2); h2 = a1/3
x1 = 2*h1-h2; x2 = -h1-h2
zz = 1; cueq.root.zz = 2*h1-h2
if x2 <> x1 then do
zz = zz+1; cueq.root.zz = -h1-h2
end
end
when s > 0 then do
h1 = 0.5*b2
f = Sqrt((h1)**2+(a2*a2*a2/27)); g = Curt(-h1+f); h = Curt(-h1-f)
h1 = g+h; h3 = a1/3
zz = 1; cueq.root.zz = h1-h3
end
end
/* Normalize */
numeric digits Digits()/2
do i = 1 to zz
cueq.root.i = Round(cueq.root.i,digits()-2)
end
/* Number of roots */
cueq.0 = zz
return zz
Lieq:
/* Roots of linear ax+b = 0 */
procedure expose glob.
arg xx; xx = xx 0 0
/* Coefficients */
parse var xx a b .
if \ Number(a) then say argument
if \ Number(b) then say argument
if a = 0 then say argument
/* Solution */
numeric digits Digits()*2
glob.Lieq.1 = -b/a; glob.lieq.0 = 1
numeric digits Digits()/2
/* Normalize */
glob.Lieq.1 = glob.lieq.1+0
/* Return number of roots */
return 1
Queq:
/* Roots of quartic ax^4+bx^3+cx^2+dx+e = 0 */
procedure expose glob.
arg xx; xx = xx 0 0 0 0 0
/* Coefficients */
parse var xx a b c d e .
if \ Number(a) then say argument
if \ Number(b) then say argument
if \ Number(c) then say argument
if \ Number(d) then say argument
if \ Number(e) then say argument
if a = 0 then say argument
/* Under construction */
numeric digits Digits()+2
zz = 0
/* ... */
glob.queq.0 = zz
numeric digits Digits()-2
/* Normalize */
do i = 1 to zz
glob.queq.i = glob.queq.i/1
end
/* Number of roots */
return zz
Sqeq:
/* Roots of square ax^2+bx+c = 0 */
procedure expose glob.
arg xx; xx = xx 0 0 0
parse var xx a b c .
if \ Number(a) then say argument
if \ Number(b) then say argument
if \ Number(c) then say argument
if a = 0 then say argument
/* Discriminant */
numeric digits Digits()+2
d = trunc(b*b-4*a*c,digits())
/* Abc formula */
h1 = -0.5/a
select
when d < 0 then do
zz = 0
end
when d = 0 then do
zz = 1
glob.sqeq.1 = h1*b
end
when d > 0 then do
zz = 2; h2 = Sqrt(d)
glob.sqeq.1 = h1*(b+h2)
glob.sqeq.2 = h1*(b-h2)
end
end
glob.sqeq.0 = zz
numeric digits Digits()-2
/* Normalize */
do i = 1 to zz
glob.sqeq.i = trunc(glob.sqeq.i,digits())/1
end
/* Number of roots */
return zz
Sequences
/* Sequences library - Build 5 Dec 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */
Abundants:
/* Abundant number sequence */
procedure expose abun.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
abun. = 0
/* Fast values */
if xx < 101 then do
a = '12 18 20 24 30 36 40 42 48 54 56 60 66 70 72 78 80 84 88 90 96 100 999'
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
abun.abundant.zz = w
end
zz = zz-1; abun.0 = zz
return zz
end
/* Check if sum(divisors) > 2*number */
zz = 0
do i = 1 to xx
s = Sigma(i)
if s > 2*i then do
zz = zz+1; abun.abundant.zz = i
end
end
abun.0 = zz
return zz
Additives:
/* Additive prime sequence */
procedure expose addi. prim.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
addi. = 0
/* Fast values */
if xx < 2 then
return 0
if xx < 101 then do
a = '2 3 5 7 11 23 29 41 43 47 61 67 83 89 999'
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
addi.prime.zz = w
end
zz = zz-1; addi.0 = zz
return zz
end
/* Get primes */
p = Primes(xx)
/* Collect additive primes */
zz = 0
do i = 1 to p
q = prim.prime.i; s = 0
do j = 1 to Length(q)
s = s+Substr(q,j,1)
end
if Prime(s) then do
zz = zz+1; addi.prime.zz = q
end
end
/* Return number of additive primes */
return zz
Amicables:
/* Amicable number sequence */
procedure expose glob. amic.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
amic. = 0
/* Scan for amicable pairs */
zz = 0
do i = 1 to xx
s = Sigma(i)-i; glob.Sigma.i = s
if i = glob.Sigma.s then do
if s = i then
iterate
zz = zz+1
amic.amicable.1.zz = s; amic.amicable.2.zz = i
end
end
amic.0 = zz
return zz
Arithmetics:
/* Arithmetic number sequence */
procedure expose arit. divi.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
arit. = 0
/* Fast values */
if xx < 101 then do
a = '1 3 5 6 7 11 13 14 15 17 19 20 21 22 23 27 29 30 31 33 ',
|| '35 37 38 39 41 42 43 44 45 46 47 49 51 53 54 55 56 57 59 60 61 62 65 66 ',
|| '67 68 69 70 71 73 77 78 79 83 85 86 87 89 91 92 93 94 95 96 97 99 999 '
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
arit.arithmetic.zz = w
end
zz = zz-1; arit.0 = zz
return zz
end
/* Check if number is Arithmetic */
zz = 0
do i = 1 to xx
if Arithmetic(i) then do
zz = zz+1; arit.arithmetic.zz = i
end
end
arit.0 = zz
return zz
Carmichaels:
/* Carmichael 3 strong prime sequence */
procedure expose carm.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
carm. = 0
/* Method Janeson */
zz = 0
do p1 = 3 by 2 to xx
if \ Prime(p1) then
iterate p1
pm = p1-1; ps = -p1*p1
do h3 = 1 to pm
t1 = (h3+p1) * pm; t2 = ps//h3
if t2 < 0 then
t2 = t2+h3
do d = 1 to h3+pm
if t1//d <> 0 then
iterate d
if t2 <> d//h3 then
iterate d
p2 = 1+t1%d
if \ Prime(p2) then
iterate d
p3 = 1+(p1*p2%h3)
if \ Prime(p3) then
iterate d
if (p2*p3)//pm <> 1 then
iterate d
zz = zz+1
carm.prime1.zz = p1; carm.prime2.zz = p2; carm.prime3.zz = p3
end d
end h3
end
carm.0 = zz
/* Return count */
return zz
Chowlas:
/* Chowla number sequence */
procedure expose chow.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
chow. = 0
/* Call aliquot */
do i = 1 to xx
chow.chowla.i = Aliquot(i)-1
end
chow.0 = xx
/* Return count */
return xx
Combinations:
/* Combination sequence */
procedure expose comb.
arg xx
if \ Whole(xx) then say argument
/* Recurring definition */
comb. = 1
if xx < 0 then do
xx = -xx; m = xx%2; a = 1
do i = 1 to m
a = a*(xx-i+1)/i; comb.xx.i = a
end
do i = m+1 to xx-1
j = xx-i; comb.xx.i = comb.xx.j
end
zz = xx+1
end
else do
do i = 1 to xx
i1 = i-1
do j = 1 to i1
j1 = j-1; comb.i.j = comb.i1.j1+comb.i1.j
end
end
zz = (xx*xx+3*xx+2)/2
end
comb.0 = zz
return zz
Composites:
/* Composite number sequence = not prime */
procedure expose comp.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
comp. = 0
/* Fast values */
if xx < 101 then do
c = '4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 33 34 35 36 38 39 ',
|| '40 42 44 45 46 48 49 50 51 52 54 55 56 57 58 60 62 63 64 65 66 68 69 70 72 ',
|| '74 75 76 77 78 80 81 82 84 85 86 87 88 90 91 92 93 94 95 96 98 99 100 999 '
do zz = 1 to Words(c)
w = Word(c,zz)
if w > xx then
leave
comp.composite.zz = w
end
zz = zz-1; comp.0 = zz
return zz
end
/* Sieve of Eratosthenes */
zz. = 0
do i = 2 to xx while i*i <= xx
do j = i+i by i to xx
zz.j = 1
end
end
zz = 0
do i = 4 to xx
if zz.i = 1 then do
zz = zz+1
comp.composite.zz = i
end
end
comp.0 = zz
return zz
Cubans:
/* Cuban prime sequence */
procedure expose cuba.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
cuba. = 0; zz = 0
/* Generate */
a = 1; b = 8; c = 2
do forever
v = b-a
if v > xx then
leave
if Prime(v) then do
zz = zz+1; cuba.prime.zz = v
end
a = b; c = c+1; b = c*c*c
end
cuba.0 = zz
return zz
Cubics:
/* Cubic prime sequence */
procedure expose cubi. prim.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get primes */
call Primes(xx)
/* Init */
cubi. = 0; cubi.prime.1 = 2; zz = 1
/* Generate */
a = 2; b = 3; k = 1
do while b <= xx
b = a + k*k*k
if prim.flag.b = 1 then do
zz = zz+1; cubi.prime.zz = b
a = b; k = 1
end
k = k+1
end
cubi.0 = zz
return zz
Cyclops:
/* Cyclop number sequence */
procedure expose cycl.
arg xx
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Init */
cycl. = 0; m = 0; zz = 1
/* Generate */
do f = 1
h = m+1; m = zz
do i = 1 to 9
do j = h to m
a = cycl.cyclop.j
do k = 1 to 9
b = i||a||k
if b > xx then
leave f
zz = zz+1; cycl.cyclop.zz = b
end
end
end
end
cycl.0 = zz
return zz
Deficients:
/* Deficient number sequence */
procedure expose defi.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
defi. = 0
/* Fast values */
if xx < 101 then do
a = '1 2 3 4 5 7 8 9 10 11 13 14 15 16 17 19 21 22 23 25 26 27 29 31 32 33 34 35 ',
|| '37 38 39 41 43 44 45 46 47 49 50 51 52 53 55 57 58 59 61 62 63 64 65 67 68 ',
|| '69 71 73 74 75 76 77 79 81 82 83 85 86 87 89 91 92 93 94 95 97 98 99 999'
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
defi.deficient.zz = w
end
zz = zz-1; defi.0 = zz
return zz
end
/* Check if sum(divisors) < 2*number */
zz = 0
do i = 1 to xx
s = Sigma(i)
if s < 2*i then do
zz = zz+1
defi.deficient.zz = i
end
end
defi.0 = zz
return zz
Divisors:
/* Divisors of an integer */
procedure expose divi.
arg xx
if \ Whole(xx) then say argument
if xx < 0 then say argument
/* Init */
divi. = 0
/* Fast values */
if xx = 1 then do
divi.divisor.1 = 1; divi.0 = 1
return 1
end
/* Euclid's method */
a = 1; divi.Left.1 = 1; b = 1; divi.Right.1 = xx
m = xx//2
do j = 2+m by 1+m while j*j < xx
if xx//j = 0 then do
a = a+1; divi.Left.a = j; b = b+1; divi.Right.b = xx%j
end
end
if j*j = xx then do
a = a+1; divi.Left.a = j
end
/* Save in table */
zz = 0
do i = 1 to a
zz = zz+1; divi.divisor.zz = divi.Left.i
end
do i = b by -1 to 1
zz = zz+1; divi.divisor.zz = divi.Right.i
end
divi.0 = zz
/* Return number of divisors */
return zz
Efactors:
/* Prime factors with exponent */
procedure expose efac.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Fast values */
if xx = 1 then do
efac.0 = 0
return 0
end
if xx < 4 then do
efac.factor.1 = xx; efac.exponent.1 = 1; efac.0 = 1
return 1
end
if Prime(xx) then do
efac.factor.1 = xx; efac.exponent.1 = 1; efac.0 = 1
return 1
end
/* Check low factors */
zz = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i); e = 0
do while xx//p = 0
e = e+1; xx = xx%p
end
if e > 0 then do
zz = zz+1; efac.factor.zz = p; efac.exponent.zz = e
end
end
/* Check higher factors */
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then do
e = 0
do while xx//j = 0
e = e+1; xx = xx%j
end
end
if e > 0 then do
zz = zz+1; efac.factor.zz = j; efac.exponent.zz = e
end
if p = 3 then
iterate
yy = j+2; e = 0
do while xx//yy = 0
e = e+1; xx = xx%yy
end
if e > 0 then do
zz = zz+1; efac.factor.zz = yy; efac.exponent.zz = e
end
end
if xx > 1 then do
if e = 0 then do
zz = zz+1; efac.factor.zz = xx; efac.exponent.zz = 1
end
end
efac.0 = zz
/* Return number of factor pairs */
return zz
Emirps:
/* Reversed primes sequence */
procedure expose emir. prim.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get primes */
call Primes(xx)
/* Init */
emir. = 0; zz = 0
/* Generate */
do i = 13 by 2 to xx
a = Reverse(i)
if a = i then
iterate i
if \ prim.flag.i then
iterate i
if \ prim.flag.a then
iterate i
zz = zz+1; emir.prime.zz = i
end
emir.0 = zz
return zz
Erdoss:
/* Erdos prime sequence */
procedure expose erdo. prim. glob.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get primes */
call Primes(xx)
/* Init */
erdo. = 0; zz = 0
/* Generate */
do i = 1 to prim.0
a = prim.prime.i
do j = 1
k = Fact(j)
if k >= a then
leave j
p = a-k
if prim.flag.p then
iterate i
end
zz = zz+1; erdo.prime.zz = a
end
erdo.0 = zz
return zz
Euclids:
/* Euclid-Mullin number sequence */
procedure expose eucl.
arg xx
/* Init */
numeric digits 100
if \ Whole(xx) then say argument
if xx < 1 then say argument
if xx > 16 then say argument
/* Init */
eucl. = 0; eucl.euclid.1 = 2; eucl.0 = 1
/* Using Ffactor = first prime factor */
p = 2
do i = 2 to xx
z = p+1; t = Ffactor(z); eucl.euclid.i = t; p = p*t
end
/* Save and return number */
eucl.0 = xx
return xx
Extras:
/* Extra prime sequence */
procedure expose extr.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
extr. = 0; extr.prime.1 = 2; zz = 1
/* Generate */
do i = 3 by 2 to xx
do j = 1 to Length(i)
if Pos(Substr(i,j,1),'014689') > 0 then
iterate i
end
if \ Prime(Digitsum(i)) then
iterate i
if \ Prime(i) then
iterate i
zz = zz+1; extr.prime.zz = i
end
extr.0 = zz
return zz
Factors:
/* Prime factors */
procedure expose fact. glob.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
fact. = 0
/* Fast values */
if xx = 1 then
return 0
if Prime(xx) then do
fact.factor.1 = xx; fact.0 = 1
return 1
end
/* Check low factors */
zz = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i)
do while xx//p = 0
zz = zz+1; fact.factor.zz = p
xx = xx%p
end
end
/* Check higher factors */
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then do
do while xx//j = 0
zz = zz+1; fact.factor.zz = j
xx = xx%j
end
end
if p = 3 then
iterate
k = j+2
do while xx//k = 0
zz = zz+1; fact.factor.zz = k
xx = xx%k
end
end
/* Last factor */
if xx > 1 then do
zz = zz+1; fact.factor.zz = xx
end
fact.0 = zz
/* Return number of factors */
return zz
Fareys:
/* Farey fraction sequence */
procedure expose fare.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
fare. = 0
/* Recurring */
fare.farey.1 = '0/1'
zz = 1
parse value 0 1 1 xx with aa bb cc dd
do while cc <= xx
k = Trunc((xx+bb)/dd)
parse value cc dd k*cc-aa k*dd-bb with aa bb cc dd
zz = zz+1
fare.farey.yy = aa'/'bb
end
fare.0 = zz
/* Number of entries */
return zz
Frobeniuss:
/* Frobenius prime sequence */
procedure expose frob. prim.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get primes */
p = Primes(xx%10)
/* Init */
frob. = 0; zz = 0
/* Generate */
do i = 1 to p-1
j = i+1; f = prim.prime.i * prim.prime.j - prim.prime.i - prim.prime.j
if f > xx then
leave
zz = zz+1; frob.prime.zz = f
end
frob.0 = zz
return zz
Giugas:
/* Giuga number sequence */
procedure expose giug.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
giug. = 0
/* Fast values */
if xx < 30 then
return 0
g = '30 858 1722 66198 2214408306 24423128562 432749205173838 14737133470010574 ',
|| '550843391309130318 244197000982499715087866346 554079914617070801288578559178 ',
|| '1910667181420507984555759916338506 999999999999999999999999999999999 '
do zz = 1 to Words(g)
w = Word(g,zz)
if w > xx then do
zz = zz-1
leave
end
giug.giuga.zz = w
end
giug.0 = zz
return zz
Hammings:
/* Hamming number sequence */
procedure expose hamm.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Ensure enough digits */
numeric digits 2**(Length(xx)+1)
/* Dijkstra */
hamm. = 1
x2 = 2; x3 = 3; x5 = 5
i2 = 1; i3 = 1; i5 = 1
do zz = 2
h = x2
if x3 < h then
h = x3
if x5 < h then
h = x5
hamm.hamming.zz = h
if zz = xx then
leave
if x2 = h then do
i2 = i2+1; x2 = hamm.hamming.i2+hamm.hamming.i2
end
if x3 = h then do
i3 = i3+1; x3 = hamm.hamming.i3+hamm.hamming.i3+hamm.hamming.i3
end
if x5 = h then do
i5 = i5+1; x5 = hamm.hamming.i5*5
end
end
hamm.0 = zz
return zz
Highcomposites:
/* Highly composite number sequence */
procedure expose high.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
high. = 0
/* Thresholds and increments */
a = '1 2 6 60 840 55440 720720 61261200 2327925600 321253732800 9999999999999'
b = '1 2 6 60 420 27720 360360 12252240 232792560 80313433200 9999999999999'
c = Words(a)-1; m = 0; zz = 0
/* Colllect cf definition */
do i = 1 to c
do j = Word(a,i) by Word(b,i) to Min(xx,Word(a,i+1)-1)
d = Divisors(j)
if d > m then do
zz = zz+1; high.highcomposite.zz = j
m = d
end
end
end
high.0 = zz
/* Return count */
return zz
Humbles:
/* Humble number sequence */
procedure expose humb.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Ensure enough digits */
numeric digits 2**(Length(xx)+1)
/* Dijkstra */
humb. = 1
x2 = 2; x3 = 3; x5 = 5; x7 = 7
i2 = 1; i3 = 1; i5 = 1; i7 = 1
do zz = 2
h = x2
if x3 < h then
h = x3
if x5 < h then
h = x5
if x7 < h then
h = x7
humb.humble.zz = h
if zz = xx then
leave
if x2 = h then do
i2 = i2+1; x2 = humb.humble.i2+humb.humble.i2
end
if x3 = h then do
i3 = i3+1; x3 = humb.humble.i3+humb.humble.i3+humb.humble.i3
end
if x5 = h then do
i5 = i5+1; x5 = humb.humble.i5*5
end
if x7 = h then do
i7 = i7+1; x7 = humb.humble.i7*7
end
end
humb.0 = zz
return zz
Kaprekars:
/* Kaprekar number sequence */
procedure expose kapr.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
kapr. = 0
/* Fast values */
if xx < 101 then do
k = '1 9 45 55'
do zz = 1 to Words(k)
w = Word(k,zz)
if w > xx then
leave
kapr.kaprekar.zz = w
end
zz = zz-1; kapr.0 = zz
return zz
end
xx = xx+1
/* Ensure enough digits */
p = Max(Digits(),2*Length(xx))
numeric digits p
zz = 1; kape.kaprekar.1 = 1
/* Method casting out nines */
do i = 9 for xx-9
a = i//9
if a > 2 then
iterate i
s = i*i
if a = s//9 then do
do j = 1 to Length(s)%2
parse var s l +(j) r
if i <> l+r then
iterate j
zz = zz+1; kapr.kaprekar.zz = i
leave
end
end
end
kapr.0 = zz
return zz
Perfects:
/* Perfect number sequence */
procedure expose perf.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
perf. = 0
/* Fast values */
if xx < 6 then
return 0
p = '6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128 ',
|| '2658455991569831744654692615953842176 ',
|| '191561942608236107294793378084303638130997321548169216 ',
|| '13164036458569648337239753460458722910223472318386943117783728128 ',
|| '14474011154664524427946373126085988481573677491474835889066354349131199152128 ',
|| '99999999999999999999999999999999999999999999999999999999999999999999999999999'
do zz = 1 to Words(p)
w = Word(p,zz)
if w > xx then
leave
perf.perfect.zz = w
end
zz = zz-1; perf.0 = zz
return zz
Practicals:
/* Practical number sequence */
procedure expose prac. fact.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
prac. = 0
/* Fast values */
if xx < 101 then do
p = '1 2 4 6 8 12 16 18 20 24 28 30 32 36 40 42 ',
|| '48 54 56 60 64 66 72 78 80 84 88 90 96 100 999 '
do zz = 1 to Words(p)
w = Word(p,zz)
if w > xx then
leave
prac.practical.zz = w
end
zz = zz-1; prac.0 = zz
return zz
end
/* Using function practical() */
zz = 1; prac.practical.1 = 1
do i = 2 by 2 to xx
if Practical(i) then do
zz = zz+1; prac.practical.zz = i
end
end
prac.0 = zz
return zz
Primes:
/* Prime sequence */
procedure expose prim. work. glob.
arg xx
numeric digits Max(Digits(),Xpon(xx)+1)
xx = Format(xx)
if \ Whole(xx) then say argument
if xx = 0 then say argument
/* Init */
prim. = 0
/* Fast values */
if xx = 1 then
return 0
/* Nth prime or threshold? */
s = (xx<0)
if s then do
xx = Abs(xx); zz = xx; xx = Primeest(xx)
end
/* Sieve of Eratosthenes */
work. = 1
do i = 3 by 2 to xx while i*i <= xx
if work.i then do
do j = i*i by i+i to xx
work.j = 0
end
end
end
/* Save results */
zz = 1; prim.prime.1 = 2; prim.flag.2 = 1
do i = 3 by 2 to xx
if work.i then do
zz = zz+1; prim.prime.zz = i; prim.flag.i = 1
end
end
if s then
zz = zz
prim.0 = zz
return zz
Primorials:
/* Primorial sequence */
procedure expose prim. prmo. work. glob.
arg xx
numeric digits Max(Digits(),Xpon(xx)+1)
xx = Format(xx)
if \ Whole(xx) then say argument
if xx = 0 then say argument
if xx < -50000 then say argument
if xx > 1e250000 then say argument
/* Init */
prmo. = 0
/* Fast values */
if Abs(xx) = 1 then
return 0
/* Xth primorial or xx threshold? */
yy = (xx<0)
if yy then do
a = Primes(xx); xx = Abs(xx)
end
else
a = Primes(600000)
/* Multiply primes */
p = prim.prime.1; zz = 1
do i = 2 to a
numeric digits Length(p)+6
p = p*prim.prime.i
if \ yy then
if p > xx then
leave
zz = zz+1
prmo.primorial.zz = p
end
prmo.0 = zz
return zz
Quadratics:
/* Quadratic prime sequence */
procedure expose quad. prim.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get primes */
call Primes(xx)
/* Init */
quad. = 0; quad.prime.1 = 2; zz = 1
/* Generate */
a = 2; b = 3; k = 1
do while b <= xx
b = a + k*k
if prim.flag.b = 1 then do
zz = zz+1; quad.prime.zz = b
a = b; k = 1
end
k = k+1
end
quad.0 = zz
return zz
Reversables:
/* Reversable prime sequence */
procedure expose reve. prim.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Get primes */
p = Primes(10**(Xpon(xx)+1))
/* Init */
reve. = 0; zz = 0
/* Generate */
do i = 1 to p
a = prim.prime.i
if a > xx then
leave i
b = Reverse(a)
if \ prim.flag.b then
iterate i
zz = zz+1; reve.prime.zz = a
end
reve.0 = zz
return zz
Totients:
/* Euler's totient number sequence */
procedure expose toti.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
toti. = 0
/* Recurring sequence */
do i = 1 to xx
toti.totient.i = i
end
do i = 2 to xx
if toti.totient.i < i then
iterate i
do j = i by i to xx
toti.totient.j = toti.totient.j-toti.totient.j/i
end
end
toti.0 = xx
return xx
Ufactors:
/* Unique prime factors */
procedure expose ufac.
arg xx
if \ Whole(xx) then say argument
if xx < 1 then say argument
/* Init */
ufac. = 0
/* Fast values */
if xx = 1 then
return 0
if xx < 4 then do
ufac.factor.1 = xx; ufac.0 = 1
return 1
end
if Prime(xx) then do
ufac.factor.1 = xx; ufac.0 = 1
return 1
end
/* Check low factors */
zz = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i); v = xx
do while xx//p = 0
xx = xx%p
end
if xx <> v then do
zz = zz+1; ufac.factor.zz = p
end
end
/* Check higher factors */
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then do
v = xx
do while xx//j = 0
xx = xx%j
end
if xx <> v then do
zz = zz+1; ufac.factor.zz = j
end
end
if p = 3 then
iterate
yy = j+2; v = xx
do while xx//yy = 0
xx = xx%yy
end
if xx <> v then do
zz = zz+1; ufac.factor.zz = yy
end
end
/* Last factor */
if xx > 1 then do
zz = zz+1; ufac.factor.zz = xx
end
ufac.0 = zz
/* Return number of factors */
return zz
Generic Quicksort
default [label]=Quicksort [lt]=< [eq]== [gt]=>
/* Sorting procedure - Build 7 Sep 2024 */
/* (C) Paul van den Eertwegh 2024 */
[label]:
/* Sort a stem on 1 or more key columns, syncing 0 or more data columns */
procedure expose (table)
arg table,keys,data
/* Collect keys */
kn = words(keys)
do x = 1 to kn
key.x = word(keys,x)
end
/* Collect data */
dn = words(data)
do x = 1 to dn
data.x = word(data,x)
end
/* Sort */
n = value(table||0); s = 1; sl.1 = 1; sr.1 = n
do until s = 0
l = sl.s; r = sr.s; s = s-1
do until l >= r
/* Up to 19 rows insertion sort is faster */
if