# 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 7 libraries. Other entries will link to this page.

See documentation for help and tips.

### Constants

```/* Constants library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Apery:
/* Apery constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 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
y = 0; v = 0
do k = 1
/* Add term */
y = y+s*((f10*f20*f21)**3*p)/(f32*f43**3)
if y = v then
leave
v = y
/* 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
y = y/24
numeric digits Digits()-2
end
return y+0

Artin:
/* Artin constant */
procedure expose glob.
/* Fast value */
y = 0.373955813619202288054728054346416415111629248606150042094742802417350182040028082344304317087250568981603
return y+0

Backhouse:
/* Backhouse constant */
procedure expose glob.
/* Fast value */
y = 1.45607494858268967139959535111654355765317837484713154027070243741400150626538989559964531940186030910992
return y+0

Bernstein:
/* Bernstein constant */
procedure expose glob.
/* Fast value */
y = 0.28016949902386913303643649123067200004248213981236
return y+0

Bloch:
/* Bloch constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.471861653452681784874468793611316149077012621739443243631460266888006805241859656902448655624594840992
else
/* Formula */
y = Sqrt((Sqrt(3)-1)*0.5)*Gamma(1/3)*Gamma(11/12)/Gamma(0.25)
return y+0

Bronze:
/* Bronze ratio constant = positive root of x^2-3x-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 3.302775637731994646559610633735247973125648286922623106355226528113583474146505222602309541009245359
else
/* Formula */
y = 0.5*(3+Sqrt(13))
return y+0

Brun:
/* Brun constant */
procedure expose glob.
/* Fast value */
y = 1.902160583104
return y+0

Cahen:
/* Cahen constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.6434105462883380261822543077575647632865878602682395059870309203074927764618326108484408955504634320
else do
numeric digits Digits()+2
/* Sylvester */
y = 0; v = 0; s = 2
do n = 0
y = y+(-1)**n/(s-1)
if y = v
then leave
v = y; s = s*s-s+1
end
numeric digits Digits()-2
end
return y+0

Capacity:
/* Logarithmic capacity of the unit disk constant */
procedure expose glob.
/* Fast value */
y = 0.5901702995080481130226689702792442936168583174407236497579321997021520903603578974892293080979039771
return y+0

Catalan1:
/* Catalan constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062
else do
numeric digits Digits()+2
s1 = 0; v = 0
do n = 0
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))
if s1 = v
then leave
v = s1
end
s2 = 0; v = 0
do n = 0
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))
if s2 = v
then leave
v = s2
end
y = 3*s1-2*s2
numeric digits Digits()-2
end
return y+0

Chaitin:
/* Chaitin constant */
procedure expose glob.
/* Fast value */
y = 0.0078749969978123844
return y+0

Champernowne:
/* Champernowne constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565
else do
p = Digits()+4
/* Concatenate integers */
y = '0.'
do n = 1 while Length(y) <= p
y = y||n
end
end
return y+0

Conway2:
/* Conway constant = real root of x^71-x^69-2x^68...-6x^2+3x-6 */
procedure expose glob.
/* Fast value or first guess */
y = 1.30357726903429639125709911215255189073070250465940487575486139062855088785246155712681576686442522555
if Digits() > 100 then do
/* Define polynomial */
f = '1 0 -1 -2 -1 2 2 1 -1 -1 -1 -1 -1 2 5 3 -2 -10 -3 -2 6 6 1 9 -3 -7 -8 -8 10 6 8 -5 -12 7 -7 7 1 -3 10 1 -6 -2 -10 -3 2 9 -3 ',
||  '14 -8 0 -7 9 3 -4 -10 -7 12 7 2 -12 -4 -2 5 0 1 -7 7 -4 12 -6 3 -6'
/* Derivatives */
f1 = Pdiff(f); f2 = Pdiff(f1)
numeric digits Digits()+2
/* Halley */
v = y
do forever
w = Peval(f,y); w1 = Peval(f1,y); w2 = Peval(f2,y)
y = y-(2*w*w1)/(2*w1*w1-w*w2)
if y = v then
leave
v = y
end
numeric digits Digits()-2
end
return y+0

Copeland:
/* Copeland-Erdos constant */
procedure expose glob.
/* Fast value */
y = 0.23571113171923293137414347535961677173798389971011031071091131271311371391491511571631671731791811911
return y+0

Devicci:
/* DeVicci constant = positive root of 4x^8-28x^6-7x^4+16x^2+16 */
procedure expose glob.
/* Fast value or first guess */
y = 1.007434756884279376098253595231099141925690114113669770234963798571152313280286777962520551474635923942
if Digits() > 100 then do
/* Define polynomial and derivative */
f = '4 0 -28 0 -7 0 16 0 16'; f1 = Pdiff(f)
numeric digits Digits()+2
/* Newton */
v = y
do forever
y = y - Peval(f,y)/Peval(f1,y)
if y = v then
leave
v = y
end
numeric digits Digits()-2
end
return y+0

Dottie:
/* Dottie constant = root of Cos(x) - x */
procedure expose glob.
/* Fast value or first guess */
y = 0.7390851332151606416553120876738734040134117589007574649656806357732846548835475945993761069317665318
if Digits() > 100 then  do
numeric digits Digits()+2
/* Halley */
v = y
do n = 1
s = Sin(y); c = Cos(y); f = c-y; f1 = -s-1; f2 = -c
y = y-(f*f1)/(f1*f1-0.5*f*f2)
if y = v then
leave
v = y
end
numeric digits Digits()-2
end
return y+0

Dubois:
/* Second Dubois-Reymond constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.194528049465325113615213730287503906590157785275923662043563911261286898039528881692156242539560897386876
else
/* Formula */
y = (E()**2-7)*0.5
return y+0

E:
/* Euler number */
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
/* Taylor */
y = 2; t = 1; v = y
do n = 2
t = t/n; y = y+t
if y = v then
leave
v = y
end
numeric digits Digits()-2
glob.e.p = y+0
end
return glob.e.p

Embree:
/* Embree-Trefethen constant */
procedure expose glob.
/* Fast value */
y = 0.70258
return y+0

Erdos1:
/* Erdos-Borwein constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.606695152415291763783301523190924580480579671505756435778079553691418420743486690565711801670155575897045429
else do
numeric digits Digits()+2
/* Series */
y = 0; v = 0
do n = 1
a = 2**n
y = y+(1/a**n)*(a+1)/(a-1)
if y = v
then leave
v = y
end
numeric digits Digits()-2
end
return y+0

Erdos2:
/* Erdos-Tenenbaum-Ford constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.08607133205593420688757309877692267776059110953033317349202023666542263581462287979938053460252876807163+0
else
/* Formula */
y = 1-(1+Ln(Ln(2)))/Ln(2)
return y+0

Euler:
/* Euler-Mascheroni constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495
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
y = s/v
numeric digits Digits()-2
end
return y+0

Feigenbaum1:
/* Feigenbaum first constant */
procedure expose glob.
/* Fast value */
y = 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716
return y+0

Feigenbaum2:
/* Feigenbaum second constant */
procedure expose glob.
/* Fast value */
y = 2.50290787509589282228390287321821578638127137672714997733619205677923546317959020670329964974643383412959
return y+0

Feller:
/* Feller-Tornier constant */
procedure expose glob.
/* Fast value */
y = 0.66131704946962233528976584627411853328547528983291635498090562622662503174312230494226174078428187
return y+0

Foias:
/* Foias constant */
procedure expose glob.
/* Fast value */
y = 1.18745235112650105459548015839651935121569268158586035301010412619878041872352540738702465760608657943378
return y+0

Fraction:
/* First continued fraction constant */
procedure expose glob.
/* Fast value */
y = 0.697774657964007982006790592551752599486658262998021232368630082816530852764641112996965654182676568723982
return y+0

Fransen:
/* Fransen-Robertson constant */
procedure expose glob.
/* Fast value */
y = 2.80777024202851936522150118655777293230808592093019829122005480959710088912190166551018530816819663814187
return y+0

Gauss:
/* Gauss constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.8346268416740731862814297327990468089939930134903470024498273701036819927095264118696911603512753241
else
/* Formula */
y = 1/AGmean(1,Sqrt(2))
return y+0

Gelfond1:
/* Gelfond constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 23.1406926327792690057290863679485473802661062426002119934450464095243423506904527835169719970675492
else
/* Formula */
y = Power(E(),Pi())
return y+0

Gelfond2:
/* Gelfond-Schneider constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.66514414269022518865029724987313984827421131371465949283597959336492044617870595486760918000519642
else
/* Formula */
y = Power(2,Sqrt(2))
return y+0

Gieseking:
/* Gieseking constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.01494160640965362502120255427452028594168930753029979201748910677659747625824402213647035422825669
else
/* Formula */
y = (Trigamma(1/3)-Trigamma(2/3))/(4*Sqrt(3))
return y+0

Glaisher:
/* Glaisher-Kinkelin constant */
procedure expose glob.
/* Fast value */
y = 1.28242712910062263687534256886979172776768892732500119206374002174040630885882646112973649195820237
return y+0

Golden:
/* Golden ratio constant = root of x^2-x-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939114
else
/* Formula */
y = (Sqrt(5)+1)*0.5
return y+0

Goldenangle:
/* Golden angle constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.39996322972865332223155550663361385312499901105811504293511275073130733823943879077996206066058396
else
/* Formula */
y = Pi()*(3-Sqrt(5))
return y+0

Golomb1:
/* Golomb constant */
procedure expose glob.
/* Fast value */
y = 0.624329988543550870992936383100837244179642620180529286973551902495638088855113254462460276195539868869
return y+0

Gompertz:
/* Gompertz constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.596347362323194074341078499369279376074177860152548781573484910482327219114874417470430497093612760344237
else do
numeric digits Digits()+2
/* Taylor */
y = 0; v = 0; s = 1; f = 1
do n = 1
s = -s; f = f*n; y = y+s/(n*f)
if y = v
then leave
v = y
end
y = -E()*(Euler()+y)
numeric digits Digits()-2
end
return y+0

Hafner:
/* Hafner-Sarnak-McCurley constant */
procedure expose glob.
/* Fast value */
y = 0.353236371854995984543516550432682011280164778566690446416085942814238325002669003483672078334335498967
return y+0

Heath:
/* Heath-Brown-Moroz constant */
procedure expose glob.
/* Fast value */
y = 0.00131764115485317810981735232251358595125073432329525167879254742178602344409610895090834
return y+0

Hermite:
/* Second Hermite constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.15470053837925152901829756100391491129520350254025375203720465296795534460586669138743079117149905
else
/* Formula */
y = 2/Sqrt(3)
return y+0

Kepler:
/* Kepler-Bouwkamp constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.1149420448532962007010401574695987428307953372008635168440233965189660128253530511779407724849858370
else do
/* Many extra digits */
p1 = Digits(); p2 = 3*p1
numeric digits p2
/* Asymptotic */
y = 0; v = 0
do k = 1 to 100
a = 2*k; b = Zeta(a); c = 2**a
y = y+((c-1)/a)*b*(b-1-1/c)
numeric digits p1
if y = v
then leave
numeric digits p2
v = y
end
numeric digits p1+2
y = Exp(-2*y)
numeric digits p1
end
return y+0

Khinchin:
/* Kinchin constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.685452001065306445309714835481795693820382293994462953051152345557218859537152002801141174931847698
else do
numeric digits Digits()+2
/* Asymptotic */
y1 = 0; v = 0
do n = 1
y2 = 0
do k = 1 to 2*n-1
y2 = y2+(-1)**(k+1)/k
end
y1 = y1+((Zeta(2*n)-1)/n)*y2
if y1 = v
then leave
v = y1
end
y = Exp(y1/Ln(2))
numeric digits Digits()-2
end
return y+0

Komornik:
/* Komornik-Loreti constant */
procedure expose glob.
/* Fast value */
y = 1.78723165018296593301327489033700838533793140296181099778478147050555749175060568699131001863407533302
return y+0

Landau1:
/* Landau-Ramanujan constant */
procedure expose glob.
/* Fast value */
y = 0.7642236535892206629906987312500923281167905413934095147216866737496146416587328588384015050131312337
return y+0

Landau2:
/* Second Landau constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.54325896534297670695272829530061323113886329375835698895573256910043468701530783731285212710225971085753+0
else
/* Formula */
y = Power(2,5/3)*Pi()**2/(3*Gamma(1/3)**3)
return y+0

Laplace:
/* Laplace limit constant */
procedure expose glob.
/* Fast value */
y = 0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168
return y+0

Lebesgue2:
/* Lebesgue constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.9894312738311469517416488090188667149242014060911110482841224326443725315845786546346329831401895579
else do
numeric digits Digits()+2
/* Series */
y = 0; v = 0
do n = 0
y = y+(Lambda(2*n+2)-1)/(2*n+1)
if y = v then
leave
v = y
end
a = Pi()**2
y = 8*y/a+4*(2*Ln(2)+Euler())/a
numeric digits Digits()-2
end
return y+0

Lemniscate:
/* Lemniscate constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.6220575542921198104648395898911194136827549514316231628168217038007905870704142502302955329614291
else
/* Formula */
y = Pi()/AGmean(1,Sqrt(2))
return y+0

Levy1:
/* Levy-Khinchin constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.18656911041562545282172297594723712056835653647205433595425429865280963205625444330034830110848687594663
else
/* Formula */
y = Pi()**2/(12*Ln(2))
return y+0

Levy2:
/* Second Levy constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 3.27582291872181115978768188245384386360847552598237414940519892419072321564496035518127754047917452949269
else
/* Formula */
y = Exp(Levy1())
return y+0

Lieb:
/* Lieb constant */
procedure expose glob.
if Digits() < 101 then
y = 1.53960071783900203869106341467188654839360467005367166938293953729060712614115558851657438822866540060055
/* Fast value */
else
/* Formula */
y = 8/(3*Sqrt(3))
return y+0

Ln2:
/* Natural log of 2 constant */
procedure expose glob.
/* Fast value */
y = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156058633269964186875
return y+0

Ln4:
/* Natural log of 4 constant */
procedure expose glob.
/* Fast value */
y = 1.386294361119890618834464242916353136151000268720510508241360018986787243939389431211726653992837375
return y+0

Ln8:
/* Natural log of 8 constant */
procedure expose glob.
/* Fast value */
y = 2.079441541679835928251696364374529704226500403080765762362040028480180865909084146817589980989256063
return y+0

Ln10:
/* Natural log of 10 constant */
procedure expose glob.
/* Fast value */
y = 2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959830
return y+0

Loch:
/* Loch constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.97027011439203392574025601921001083378128470478516128661035052993125419989173704803621267490802902646924
else do
/* Formula */
y = 6*Ln(2)*Ln(10)/Pi()**2
end
return y+0

Magicangle:
/* Magic angle constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.9553166181245092781638571025157577542434146950100054909596981293219120459039764553873916025856280734
else
/* Formula */
y = ACos(1/Sqrt(3))
return y+0

Meissel:
/* Meissel-Mertens constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.261497212847642783755426838608695859051566648261199206192064213924924510897368209714142631434246651051617
else do
numeric digits digits()+2
y = 0; v = 0
do n = 2 to 1000
t = Moebius(n)*Ln(Zeta(n))/n
if t <> 0 then do
y = y+t
if y = v then
leave
v = y
end
end
y = Euler()+y
end
return y+0

Mills:
/* Mills constant */
procedure expose glob.
/* Fast value */
y = 1.30637788386308069046861449260260571291678458515671364436805375996643405376682659882150140370119739570729
return y+0

Mrb:
/* MRB constant */
procedure expose glob.
/* Fast value */
y = 0.187859642462067120248517934054273230055903094900138786172004684089477231564660213703296654433107496903842
return y+0

Nielsen:
/* Nielsen-Ramanujan constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.82246703342411321823620758332301259460947495060339921886777911468500373520160043691681445030987935
else
/* Formula */
y = Pi()**2/12
return y+0

Niven:
/* Niven constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.70521114010536776428855145343450816076202765165346909999428490654731319216812249193424413210087100179
else do
numeric digits Digits()+2
/* Series */
y = 0; v = 0
do k = 2 to 100
y = y+1-1/Zeta(k)
if y = v
then leave
v = y
end
y = y+1
numeric digits Digits()-2
end
return y+0

Omega:
/* Omega constant = Lambertw(1) = root of x*e^x-1 */
procedure expose glob.
/* Fast value or first guess */
y = 0.5671432904097838729999686622103555497538157871865125081351310792230457930866845666932194469617522946
if Digits() > 100 then do
numeric digits Digits()+2
/* Halley */
v = y
do n = 1
a = Exp(y); b = a*y-1
y = y-b/(a*(y+1)-((y+2)*b/(2*y+2)))
if y = v
then leave
v = y
end
numeric digits Digits()-2
end
return y+0

Paperfolding:
/* Paperfolding constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.850736188201867260367797760532066604411399493082710647728168261031830158459319445348545982642193923996091
else do
numeric digits Digits()+2
/* Definition */
y = 0; v = 0
do n = 0
y = y+8**(2**n)/(2**(2**(n+2))-1)
if y = v
then leave
v = y
end
numeric digits Digits()-2
end
return y+0

Parabolic:
/* Universal parabolic constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.29558714939263807403429804918949038759783220363858348392997534664410966268413312668409442623789762
else do
/* Formula */
a = Sqrt(2); y = Ln(a+1)+a
end
return y+0

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 */
y = 0
do n = 0
v = y; y = y + Fact(6*n)*(13591409+545140134*n)/(Fact(3*n)*Fact(n)**3*-640320**(3*n))
if y = v then
leave
end
y = 4270934400/(Sqrt(10005)*y)
end
else do
/* Agmean */
y = 0.25; a = 1; g = Sqrt(0.5); n = 1
do until a = v
v = a
x = (a+g)*0.5; g = Sqrt(a*g)
y = y-n*(x-a)**2; n = n+n; a = x
end
y = a*a/y
end
numeric digits Digits()-2
glob.pi.p = y+0
end
return glob.pi.p

Plastic:
/* Plastic ratio constant = real root of x^3-x-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.32471795724474602596090885447809734073440405690173336453401505030282785124554759405469934798178728
else do
/* Formula */
a = Sqrt(69)/18
y = Cbrt(0.5+a)+Cbrt(0.5-a)
end
return y+0

Porter:
/* Porter constant */
procedure expose glob.
/* Fast value */
y = 1.46707807943397547289779848470722995344990332241488777739968581761660674432904480843036932751117401521266
return y+0

Prouhet:
/* Prouhet-Thue-Morse constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.4124540336401075977833613682584552830894783744557695575733794153487935923657825889638045404862121334
else do
numeric digits Digits()+2
/* Infinite product */
y = 1; v = 0
do n = 0
y = y*(1-1/2**(2**n))
if y = v
then leave
v = y
end
y = (2-y)*0.25
numeric digits Digits()-2
end
return y+0

Ramanujan:
/* Ramanujan constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 262537412640768743.999999999999250072597198185688879353856337336990862707537410378210647910118607313
else
/* Formula */
y = Exp(Pi()*Sqrt(163))
return y+0

Recifibo:
/* Reciprocal Fibonacci constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 3.35988566624317755317201130291892717968890513373196848649555381532513031899668338361541621645679008729704
else do
numeric digits Digits()+2
/* Recurring */
y = 1; v = 0; a = 0; b = 1
do n = 2
c = a+b; y = y+1/c
if y = v
then leave
a = b; b = c; v = y
end
numeric digits Digits()-2
end
return y+0

Robbins:
/* Robbins constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.661707182267176235155831133248413581746400135790953604808944229479584646138597631306652480768107120151709
else
/* Formula */
y = (4+17*Sqrt(2)-6*Sqrt(3)-7*Pi())/105+(Ln(1+Sqrt(2)))*0.2+(Ln(2+Sqrt(3)))*0.4
return y+0

Salem:
/* Salem constant = largest positive root of x^10+x^9-x^7-x^6-x^5-x^4-x^3+x+1 */
procedure expose glob.
/* Fast value or first guess */
y = 1.176280818259917506544070338474035050693415806564695259830106347029688376548549962096830115581815395
if Digits() > 100 then  do
/* Define polynomial and derivative */
f = '1 1 0 -1 -1 -1 -1 -1 0 1 1'; f1 = Pdiff(f)
numeric digits Digits()+2
/* Newton */
v = 0
do n = 1 to 100
y = y - Peval(f,y)/Peval(f1,y)
if y = v then
leave
v = y
end
numeric digits Digits()-2
end
return y+0

Sierpinski:
/* Sierpinski constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.584981759579253217065893587383171160088051651852630917321544987971932044001157120211117724527064
else
/* Formula */
y = Pi()*Ln(Exp(2*Euler())/(2*Gauss()**2))
return y+0

Silver:
/* Silver ratio constant = positive root of x^2-2x-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157
else
/* Formula */
y = Sqrt(2)+1
return y+0

Soldner:
/* Soldner constant */
procedure expose glob.
/* Fast value */
y = 1.45136923488338105028396848589202744949303228364801586309300455766242559575451783565953135771108682884
return y+0

Somos:
/* Somos quadratic recurrence constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.66168794963359412129581892274995074996441863502506820818971116802560902982638372790836917641146116715528
else do
numeric digits Digits()+2
/* Definition */
y = 1; v = 0; a = 2
do n = 1
y = y*Power(n,1/a)
if y = v
then leave
v = y; a = 2*a
end
numeric digits Digits()-2
end
return y+0

Sqrt2:
/* Square root of 2 constant */
procedure expose glob.
/* Fast value */
y = 1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573
return y+0

Sqrt3:
/* Square root of 3 constant */
procedure expose glob.
/* Fast value */
y = 1.732050807568877293527446341505872366942805253810380628055806979451933016908800037081146186757248576
return y+0

Sqrt5:
/* Square root of 5 constant */
procedure expose glob.
/* Fast value */
y = 2.236067977499789696409173668731276235440618359611525724270897245410520925637804899414414408378782275
return y+0

Stephens:
/* Stephens constant */
procedure expose glob.
/* Fast value */
y = 0.57595996889294543964316337549249669250651396717649236360064079866537255169886852843640987209172618
return y+0

Supergolden:
/* Super golden ratio constant = real solution of x^3-x^2-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.46557123187676802665673122521993910802557756847228570164318311124926299668501784047812580119490927
else do
/* Formula */
a = 3*Sqrt(93)
y = (1+Cbrt((29+a)*0.5)+Cbrt((29-a)*0.5))/3
end
return y+0

Taniguchi:
/* Taniguchi constant */
procedure expose glob.
/* Fast value */
y = 0.678234491917391978035538279482894814096332239189440103036460415964983370740123233213762122933484616888328
return y+0

Tau1:
/* Tau constant = 2*pi */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 6.28318530717958647692528676655900576839433879875021164194988918461563281257241799725606965068423414
else
/* Formula */
y = 2*Pi()
return y+0

Tetranacci2:
/* Tetranacci constant = root of x^4-x^3-x^2-x-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.92756197548292530426190586173662216869855425516338472714664703800966606229781555914981825346189065325
else do
/* Formula */
a = Cbrt((Sqrt(177304464)+7020))
y = 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 y+0

Tribonacci2:
/* Tribonacci constant = root of x^3-x^2-x-1 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 1.839286755214161132551852564653286600424178746097592246778758639404203222081966425738435419428307014
else do
/* Formula */
a = 3*Sqrt(33)
y = (1+Cbrt(19+a)+Cbrt(19-a))/3
end
return y+0

Twinprimes:
/* Twin primes constant */
procedure expose glob.
/* Fast value */
y = 0.660161815846869573927812110014555778432623360284733413319448423335405642304495277143760031413839867911779
return y+0

Vanderpauw:
/* Van der Pauw constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 4.532360141827193809627682945716666810171861467723795584186016547940600953721305102259083879604016089653
else
/* Formula */
y = Pi()/Ln(2)
return y+0

Viswanath:
/* Viswanath constant */
procedure expose glob.
/* Fast value */
return 1.1319882487943+0

Wallis:
/* Wallis constant = real root of x^3-2x-5 */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 2.09455148154232659148238654057930296385730610562823918030412852904531218998348366714626728177715776
else do
/* Formula */
a = Sqrt(1929)
y = Cbrt((45-a)/18)+Cbrt((45+a)/18)
end
return y+0

Weierstrass:
/* Weierstrass constant */
procedure expose glob.
if Digits() < 101 then
/* Fast value */
y = 0.474949379987920650332504636327982968559549373217202982283331024864557929174883860274275641250502144418903
else
/* Formula */
y = Power(2,1.25)*Sqrt(Pi())*Exp(Pi()*0.125)/Gamma(0.25)**2
return y+0

Zscore:
/* Z-score 97.5 percentile constant */
procedure expose glob.
/* Fast value */
y = 1.95996398454005423552459443052055152795555007786954839847695264636163527414488266779825470949281420601799
return y+0
```

### Functions

```/* Functions library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Acos:
/* Arc cosine function */
procedure expose glob.
arg x
/* Validity */
if Abs(x) > 1 then
return 'X'
/* Formula */
return Pi()/2-ASin(x)

Acosh:
/* Arc cosine hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if x < 1 then
return 'X'
/* Formula */
return Ln(x+Sqrt(x*x-1))

Acot:
/* Arc cotangent function */
procedure expose glob.
arg x
/* Formula */
return ACos(x/Sqrt(x*x+1))

Acoth:
/* Arc cotangent hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if Abs(x) <= 1 then
return 'X'
/* Formula */
return Ln((x+1)/(x-1))*0.5

Acsc:
/* Arc cosecant function */
procedure expose glob.
arg x
/* Validity */
if Abs(x) < 1 then
return 'X'
/* Formulas */
if x > 0 then
return ASin(1/x)
else
return -(Pi()+ASin(1/x))

Acsch:
/* Arc cosecant hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if x = 0 then
return 'X'
/* Formula */
a = 1/x
return Ln(a+Sqrt(a*a+1))

Agmean:
/* Arithmetic geometric mean function */
procedure expose glob.
arg x,y
/* Validity */
if x < 0 then
return 'X'
if y < 0 then
return 'X'
if x = 0 then
return y/2
/* Fast values */
if y = 0 then
return 0
if x = y then
return x
numeric digits Digits()+2
/* Series */
vx = x+1
do while x <> vx
vx = x; vy = y
x = (vx+vy)*0.5; y = Sqrt(vx*vy)
end
numeric digits Digits()-2
return x+0

Alog:
/* Antilog function = 10^x */
procedure expose glob.
arg x
/* Formula */
return Power(10,x)

Amean:
/* Arithmetic mean function */
procedure expose glob.
arg x,y
/* Formula */
return (x+y)/2

Asec:
/* Arc secant function */
procedure expose glob.
arg x
/* Validity */
if Abs(x) < 1 then
return 'X'
/* Formula */
return Sign(x)*ACos(1/x)

Asech:
/* Arc secant hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if x <= 0 then
return 'X'
if x > 1 then
return 'X'
/* Formula */
a = 1/x
return Ln(a+Sqrt(a*a-1))

Asin:
/* Arc sine function */
procedure expose glob.
arg x
/* Validity */
if Abs(x) > 1 then
return 'X'
/* Acos is faster for x in [3/4,1) */
if Abs(x) >= 0.75 then
return Sign(x)*ACos(Sqrt(1-x*x))
numeric digits Digits()+2; numeric fuzz 2
/* Taylor series */
t = x; y = x; x = x*x; v = y
do n = 2 by 2
t = t*x*(n-1)/n; y=y+t/(n+1)
if y = v then
leave
v = y
end
numeric digits Digits()-2
return y+0

Asinh:
/* Arc sine hyperbolic function */
procedure expose glob.
arg x
/* Formula */
return Ln(x+Sqrt(x*x+1))

Atan:
/* Arc tangent function */
procedure expose glob.
arg x
numeric digits Digits()+2
/* Twice argument reduction */
x = x/(1+Sqrt(x*x+1)); x = x/(1+Sqrt(x*x+1))
/* Taylor series */
numeric fuzz 2
t = x; y = x; v = y; x = x*x
do n = 3 by 2
t = -t*x; y = y+t/n
if y = v then
leave
v = y
end
y = 4*y
numeric digits Digits()-2
return y+0

Atanh:
/* Arc tangent hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if Abs(x) >= 1 then
return 'X'
/* Formula */
return Ln((1+x)/(1-x))*0.5

Atan2:
/* Arc tangent correct quadrant function */
procedure expose glob.
arg y,x
/* Validity */
if x = 0 & y = 0 then
return 'X'
/* Formulas */
if x > 0 then
return ATan(y/x)
if x < 0 then
if y < 0 then
return ATan(y/x)-Pi()
else
return ATan(y/x)+Pi()
if y < 0 then
return -Pi()*0.5
else
return Pi()*0.5

Beta:
/* Beta function */
procedure expose glob.
arg x,y
/* Formula */
return (Gamma(x)*Gamma(y))/Gamma(x+y)

Cbrt:
/* Cubic root function = x^(1/3) */
procedure expose glob.
arg x
/* Fast values */
if x = 0 then
return 0
if x = 1 then
return 1
/* Negative argument allowed */
if x < 0 then do
s = -1; x = Abs(x)
end
else
s = 1
p = Digits()
numeric digits p+2
/* Argument reduction to [0,1000) */
i = Xpon(x); i = (i-(2*(i<0)))%3; x = x/1000**i
/* First guess 1 digit accurate */
t = '4.5 17.5 45.5 94.5 170.5 279.5 427.6 620.5 864.5 1000'
do y = 1 until Word(t,y) > x
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
y = (x+2*y*y*y)/(3*y*y)
end
/* Inverse reduction */
y = y*10**i
numeric digits p
return s*y

Ceil:
/* Ceiling function */
procedure expose glob.
arg x
/* Formulas */
if Whole(x) then
return x
else
return Trunc(x)+(x>=0)

Cos:
/* Cosine function */
procedure expose glob.
arg x
numeric digits Digits()+2
/* Argument reduction */
x = x//(2*Pi())
if Abs(x) > Pi() then
x = x-Sign(x)*2*Pi()
/* Series Taylor */
t = 1; y = 1; x = x*x; v = y
do n = 2 by 2
t = -t*x/(n*(n-1)); y = y+t
if y = v then
leave
v = y
end
numeric digits Digits()-2
return y+0

Cosh:
/* Cosine hyperbolic function */
procedure expose glob.
arg x
y = Exp(x)
/* Formula */
return (y+1/y)*0.5

Cot:
/* Cotangent function */
procedure expose glob.
arg x
y = Sin(x)
/* Validity */
if y = 0 then
return 'X'
/* Formula */
return Cos(x)/y

Coth:
/* Cotangent hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if x = 0 then
return 'X'
/* Formula */
return 1/Tanh(x)

Csc:
/* Cosecant function */
procedure expose glob.
arg x
y = Sin(x)
/* Validity */
if y = 0 then
return 'X'
/* Formula */
return 1/y

Csch:
/* Cosecant hyperbolic function */
procedure expose glob.
arg x
/* Validity */
if x = 0 then
return 'X'
/* Formula */
return 1/Sinh(x)

Deg:
/* Radians->degrees function */
procedure expose glob.
arg x
/* Formula */
return x*180/Pi()

Digamma:
/* Digamma function */
procedure expose glob.
arg x
/* Validity */
if x <= 0 & Whole(x) then
return 'X'
/* Formula */
if x = 0.125 then
return -0.5*Pi()-4*Ln(2)-( Pi()+Ln(Sqrt(2)+1)-Ln(Sqrt(2)-1))/Sqrt(2)-Euler()
if x = 0.25 then
return -0.5*Pi()-3*Ln(2)-Euler()
if x = 0.5 then
return -2*Ln(2)-Euler()
if x = 1 then
return -Euler()
if x > 0 & x < 500 & Whole(x) then
return Harmonic(x-1)-Euler()
/* Fast series */
if x > 0 & Half(x) then do
x = Trunc(x)
y = 0
do n = 1 to x
y = y+2/(2*n-1)
end
y = y-2*Ln(2)-Euler()
return y
end
numeric digits Digits()+2
/* Argument reduction */
p = Digits()
a = 0
do while x < p
a = a+1/x; x = x+1
end
/* Series using Zeta function */
y = -0.5/x; v = 0
do n = 2 by 2
y =  y+Zeta(1-n)/(x**n)
numeric fuzz 2
if y = v then
leave
numeric fuzz
v = y
end
/* Inverse reduction */
y = Ln(x)+y-a
numeric digits Digits()-2
return y+0

Dxgamma:
/* First derivative of Gamma function */
procedure expose glob.
arg x
/* Validity */
if x <= 0 & Whole(x) then
return 'X'
p = Digits()
numeric digits p*2+2
/* From derivative definition */
d = 10**(-p)
y = (Gamma(x+d)-Gamma(x))/d
numeric digits p
return y+0

Dxzeta:
/* First derivative of Zeta function */
procedure expose glob.
arg x
/* Validity */
if x = 1 then
return 'X'
/* Formula */
if x < 0 then
if Whole(x) then
if Even(x) then
return (-1)**(-0.5*x)*Fact(-x)*Zeta(1-x)/(2*(2*Pi())**-x)
if x = -2 then
return -Zeta(3)/(4*Pi()**2)
if x = -1 then
return 1/12-Ln(Glaisher())
if x = 0 then
return -0.5*Ln(2*Pi())
if x = 0.5 then
return 0.25*(Pi()+2*Euler()+6*Ln(2)+2*Ln(Pi()))*Zeta(0.5)
if x = 2 then
return Pi()*Pi()*(Euler()+Ln(2)-12*Ln(Glaisher())+Ln(Pi()))/6
if x > 5 then do
/* Taylor series for the Gamma derivative */
numeric digits Digits()+2
y = 0; v = 1
a = Power(2,x); b = Ln(4); c = a-2
do n = 1 to 10000
if n//2 = 0 then
t1 = 1
else
t1 = -1
y = y+(t1*a*Power(n,-x))*(b+c*Ln(n))/(c*c)
if y = v
then leave
v = y
end
numeric digits Digits()-2
end
else do
/* Many extra digits needed */
p = Digits()
numeric digits p*4
/* From derivative definition */
d = 10**(-p+p%5); y = (Zeta(x+d)-Zeta(x))/d
numeric digits p
end
return y+0

Erf:
/* Error function */
procedure expose glob.
arg x
numeric digits Digits()+2
/* Asymptotic series */
a = (2*x*Exp(-x*x))/Sqrt(Pi()); b = 2*x*x
y = 0; v = y; f = 1; g = 1
do n = 0
y = y+f/g
if y = v then
leave
v = y; f = f*b; g = g*(2*n+3)
end
y = a*y
numeric digits Digits()-2
return y+0

Erfc:
/* Error complementary function */
procedure expose glob.
arg x
/* Formula */
return 1-Erf(x)

Exp:
/* Exponential function = e^x */
procedure expose glob.
arg x
numeric digits Digits()+2
/* Fast values */
if Whole(x) then
return E()**x
/* Argument reduction */
i = x%1
if Abs(x-i) > 0.5 then
i = i+Sign(x)
/* Taylor series */
x = x-i; y = 1; t = 1; v = y
do n = 1
t = (t*x)/n; y = y+t
if y = v then
leave
v = y
end
/* Inverse reduction */
y = y*e()**i
numeric digits Digits()-2
return y+0

Factorial:
/* Factorial function = x! */
procedure expose glob.
arg x
/* Validity */
if x < 0 & Whole(x) then
return 'X'
/* Fast values */
if x = 0 then
return 1
/* Formula */
return Gamma(x+1)

Floor:
/* Floor function */
procedure expose glob.
arg x
/* Formula */
if Whole(x) then
return x
else
return Trunc(x)-(x<0)

Frac:
/* Fractional part function */
procedure expose glob.
arg x
/* Formula */
return x-x%1

Gamma:
/* Gamma function */
procedure expose glob.
arg x
/* Validity */
if x < 1 & Whole(x) then
return 'X'
/* Formulas for negative and positive (half)integers */
if x < 0 then do
if Half(x) then do
numeric digits Digits()+2
i = Abs(Floor(x)); y = (-1)**i*2**(2*i)*Fact(i)*Sqrt(Pi())/Fact(2*i)
numeric digits Digits()-2
return y+0
end
end
if x > 0 then do
if Whole(x) then
return Fact(x-1)
if Half(x) then do
numeric digits Digits()+2
i = Floor(x); y = Fact(2*i)*Sqrt(Pi())/(2**(2*i)*Fact(i))
numeric digits Digits()-2
return y+0
end
end
p = Digits()
if p < 61 then do
/* Lanczos with predefined coefficients */
/* Map negative x to positive x */
if x < 0 then
return Pi()/(Gamma(1-x)*Sin(Pi()*x))
/* Argument reduction to interval (0.5,1.5) */
numeric digits p+2
c = Trunc(x); x = x-c
if x < 0.5 then do
x = x+1; c = c-1
end
/* Series coefficients 1/Gamma(x) 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 */
x = x-1; s = 0
do k = 52 by -1 to 1
s = s*x+c.k
end
y = 1/s
/* Undo reduction */
if c = -1 then
y = y/x
else do
do i = 1 to c
y = (x+i)*y
end
end
end
else do
x = x-1
/* Spouge */
/* Estimate digits and iterations */
q = Floor(p*1.5); a = Floor(p*1.3)
numeric digits q
/* Series */
s = 0
do k = 1 to a-1
s = s+((-1)**(k-1)*Power(a-k,k-0.5)*Exp(a-k))/(Fact(k-1)*(x+k))
end
s = s+Sqrt(2*Pi()); y = Power(x+a,x+0.5)*Exp(-a-x)*s
end
/* Normalize */
numeric digits p
return y+0

Gmean:
/* Geometric mean function */
procedure expose glob.
arg x,y
/* Formula */
return Sqrt(x*y)

Isqrt:
/* Integer square root function = Floor(x^(1/2)) */
procedure expose glob.
arg x
/* Validity */
if x < 0 then
return 'X'
/* Fast values */
if x < 1 then
return 0
if x < 4 then
return 1
if x < 9 then
return 2
/* Uses only integers */
x = x%1; q = 1
do until q > x
q = q*4
end
y = 0
do while q > 1
q = q%4; t = x-y-q; y = y%2
if t >= 0  then do
x = t; y = y+q
end
end
return y

Lambda:
/* Dirichlet lambda function */
procedure expose glob.
arg x
/* Formula */
return (1-Power(2,-x))*Zeta(x)

Lambertw0:
/* Lambert W0 function */
procedure expose glob.
arg x
/* Fast values */
if x = 0 then
return 0
numeric digits Digits()+2
/* Argument reduction */
e = E()
select
when x > E() then do
a = Ln(x); y = a-Ln(a)
end
when x > 0 then
y = x/E()
when x >= -1/E() then do
a = E()*x; b = a+1; c = Sqrt(b); y = (a*Ln(1+c))/(b+c)
end
/* Validity */
otherwise
return 'X'
end
/* Halley */
numeric fuzz 2
v = y
do forever
a = Exp(y); b = a*y; y = y-(b-x)/(a+b-((y+2)*(b-x))/(2*y+2))
if y = v then
leave
v = y
end
numeric fuzz
numeric digits Digits()-2
return y+0

Lambertwm1:
/* Lambert W-1 function */
procedure expose glob.
arg x
/* Validity */
if x < -1/E() then
return 'X'
if x >= 0 then
return 'X'
numeric digits Digits()+2
/* Argument reduction */
select
when x > -0.25 then do
a = Ln(-x); y = a-Ln(-a)
end
when x > -1/E() then do
y = -1-Sqrt(2)*Sqrt(E()*x+1)
end
end
/* Halley */
v = y
numeric fuzz 2
do forever
a = Exp(y); b = a*y; y = y-(b-x)/(a+b-((y+2)*(b-x))/(2*y+2))
if y = v then
leave
v = y
end
numeric fuzz
numeric digits Digits()-2
return y+0

Lcm:
/* Least common multiple function */
procedure expose glob.
arg x,y
x = Abs(x); y = Abs(y)
/* Special values */
if y = 0 then
return 0
/* Euclid */
a = x*y
do until y = 0
parse value x//y y with y x
end
return a%x

Ln:
/* Natural logarithm function */
procedure expose glob.
arg x
/* Validity */
if x <= 0 then
return 'X'
/* Fast values */
if x = 1 then
return 0
p = Digits()
/* In memory? */
if glob.ln.x.p <> '' then
return glob.ln.x.p
/* Precalculated values */
if x = 2 & p < 101 then do
glob.ln.x.p = Ln2()
return glob.ln.x.p
end
if x = 4 & p < 101 then do
glob.ln.x.p = Ln4()
return glob.ln.x.p
end
if x = 8 & p < 101 then do
glob.ln.x.p = Ln8()
return glob.ln.x.p
end
if x = 10 & p < 101 then do
glob.ln.x.p = Ln10()
return glob.ln.x.p
end
numeric digits p+2
/* Argument reduction */
z = x; i = 0; e = 1/E()
if z < 0.5 then do
y = 1/z
do while y > 1.5
y = y*e; i = i-1
end
z = 1/y
end
if z > 1.5 then do
do while z > 1.5
z = z*e; i = i+1
end
end
/* Taylor series */
q = (z-1)/(z+1); f = q; y = q; v = q; q = q*q
do n = 3 by 2
f = f*q; y = y+f/n
if y = v then
leave
v = y
end
numeric digits p
/* Inverse reduction */
glob.ln.x.p = 2*y+i
return glob.ln.x.p

Log:
/* Briggs logarithm function */
procedure expose glob.
arg x
/* Formula */
return Ln(x)/Ln(10)

Logxy:
/* Logarithm x base y function */
procedure expose glob.
arg x,y
/* Validity */
if y = 1 then
return 'X'
/* Formula */
return Ln(x)/Ln(y)

Mant:
/* Mantissa function */
procedure expose glob.
arg x
/* Fast values */
if x = 0 then
return 0
/* Formula */
y = x*1E+99999
return Left(y,Pos('E',y)-1)

Nroot:
/* Nth root function = x^(1/n) */
procedure expose glob.
arg x,n
/* Validity */
if \ Whole(n) then
return 'X'
if n < 1 then
return 'X'
if x < 0 then
if Even(n) then
return 'X'
/* Fast values */
if x = 0 then
return 0
if x = 1 then
return 1
/* Formulas using faster methods */
if n = 2 then
return Sqrt(x)
if n = 3 then
return Cbrt(x)
if n = 4 then
return Qtrt(x)
/* Calculate */
sx = Sign(x); x = Abs(x)
p1 = Digits(); p2 = p1+2
numeric digits 3
/* First guess low accuracy */
y = 1/Exp(Ln(x)/n)
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/n; b = n+1
do j = k to 1 by -1
numeric digits d.j
y = y*a*(b-x*y**n)
end
y = 1/y
if sx < 0 then
y = -y
numeric digits p1
return y+0

Power:
/* Power function function = x^y */
procedure expose glob.
arg x,y
/* Validity */
if x < 0 then
return 'X'
/* Fast values */
if x = 0 then
return 0
if y = 0 then
return 1
/* Fast formula */
if Whole(y) then
return x**y
/* Formulas */
if Abs(y//1) = 0.5 then
return Sqrt(x)**Sign(y)*x**(y%1)
else
return Exp(y*Ln(x))

Powermod:
/* Power modulus function = x^(y mod z) */
procedure expose glob.
arg x,y,z
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if \ Whole(z) then
return 'X'
if x < 0 then
return 'X'
if y < 0 then
return 'X'
if z < 1 then
return 'X'
/* Fast values */
if z = 1 then
return 0
/* Binary method */
numeric digits Max(Length(Format(x,,,0)),Length(Format(y,,,0)),2*Length(Format(z,,,0)))
b = x//z; r = 1
do while y > 0
if y//2 then
r = r*x//z
x = x*x//z; y = y%2
end
return r

Qtrt:
/* Quartic root function = x^(1/4) */
procedure expose glob.
arg x
/* Validity */
if x < 0 then
return 'X'
/* Fast values */
if x = 0 then
return 0
if x = 1 then
return 1
p = Digits()
numeric digits p+2
/* Argument reduction to [0,10000) */
i = Xpon(x); i = (i+(i<0))%4-(i<0); x = x/10000**i
/* First guess 1 digit accurate */
t = '8.5 48.5 168.5 440.5 960.5 1848.5 3248.5 5328.5 8280.5 10000'
do y = 1 until Word(t,y) > x
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*x
do k = n to 1 by -1
numeric digits d.k
y = f1*y+f2/(y*y*y)
end
numeric digits p
/* Inverse reduction */
return y*10**i

/* Degrees->radians function */
procedure expose glob.
arg x
/* Formula */
return x*Pi()/180

Rand:
/* Random number between x and y function in max 12 digits precision */
procedure expose glob.
arg x,y
/* Validity */
if x <> '' then do
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if x >= y then
return 'X'
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) */
z = '0.'Left(glob.rand,Length(glob.rand)-3)
numeric digits 12
if x = '' then
return z/1
else
return Floor(z*(y+1-x)+x)

Round:
/* Round to nearest integer function */
procedure expose glob.
arg x,y
if y = '' then
y = 0
/* Validity */
if \ Whole(y) then
return 'X'
if y < 0 then
return 'X'
/* Formula */
return Format(x,,y)/1

Sec:
/* Secant function */
procedure expose glob.
arg x
y = Cos(x)
/* Validity */
if y = 0 then
return 'X'
/* Formula */
return 1/y

Sech:
/* Secant hyperbolic function */
procedure expose glob.
arg x
/* Formula */
return 1/Cosh(x)

Sin:
/* Sine function */
procedure expose glob.
arg x
numeric digits Digits()+2
/* Argument reduction */
u = Pi(); x = x//(2*u)
if Abs(x) > u then
x = x-Sign(x)*2*u
/* Taylor series */
t = x; y = x; x = x*x; v = y
do n = 2 by 2
t = -t*x/(n*(n+1)); y = y+t
if y = v then
leave
v = y
end
numeric digits Digits()-2
return y+0

Sinh:
/* Sine hyperbolic function */
procedure expose glob.
arg x
y = Exp(x)
/* Formula */
return (y-1/y)*0.5

Sqrt:
/* Square root function = x^(1/2) */
procedure expose glob.
arg x
/* Validity */
if x < 0 then
return 'X'
/* Fast values */
if x = 0 then
return 0
if x = 1 then
return 1
p = Digits()
/* Predefined values */
if x = 2 & p < 101 then
return Sqrt2()
if x = 3 & p < 101 then
return Sqrt3()
if x = 5 & p < 101 then
return Sqrt5()
numeric digits p+2
/* Argument reduction to [0,100) */
i = Xpon(x); i = (i-(i<0))%2; x = x/100**i
/* First guess 1 digit accurate */
t = '2.5 6.5 12.5 20.5 30.5 42.5 56.5 72.5 90.5 100'
do y = 1 until Word(t,y) > x
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
y = (y+x/y)*0.5
end
numeric digits p
return y*10**i

Tan:
/* Tangent function */
procedure expose glob.
arg x
y = Cos(x)
/* Validity */
if y = 0 then
return 'X'
/* Formula */
return Sin(x)/y

Tanh:
/* Tangent hyperbolic function */
procedure expose glob.
arg x
/* Formula */
y = Exp(2*x)
return (y-1)/(y+1)

Tau2:
/* Ramanujan tau function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Niebur series */
numeric digits Digits()+2
s = 0; a = 18*x*x
do n = 1 to x-1
b = n*n
s = s+b*(35*b-52*n*x+a)*Sigma(n)*Sigma(x-n)
end
y = x**4*Sigma(x)-24*s
numeric digits Digits()-2
return y+0

Trigamma:
/* Trigamma function */
procedure expose glob.
arg x
/* Validity */
if x <= 0 & Whole(x) then
return 'X'
/* Fast values */
if x = 1 then
return Zeta(2)
numeric digits Digits()+2
p = Digits()
/* Argument reduction */
a = 0
do while x < p
a = a+1/(x*x); x = x+1
end
/* Series using Bernoulli numbers */
y = 0; v = 0
do n = 2 by 2
y =  y+Bernoulli(n)/x**(n+1)
if y = v then
leave
v = y
end
/* Inverse reduction and result */
y = y+1/x+1/(2*x*x)+a
numeric digits Digits()-2
return y+0

Xpon:
/* Exponent function */
procedure expose glob.
arg x
/* Formula */
if x = 0 then
return 0
else
return Right(x*1E+99999,6)-99999

Zeta:
/* Zeta function */
procedure expose glob.
arg x
/* Validity */
if x = 1 then
return 'X'
/* Fast values */
if x < 0 & Whole(x) then
if Even(x) then
return 0
if x = 0 then
return -0.5
if x = 3 then
return Apery()
/* Formulas */
if x < 0 & Whole(x) then
if Odd(x) then
return ((-1)**-x*Bernoulli(1-x))/(1-x)
if x > 0 & Whole(x) then
if Even(x) then
return ((-1)**(x/2+1)*Bernoulli(x)*(2*Pi())**x)/(2*Fact(x))
if x < 0 then
return -2*Factorial(-x)*Sin(Pi()*-x/2)*Zeta(1-x)/Power(Pi()*2,1-x)
/* Selected odd integers fast cf Plouffe */
p = Digits()
numeric digits p+2
if x > 2 & x < 22 & Whole(x) then do
if Odd(x) then do
select
when x = 3 then do
a = 180; b = 7; c = 360; d = 0
end
when x = 5 then do
a = 1470; b = 5; c = 3024; d = 84
end
when x = 7 then do
a = 56700; b = 19; c = 113400; d = 0
end
when x = 9 then do
a = 18523890; b = 625; c = 37122624; d = 74844
end
when x = 11 then do
a = 425675250; b = 1453; c = 851350500; d = 0
end
when x = 13 then do
a = 257432175; b = 89; c = 514926720; d = 62370
end
when x = 15 then do
a = 390769879500; b = 13687; c = 781539759000; d = 0
end
when x = 17 then do
a = 1904417007743250; b = 6758333; c = 3808863131673600; d = 29116187100
end
when x = 19 then do
a = 21438612514068750; b = 7708537; c = 42877225028137500; d = 0
end
when x = 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**x; 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
y = (b*Pi()**x-c*s1-d*s2)/a
numeric digits p
return y+0
end
end
/* Standard for all other x */
/* 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-Power(2,1-x)))
/* Series */
s = 0; y = 0
do j = 0 to 2*n-1
b = (-1)**j; c = Power(j+1,x)
if j < n then
s = 0
else
s = s+Comb(n,j-n)
y = y+(s-2**n)*b/c
end
/* Inverse reduction */
y = a*y
numeric digits Digits()-2
return y+0
```

### Numbers

```/* Numbers library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Abundant:
/* Is a number abundant? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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(x,a) = 0 then
return 0
else
return 1
end
/* Cf definition */
if Sigma(x) > 2*x then
return 1
else
return 0

Abundants:
/* Abundant numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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 n = 1 to Words(a)
w = Word(a,n)
if w > x then
leave
glob.abundant.n = w
end
n = n-1; glob.abundant.0 = n
return n
end
/* Check if sum(divisors) > 2*number */
n = 0
do i = 1 to x
s = Sigma(i)
if s > 2*i then do
n = n+1; glob.abundant.n = i
end
end
glob.abundant.0 = n
return n

Ackerman:
/* Ackerman's function */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if x < 1 then
return 'X'
if y < 1 then
return 'X'
/* Fast values */
if x = 0  then
return y+1
if x = 1  then
return y+2
if x = 2  then
return y+y+3
if x = 3  then
return 2**(y+3)-3
if x = 4  then do
if y = 0 then
return 2**(2**2)-3
if y = 1 then
return 2**(2**(2**2))-3
if y = 2 then
return 2**(2**(2**(2**2)))-3
end
/* Cf definition */
if n=0  then
return Ackermann(m-1,1)
else
return Ackermann(m-1,Ackermann(m,n-1))

/* Is a number additive prime? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 101 then do
a = '2 3 5 7 11 23 29 41 43 47 61 67 83 89'
if Wordpos(x,a) = 0 then
return 0
else
return 1
end
/* Cf definition */
return 0

/* Additive prime numbers */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 2 then do
return 0
end
if x < 101 then do
a = '2 3 5 7 11 23 29 41 43 47 61 67 83 89 999'
do n = 1 to Words(a)
w = Word(a,n)
if w > x then
leave
end
n = n-1; glob.additiveprime.0 = n
return n
end
/* Get primes */
p = Primes(x)
/* Collect additive primes */
n = 0
do i = 1 to p
q = glob.prime.i; s = 0
do j = 1 to Length(q)
s = s+Substr(q,j,1)
end
if Prime(s) then do
n = n+1; glob.additiveprime.n = q
end
end
/* Return number of additive primes */
return n

Almostprime:
/* Is a number almost k-prime? */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if x < 1 then
return 'X'
if y < 1 then
return 'X'
/* Fast values */
if y = 1 then
return Prime(x)
if y = 2 then
return Semiprime(x)
/* Cf definition */
if Factors(x) = y then
return 1
else
return 0

Amicable:
/* Amicable number pairs */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Scan for amicable pairs */
n = 0
do i = 1 to x
s = Sigma(i)-i; glob.sigma.i = s
if i = glob.sigma.s then do
if s = i then
iterate
n = n+1
glob.amicable.1.n = s; glob.amicable.2.n = i
end
end
glob.amicable.0 = n
return n

Arithmetic:
/* Is a number arithmetic? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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(x,a) = 0 then
return 0
else
return 1
end
/* Cf definition */
s = Sigma(x)
if Whole(s/glob.divisor.0) then
return 1
else
return 0

Arithmetics:
/* Arithmetic numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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 n = 1 to Words(a)
w = Word(a,n)
if w > x then
leave
glob.Arithmetic.n = w
end
n = n-1; glob.Arithmetic.0 = n
return n
end
/* Check if numer is Arithmetic */
n = 0
do i = 1 to x
if Arithmetic(i) then do
n = n+1; glob.arithmetic.n = i
end
end
glob.arithmetic.0 = n
return n

Comb:
/* Combinations function = binomial coefficients */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if x < 0 then
return 'X'
if y < 0 then
return 'X'
if x < y then
return 'X'
/* Fast values */
if x = y then
return 1
if y = 0 then
return 1
/* Argument reduction */
if x-y < y then
y = x-y
/* In memory? */
if glob.comb.x.y <> '' then
return glob.comb.x.y
numeric digits Digits()+2
/* Calculate combination */
z = 1
do n = (x-y)+1 to x
z = z*n
end
do n = 2 to y
z = z/n
end
numeric digits Digits()-2
glob.comb.x.y = z+0
return glob.comb.x.y

Composite:
/* Is a number composite? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 2 then
return 'X'
/* Fast values */
if x < 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(x,c) = 0 then
return 0
else
return 1
end
if x//2 = 0 then
return 1
if x//3 = 0 then
return 1
if Right(x,1) = 5 then
return 1
/* Method trial division */
do n = 5 by 6 while n*n <= x
if x//n = 0 then
return 1
if x//(n+2) = 0 then
return 1
end
return 0

Composites:
/* Composite numbers = not prime */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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 n = 1 to Words(p)
w = Word(p,n)
if w > x then
leave
glob.composite.n = w
end
n = n-1; glob.composite.0 = n
return n
end
/* Sieve of Eratosthenes */
z. = 0
do i = 2 to x while i*i <= x
do j = i+i by i to x
z.j = 1
end
end
n = 0
do i = 4 to x
if z.i = 1 then do
n = n+1
glob.composite.n = i
end
end
glob.composite.0 = n
return n

Deficient:
/* Is a number deficient? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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(x,d) = 0 then
return 0
else
return 1
end
/* Cf definition */
if Sigma(x) < 2*x then
return 1
else
return 0

Deficients:
/* Deficient numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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 n = 1 to Words(a)
w = Word(a,n)
if w > x then
leave
glob.Deficient.n = w
end
n  = n-1; glob.Deficient.0 = n
return n
end
/* Check if sum(divisors) < 2*number */
n = 0
do i = 1 to x
s = Sigma(i)
if s < 2*i then do
n = n+1
glob.Deficient.n = i
end
end
glob.Deficient.0 = n
return n

Dfact:
/* Double factorial function = n!! */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
/* Fast values */
if x < 2 then
return 1
/* Loops cf definition */
if x//2 = 1 then do
y = 1
do n = 3 to x by 2
y = y*n
end
end
else do
y = 2
do n = 4 to x by 2
y = y*n
end
end
return y

Digitsum:
/* Digitsum function = sum(digits of a number) */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Sum digits */
y = 0
do n = 1 to Length(x)
y = y+Substr(x,n,1)
end
return y

Divisor:
/* Divisor function = sum(Divisors(x)^y) */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if y = '' then
y = 1
numeric digits Digits()+2
/* Get divisors */
d = Divisors(x)
/* Loop cf definition */
select
when y = 0 then
z = d
when y = 1 then do
z = 0
do n = 1 to d
z = z+glob.divisor.n
end
end
when Whole(y) then do
z = 0
do n = 1 to d
z = z+glob.divisor.n**y
end
end
otherwise do
z = 0
do n = 1 to d
z = z+Power(glob.divisor.n,y)
end
end
end
numeric digits Digits()-2
return z+0

Divisors:
/* Divisors of an integer */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Fast values */
if x = 1 then do
glob.divisor.1 = 1; glob.divisor.0 = 1
return 1
end
/* Euclid's method */
a = 1; glob.links.1 = 1; b = 1; glob.rechts.1 = x
m = x//2
do j = 2+m by 1+m to ISqrt(x)
if x//j = 0 then do
a = a+1; glob.links.a = j; b = b+1; glob.rechts.b = x%j
end
end
if j*j = x then do
a = a+1; glob.links.a = j
end
/* Save in globals */
n = 0
do i = 1 to a
n = n+1; glob.divisor.n = glob.links.i
end
do i = b by -1 to 1
n = n+1; glob.divisor.n = glob.rechts.i
end
glob.divisor.0 = n
/* Return number of divisors */
return n

Euclids:
/* Euclid-Mullin numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
if x > 16 then
return 'X'
/* Fast values */
if x = 1 then do
glob.euclid.1 = 2
glob.euclid.0 = 1
return x
end
/* First divisor = lowest prime factor */
p = 2
do i = 2 to x
q = p+1
select
when Even(q) then
f = 2
otherwise do
f = q; j = 3
do while j*j <= q
if q//j = 0 then do
f = j
leave
end
j = j+2
end
end
end
glob.euclid.i = f
p = p*f
end
glob.euclid.0 = x
return x

Even:
/* Is a number even? function */
procedure expose glob.
arg x
/* Validity */
if  \ Whole(x) then
return 'X'
/* Formula */
return \ Abs(x//2)

Fact:
/* Factorial function = n! */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Current in memory? */
if glob.fact.x <> '' then
return glob.fact.x
w = x-1
/* Previous in memory? */
if glob.fact.w = '' then do
/* Loop cf definition */
y = 1
do n = 2 to x
y = y*n
end
glob.fact.x = y
end
else
/* Multiply */
glob.fact.x = glob.fact.w*x
return glob.fact.x

Factors:
/* Prime factors of an integer */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x = 1 then do
glob.factor.0 = 0
return 0
end
if x < 4 then do
glob.factor.1 = x; glob.factor.0 = 1
return 1
end
if Prime(x) then do
glob.factor.1 = x; glob.factor.0 = 1
return 1
end
/* 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)
do while x//p = 0
n = n+1; glob.factor.n = p
x = x%p
end
end
/* Check higher factors */
do j = 29 by 6 while j*j <= x
p = Right(j,1)
if p <> 5 then do
do while x//j = 0
n = n+1; glob.factor.n = j
x = x%j
end
end
if p = 3 then
iterate
y = j+2
do while x//y = 0
n = n+1; glob.factor.n = y
x = x%y
end
end
/* Last factor */
if x > 1 then  do
n = n+1; glob.factor.n = x
end
glob.factor.0 = n
/* Return number of factors */
return n

Fareys:
/* Farey fraction numbers */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Recurring */
glob.farey.1 = '0/1'
n = 1
parse value 0 1 1 x with aa bb cc dd
do while cc <= x
k = Trunc((x+bb)/dd)
parse value cc dd k*cc-aa k*dd-bb with aa bb cc dd
n = n+1
glob.farey.y = aa'/'bb
end
glob.farey.0 = n
/* Number of entries */
return n

Ffact:
/* Falling factorial function */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if x <= 0 then
return 'X'
/* Fast values */
if x = 1 then
return 1
if y = 0 then
return 1
/* Calculate cf definition */
u = 1
do k = 0 to y-1
u = u*(x-k)
end
return u

Gcd:
/* Greatest common divisor function */
procedure expose glob.
arg x,y
x = Abs(x); y = Abs(y)
/* Fast values */
if y = 0 then
return x
/* Euclidian division */
do until y = 0
parse value x//y y with y x
end
return x

Giugas:
/* Giuga numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
if x < 30 then do
glob.giuga.0 = 0
return 0
end
/* Fast values */
g = '30 858 1722 66198 2214408306 24423128562 432749205173838 14737133470010574 ',
||  '550843391309130318 244197000982499715087866346 554079914617070801288578559178 ',
||  '1910667181420507984555759916338506 999999999999999999999999999999999 '
do n = 1 to Words(g)
w = Word(g,n)
if w > x then do
n = n-1
leave
end
glob.giuga.n = w
end
glob.giuga.0 = n
return n

Half:
/* Is a number half integer? function */
procedure expose glob.
arg x
/* Formula */
return (Frac(Abs(x))=0.5)

Kaprekar:
/* Is a number Kaprekar? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 101 then do
k = '1 9 45 55 99'
if Wordpos(x,k) = 0 then
return 0
else
return 1
end
/* Method casting out nines */
p = Max(Digits(),2*Length(x))
numeric digits p
a = x//9
if a > 2 then
return 0
s = x*x
if a = s//9 then do
do i = 1 to Length(s)%2
parse var s l +(i) r
if x = l+r then
return 1
end
end
return 0

Kaprekars:
/* Kaprekar numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 101 then do
k = '1 9 45 55'
do n = 1 to Words(k)
w = Word(k,n)
if w > x then
leave
glob.kaprekar.n = w
end
n = n-1; glob.kaprekar.0 = n
return n
end
x = x+1
/* Ensure enough digits */
p = Max(Digits(),2*Length(x))
numeric digits p
n = 1; glob.kaprekar.n = 1
/* Method casting out nines */
do i = 9 for x-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
n = n+1; glob.kaprekar.n = i
leave
end
end
end
glob.kaprekar.0 = n
return n

Lcm:
/* Least common multiple function */
procedure expose glob.
arg x,y
x = Abs(x); y = Abs(y)
/* Fast values */
if y = 0 then
return 0
/* Euclid */
a = x*y
do until y = 0
parse value x//y y with y x
end
return a%x

Odd:
/* Is a number odd? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
/* Formula */
return Abs(x//2)

Perfect:
/* Is a number perfect? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values (all known) */
p = '6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128 ',
||  '2658455991569831744654692615953842176 ',
||  '191561942608236107294793378084303638130997321548169216 ',
||  '13164036458569648337239753460458722910223472318386943117783728128 ',
||  '14474011154664524427946373126085988481573677491474835889066354349131199152128 '
if Wordpos(x,p) = 0 then
return 0
else
return 1

Perfects:
/* Perfect numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
if x < 6 then do
glob.perfect.0 = 0
return 0
end
/* Fast values */
p = '6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128 ',
||  '2658455991569831744654692615953842176 ',
||  '191561942608236107294793378084303638130997321548169216 ',
||  '13164036458569648337239753460458722910223472318386943117783728128 ',
||  '14474011154664524427946373126085988481573677491474835889066354349131199152128 ',
||  '99999999999999999999999999999999999999999999999999999999999999999999999999999'
do n = 1 to Words(p)
w = Word(p,n)
if w > x then
leave
glob.perfect.n = w
end
n = n-1; glob.perfect.0 = n
return n

Perm:
/* Permutations function */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(y) then
return 'X'
if x < 0 then
return 'X'
if y < 0 then
return 'X'
if x < y then
return 'X'
/* Cf definition */
z = 1
do n = (x-y)+1 to x
z = z*n
end
return z

Phi:
/* Euler's totient function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 3 then
return 1
if x < 5 then
return 2
/* Multiplicative property using factors */
f = Factors(x)+1; glob.factor.f = 0
y = 1; v = glob.factor.1; n = 1
do i = 2 to f
a = glob.factor.i
if a = v then
n = n+1
else do
y = y*v**(n-1)*(v-1)
v = a; n = 1
end
end
return y

Practical:
/* Is a number practical? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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(x,p) = 0 then
return 0
else
return 1
end
if Odd(x) then
return 0
/* Method Srinivasan-Stewart-Sierpinski */
n = Factors(x)
f = 2; s = 2
do i = 2 to n
if glob.factor.i > s+1
then return 0
f = f*glob.factor.i; s = Sigma(f)
end
return 1

Practicals:
/* Practical numbers */
procedure expose glob.
arg x
/* Validity */
if Whole(x) = 0 then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x < 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 n = 1 to Words(p)
w = Word(p,n)
if w > x then
leave
glob.prime.n = w
end
n = n-1; glob.prime.0 = n
return n
end
/* Using function Practical() */
n = 1; glob.practical.1 = 1
do i = 2 by 2 to x
if Practical(i) then do
n = n+1; glob.practical.n = i
end
end
glob.practical.0 = n
return n

Prim:
/* Primorial function = n# */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Fast values */
if x < 2 then
return 1
/* Get primes <= x */
p = Primes(x)
/* Calculate product */
y = 1
numeric digits Digits()+2
do n = 1 to p
y = y*glob.prime.n
end
numeric digits Digits()-2
return y+0

Prime:
/* Is a number prime? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 2 then
return 'X'
/* Low primes also used as deterministic witnesses */
lp = '2 3 5 7 11 13 17 19 23 29 31 37 41'
/* Fast values */
if x < 42 then do
if Wordpos(x,lp) = 0 then
return 0
else
return 1
end
if x//2 = 0 then
return 0
if x//3 = 0 then
return 0
if Right(x,1) = 5 then
return 0
/* Miller-Rabin primality test */
numeric digits 2*Length(x)
d = x-1; e = d
/* Reduce x-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 x < Word(th,w) then do
do k = 1 to w
if x < 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 k = 1 to k
a = Word(lp,k); y = Powermod(a,d,x)
if y = 1 then
iterate
if y = e then
iterate
do s
y = (y*y)//x
if y = 1 then
return 0
if y = e then
leave
end
if y <> e then
return 0
end
return 1

Primes:
/* Prime numbers */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x = 1 then do
glob.prime.0 = 0
return 0
end
if x < 101 then do
p = '2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 999'
do n = 1 to Words(p)
w = Word(p,n)
if w > x then
leave
glob.prime.n = w
end
n = n-1; glob.prime.0 = n
return n
end
/* Wheeled sieve of Eratosthenes */
do i = 3 by 2 to x while i*i <= x
if glob.priem.i = '' then do
do j = i*3 by i+i to x
glob.priem.j = 0
end
end
end
/* Save results */
n = 1; glob.prime.1 = 2
do i = 3 by 2 to x
if glob.priem.i = '' then do
n = n+1; glob.prime.n = i
end
end
glob.prime.0 = n
return n

/* Radical function = product of unique prime factors */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Get unique factors */
n = Ufactors(x)
/* Calculate product */
y = 1
do i = 1 to n
y = y*glob.ufactor.i
end
return y

Rfact:
/* Rising factorials function */
procedure expose glob.
arg x,y
/* Validity */
if \ Whole(y) then
return 'X'
if x = 0 then
return 'X'
if y < 0 then
return 'X'
/* Fast values */
if x = 1 then
return 1
if y = 0 then
return 1
/* Cf definition */
u = 1
do k = 0 to y-1
u = u*(x+k)
end
return u

Semiprime:
/* Is a number semi prime? function */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
if x < 4 then
return 0
/* 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 x < 101 then do
if Wordpos(x,s) = 0 then
return 0
else
return 1
end
/* Wheeled scan */
do i = 2 for 2
if x//i = 0 then
if Prime(x%i) then
return 1
else
return 0
end
do i = 5 by 6 until i*i > x
do j = i by 2 for 2
if x//j==0 then
if Prime(x%j) then
return 1
else
return 0
end
end
return 0

Sigma:
/* Sigma function = Sum of all divisors of x including 1 and x */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x = 1 then do
glob.divisor.0 = 1
return 1
end
/* Euclid's method */
m = x//2; y = 1+x; n = 2
do j = 2+m by 1+m to ISqrt(x)
if x//j = 0 then do
y = y+j+x%j; n = n+2
end
end
if j*j = x then do
y = y+j; n = n+1
end
/* Store number of divisors */
glob.divisor.0 = n
/* Return sum */
return y

Ufactors:
/* Unique prime factors of an integer */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Fast values */
if x = 1 then do
glob.ufactor.0 = 0
return 0
end
if x < 4 then do
glob.ufactor.1 = x; glob.ufactor.0 = 1
return 1
end
if Prime(x) then do
glob.ufactor.1 = x; glob.ufactor.0 = 1
return 1
end
/* Factors has been called */
if glob.factor.0 <> '' then do
n = 0; v = 0
do i = 1 to glob.factor.0
if glob.factor.i <> v then do
n = n+1; glob.ufactor.n = glob.factor.i
v = glob.factor.i
end
end
glob.ufactor.0 = n
return n
end
/* Check low factors */
n = 0; v = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i)
do while x//p = 0
if p <> v then do
n = n+1; glob.ufactor.n = p
v = p
end
x = x%p
end
end
/* Check higher factors */
do j = 29 by 6 while j*j <= x
p = Right(j,1)
if p <> 5 then do
do while x//j = 0
if j <> v then do
n = n+1; glob.ufactor.n = j
v = j
end
x = x%j
end
end
if p = 3 then
iterate
y = j+2
do while x//y = 0
if y <> v then do
n = n+1; glob.ufactor.n = y
v = y
end
x = x%y
end
end
/* Last factor */
if x > 1 then  do
if x <> v then do
n = n+1; glob.ufactor.n = x
end
end
glob.ufactor.0 = n
/* Return number of factors */
return n

Whole:
/* Is a number integer? function */
procedure expose glob.
arg x
/* Formula */
return Datatype(x,'w')
```

### Polynomial

```/* Polynomial library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Pdiff:
/* Polynomial differentiation function */
procedure expose glob.
arg p
d = ''; w = Words(p)
do n = 1 to w-1
d = d (w-n)*Word(p,n)
end
return d

Peval:
/* Polynomial evaluation function */
procedure expose glob.
arg p,x
y = 0; w = Words(p)
numeric digits Digits()+2
do n = 1 to w
y = y*x+Word(p,n)
end
numeric digits Digits()-2
return y+0
```

### Rational

```/* Rational library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Rabs:
/* Rational absolute value function */
procedure expose glob.
arg x
parse var x nx '/' dx
nx = abs(nx); dx = abs(dx)
return nx'/'dx

/* Rational add function */
procedure expose glob.
arg x,y
parse var x nx '/' dx; parse var y ny '/' dy
return nx*dy+ny*dx'/'dx*dy

Rdivide:
/* Rational divide function */
procedure expose glob.
arg x,y
parse var x nx '/' dx; parse var y ny '/' dy
return nx*dy'/'dx*ny

Rinvert:
/* Rational invert function */
procedure expose glob.
arg x
parse var x nx '/' dx
return dx'/'nx

Rmultiply:
/* Rational multiply function */
procedure expose glob.
arg x,y
parse var x nx '/' dx; parse var y ny '/' dy
return nx*ny'/'dx*dy

Rnegate:
/* Rational negate function */
procedure expose glob.
arg x
parse var x nx '/' dx
return -nx'/'dx

Rpower:
/* Rational power function */
procedure expose glob.
arg x,y
parse var x nx '/' dx
return nx**y'/'dx**y

Rrational:
/* Integer->Rational function */
procedure expose glob.
arg x
in = x; 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(x)
/* Calculate numerator and denominator */
z = tnum; tnum = n*tnum+num; num = z
z = tden; tden = n*tden+den; den = z
/* Check precision */
if n = x | tnum/tden = in then do
if tnum > h | tden > h then
leave
/* Possible final values */
num = tnum; den = tden
leave
end
/* Next round */
x = 1/(x-n)
end
return num'/'den

Rsimplify:
/* Rational simplify function */
procedure expose glob.
arg x
parse var x nx '/' dx
g = gcd(nx,dx)
return nx/g'/'dx/g

Rsubtract:
/* Rational subtract function */
procedure expose glob.
arg x,y
parse var x nx '/' dx; parse var y ny '/' dy
return nx*dy-ny*dx'/'dx*dy
```

### Roots

```/* Roots library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Cbeq:
/* Cubic roots of ax^3+bx^2+cx+d = 0 */
procedure expose glob.
arg x
/* Validity */
if Words(x) <> 4 then
return 'X'
if Word(x,1) = 0 then
return 'X'
numeric digits Digits()*2
/* Coefficients */
a = Word(x,1); b = Word(x,2); c = Word(x,3); d = Word(x,4)
/* Method Cardano */
/* Reduce to x^3+bx+c = 0 */
t = 1/a; p = b*t; q = c*t; r = d*t
/* Discriminant */
a = q-p*p/3; b = r-p*q/3+2*p*p*p/27
d = a*a*a/27+b*b/4
if Abs(d) < 10**-(Digits()-2) then
d = 0
/* Roots */
select
when d < 0 then do
t = Sqrt(-a/3); u = ACos(Sqrt(-27/(a*a*a))*-b/2)
n = 3; glob.cbeq.1 = 2*t*Cos(u/3); glob.cbeq.2 = 2*t*Cos((u+2*Pi())/3); glob.cbeq.3 = 2*t*Cos((u+4*Pi())/3)
end
when d = 0 then do
t = Cbrt(-b/2)
n = 2; glob.cbeq.1 = 2*t; glob.cbeq.2 = -t
end
when d > 0 then do
t = -b/2; u = Sqrt(d)
n = 1; glob.cbeq.1 = Cbrt(t+u)+Cbrt(t-u)
end
end
glob.cbeq.0 = n
/* Undo reduction */
do i = 1 to n
glob.cbeq.i = glob.cbeq.i-p/3
end
numeric digits Digits()/2
/* Normalize */
do i = 1 to n
glob.cbeq.i = glob.cbeq.i+0
end
/* Return number of roots */
return n

Qteq:
/* Quartic roots of ax^4+bx^3+cx^2+dx+e = 0 */
procedure expose glob.
arg x
/* Validity */
if Words(x) <> 5 then
return 'X'
if Word(x,1) = 0 then
return 'X'
numeric digits Digits()+2
/* Method Ferrari */
/* Reduce to x^4+ax^3+bx^2+cx+d = 0 */
t = Word(x,1); a1 = Word(x,2)/t; b1 = Word(x,3)/t; c1 = Word(x,4)/t; d1 = Word(x,5)/t
/* Build and solve cubic equation x^3+ax^2+bx+c = 0 */
a2 = -b1; b2 = a1*c1-4*d1; c2 = d1*(4*b1-a1*a1)-c1*c1; p = 1 a2 b2 c2; call Cbeq(p)
/* Get root with highest absolute value */
y = Abs(glob.cbeq.1)
do n = 2 to glob.cbeq.0
a  = Abs(glob.cbeq.n)
if a > y then
y = a
end
/* Calculate additional values */
a = a1/2; b = y/2; d = Sqrt(b*b-d1)
if d = 0 then
c = Sqrt(a*a-b1+y)
else
c = (a*b-c1/2)/d
/* Build and solve first quadratic equation */
e = a+c; f = b+d; p = 1 e f; call Sqeq(p)
n = 0
do i = 1 to glob.sqeq.0
n = n+1; glob.qteq.n = glob.sqeq.i
end
/* Build and solve second quadratic equation */
e = a-c; f = b-d; p = 1 e f; call Sqeq(p)
do i = 1 to glob.sqeq.0
n = n+1; glob.qteq.n = glob.sqeq.i
end
glob.qteq.0 = n
numeric digits Digits()-2
/* Normalize */
do i = 1 to n
glob.qteq.n = glob.qteq.n+0
end
/* Return number of roots */
return n

Sqeq:
/* Square roots of ax^2+bx+c = 0 */
procedure expose glob.
arg x
/* Validity */
if Words(x) <> 3 then
return 'X'
if Word(x,1) = 0 then
return 'X'
numeric digits Digits()+2
/* Coefficients */
a = Word(x,1); b = Word(x,2); c = Word(x,3)
/* Discriminant */
d = b*b-4*a*c
/* Abc formula */
select
when d < 0 then do
n = 0
end
when d = 0 then do
n = 1; glob.sqeq.1 = -0.5*b/a
end
when d > 0 then do
if b = 0 then
s = 1
else
s = Sign(b)
q = -0.5*(b+s*Sqrt(d))
n = 2; glob.sqeq.1 = q/a; glob.sqeq.2 = c/q
end
end
glob.sqeq.0 = n
numeric digits Digits()-2
/* Normalize */
do i = 1 to n
glob.sqeq.i = glob.sqeq.i+0
end
/* Number of roots */
return n
```

### Sequences

```/* Sequences library 28 Jul 2024 */
/* (C) Paul van den Eertwegh 2023-2024 */

Alcuin:
/* Alcuin sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Special values */
if x < 4 then
return 0
/* Formula */
if Even(x) then
return Round(x*x/48)
else
return Round((x+3)**2/48)

Bell:
/* Bell sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
if x = 1 then
return 1
/* Recurring */
numeric digits Digits()+2
if x < 36 then do
y = 0
do k = 0 to x
y = y+Stirling2(x,k)
end
end
else do
y. = 0; y.0 = 1; y.1 = 1
do n = 2 to x
s = 0
do k = 0 to n-1
s = s+Comb(n-1,k)*y.k
end
y.n = s
end
y = y.x
end
numeric digits Digits()-2
return y+0

Beraha:
/* Beraha sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 1 then
return 4
if x = 2 then
return 0
if x = 3 then
return 1
if x = 4 then
return 2
if x = 6 then
return 3
if x = 8 then
return 2+Sqrt(2)
if x = 10 then
return 0.5*(5+Sqrt(5))
/* Formula */
return 2+2*Cos(2*Pi()/x)

Bernoulli:
/* Bernoulli sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
if x = 1 then
return 0.5
if x//2 = 1 then
return 0
p = Digits()
/* In memory? */
if glob.bernoulli.x.p <> '' then
return glob.bernoulli.x.p
/* Take extra digits */
numeric digits p+x
/* Recurring method Akiyama-Tanigawa */
do m = 0 by 1 to x
b.m = 1/(m+1)
do j = m by -1 to 1
k = j-1; b.k = j*(b.k-b.j)
end
end
numeric digits p
/* Result */
glob.bernoulli.x.p = b.0+0
return glob.bernoulli.x.p

Cake:
/* Cake sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Formula */
return (x*x*x+5*x+6)/6

Catalan2:
/* Catalan sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
numeric digits digits()+2
y = 1
do n = 1 to x
y = y*(4*n-2)/(n+1)
end
numeric digits digits()-2
return y+0

Caterer:
/* Lazy caterer's sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Formula */
return (x*x+x+2)/2

Central:
/* Central binomial coefficient sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Formula */
return Comb(2*x,x)

Conway1:
/* Conway look-and-say sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Series cf definition */
s = '1'
do n = 1 to x
ns = ''; c = 1; m = Length(s)
do j = 1 to m
if j = m | Substr(s,j,1) <> Substr(s,j+1,1) then do
ns = ns||c||Substr(s,j,1); c = 1
end
else
c = c+1
end
s = ns
end
return s

Cullen:
/* Cullen sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Formula */
return x*2**x+1

Derangement:
/* Derangement sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Recurring cf definition */
numeric digits Digits()+2
a = 1; b = 0; y = 0
do n = 2 to x
y = (a+b)*(n-1)
a = b; b = y
end
numeric digits Digits()-2
return y+0

Favard:
/* Favard sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
if x = 1 then
return Pi()*0.5
if x = 2 then
return Pi()**2*0.125
if x = 3 then
return Pi()**3/24
if x = 4 then
return 5*Pi()**4/384
if x = 5 then
return Pi()**5/240
if x = 6 then
return 61*Pi()**6/46080
if x = 7 then
return 17*Pi()**7/40320
if x = 8 then
return 277*Pi()**8/2064384
if x = 9 then
return 31*Pi()**9/725760
if x = 10 then
return 50521*Pi()**10/3715891200
if x = 11 then
return 691*Pi()**11/159667200
if x = 12 then
return 540553*Pi()**12/392398110720
numeric digits Digits()+2
/* Recurring */
f.0 = 1; f.1 = Pi()*0.5
do j = 2 to x
s = 0
do k = 1 to Floor(j*0.5)
a = 2*k; b = a-1; c = j-a
s = s+f.b*f.c
end
f.j = s*Pi()/(2*j)
end
numeric digits Digits()-2
return f.x+0

Fermat:
/* Fermat sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
numeric digits Digits()+2
/* Formula */
y = 2**(2**x)
numeric digits Digits()-2
return y+0

Fibonacci:
/* Fibonacci sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 2 then
return x
/* Recurring */
a = 0; b = 1
do n = 2 to x
y = a+b; a = b; b = y
end
return y+0

Golomb2:
/* Golomb-Silverman sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Special values */
if x < 1 then
return 1
/* Recurring */
numeric digits Digits()+2
y. = 0; y.1 = 1
do n = 2 to x
a = n-1; b = y.a; c = y.b; d = n-c
y.n = 1+y.d
end
numeric digits Digits()-2
return y.x+0

Gould:
/* Gould sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
/* Recurring */
numeric digits Digits()+2
y. = 0; y.0 = 1
do n = 1 to x
if Even(n) then do
a = n%2; y = y.a
end
else do
a = n%2; y = 2*y.a
end
y.n = y
end
numeric digits Digits()-2
return y.x+0

Gregory1:
/* Gregory sequence 1 */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Special values */
if x = 1 then
return 0.5
/* Recurring */
numeric digits Digits()+2
y.1 = 0.5
if x = 1 then
return y.1
do n = 2 to x
s = 0
do k = 1 to n-1
s = s+y.k/(n+1-k)
end
y.n = -s+1/(n+1)
end
numeric digits Digits()-2
return (-1)**(x-1)*y.x

Gregory2:
/* Gregory sequence 2 */
procedure expose glob.
arg x
/* Validity */
if x < 1 then
return 'X'
/* Special values */
if x = 1 then
return 0.25*Pi()
/* Formula */
return ATan(1/x)

Harmonic:
/* Harmonic sequence */
procedure expose glob.
arg x,y
if y = '' then
y = 1
/* Validity */
if \ Whole(x) then
return 'X'
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
if y < 0 then
return 'X'
/* Special values */
if x = 0 then
return 0
if x < 2 then
return 1
numeric digits Digits()+2
/* From definition */
z = 0
if y = 1 then do
do n = 1 to x
z = z+1/n
end
end
else do
do n = 1 to x
z = z+1/n**y
end
end
numeric digits Digits()-2
return z+0

Jacobsthal:
/* Jacobsthal sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 3 then
return 1
numeric digits Digits()+2
/* Formula */
y = (2**x-(-1)**x)/3
numeric digits Digits()-2
return y+0

Lah:
/* Lah sequence */
procedure expose glob.
arg n,k
/* Validity */
if \ Whole(n) then
return 'X'
if \ Whole(k) then
return 'X'
if n < 1 then
return 'X'
if k < 1 then
return 'X'
if n < k then
return 'X'
/* Special values */
if n = k then
return 1
/* Formulas */
s = (-1)**n
if k = 1 then
return s*Fact(n)
else
return s*Comb(n-1,k-1)*Fact(n)/Fact(k)

Lebesgue1:
/* Lebesgue sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
if x = 1 then
return 1/3+2*Sqrt(3)/Pi()
if x = 2 then
return 0.2+Sqrt(25-2*Sqrt(5))/Pi()
if x = 3 then
return 1/7+(22*Sin(Pi()/7)-2*Cos(Pi()/14)+10*Cos(3*Pi()/14))/(3*Pi())
if x = 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 */
y = 0
do n = 1 to x
y = y+Tan(n*Pi()/(2*x+1))/n
end
y = 1/(2*x+1)+2*y/Pi()
numeric digits Digits()-2
return y+0

Leonardo:
/* Leonardo sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Special values */
if x < 3 then
return 1
/* From definition */
a = 1; b = 1
do n = 3 to x
y = a+b+1
a = b; b = y
end
return y+0

Lucas:
/* Lucas sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 1 then
return 2
if x = 2 then
return 1
/* Recurring */
numeric digits Digits()+2
a = 2; b = 1; y = 0
do n = 2 to x
y = a+b
a = b; b = y
end
numeric digits Digits()-2
return y+0

Magic:
/* Magic sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Formula */
return (x*x*x+x)/2

Metallic:
/* Metallic ratio sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
if x = 1 then
return Golden()
if x = 2 then
return Silver()
if x = 3 then
return Bronze()
/* Formula */
return 0.5*(x+Sqrt(x*x+4))

Motzkin:
/* Motzkin sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 3 then
return 1
numeric digits Digits()+2
/* Recurring */
a = 1; b = 1; y = 0
do n = 2 to x
y = ((3*n-3)*a+(2*n+1)*b)/(n+2)
a = b; b = y
end
numeric digits Digits()-2
return y+0

Moebius:
/* Moebius sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Using # of (unique) prime factors */
call Factors(x)
call Ufactors(x)
if glob.factor.0 = glob.ufactor.0 then
if Even(glob.factor.0) then
return 1
else
return -1
else
return 0

Narayana:
/* Narayana cows sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 3 then
return 1
/* Recurring */
numeric digits Digits()+2
a = 1; b = 1; c = 1; y = 0
do n = 3 to x
y = a+c
a = b; b = c; c = y
end
numeric digits Digits()-2
return y+0

/* Padovan sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 3 then
return 1
numeric digits Digits()+2
/* Recurring */
a = 1; b = 1; c = 1; y = 0
do n = 3 to x
y = a+b
a = b; b = c; c = y
end
numeric digits Digits()-2
return y+0

Partition:
/* Partition sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 1
numeric digits Digits()+2
/* Recurring */
p. = 0; p.0 = 1
do n = 1 to x
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)*(p.a+p.b)
end
p.n = s
end
numeric digits Digits()-2
return p.x+0

Pell1:
/* Pell sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 2 then
return x
numeric digits Digits()+2
/* Recurring */
a = 0; b = 1; y = 0
do n = 2 to x
y = a+2*b
a = b; b = y
end
numeric digits Digits()-2
return y+0

Pell2:
/* Pell-Lucas companion sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 2 then
return 2
/* Recurring */
numeric digits Digits()+2
a = 2; b = 2; y = 0
do n = 2 to x
y = a+2*b
a = b; b = y
end
numeric digits Digits()-2
return y+0

Perrin:
/* Perrin sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 3
if x = 1 then
return 0
if x = 2 then
return 2
numeric digits Digits()+2
/* Recurring */
a = 3; b = 0; c = 2; y = 0
do n = 3 to x
y = a+b
a = b; b = c; c = y
end
numeric digits Digits()-2
return y+0

Pronic:
/* Pronic sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Formula */
numeric digits Digits()+2
y = x*x+x
numeric digits Digits()-2
return y+0

Sorting:
/* Sorting sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 3
if x = 1 then
return 0
if x = 2 then
return 2
/* Recurring */
y = 0; i = x-1; z = 1
do while i >= 0
y  = y+i; i = i-z; z = z+z
end
return y+0

Recaman:
/* Recaman sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x = 0 then
return 0
/* Recurring cf definition */
r. = 0; r.0 = 0; s. = 0; s.0 = 1
do n = 1 to x
m = n-1; a = r.m-n
if a > 0 & s.a = 0 then do
r.n = a; s.a = 1
end
else do
b = r.m+n
r.n = b; s.b = 1
end
end
return r.x

Stirling1:
/* Stirling first kind sequence */
procedure expose glob.
arg n,k
/* Validity */
if \ Whole(n) then
return 'X'
if \ Whole(k) then
return 'X'
if n < 0 then
return 'X'
if k < 0 then
return 'X'
/* Special values */
if k = n then
return 1
if k > n then
return 0
if k = 0 then
return 0
/* Formulas */
s = (-1)**(n-k)
if k = 1 then
return s*Fact(n-1)
if k = 2 then
return s*Round(Fact(n-1)*Harmonic(n-1))
if k = 3 then
return s*Round(Fact(n-1)*(Harmonic(n-1)**2-Harmonic(n-1,2))/2)
if k = n-3 then
return s*Comb(n,2)*Comb(n,4)
if k = n-2 then
return s*Round((3*n-1)*Comb(n,3)/4)
if k = n-1 then
return s*Comb(n,2)
numeric digits Digits()+2
/* Definition */
c. = 0; c.1 = 1
do i = 1 to n-1
a = i
do j = i+1 to 1 by -1
b = j-1; c.j = c.b+a*c.j
end
end
numeric digits Digits()-2
return s*c.k

Stirling2:
/* Stirling second kind sequence */
procedure expose glob.
arg n,k
/* Validity */
if \ Whole(n) then
return 'X'
if \ Whole(k) then
return 'X'
if n < 0 then
return 'X'
if k < 0 then
return 'X'
if n < k then
return 'X'
/* Special values */
if n = 0 then
return 0
if k = 0 then
return 0
if n = k then
return 1
if k = 1 then
return 1
if k = 2 then
return 2**(n-1)-1
if k = 3 then
return (3**n-3*2**n+3)/6
if k = n-1 then
return Comb(n,2)
numeric digits Digits()+2
/* Definition */
s. = 1; x = n-k
do i = 2 to k
do j = 1 to x
a = j-1; s.j = s.j+i*s.a
end
end
numeric digits Digits()-2
return s.x+0

Sylvester:
/* Sylvester sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
if x > 33 then
return 'X'
numeric digits Digits()+2
/* Definition */
y = 2
do n = 2 to x
y = y*y-y+1
end
numeric digits Digits()-2
return y+0

Tetranacci1:
/* Tetranacci sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 4 then
return 0
/* Recurring */
a = 0; b = 0; c = 0; d = 1
do n = 3 to x
y = a+b+c; a = b; b = c; c = d; d = y
end
return y+0

Thue:
/* Thue-Morse numbers */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Concatenate 2s complements */
y = 0
do i = 1 to x
y = y||translate(y,10,01)
end
return y

Tribonacci1:
/* Tribonacci sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 3 then
return 0
/* Recurring */
a = 0; b = 0; c = 1
do n = 3 to x
y = a+b+c; a = b; b = c; c = y
end
return y+0

Ulam:
/* Ulam sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Special values */
if x < 3 then
return x
/* Series */
numeric digits Digits()+2
c.1 = 1; c.2 = 2; d = 2; e.= 0; e.1 = 1; e.2 = 1; z = 3
do until d = x
n = 0
do j = 1 to d
b = z-c.j
if e.b then do
if c.j <> b then do
n = n+1
if n > 2 then
leave j
end
end
end
if n = 2 then do
d = d+1; c.d = z; e.z = 1
end
z = z+1
end
y = c.d
numeric digits Digits()-2
return y+0

Wedderburn:
/* Wdderburn sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 0 then
return 'X'
/* Special values */
if x < 2 then
return x
/* Recurring */
numeric digits Digits()+2
y. = 0; y.1 = 1
do n = 2 to x
if Even(n) then do
a = n%2; s = y.a*(y.a+1)/2
end
else do
a = n%2+1; s = 0
end
do i = 1 to a-1
j = n-i; s = s+y.i*y.j
end
y.n = s
end
numeric digits Digits()-2
return y.x+0

Woodall:
/* Woodall sequence */
procedure expose glob.
arg x
/* Validity */
if \ Whole(x) then
return 'X'
if x < 1 then
return 'X'
/* Formula */
numeric digits Digits()+2
y = x*2**x-1
numeric digits Digits()-2
return y+0

Zigzag:
/* Zigzag or alternating permutations sequence */
procedure expose glob.
arg x
/* Validity */
if x < 0 then
return 'X'
if \ Whole(x) then
return 'X'
/* Special values */
if x = 0 then
return 1
if x = 1 then
return 1
/* Recurring */
numeric digits Digits()+2
z. = 0; z.0 = 1; z.1 = 1
do n = 2 to x
s = 0
do k = 0 to n-1
a = Comb(n-1,k); b = n-k-1
s = s+a*z.k*z.b
end
z.n = s/2
end
y = z.x
numeric digits Digits()-2
return y+0
```