Mathematics.rex
Introduction
In 2024 I started posting on RosettaCode and found the copy-and-paste way to get a working program annoying. Therefore I'll gather all procedures I have until now and present them here in a set modules, stand-alone files and documentation. All files are contained in below self-extractor. Follow the instructions as given in the code and obtain the full rexx math package.
After running the self-extractor, you'll find some useful text files (among others):
Modules.txt lists all modules.
Labels.txt lists all procedures per module.
See documentation for help and tips.
Self-extractor
-- Math.rex - Mathemetics toolbox self-extractor
-- REXX Modules, snippets, programs and documentation
-- Copy and save this self-extractor in any folder. Name it as
-- you like. Run it with Regina or ooRexx. Except for the
-- version info, it will run silently unless you specify
-- parameter V (= verbose). It will create bunch of programs
-- and includes. Existing ones will be overwritten without
-- warning. Run.rex is the preprocessor and Math.inc contains
-- all math modules.
-- Modules.txt (module list)
-- Labels.txt (procedure and routine list)
-- Comment.txt (extract of all program comments)
-- Stems.txt (global stem list)
-- Crossref.txt (cross reference global stems and procedures)
-- [regina|rexx] run demo1: shows include only required modules
-- [regina|rexx] run demo2: shows include all modules in 1 file
-- Model1.rex: may be used for your own playground (recommended)
Main:
-- Main program flow
say 'MATH SELF-EXTRACTOR - 11 Apr 2025'; say
call Time('r')
parse upper arg v
verbose = (v = 'V')
call Extract
call Documentation
call Stems
call Crossreference
if verbose then
say Time('e')/1 'seconds'
exit
Extract:
-- Source generator
if verbose then
say 'Generate source code...'
parse source . . program
call Stream program,'c','open read'
math = 'Math.inc'
call Stream math,'c','open write replace'
tot = 0; module = ''
do while Lines(program) > 0
prog = Linein(program)
if Left(prog,10) = '-- Module ' then do
if module <> '' then
call Stream module,'c','close'
parse var prog . . module .
if verbose then
say Substr(prog,4,length(prog)-3)
call Stream module,'c','open write replace'
end
if module <> '' then do
if Pos('.inc',module) > 0 then
call Lineout math,prog
call Lineout module,prog
tot = tot+1
end
end
call Stream module,'c','close'
call Stream math,'c','close'
call Stream program,'c','close'
if verbose then do
say tot 'lines'
say 'Done'
say
end
return
Documentation:
-- Documentation generator
if verbose then
say 'Generate source documentation...'
modules = 'Modules.txt'
call Stream modules,'c','open write replace'
labels = 'Labels.txt'
call Stream labels,'c','open write replace'
comment = 'Comment.txt'
call Stream comment,'c','open write replace'
tot = 0
call lineout modules,'LIST OF MODULES'
call lineout modules,''
call lineout labels,'LIST OF PROCEDURES PER MODULE'
call lineout labels,''
call lineout comment,'LIST OF COMMENTS PER PROCEDURE'
call lineout comment,''
tot = tot+6
call Stream program,'c','open read'
do while Lines(program) > 0
prog = Linein(program)
select
when left(prog,10) = '-- Module ' then do
s = substr(prog,4)
parse var s . module .
prog = Linein(program)
s = substr(prog,4)
module = left('Module' module,22) s
call lineout modules,module
call lineout labels,''
call lineout labels,module
call lineout comment,''
call lineout comment,''
call lineout comment,module
tot = tot+6
end
when left(prog,2) = '--' then do
s = substr(prog,4)
call lineout comment,s
tot = tot+1
end
when right(prog,1) = ':' then do
label1 = left(prog,length(prog)-1)
prog = Linein(program)
s = substr(prog,4)
label2 = left('Label' label1,22) s
call lineout labels,label2
call lineout comment,''
call lineout comment,label2
tot = tot+3
end
otherwise
nop
end
end
call Stream program,'c','close'
call Stream modules,'c','close'
call Stream labels,'c','close'
call Stream comment,'c','close'
if verbose then do
say tot 'lines'
say 'Done'
say
end
return
Stems:
-- Global stems list
if verbose then
say 'Generate global stem list...'
stems = 'Stems.txt'
call Stream stems,'c','open write replace'
tot = 0
call lineout stems,'LIST OF GLOBAL STEMS'
call lineout stems,''
call lineout stems,'abun. Abundant numbers'
call lineout stems,'addi. Additive primes'
call lineout stems,'amic. Amicable numbers'
call lineout stems,'arit. Arithmetic numbers'
call lineout stems,'bern. Bernoulli numbers'
call lineout stems,'carm. Carmichael strong 3 primes'
call lineout stems,'cata. Catalan numbers'
call lineout stems,'cbrt. Cubic roots'
call lineout stems,'chow. Chowla numbers'
call lineout stems,'comb. Combinations'
call lineout stems,'comp. Composite numbers'
call lineout stems,'cuba. Cuban primes'
call lineout stems,'cubi. Cubic primes'
call lineout stems,'cycl. Cyclop numbers'
call lineout stems,'defi. Deficient numbers'
call lineout stems,'divi. Divisors'
call lineout stems,'efac. Exponent factors'
call lineout stems,'emir. Emirp primes'
call lineout stems,'erdo. Erdos numbers'
call lineout stems,'eucl. Euclid numbers'
call lineout stems,'extr. Extra primes'
call lineout stems,'fact. Factors'
call lineout stems,'fibo. Fibonacci numbers'
call lineout stems,'flag. Flags'
call lineout stems,'frob. Frobenius primes'
call lineout stems,'giug. Giuga numbers'
call lineout stems,'glob. Global variables'
call lineout stems,'hamm. Hamming numbers'
call lineout stems,'humb. Humble numbers'
call lineout stems,'high. High composite numbers'
call lineout stems,'kapr. Kaprekar numbers'
call lineout stems,'labl. Labels'
call lineout stems,'lirt. Linear root'
call lineout stems,'ludi. Ludic numbers'
call lineout stems,'magi. Magic numbers'
call lineout stems,'magn. Magmanimous numbers'
call lineout stems,'motz. Motzkin numbers'
call lineout stems,'perf. Perfect numbers'
call lineout stems,'poly. Polynomial coefficients'
call lineout stems,'prac. Practical numbers'
call lineout stems,'prim. Prime numbers'
call lineout stems,'prmo. Primorial numbers'
call lineout stems,'qtrt. Quartic roots'
call lineout stems,'quad. Quadratic primes'
call lineout stems,'reve. Reversable primes'
call lineout stems,'smoo. Smooth numbers'
call lineout stems,'sqrt. Square roots'
call lineout stems,'squa. Square free numbers'
call lineout stems,'toti. Totient numbers'
call lineout stems,'ufac. Unique factors'
tot = tot+43
call Stream stems,'c','close'
if verbose then do
say tot 'lines'
say 'Done'
say
end
return
Crossreference:
-- Stem cross reference lists
if verbose then
say 'Generate stem cross references...'
call Stream program,'c','open read'
labl. = 0; n = 0
do while Lines(program) > 0
prog = Linein(program)
select
when left(prog,2) = '--' then
nop
when right(prog,1) = ':' then do
label = left(prog,length(prog)-1)
end
when left(prog,10) = 'procedure ' then do
if pos('expose',prog) > 0 then do
x = substr(prog,18)
do while x <> ''
parse var x expo x
n = n+1; labl.proc.n = label; labl.stem.n = expo
end
end
end
otherwise
nop
end
end
labl.0 = n
call Stream program,'c','close'
tot = 0
crossref = 'Crossref.txt'
call Stream crossref,'c','open write replace'
call Sortproc
call lineout crossref,'STEMS BY PROCEDURE'
head1 = left('Procedure',15) left('Stem',5)
head2 = copies('-',length(head1))
call lineout crossref,head2
call lineout crossref,head1
call lineout crossref,head2
do i = 1 to n
call lineout crossref,left(labl.proc.i,15) labl.stem.i
tot = tot+1
end
call lineout crossref,head2
call lineout crossref,''
tot = tot+6
call Sortstem
call lineout crossref,'PROCEDURES BY STEM'
head1 = left('Stem',5) left('Procedure',15)
call lineout crossref,head2
call lineout crossref,head1
call lineout crossref,head2
do i = 1 to n
call lineout crossref,labl.stem.i left(labl.proc.i,15)
tot = tot+1
end
call lineout crossref,head2
tot = tot+5
call Stream crossref,'c','close'
if verbose then do
say tot 'lines'
say 'Done'
say
end
return
Sortproc:
-- Sort by procedure and stem
procedure expose labl.
n = labl.0; s = 1; sl.1 = 1; sr.1 = n
do until s = 0
l = sl.s; r = sr.s; s = s-1
do until l >= r
m = (l+r)%2; i = l; j = r
a = labl.proc.m; b = labl.stem.m
do until i > j
do i = i
if labl.proc.i > a then
leave
if labl.proc.i = a then
if labl.stem.i >= b then
leave
end
do j = j by -1
if labl.proc.j < a then
leave
if labl.proc.j = a then
if labl.stem.j <= b then
leave
end
if i <= j then do
t = labl.proc.i; labl.proc.i = labl.proc.j; labl.proc.j = t
t = labl.stem.i; labl.stem.i = labl.stem.j; labl.stem.j = t
i = i+1; j = j-1
end
end
if j-l < r-i then do
if i < r then do
s = s+1; sl.s = i; sr.s = r
end
r = j
end
else do
if l < j then do
s = s+1; sl.s = l; sr.s = j
end
l = i
end
end
end
return
Sortstem:
-- Sort by stem and procedure
procedure expose labl.
n = labl.0; s = 1; sl.1 = 1; sr.1 = n
do until s = 0
l = sl.s; r = sr.s; s = s-1
do until l >= r
m = (l+r)%2; i = l; j = r
a = labl.stem.m; b = labl.proc.m
do until i > j
do i = i
if labl.stem.i > a then
leave
if labl.stem.i = a then
if labl.proc.i >= b then
leave
end
do j = j by -1
if labl.stem.j < a then
leave
if labl.stem.j = a then
if labl.proc.j <= b then
leave
end
if i <= j then do
t = labl.stem.i; labl.stem.i = labl.stem.j; labl.stem.j = t
t = labl.proc.i; labl.proc.i = labl.proc.j; labl.proc.j = t
i = i+1; j = j-1
end
end
if j-l < r-i then do
if i < r then do
s = s+1; sl.s = i; sr.s = r
end
r = j
end
else do
if l < j then do
s = s+1; sl.s = l; sr.s = j
end
l = i
end
end
end
return
-- Module Run.rex - 1 Apr 2025
-- REXX Preprocessor and Executor
-- Usage: rexx|regina run <your program> <your arguments>
Main:
-- Main program flow
-- Trap conditions
signal on halt name Abend
signal on notready name Abend
signal on syntax name Abend
call Time('r')
-- Get various parameters
parse arg program parms
parse version version
-- Program to process, suffix 'rex' assumed if not specified
if Pos('.',program) = 0 then
program = program'.rex'
-- Get search path from Windows environment
-- Current folder and main folder '\Rex' assumed
p = Value(path,,environment)
n = 1; srce.n = ''; search = ''
do while p <> ''
parse var p a';'p
if Translate(Left(a,5)) = '\REX\' then do
n = n+1; srce.n = a'\'; search = search a
end
end
srce.0 = n; search = Strip(search)
-- Find program to process
do n = 1 to srce.0 until s <> ''
rex = srce.n||program
s = Stream(rex,'c','query exists')
end
if s = '' then
call Abend 'P'
n1 = 0; nr = 0; stat. = 0; sn = 0
-- Open output program to generate
runner = 'Runner.rex'
stampr = Stream(rex,'c','query timestamp')
stampr = Date('n',Left(stampr,10),'i') Right(stampr,8)
call Stream runner,'c','open write replace'
-- Write source info
call Lineout runner,'--' rex
nr = nr+1
-- Process source
le = Length(rex)
call Stream rex,'c','open read'
do while Lines(rex) = 1
n1 = n1+1; rexline = Linein(rex)
-- Something to include? Verb 'include' assumed
if Translate(Word(rexline,1)) = 'INCLUDE' then do
-- Find source to include. Suffix '.inc' assumed
do n = 1 to srce.0 until s <> ''
inc = srce.n||Word(rexline,2)'.inc'
s = Stream(inc,'c','query exists')
end
if s = "" then
call Abend 'I'
-- Saved on
stampi = Stream(inc,'c','query timestamp')
stampi = Date('n',Left(stampi,10),'i') Right(stampi,8)
-- Write header
call Lineout runner,'--' inc
nr = nr+1
-- Save change parameters
parm = 0
do wrd = 3 to Words(rexline)
parm = parm+1; parm.parm = Word(rexline,wrd)
end
parm.0 = parm; eerste = 1; n2 = 0
-- Process file to include
le = Max(Length(inc),le)
call Stream inc,'c','open read'
do while Lines(inc) = 1
n2 = n2+1
incline = Linein(inc)
-- Verb 'default' assumed for text replacement
if Translate(Word(incline,1)) = 'DEFAULT' then do
-- Add default change parameters
do wrd = 2 to Words(incline)
tmp = Word(incline,wrd)
parse var tmp change '=' .
if Pos(change,rexline) = 0 then do
parm = parm+1; parm.parm = Word(incline,wrd)
end
end
parm.0 = parm
end
else do
-- Write change parameters
if eerste = 1 then do
if parm.0 > 0 then do
changes = ''
do prm = 1 to parm.0
changes = changes parm.prm
end
call Lineout runner,'-- Replace 'changes
nr = nr+1
end
eerste = 0
end
-- Process changes
do prm = 1 to parm.0
parse var parm.prm change '=' in
incline = Changestr(change,incline,in)
end
-- Write line from include
if Pos('ooRexx',version) > 0 then do
call Lineout runner,incline
nr = nr+1
end
else do
if Left(incline,2) <> '::' then do
call Lineout runner,incline
nr = nr+1
end
end
end
end
-- Write some footers
call Stream inc,'c','close'
call Lineout runner,'--' inc
nr = nr+1
-- Save stats
sn = sn+1; stat.include.sn = inc; stat.stamp.sn = stampi; stat.Lines.sn = n2
end
else do
-- Write original source
if Pos('ooRexx',version) > 0 then do
call Lineout runner,rexline
nr = nr+1
end
else do
if Left(rexline,2) <> '::' then do
call Lineout runner,rexline
nr = nr+1
end
end
end
end
-- Source processing done
call Stream rex,'c','close'
-- Write stats
se = '--' Copies('-',le+27)
call Lineout runner,''
call Lineout runner,se
call Lineout runner,'--' Left('Summary',le) 'Lines' Left('Timestamp',20)
call Lineout runner,se
call Lineout runner,'--' Left(rex,le) Right(n1,5) Right(stampr,20)
nr = nr+5
do i = 1 to sn
call Lineout runner,'--' Left(stat.include.i,le) Right(stat.Lines.i,5) Right(stat.stamp.i,20)
nr = nr+1
end
nr = nr+3
call Lineout runner,se
call Lineout runner,'--' Left('Runner.rex',le) Right(nr,5) Right(Date() Time(),20)
call Lineout runner,se
call Stream runner,'c','close'
-- Run the generated program. Regina and ooRexx precoded.
select
when Pos('Regina',version) > 0 then
'regina' runner parms
when Pos('ooRexx',version) > 0 then
'rexx' runner parms
end
exit
Abend:
-- Condition handling specific for the preprocessor
arg type
say
signal off halt
signal off notready
signal off syntax
parse version version
say 'Rexx :' version
parse source . . source
say 'Source :' source
say 'Progr :' program
select
when type = 'P' then do
say 'Search :' search
say 'Error : Program not found'
end
when type = 'I' then do
say 'Incl :' Word(rexline,2)
say 'Search :' search
say 'Error : Include not found'
end
otherwise do
cc = Condition('c'); dd = Condition('d')
select
when cc = 'HALT' then
cc = 'Program interrupted'
when cc = 'NOTREADY' then
cc = 'Program or Include is in use'
when cc = 'NOVALUE' then
cc = 'Variable has no value'
when cc = 'SYNTAX' then
cc = 'Runtime error'
otherwise
cc = ''
end
if cc <> '' then
say 'Signal :' cc
if Datatype(rc) = 'NUM' then
if rc > 0 then
say 'Error :' rc Errortext(rc)
if dd <> '' then
if dd <> 'SIGINT' then
say 'Reason :' dd
end
end
exit
-- Module Calc.rex - 1 Apr 2025
-- REXX calculator
-- Regina|Rexx calc and play at the REXX prompt
IgnoreError16:
include Settings
do forever
call Charout, 'REXX '
pull command
if command = '' then
leave
interpret command
say
end
return
include Abend
include Complex
include Constants
include Functions
include Numbers
include Polynomial
include Rational
include Roots
include Sequences
-- Module Demo1.rex - 22 Feb 2025
-- Demo having only required modules included
Main:
-- Main program flow
include Settings
say version; say 'Complex arithmetic'; say
a = '1 2'; b = '3 4'; c = '5 6'; d = '7 8'
say 'VALUES'
say 'a =' Crect2form(a)
say 'b =' Crect2form(b)
say 'c =' Crect2form(c)
say 'd =' Crect2form(d)
say
say 'BASICS'
say 'a+b =' Crect2form(Cadd(a,b))
say 'a+b+c+d =' Crect2form(Cadd(a,b,c,d))
say 'a-b =' Crect2form(Csub(a,b))
say 'a-b-c-d =' Crect2form(Csub(a,b,c,d))
say 'a*b =' Crect2form(Cmul(a,b))
say 'a*b*c*d =' Crect2form(Cmul(a,b,c,d))
say 'a/b =' Crect2form(Cdiv(a,b))
say 'a/b/c/d =' Crect2form(Cdiv(a,b,c,d))
say '-a =' Crect2form(Cneg(a))
say '1/a =' Crect2form(Cinv(a))
say
say 'BONUS'
say 'Argument(a) =' Carg(a)
say 'Conjugate(a) =' Crect2form(Cconj(a))
say 'ImagPart(a) =' Cim(a)
say 'Modulus(a) =' Cmod(a)
say 'Polar(a) =' Cpol2form(Cpol(a))
say 'RealPart(a) =' Cre(a)
say
say 'FORMULA'
say 'a^2-2ab+3c-4ad^4+5 = ',
Crect2form(Cadd(Cpow(a,2),Cmul(-2,a,b),Cmul(3,c),Cmul(-4,a,Cpow(d,4)),5))
exit
include Complex
include Functions
include Constants
include Abend
-- Module Demo2.rex - 22 Feb 2025
-- Demo having the whole package included
Main:
-- Main program flow
include Settings
say version; say 'Rational arithmetic'; say
a = '1 2'; b = '-3 4'; c = '5 -6'; d = '-7 -8'; e = 3; f = 1.666666666
say 'VALUES'
say 'a =' Rlst2form(a)
say 'b =' Rlst2form(b)
say 'c =' Rlst2form(c)
say 'd =' Rlst2form(d)
say 'e =' e
say 'f =' f
say
say 'BASICS'
say 'a+b =' Rlst2form(Radd(a,b))
say 'a+b+c+d =' Rlst2form(Radd(a,b,c,d))
say 'a-b =' Rlst2form(Rsub(a,b))
say 'a-b-c-d =' Rlst2form(Rsub(a,b,c,d))
say 'a*b =' Rlst2form(Rmul(a,b))
say 'a*b*c*d =' Rlst2form(Rmul(a,b,c,d))
say 'a/b =' Rlst2form(Rdiv(a,b))
say 'a/b/c/d =' Rlst2form(Rdiv(a,b,c,d))
say '-a =' Rlst2form(Rneg(a))
say '1/a =' Rlst2form(Rinv(a))
say
say 'COMPARE'
say 'a<b =' Rlt(a,b)
say 'a<=b =' Rle(a,b)
say 'a=b =' Req(a,b)
say 'a>=b =' Rge(a,b)
say 'a>b =' Rgt(a,b)
say 'a<>b =' Rne(a,b)
say
say 'BONUS'
say 'Abs(c) =' Rlst2form(Rabs(c))
say 'Equal(a,b) =' Req(a,b)
say 'Float(b) =' Rfloat(b)
say 'Neg(d) =' Rlst2form(Rneg(d))
say 'Power(a,e) =' Rlst2form(Rpower(a,e))
say 'Rational(f) =' Rlst2form(Rrat(f))
say
say 'FORMULA'
say 'a^2-2ab+3c-4ad^4+5 = ',
Rlst2form(Radd(Rpower(a,2),Rmul(-2,a,b),Rmul(3,c),Rmul(-4,a,Rpower(d,4)),5))
exit
include Math
-- Module Model1.rex - 23 Feb 2025
-- Template for your own math playground
Main:
-- Main program flow
include Settings
-- Your code...
exit
include Math
-- Module Abend.inc - 11 Apr 2025
-- Condition/abend handler and debugging tool
-- Also used for argument error trapping
-- and display of the program call stack
Abend:
-- Condition and abend handler
-- Disable signals
Trace off
signal off halt
signal off notready
signal off novalue
signal off syntax
-- Condition info
parse version version!; parse source . . source!
s! = Strip(Sourceline(sigl)); c! = Condition('c'); d! = Condition('d')
-- Source info
if rc <> 16 then do
say
say 'Rexx :' version!
say 'Source :' source!
end
say 'Line :' sigl s!
-- Variables
call GetVariables
call ShowVariables
-- Condition message
select
when d! = 'ARGUMENT' then
c! = 'Argument outside domain, invalid or missing'
when d! = 'STACK' then
c! = 'Program call stack display forced by user'
when c! = 'HALT' then
c! = 'Program interrupted by user'
when c! = 'NOTREADY' then
c! = 'File does not exist or is in use by another app'
when c! = 'NOVALUE' then
c! = 'Variable has no value'
when c! = 'SYNTAX' then
c! = 'Runtime error'
otherwise
c! = ''
end
-- Additional info
if rc <> 16 then do
say 'Digits :' Digits()
if c! <> '' then
say 'Signal :' c!
if d! <> 'ARGUMENT' & d! <> 'STACK' then do
if Datatype(rc) = 'NUM' then
if rc > 0 then
say 'Error :' rc Errortext(rc)
if d! <> '' then
if d! <> 'SIGINT' then
say 'Reason :' d!
end
end
say
-- Force program call stack display
signal on novalue name IgnoreError16
say IgnoreError16
return
Debug:
-- Debugging info
-- Show relevant info
s! = Strip(Sourceline(sigl)); d! = Condition('d')
say
say 'Line :' sigl s!
call GetVariables
call ShowVariables
return
GetVariables:
-- Variable collector
-- Collect variables from source line
v!. = 0
s! = Upper(s!)
s! = Translate(s!," ","+-*|<>,;=/\%()'")
do i! = 1 to Words(s!)
w! = Word(s!,i!)
if Symbol(w!) = 'VAR' then do
call Variable
if Pos('.',w!) > 0 then do
t! = Translate(w!," ",".")
do j! = 2 to Words(t!)
w! = Word(t!,j!)
if Symbol(w!) = 'VAR' then
call Variable
end
end
end
end
return
Variable:
-- Variable saver
-- Save variable
do k! = 1 to v!.0
if w! = v!.k! then
leave
end
if k! > v!.0 then do
v!.k! = w!; v!.0 = k!
end
return
ShowVariables:
-- Variable displayer
-- Show values
do i! = 1 to v!.0
v! = v!.i!; a! = Value(v!)
if d! = 'ARGUMENT' then
say 'Arg :' v! '=' '"'a!'"'
else
say 'Var :' v! '=' '"'a!'"'
end
return
-- Module Complex.inc - 22 Mar 2025
-- Complex numbers arithmetic, functions and tools
-- Numbers are coded as string 'real imaginary'
-- Examples '1 2' = 1+2i, '0 1' = i, '1' = 1
Cabs:
-- Absolute value (modulus, magnitude)
-- |z|
procedure expose glob.
arg xx
-- Formula
return Cmod(xx)
Cadd:
-- Addition
-- z1+z2+z3+...
procedure
if Arg() = 0 then say argument
-- Add all arguments
do n = 1 to Arg()
xx = Arg(n) 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument n
if \ Datatype(xi,'n') then say argument n
if n = 1 then do
zr = xr; zi = xi
end
else do
zr = zr+xr; zi = zi+xi
end
end
return (zr) (zi)
Carccos:
-- Inverse cosine
-- -i*Ln(z+Sqrt(z^2-1))
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return Pi()
if Ceq(xx,-0.5) then
return Tau()/3
if Ceq(xx,0) then
return Halfpi()
if Ceq(xx,0.5) then
return Pi()/3
if Ceq(xx,1) then
return 0
-- Formula
if Ceq(xx,Cneg(I())) then
return (Halfpi()) (Ln(Sqrt(2)+1))
if Ceq(xx,I()) then
return (Halfpi()) (Ln(Sqrt(2)-1))
return Cmul(-1,I(),Cln(Cadd(xx,Csqrt(Csub(Csquare(xx),1)))))
Carccot:
-- Inverse cotangent
-- -i*Ln((z+i)/(z-i))/2
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return -0.25*Pi()
if Ceq(xx,0) then
return Halfpi()
if Ceq(xx,1) then
return 0.25*Pi()
-- Formula
return Cmul(-0.5,I(),Cln(Cdiv(Cadd(xx,I()),Csub(xx,I()))))
Carccsc:
-- Inverse cosecant
-- i*Ln(Sqrt(1-1/z^2)-i/z)
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-2) then
return -Pi()/6
if Ceq(xx,-1) then
return -Halfpi()
if Ceq(xx,1) then
return Halfpi()
if Ceq(xx,2) then
return Pi()/6
-- Formula
return Cmul(I(),Cln(Csub(Csqrt(Csub(1,Cinv(Csquare(xx)))),Cdiv(I(),xx))))
Carcosh:
-- Inverse cosine hyperbolic
-- Ln(z+Sqrt(z^2-1))
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return (0) (Pi())
if Ceq(xx,-0.5) then
return (0) (Tau()/3)
if Ceq(xx,0) then
return (0) (Halfpi())
if Ceq(xx,0.5) then
return (0) (Pi()/3)
if Ceq(xx,1) then
return 0
-- Formula
if Ceq(xx,Cneg(I())) then
return Cln(Cmul(-1,I(),1+Sqrt(2)))
if Ceq(xx,I()) then
return Cln(Cmul(I(),1+Sqrt(2)))
return Cln(Cadd(xx,Csqrt(Csub(Csquare(xx),1))))
Carcoth:
-- Inverse cotangent hyperbolic
-- Ln((1+z)/(1-z))/2
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,0) then
return (0) (Halfpi())
if Ceq(xx,Cneg(I())) then
return (0) (0.25*Pi())
if Ceq(xx,I()) then
return (0) (-0.25*Pi())
-- Formula
return Cmul(0.5,Cln(Cdiv(Cadd(xx,1),Csub(xx,1))))
Carcsch:
-- Inverse cosecant hyperbolic
-- Ln(1/z+Sqrt(1+1/z^2))
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return Ln(Sqrt(2)-1)
if Ceq(xx,1) then
return Ln(Sqrt(2)+1)
if Ceq(xx,Cmul(-2,I())) then
return (0) (Pi()/6)
if Ceq(xx,Cneg(I())) then
return (0) (Halfpi())
if Ceq(xx,I()) then
return (0) (-Halfpi())
if Ceq(xx,Cmul(2,I())) then
return (0) (-Pi()/6)
-- Formula
a = Cinv(xx); b = Csquare(a)
return Cln(Cadd(a,Csqrt(Cadd(b,1))))
Carcsec:
-- Inverse secant
-- -i*Ln(i*Sqrt(1-1/z^2)+1/z)
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-2) then
return Tau()/3
if Ceq(xx,-1) then
return 0
if Ceq(xx,1) then
return 0
if Ceq(xx,2) then
return Pi()/3
-- Formula
if Ceq(xx,Cneg(I())) then
return (Halfpi()) (Ln(Sqrt(2)-1))
if Ceq(xx,I()) then
return (Halfpi()) (Ln(Sqrt(2)+1))
return Cmul(-1,I(),Cln(Cadd(Cmul(I(),Csqrt(Csub(1,Cinv(Csquare(xx))))),Cinv(xx))))
Carcsin:
-- Inverse sine
-- i*Ln(Sqrt(1-z^2)-i*z)
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return -Halfpi()
if Ceq(xx,-0.5) then
return -Pi()/6
if Ceq(xx,0) then
return 0
if Ceq(xx,0.5) then
return Pi()/6
if Ceq(xx,1) then
return Halfpi()
-- Formula
if Ceq(xx,Cneg(I())) then
return (0) (Ln(Sqrt(2)-1))
if Ceq(xx,I()) then
return (0) (Ln(Sqrt(2)+1))
return Cmul(I(),Cln(Csub(Csqrt(Csub(1,Csquare(xx))),Cmul(I(),xx))))
Carctan:
-- Inverse tangent
-- Ln((1+i*z)/(1-i*z))/(2i)
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return -0.25*Pi()
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 0.25*Pi()
-- Formula
a = Cmul(I(),xx)
return Cdiv(Cln(Cdiv(Cadd(1,a),Csub(1,a))),I(),2)
Carg:
-- Argument (phase)
-- Arctan2(Re(z),Im(z))
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Formula
numeric digits Digits()+1
return Arctan2(xr,xi)
Carsech:
-- Inverse secant hyperbolic
-- Ln(1/z+Sqrt(1-1/z^2))
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-2) then
return (0) (Tau()/3)
if Ceq(xx,-1) then
return (0) (Pi())
if Ceq(xx,1) then
return 0
if Ceq(xx,2) then
return (0) (Pi()/3)
-- Formula
if Ceq(xx,Cneg(I())) then
return (Ln(Sqrt(2)+1)) (Cmul(0.5,Pi()))
if Ceq(xx,I()) then
return (Ln(Sqrt(2)+1)) (Cmul(-0.5,Pi()))
a = Cinv(xx); b = Csquare(a)
return Cln(Cadd(a,Csqrt(Csub(b,1))))
Carsinh:
-- Inverse sine hyperbolic
-- Ln(z+Sqrt(z^2+1))
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,-1) then
return Ln(Sqrt(2)-1)
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return Ln(Sqrt(2)+1)
if Ceq(xx,Cneg(I())) then
return (0) (-Halfpi())
if Ceq(xx,Cmul(-0.5,I())) then
return (0) (-Pi()/6)
if Ceq(xx,Cmul(0.5,I())) then
return (0) (Pi()/6)
if Ceq(xx,I()) then
return (0) (Halfpi())
-- Formula
return Cln(Cadd(xx,Csqrt(Cadd(Csquare(xx),1))))
Cartanh:
-- Inverse tangent hyperbolic
-- Ln((1+z)/(1-z))/2
procedure expose glob.
arg xx
-- Fast
numeric digits Digits()+1
if Ceq(xx,0) then
return 0
if Ceq(xx,Cneg(I())) then
return (0) (-0.25*Pi())
if Ceq(xx,I()) then
return (0) (0.25*Pi())
-- Formula
return Cmul(0.5,Cln(Cdiv(Cadd(1,xx),Csub(1,xx))))
Ccbrt:
-- Cubic root
-- z^(1/3)
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 1
-- Formula
numeric digits Digits()+2
zz = Cmul(Cbrt(Cmod(xx)),Cexp(Cmul(I(),Carg(xx),1/3)))
return zz
Cconj:
-- Conjugate
-- Re(z) -Im(z)
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Formula
return (xr) (-xi)
Ccos:
-- Cosine
-- Cos(Re(z))Cosh(Im(z)) -Sin(Re(z))Sinh(Im(z))
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 1
-- Formula
numeric digits Digits()+1
return (Cos(xr)*Cosh(xi)) (-Sin(xr)*Sinh(xi))
Ccosh:
-- Cosine hyperbolic
-- Cosh(Re(z))Cos(Im(z)) Sinh(Re(z))Sin(Im(z))
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 1
-- Formula
numeric digits Digits()+1
return (Cosh(xr)*Cos(xi)) (Sinh(xr)*Sin(xi))
Ccot:
-- Cotangent
-- Cos(z)/Sin(z)
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return Cdiv(Ccos(xx),Csin(xx))
Ccoth:
-- Cotangent hyperbolic
-- Cosh(z)/Sinh(z)
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return Cdiv(Ccosh(xx),Csinh(xx))
Ccsc:
-- Cosecant
-- 1/Sin(z)
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return Cinv(Csin(xx))
Ccsch:
-- Cosecant hyperbolic
-- 1/Sinh(z)
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return Cinv(Csinh(xx))
Ccube:
-- Cube
-- z^3
procedure
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 1
-- Formula
numeric digits Digits()+1
r2 = xr*xr; r3 = r2*xr; i2 = xi*xi; i3 = i2*xi
return (r3-3*xr*i2) (3*r2*xi-i3)
Cdiv:
-- Division
-- z1/z2/z3/...
procedure
if Arg() = 0 then say argument
-- Divide first argument by all other ones
numeric digits Digits()+1
do n = 1 to Arg()
xx = Arg(n) 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument n
if \ Datatype(xi,'n') then say argument n
if n = 1 then do
zr = xr; zi = xi
end
else do
if Ceq(xx,0) then say argument n
d = xr*xr+xi*xi
r = (zr*xr+zi*xi)/d; i = (zi*xr-zr*xi)/d
zr = r; zi = i
end
end
return (zr) (zi)
Ceq:
-- Equal?
procedure
arg xx,yy; xx = xx 0; yy = yy 0; parse var xx xr xi .; parse var yy yr yi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
if \ Datatype(yr,'n') then say argument
if \ Datatype(yi,'n') then say argument
-- z1 = z2?
return (xr=yr & xi=yi)
Ceval:
-- Evaluation polynomial with real coefficients
procedure
arg xx,yy; yy = yy 0; parse var yy yr yi .
if xx = '' then say argument
-- Horner's scheme
zz = 0 0
numeric digits Digits()+2
do n = 1 to Words(xx)
w = Word(xx,n)
zz = Cadd(Cmul(zz,yy),w)
end
return zz
Cexp:
-- Exponential
-- e^Re(z)*Cos(Im(z)) e^Re(z)*Sin(Im(z))
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 1
if Ceq(xx,1) then
return E()
-- Formula
numeric digits Digits()+1
a = Exp(xr)
return (a*Cos(xi)) (a*Sin(xi))
Cim:
-- Imaginary part
-- Im(z)
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Formula
return (xi)
Cinv:
-- Inversion
-- 1/z
procedure
arg xx
-- Formula
numeric digits Digits()+1
return Cdiv(1,xx)
Cln:
-- Natural logarithm
-- Ln|z| Arg(z)
procedure expose glob.
arg xx
if Ceq(xx,0) then say argument
-- Formula
numeric digits Digits()+1
return (Ln(Cmod(xx))) (Carg(xx))
Cmod:
-- Modulus (absolute value, magnitude)
-- |z|
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Formula
numeric digits Digits()+1
return Sqrt(xr*xr+xi*xi)
Cmul:
-- Multiplication
-- z1*z2*z3*...
procedure
if Arg() = 0 then say argument
-- Multiply all arguments
numeric digits Digits()+1
do n = 1 to Arg()
xx = Arg(n) 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument n
if \ Datatype(xi,'n') then say argument n
if n = 1 then do
zr = xr; zi = xi
end
else do
r = zr*xr-zi*xi; i = zr*xi+zi*xr
zr = r; zi = i
end
end
return (zr) (zi)
Cne:
-- Not equal
-- z1 <> z2?
procedure
arg xx,yy
-- Formula
return \ Ceq(xx,yy)
Cneg:
-- Negation
-- -Re(z) -Im(z)
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Formula
return (-xr) (-xi)
Cnorm:
-- Normalization
-- Re(z)/1 Im(z)/1
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Formula
return (xr/1) (xi/1)
Complex:
-- Are all arguments in C?
procedure
do n = 1 to Arg()
xx = Arg(n) 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument n
if \ Datatype(xi,'n') then say argument n
end
return 1
Cpol:
-- Rectangular -> Polar
-- Re(z) Im(z) -> |z| Arg(z)
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return (Cmod(xx)) (Carg(xx))
Cpol2form:
-- Polar -> Formula
procedure
arg xx; xx = xx 0; parse var xx xm xa .
if \ Datatype(xm,'n') then say argument
if \ Datatype(xa,'n') then say argument
-- Compose formula
if xm = 0 then
return 0
za = xm'r'
select
when xa < 0 then
zb = '-'Abs(xa)'th'
when xa = 0 then
zb = ''
when xa > 0 then
zb = ','xa'th'
end
return za||zb
Cpow:
-- Integer power
-- z^n, n in Z
procedure
arg xx,yy
if \ Datatype(yy,'w') then say argument
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 0
if Ceq(yy,0) then
return 1
if Ceq(yy,1) then
return xx
-- Calculate power by repeated divide or multiply
numeric digits Digits()+2
if yy < 0 then do
zz = Cinv(xx)
do Abs(yy)-1
zz = Cdiv(zz,xx)
end
end
else do
zz = xx
do yy-1
zz = Cmul(zz,xx)
end
end
return zz
Cpower:
-- Complex power
-- e^(w*Ln(z)), w in R or C
procedure expose glob.
arg xx,yy
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 0
if Ceq(yy,0) then
return 1
if Ceq(yy,1) then
return xx
-- Formula
numeric digits Digits()+1
return Cexp(Cmul(yy,Cln(xx)))
Cqtrt:
-- Quartic root
-- z^(1/4)
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 1
-- Formula
numeric digits Digits()+2
zz = Cmul(Cbrt(Cmod(xx)),Cexp(Cmul(I(),Carg(xx),0.25)))
return zz
Cquartic:
-- Quartic
-- z^4
procedure
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 1
-- Formula
numeric digits Digits()+1
-- Formula
numeric digits Digits()+1
r2 = xr*xr; r3 = r2*xr; r4 = r2*r2; i2 = xi*xi; i3 = i2*xi; i4 = i2*i2
return (r4-6*r2*i2+i4) (4*r3*xi-4*xr*i3)
Cre:
-- Real part
-- Re(z)
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Formula
return xr
Crect:
-- Polar -> Rectangular
-- |z| Arg(z) -> Re(z) Im(z)
procedure expose glob.
arg xx; xx = xx 0; parse var xx xm xa .
if \ Datatype(xm,'n') then say argument
-- Formula
numeric digits Digits()+1
return (xm*Cos(xa)) (xm*Sin(xa))
Crect2form:
-- Rectangular -> Formula
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- Compose formula
if xr = 0 then
za = ''
else
za = xr
select
when xi < 0 then
if xi = -1 then
zb = '-i'
else
zb = xi'i'
when xi = 0 then
if xr = 0 then
zb = 0
else
zb = ''
when xi > 0 then
if xi = 1 then
zb = '+i'
else
zb = '+'xi'i'
end
return za||zb
Cround:
-- Round to specified decimals
procedure
arg xx,yy; xx = xx 0; parse var xx xr xi .
if yy = '' then
yy = 0
-- Formula
return (Round(xr,yy)) (Round(xi,yy))
Csci:
-- Convert to scientific format
procedure
arg xx; xx = xx 0; parse var xx xr xi .
-- Formula
return Sci(xr) Sci(xi)
Csec:
-- Secant
-- 1/Cos(z)
procedure expose glob.
arg xx
-- Fast
if Ceq(xx,0) then
return 1
-- Formula
numeric digits Digits()+1
return Cinv(Ccos(xx))
Csech:
-- Secant hyperbolic
-- 1/Cosh(z)
procedure expose glob.
arg xx
-- Fast
if Ceq(xx,0) then
return 1
-- Formula
numeric digits Digits()+1
return Cinv(Ccosh(xx))
Csin:
-- Sine
-- Sin(Re(z))Cosh(Im(z)) Cos(Re(z))Sinh(Im(z))
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
-- Formula
numeric digits Digits()+1
return (Sin(xr)*Cosh(xi)) (Cos(xr)*Sinh(xi))
Csinh:
-- Sine hyperbolic
-- Sinh(Re(z))Cos(Im(z)) Cosh(Re(z))Sin(Im(z))
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
-- Formula
numeric digits Digits()+1
return (Sinh(xr)*Cos(xi)) (Cosh(xr)*Sin(xi))
Csqrt:
-- Square root
-- z^(1/2)
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 1
-- Formula
numeric digits Digits()+2
m = Cmod(xx)
a = Round(m+xr,Digits()-2); b = Round(m-xr,Digits()-2)
if xi < 0 then
s = -1
else
s = 1
r = Sqrt(0.5*a); i = Sqrt(0.5*b)*s
return (r) (i)
Csquare:
-- Square
-- z^2
procedure
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
if Ceq(xx,1) then
return 1
-- Formula
xr2 = xr*xr; xi2 = xi*xi
return (xr2-xi2) (2*xr*xi)
Cstd:
-- Convert to standard format
procedure
arg xx; xx = xx 0; parse var xx xr xi .
-- Formula
return Std(xr) Std(xi)
Csub:
-- Subtraction
-- z1-z2-z3-...
procedure
if Arg() = 0 then say argument
-- Subtract second and following arguments from the first one
numeric digits Digits()+1
do n = 1 to Arg()
xx = Arg(n) 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument n
if \ Datatype(xi,'n') then say argument n
if n = 1 then do
zr = xr; zi = xi
end
else do
zr = zr-xr; zi = zi-xi
end
end
return (zr) (zi)
Ctan:
-- Tangent
-- Sin(z)/Cos(z)
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
-- Formula
numeric digits Digits()+1
return Cdiv(Csin(xx),Ccos(xx))
Ctanh:
-- Tangent hyperbolic
-- Sinh(z)/Cosh(z)
procedure expose glob.
arg xx; xx = xx 0; parse var xx xr xi .
-- Fast
if Ceq(xx,0) then
return 0
-- Formula
numeric digits Digits()+1
return Cdiv(Csinh(xx),Ccosh(xx))
Czero:
-- Zero?
procedure
arg xx; xx = xx 0; parse var xx xr xi .
if \ Datatype(xr,'n') then say argument
if \ Datatype(xi,'n') then say argument
-- z1 = 0?
return (xr=0 & xi=0)
-- Module Constants.inc - 22 Mar 2025
-- Mathematical constants like Pi, E, Phi, i and many more
Apery:
-- Apery constant
procedure
-- Fast
if Digits() < 111 then
return 1.2020569031595942853997381615114499907649862923404988817922715553418382057863130901864558736093352581461991578
-- Wedeniwski recurring
numeric digits Digits()+2
-- Term k = 0
f10 = 1; f20 = 1; f21 = 1; f32 = 2; f43 = 6; p = 12463; s = 1
zz = 0; v = 0
numeric fuzz 2
do k = 1
-- Add term
zz = zz+s*((f10*f20*f21)**3*p)/(f32*f43**3)
if zz = v then
leave
v = zz
-- Prepare next term
k2 = k+k; k3 = k2+k; k4 = k3+k
-- Factorials
f10 = f10*k; f20 = f20*k2*(k2-1); f21 = f21*k2*(k2+1); f32 = f32*k3*(k3+1)*(k3+2); f43 = f43*k4*(k4+1)*(k4+2)*(k4+3)
-- Polynomial
p = ((((126392*k+412708)*k+531578)*k+336367)*k+104000)*k+12463
-- Sign
s = -s
end
return zz/24
Artin:
-- Artin constant
procedure
-- Fast
return 0.373955813619202288054728054346416415111629248606150042094742802417350182040028082344304317087250568981603
Backhouse:
-- Backhouse constant
procedure
-- Fast
return 1.45607494858268967139959535111654355765317837484713154027070243741400150626538989559964531940186030910992
Bernstein:
-- Bernstein constant
procedure
-- Fast
return 0.28016949902386913303643649123067200004248213981236
Bloch:
-- Bloch constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.471861653452681784874468793611316149077012621739443243631460266888006805241859656902448655624594840993
-- Formula
return Sqrt((Sqrt(3)-1)*0.5)*Gamma(1/3)*Gamma(11/12)/Gamma(0.25)
Bronze:
-- Bronze ratio constant
-- Positive root of x^2-3x-1
procedure
-- Fast
if Digits() < 111 then
return 3.3027756377319946465596106337352479731256482869226231063552265281135834741465052226023095410092453588367570910
-- Formula
return 0.5*(3+Sqrt(13))
Brun:
-- Brun constant
procedure
-- Fast
return 1.902160583104
Cahen:
-- Cahen constant
procedure
-- Fast
if Digits() < 111 then
return 0.64341054628833802618225430775756476328658786026823950598703092030749277646183261084844089555046343195405372901
-- Sylvester
numeric digits Digits()+2
zz = 0; v = 0; s = 2
do n = 0
zz = zz+(-1)**n/(s-1)
if zz = v
then leave
v = zz; s = s*s-s+1
end
return zz
Capacity:
-- Logarithmic capacity of the unit disk constant
procedure
-- Fast
return 0.5901702995080481130226689702792442936168583174407236497579321997021520903603578974892293080979039771
Catalan:
-- Catalan constant
procedure
-- Fast
if Digits() < 111 then
return 0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442262
-- Broadhurst
numeric digits Digits()+2
numeric fuzz 2
s1 = 0; v = 0
do n = 0 until s1 = v
v = s1
a = 1/2**(4*n); b = 8*n
c1 = (b+1)*(b+1); c2 = (b+2)*(b+2); c3 = (b+3)*(b+3); c5 = (b+5)*(b+5); c6 = (b+6)*(b+6); c7 = (b+7)*(b+7)
s1 = s1+a*(1/(2*c1)-1/(2*c2)+1/(4*c3)-1/(8*c5)+1/(8*c6)-1/(16*c7))
end
s2 = 0; v = 0
do n = 0 until s2 = v
v = s2
a = 1/2**(12*n); b = 8*n
c1 = (b+1)*(b+1); c2 = (b+2)*(b+2); c3 = (b+3)*(b+3); c5 = (b+5)*(b+5); c6 = (b+6)*(b+6); c7 = (b+7)*(b+7)
s2 = s2+a*(1/(8*c1)+1/(16*c2)+1/(64*c3)-1/(512*c5)-1/(1024*c6)-1/(4096*c7))
end
return 3*s1-2*s2
Chaitin:
-- Chaitin constant
procedure
-- Fast
return 0.0078749969978123844
Champernowne:
-- Champernowne constant
procedure
-- Fast
if Digits() < 111 then
return 0.123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616
-- Concatenate integers
p = Digits()+4
zz = '0.'
do n = 1 while Length(zz) <= p
zz = zz||n
end
return zz
Conway:
-- Conway constant
-- Real root of x^71-x^69-2x^68...-6x^2+3x-6
procedure
-- Fast
if Digits() < 111 then
return 1.3035772690342963912570991121525518907307025046594048757548613906285508878524615571268157668644252255534713930471
-- Polynomial and derivatives
p =,
'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'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.44*d)
end
d.n = 4
-- Halley
zz = 1.3
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz - (2*f*f1) / (2*f1*f1-f*f2)
end
return zz
Copeland:
-- Copeland-Erdos constant
procedure
-- Fast
return 0.23571113171923293137414347535961677173798389971011031071091131271311371391491511571631671731791811911
Devicci:
-- DeVicci constant
-- Positive root of 4x^8-28x^6-7x^4+16x^2+16
procedure
-- Fast
if Digits() < 111 then
return 1.0074347568842793760982535952310991419256901141136697702349637985711523132802867779625205514746359239420763739328
-- Polynomial and derivatives
p = '4 0 -28 0 -7 0 16 0 16'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.38*d)
end
d.n = 4
-- Halley
zz = 1.0
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Dottie:
-- Dottie constant
-- Root of Cos(x)-x
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.73908513321516064165531208767387340401341175890075746496568063577328465488354759459937610693176653184980124664
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.36*d)
end
d.n = 4
-- Halley
zz = 0.74
do k = n by -1 to 1
numeric digits d.k
c = Cos(zz); s = Sin(zz); f = c-zz; f1 = -s-1; f2 = -c
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Dubois:
-- Second Dubois-Reymond constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.1945280494653251136152137302875039065901577852759236620435639112612868980395288816921562425395608973868765806327394
-- Formula
return (E()**2-7)*0.5
E:
-- Euler constant
procedure expose glob.
-- Fast
p = Digits()
if p < 111 then
return 2.7182818284590452353602874713526624977572470936999595749669676277240766303535475945713821785251664274274663919
if glob.e.p <> '' then
return glob.e.p
-- Taylor
numeric digits Digits()+2
numeric fuzz 2
zz = 2; t = 1
do n = 2 until zz = v
v = zz
t = t/n; zz = zz+t
end
glob.e.p = zz
return zz
Embree:
-- Embree-Trefethen constant
procedure
-- Fast
return 0.70258
Erdos1:
-- Erdos-Borwein constant
procedure
-- Fast
if Digits() < 111 then
return 1.6066951524152917637833015231909245804805796715057564357780795536914184207434866905657118016701555758970454290631544131
-- Series
numeric digits Digits()+2
zz = 0
numeric fuzz 2
do n = 1 until zz = v
v = zz
a = 2**n; zz = zz+(1/a**n)*(a+1)/(a-1)
end
return zz
Erdos2:
-- Erdos-Tenenbaum-Ford constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.086071332055934206887573098776922677760591109530333173492020236665422635814622879799380534602528768071633655969112
-- Formula
return 1-(1+Ln(Ln(2)))/Ln(2)
Euler:
-- Euler-Mascheroni constant
procedure expose glob.
-- Fast
p = Digits()
if p < 111 then
return 0.57721566490153286060651209008240243104215933593992359880576723488486772677766467093694706329174674951463144725
if glob.euler.p <> '' then
return glob.euler.p
-- Brent McMillan
numeric digits Digits()+2
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
glob.euler.p = s/v
return glob.euler.p
Feigenbaum1:
-- Feigenbaum first constant
procedure
-- Fast
return 4.66920160910299067185320382046620161725818557747576863274565134300413433021131473713868974402394801381716
Feigenbaum2:
-- Feigenbaum second constant
procedure
-- Fast
return 2.50290787509589282228390287321821578638127137672714997733619205677923546317959020670329964974643383412959
Feller:
-- Feller-Tornier constant
procedure
-- Fast
return 0.66131704946962233528976584627411853328547528983291635498090562622662503174312230494226174078428187
Foias:
-- Foias constant
procedure
-- Fast
return 1.18745235112650105459548015839651935121569268158586035301010412619878041872352540738702465760608657943378
Fraction:
-- First continued fraction constant
procedure
-- Fast
return 0.697774657964007982006790592551752599486658262998021232368630082816530852764641112996965654182676568723982
Fransen:
-- Fransen-Robertson constant
procedure
-- Fast
return 2.80777024202851936522150118655777293230808592093019829122005480959710088912190166551018530816819663814187
Gauss:
-- Gauss constant
procedure
-- Fast
if Digits() < 111 then
return 0.83462684167407318628142973279904680899399301349034700244982737010368199270952641186969116035127532412906785035
-- Formula
return 1/Agmean(1,Sqrt(2))
Gelfond1:
-- Gelfond constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 23.14069263277926900572908636794854738026610624260021199344504640952434235069045278351697199706754921967595270
-- Formula
return Power(E(),Pi())
Gelfond2:
-- Gelfond-Schneider constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 2.665144142690225188650297249873139848274211313714659492835979593364920446178705954867609180005196416941989364
-- Formula
return Power(2,Sqrt(2))
Gieseking:
-- Gieseking constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 1.014941606409653625021202554274520285941689307530299792017489106776597476258244022136470354228256694945861828
-- Formula
return (Trigamma(1/3)-Trigamma(2/3))/(4*Sqrt(3))
Glaisher:
-- Glaisher-Kinkelin constant
procedure
-- Fast
return 1.28242712910062263687534256886979172776768892732500119206374002174040630885882646112973649195820237
Golden:
-- Golden ratio constant
-- Root of x^2-x-1
procedure
-- Fast
if Digits() < 111 then
return 1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137484754088
-- Formula
return (Sqrt(5)+1)*0.5
Goldenangle:
-- Golden angle constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 2.399963229728653322231555506633613853124999011058115042935112750731307338239438790779962060660583963731372963
-- Formula
return Pi()*(3-Sqrt(5))
Golomb:
-- Golomb constant
procedure
-- Fast
return 0.624329988543550870992936383100837244179642620180529286973551902495638088855113254462460276195539868869
Gompertz:
-- Gompertz constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.5963473623231940743410784993692793760741778601525487815734849104823272191148744174704304970936127603442370347484286
-- Taylor
numeric digits Digits()+2
zz = 0; s = 1; f = 1
numeric fuzz 2
do n = 1 until zz = v
v = zz
s = -s; f = f*n; zz = zz+s/(n*f)
end
numeric fuzz
return -E()*(Euler()+zz)
Hafner:
-- Hafner-Sarnak-McCurley constant
procedure
-- Fast
return 0.353236371854995984543516550432682011280164778566690446416085942814238325002669003483672078334335498967
Halfpi:
-- Pi/2
procedure expose glob.
-- Fast
if Digits() < 111 then
return 1.570796326794896619231321691639751442098584699687552910487472296153908203143104499314017412671058533991074043
-- Formula
return 0.5*Pi()
Heath:
-- Heath-Brown-Moroz constant
procedure
-- Fast
return 0.00131764115485317810981735232251358595125073432329525167879254742178602344409610895090834
Hermite:
-- Second Hermite constant
procedure
-- Fast
if Digits() < 111 then
return 1.154700538379251529018297561003914911295203502540253752037204652967955344605866691387430791171499050450417428
-- Formula
return 2/Sqrt(3)
I:
-- Imaginary unit
procedure
-- Fast
return (0 1)
Kepler:
-- Kepler-Bouwkamp constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.11494204485329620070104015746959874283079533720086351684402339651896601282535305117794077248498583699376347675
d = Digits()
numeric digits d*3
-- Optimized asymptotic
zz = 0; v = 0
do k = 1 to 1000 until zz = v
v = zz
numeric fuzz
a = 2*k; b = Zeta(a); c = 2**a
numeric fuzz Digits()-d-2
zz = zz+((c-1)/a)*b*(b-1-1/c)
end
numeric fuzz
return Exp(-2*zz)
Khinchin:
-- Khinchin constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 2.6854520010653064453097148354817956938203822939944629530511523455572188595371520028011411749318476979951534659
-- Asymptotic
numeric digits Digits()+2
y1 = 0
do n = 1 until y1 = v
v = y1
y2 = 0
do k = 1 to 2*n-1
y2 = y2+(-1)**(k+1)/k
end
numeric fuzz
y1 = y1+((Zeta(2*n)-1)/n)*y2
numeric fuzz 2
end
numeric fuzz
return Exp(y1/Ln(2))
Komornik:
-- Komornik-Loreti constant
procedure
-- Fast
return 1.78723165018296593301327489033700838533793140296181099778478147050555749175060568699131001863407533302
Landau1:
-- Landau-Ramanujan constant
procedure
-- Fast
return 0.7642236535892206629906987312500923281167905413934095147216866737496146416587328588384015050131312337
Landau2:
-- Landau constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.543258965342976706952728295300613231138863293758356988955732569100434687015307837312852127102259710857533437902199
-- Formula
return Power(2,5/3)*Pi()**2/(3*Gamma(1/3)**3)
Laplace:
-- Laplace limit constant
procedure
-- Fast
return 0.66274341934918158097474209710925290705623354911502241752039253499097185308651127724965480259895818168
Lebesgue:
-- Lebesgue constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.98943127383114695174164880901886671492420140609111104828412243264437253158457865463463298314018955793304417328
-- Series
numeric digits Digits()+2
zz = 0
do n = 0 until zz = v
v = zz
numeric fuzz
zz = zz+(Lambda(2*n+2)-1)/(2*n+1)
numeric fuzz 2
end
numeric fuzz
a = Pi()**2
return 8*zz/a+4*(2*Ln(2)+Euler())/a
Lemniscate:
-- Lemniscate constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 2.62205755429211981046483958989111941368275495143162316281682170380079058707041425023029553296142909344613575
-- Formula
return Pi()/Agmean(1,Sqrt(2))
Levy1:
-- Levy-Khinchin constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 1.186569110415625452821722975947237120568356536472054335954254298652809632056254443300348301108486875946639263747175
-- Formula
return Pi()**2/(12*Ln(2))
Levy2:
-- Second Levy constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 3.275822918721811159787681882453843863608475525982374149405198924190723215644960355181277540479174529492698526243402
-- Formula
return Exp(Levy1())
Lieb:
-- Lieb constant
procedure
-- Fast
if Digits() < 111 then
return 1.539600717839002038691063414671886548393604670053671669382939537290607126141155588516574388228665400600556570147028
-- Formula
return 8/(3*Sqrt(3))
Ln10:
-- Natural log of 10 constant
procedure expose glob.
-- Fast
p = Digits()
if p < 111 then
return 2.3025850929940456840179914546843642076011014886287729760333279009675726096773524802359972050895982983419677840
if glob.ln10.p <> '' then
return glob.ln10.p
-- Weisstein
numeric digits Digits()+2
zz = 0
do k = 0 until zz = v
v = zz
k8 = 8*k; zz = zz+(729/(k8+1)+162/(k8+2)+81/(k8+3)+9/(k8+5)+2/(k8+6)+1/(k8+7))/6561**k
end
zz = 2*zz/729; glob.ln10.p = zz
return zz
Ln2:
-- Natural log of 2 constant
procedure expose glob.
-- Fast
p = Digits()
if p < 111 then
return 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754200148102
if glob.Ln2.p <> '' then
return glob.Ln2.p
-- Weisstein
numeric digits Digits()+2
zz = 0
do k = 0 until zz = v
v = zz
k18 = 18*k; zz = zz+(2187/(k18+1)-1458/(k18+3)-729/(k18+4)-81/(k18+7)+54/(k18+9)+27/(k18+10)+3/(k18+13)-2/(k18+15)-1/(k18+16))/(-19683)**k
end
zz = zz/2187; glob.Ln2.p = zz
return zz
Ln4:
-- Natural log of 4 constant
procedure expose glob.
-- Formula
return 2*Ln2()
Loch:
-- Loch constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.970270114392033925740256019210010833781284704785161286610350529931254199891737048036212674908029026469241589522939
-- Formula
return 6*Ln(2)*Ln(10)/Pi()**2
Magicangle:
-- Magic angle constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.95531661812450927816385710251575775424341469501000549095969812932191204590397645538739160258562807342160190353
-- Formula
return Arccos(1/Sqrt(3))
Maxnumber:
-- Maximum number allowed
procedure
-- Value
return 9E999999999
Meissel:
-- Meissel-Mertens constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.2614972128476427837554268386086958590515666482611992061920642139249245108973682097141426314342466510516177288764860
-- Series
numeric digits Digits()+2
zz = 0; v = 0
do n = 2 to 1000
t = Moebius(n)*Ln(Zeta(n))/n
if t <> 0 then do
zz = zz+t
if zz = v then
leave
v = zz
end
end
return Euler()+zz
Mills:
-- Mills constant
procedure
-- Fast
return 1.30637788386308069046861449260260571291678458515671364436805375996643405376682659882150140370119739570729
Mrb:
-- MRB constant
procedure
-- Fast
return 0.187859642462067120248517934054273230055903094900138786172004684089477231564660213703296654433107496903842
Nielsen:
-- Nielsen-Ramanujan constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.822467033424113218236207583323012594609474950603399218867779114685003735201600436916814450309879352652002159
-- Formula
return Pi()**2/12
Niven:
-- Niven constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 1.705211140105367764288551453434508160762027651653469099994284906547313192168122491934244132100871001797653827048
-- Series
numeric digits Digits()+2
zz = 0
do k = 2 to 1000 until zz = v
v = zz
numeric fuzz
zz = zz+1-1/Zeta(k)
numeric fuzz 2
end
numeric fuzz
return zz+1
Omega:
-- Omega constant
-- Root of x*e^x-1
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.56714329040978387299996866221035554975381578718651250813513107922304579308668456669321944696175229455763802497
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.40*d)
end
d.n = 4
-- Halley
zz = 0.57
do k = n by -1 to 1
numeric digits d.k
a = Exp(zz); b = a*zz-1
zz = zz-b/(a*(zz+1)-((zz+2)*b/(2*zz+2)))
end
return zz
Paperfolding:
-- Paper folding constant
procedure
-- Fast
if Digits() < 111 then
return 0.8507361882018672603677977605320666044113994930827106477281682610318301584593194453485459826421939239960918042677191
-- Definition
numeric digits Digits()+2
numeric fuzz 2
zz = 0
do n = 0 until zz = v
v = zz
zz = zz+8**(2**n)/(2**(2**(n+2))-1)
end
return zz
Parabolic:
-- Universal parabolic constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 2.295587149392638074034298049189490387597832203638583483929975346644109662684133126684094426237897615591753662
-- Formula
return Ln(Sqrt(2)+1)+Sqrt(2)
Pi:
-- Pi constant
procedure expose glob.
-- Fast
p = Digits()
if p < 111 then
return 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148087
if glob.pi.p <> '' then
return glob.pi.p
-- Chudnovsky
numeric digits Digits()+2
if p < 201 then do
zz = 0
do n = 0 until zz = v
numeric fuzz
v = zz; zz = zz + Fact(6*n)*(13591409+545140134*n)/(Fact(3*n)*Fact(n)**3*-640320**(3*n))
numeric fuzz 2
end
numeric fuzz
zz = 4270934400/(Sqrt(10005)*zz)
end
else do
-- Agmean
zz = 0.25; z = 1; g = Sqrt(0.5); n = 1
do until z = v
v = z
numeric fuzz
x = (z+g)*0.5; g = Sqrt(z*g)
numeric fuzz 2
zz = zz-n*(x-z)**2; n = n+n; z = x
end
numeric fuzz
zz = z*z/zz
end
glob.pi.p = zz
return zz
Plastic:
-- Plastic ratio constant
-- Real root of x^3-x-1
procedure
-- Fast
if Digits() < 111 then
return 1.324717957244746025960908854478097340734404056901733364534015050302827851245547594054699347981787280329910921
-- Polynomial and derivatives
p = '1 0 -1 -1'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.36*d)
end
d.n = 4
-- Halley
zz = 1.3
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Porter:
-- Porter constant
procedure
-- Fast
return 1.46707807943397547289779848470722995344990332241488777739968581761660674432904480843036932751117401521266
Prouhet:
-- Prouhet-Thue-Morse constant
procedure
-- Fast
if Digits() < 111 then
return 0.41245403364010759778336136825845528308947837445576955757337941534879359236578258896380454048621213339625636645
-- Infinite product
numeric digits Digits()+2
zz = 1; v = 1
do n = 0
zz = zz*(1-1/2**(2**n))
if zz = v
then leave
v = zz
end
return (2-zz)*0.25
Ramanujan:
-- Ramanujan constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 262537412640768743.9999999999992500725971981856888793538563373369908627075374103782106479101186073129511813462
-- Formula
return Exp(Pi()*Sqrt(163))
Recifibo:
-- Reciprocal Fibonacci constant
procedure
-- Fast
if Digits() < 111 then
return 3.359885666243177553172011302918927179688905133731968486495553815325130318996683383615416216456790087297045342928854
-- Recurring
numeric digits Digits()+2
zz = 1; a = 0; b = 1
numeric fuzz 2
do n = 2 until zz = v
v = zz
c = a+b; zz = zz+1/c
a = b; b = c
end
return zz
Robbins:
-- Robbins constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.6617071822671762351558311332484135817464001357909536048089442294795846461385976313066524807681071201517097753107594
-- Formula
return (4+17*Sqrt(2)-6*Sqrt(3)-7*Pi())/105+(Ln(1+Sqrt(2)))*0.2+(Ln(2+Sqrt(3)))*0.4
Salem:
-- Salem constant
-- Largest positive root of x^10+x^9-x^7-x^6-x^5-x^4-x^3+x+1
procedure
-- Fast
if Digits() < 111 then
return 1.1762808182599175065440703384740350506934158065646952598301063470296883765485499620968301155818153946592071814
-- Polynomial and derivatives
p = '1 1 0 -1 -1 -1 -1 -1 0 1 1'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.39*d)
end
d.n = 4
-- Halley
zz = 1.2
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Sierpinski:
-- Sierpinski constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 2.5849817595792532170658935873831711600880516518526309173215449879719320440011571202111177245270642830313440
-- Formula
return Pi()*Ln(Exp(2*Euler())/(2*Gauss()**2))
Silver:
-- Silver ratio constant
-- Positive root of x^2-2x-1
procedure
-- Fast
if Digits() < 111 then
return 2.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641572735013846
-- Formula
return Sqrt(2)+1
Soldner:
-- Soldner constant
procedure
-- Fast
return 1.45136923488338105028396848589202744949303228364801586309300455766242559575451783565953135771108682884
Somos:
-- Somos quadratic recurrence constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 1.661687949633594121295818922749950749964418635025068208189711168025609029826383727908369176411461167155281345509912
-- Definition
numeric digits Digits()+2
zz = 1; v = 0; a = 2
do n = 1
zz = zz*Power(n,1/a)
if zz = v
then leave
v = zz; a = 2*a
end
return zz
Sqrt2:
-- Square root of 2 constant
-- Real root of x^2-2
procedure
-- Fast
if Digits() < 111 then
return 1.4142135623730950488016887242096980785696718753769480731766797379907324784621070388503875343276415727350138462
-- Dynamic precision
numeric digits Digits()+2
d = Digits()
do n = 1 while d > 4
d.n = d; d = Trunc(0.36*d)
end
d.n = 4
-- Halley
zz = 1.4
do k = n by -1 to 1
numeric digits d.k
z2 = zz*zz; zz = zz-2*zz*(z2-2)/(3*z2+2)
end
return zz
Sqrt3:
-- Square root of 3 constant
-- Real root of x^2-3
procedure
-- Fast
if Digits() < 111 then
return 1.7320508075688772935274463415058723669428052538103806280558069794519330169088000370811461867572485756756261414
-- Dynamic precision
numeric digits Digits()+2
d = Digits()
do n = 1 while d > 4
d.n = d; d = Trunc(0.36*d)
end
d.n = 4
-- Halley
zz = 1.7
do k = n by -1 to 1
numeric digits d.k
z2 = zz*zz; zz = zz-2*zz*(z2-3)/(3*z2+3)
end
return zz
Sqrt5:
-- Square root of 5 constant
-- Real root of x^2-5
procedure
-- Fast
if Digits() < 111 then
return 2.2360679774997896964091736687312762354406183596115257242708972454105209256378048994144144083787822749695081762
-- Dynamic precision
numeric digits Digits()+2
d = Digits()
do n = 1 while d > 4
d.n = d; d = Trunc(0.36*d)
end
d.n = 4
-- Halley
zz = 2.2
do k = n by -1 to 1
numeric digits d.k
z2 = zz*zz; zz = zz-2*zz*(z2-5)/(3*z2+5)
end
return zz
Stephens:
-- Stephens constant
procedure
-- Fast
return 0.57595996889294543964316337549249669250651396717649236360064079866537255169886852843640987209172618
Supergolden:
-- Super golden ratio constant
-- Real root of x^3-x^2-1
procedure
-- Fast
if Digits() < 111 then
return 1.465571231876768026656731225219939108025577568472285701643183111249262996685017840478125801194909270064382456
-- Polynomial and derivatives
p = '1 -1 0 -1'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.36*d)
end
d.n = 4
-- Halley
zz = 1.5
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Taniguchi:
-- Taniguchi constant
procedure
-- Fast
return 0.678234491917391978035538279482894814096332239189440103036460415964983370740123233213762122933484616888328
Tau:
-- Tau constant
-- 2pi
procedure expose glob.
-- Fast
if Digits() < 111 then
return 6.283185307179586476925286766559005768394338798750211641949889184615632812572417997256069650684234135964296173
-- Formula
return 2*Pi()
Tetranacci:
-- Tetranacci constant
-- Root of x^4-x^3-x^2-x-1
procedure
-- Fast
if Digits() < 111 then
return 1.927561975482925304261905861736622168698554255163384727146647038009666062297815559149818253461890653254266934313
-- Polynomial and derivatives
p = '1 -1 -1 -1 -1'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.38*d)
end
d.n = 4
-- Halley
zz = 1.9
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Tribonacci:
-- Tribonacci constant
-- Root of x^3-x^2-x-1
procedure
-- Fast
if Digits() < 111 then
return 1.8392867552141611325518525646532866004241787460975922467787586394042032220819664257384354194283070141419798269
-- Polynomial and derivatives
p = '1 -1 -1 -1'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.38*d)
end
d.n = 4
-- Halley
zz = 1.8
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Twinprimes:
-- Twin primes constant
procedure
-- Fast
return 0.660161815846869573927812110014555778432623360284733413319448423335405642304495277143760031413839867911779
Vanderpauw:
-- Van der Pauw constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 4.5323601418271938096276829457166668101718614677237955841860165479406009537213051022590838796040160896530664627078
-- Formula
return Pi()/Ln(2)
Viswanath:
-- Viswanath constant
procedure
-- Fast
return 1.1319882487943
Wallis:
-- Wallis constant
-- Real root of x^3-2x-5
procedure
-- Fast
if Digits() < 111 then
return 2.094551481542326591482386540579302963857306105628239180304128529045312189983483667146267281777157757860839521
-- Polynomial and derivatives
p = '1 0 -2 -5'
p1 = Pdif(p); p2 = Pdif(p1)
-- Dynamic precision
d = Digits()+2
do n = 1 while d > 4
d.n = d; d = Trunc(0.38*d)
end
d.n = 4
-- Halley
zz = 2.1
do k = n by -1 to 1
numeric digits d.k
f = Peval(p,zz); f1 = Peval(p1,zz); f2 = Peval(p2,zz)
zz = zz-(2*f*f1)/(2*f1*f1-f*f2)
end
return zz
Weierstrass:
-- Weierstrass constant
procedure expose glob.
-- Fast
if Digits() < 111 then
return 0.4749493799879206503325046363279829685595493732172029822833310248645579291748838602742756412505021444189037849426240
-- Formula
return Power(2,1.25)*Sqrt(Pi())*Exp(Pi()*0.125)/Gamma(0.25)**2
Zscore:
-- Z-score 97.5 percentile constant
procedure
-- Fast
return 1.95996398454005423552459443052055152795555007786954839847695264636163527414488266779825470949281420601799
-- Module Functions.inc - 11 Apr 2025
-- Standard and advanced mathematical functions
-- This module supports only real arguments and results
Ackermann:
-- Ackermann's function
procedure
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
-- Fast
if xx = 0 then
return yy+1
if xx = 1 then
return yy+2
if xx = 2 then
return yy+yy+3
if xx = 3 then
return 2**(yy+3)-3
if xx = 4 then do
if yy = 0 then
return 2**(2**2)-3
if yy = 1 then
return 2**(2**(2**2))-3
if yy = 2 then
return 2**(2**(2**(2**2)))-3
end
-- Definition
if yy = 0 then
return Ackermann(xx-1,1)
else
return Ackermann(xx-1,Ackermann(xx,yy-1))
Agmean:
-- Arithmetic geometric mean
-- AGM
procedure
arg xx,yy
if \ Datatype(xx,'n') then say argument
if xx <= 0 then say argument
if \ Datatype(yy,'n') then say argument
if yy < 0 then say argument
-- Fast
if yy = 0 then
return 0
if xx = yy then
return xx
-- Series
numeric digits Digits()+1
vx = xx+1
do until xx = vx
vx = xx; vy = yy
xx = (vx+vy)*0.5; yy = Sqrt(vx*vy)
end
return xx
Aliquot:
-- Aliquot function
-- Sum of proper divisors
procedure
arg xx
-- Formula
return Divisor(xx,1)-xx
Alog10:
-- Antilog base 10
-- 10^x
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return Power(10,xx)
Amean:
-- Arithmetic mean
-- AM
procedure
nn = Arg()
if nn = 0 then say argument
-- Cf definition
numeric digits Digits()+1
zz = 0
do n = 1 to nn
xx = Arg(n)
if \ Datatype(xx,'n') then say argument
zz = zz+xx
end
return zz/nn
Arccos:
-- Inverse cosine
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if Abs(xx) > 1 then say argument
-- Reflect
numeric digits Digits()+1
if xx < 0 then
return Pi()-Arccos(-xx)
-- Fast
if xx = 0 then
return Halfpi()
if xx = 0.5 then
return Pi()/3
if xx = 1 then
return 0
-- Formula
if xx < 0.71 then
return Halfpi()-Arcsin(xx)
if xx < 0.81 then do
numeric digits Digits()+1
return Arctan(Sqrt(1-xx*xx)/xx)
end
else do
numeric digits Digits()-Xpon(1-xx)
return Arcsin(Sqrt(1-xx*xx))
end
Arccot:
-- Inverse cotangent
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Reflect
numeric digits Digits()+1
if xx < 0 then
return Pi()-Arccot(-xx)
-- Fast
if xx = 0 then
return Halfpi()
if xx = 1 then
return 0.25*Pi()
-- Formula
if xx < 0.40 then
return Arccos(xx/Sqrt(xx*xx+1))
else
return Arctan(1/xx)
Arccsc:
-- Inverse cosecant
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if Abs(xx) < 1 then say argument
-- Reflect
if xx < 0 then
return -Arccsc(-xx)
-- Fast
numeric digits Digits()+1
if xx = 1 then
return Halfpi()
if xx = 2 then
return Pi()/6
-- Formula
return Arcsin(1/xx)
Arcosh:
-- Inverse cosine hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 1 then say argument
-- Fast
if xx = 1 then
return 0
-- Formula
numeric digits Digits()-(xx<1.1)*Xpon(xx-1)+2
return Ln(xx+Sqrt(xx*xx-1))
Arcoth:
-- Inverse cotangent hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx <= 1 then say argument
-- Formula
numeric digits Digits()+Xpon(xx)+1
return 0.5*Ln((xx+1)/(xx-1))
Arcsch:
-- Inverse cosecant hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx = 0 then say argument
-- Formula
numeric digits Digits()+Max(Xpon(xx),0)+3
zz = 1/xx
return Ln(zz+Sqrt(zz*zz+1))
Arcsec:
-- Inverse secant
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if Abs(xx) < 1 then say argument
-- Reflect
numeric digits Digits()+1
if xx < 0 then
return Pi()-Arcsec(-xx)
-- Fast
if xx = 1 then
return 0
if xx = 2 then
return Pi()/3
-- Formula
return Arccos(1/xx)
Arcsin:
-- Inverse sine
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if Abs(xx) > 1 then say argument
-- Reflect
if xx < 0 then
return -Arcsin(-xx)
-- Fast
numeric digits Digits()+2
if xx = 0 then
return 0
if xx = 0.5 then
return Pi()/6
if xx = 1 then
return Halfpi()
-- Formula
if xx > 0.80 then do
if xx > 0.90 then
numeric digits Digits()-Xpon(1-xx)
return Halfpi()-Arcsin(Sqrt(1-xx*xx))
end
-- Taylor
t = xx; zz = xx; xx = xx*xx
do n = 2 by 2 until zz = v
v = zz
t = t*xx*(n-1)/n; zz = zz+t/(n+1)
end
return zz
Arctan:
-- Inverse tangent
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Reflect
if xx < 0 then
return -Arctan(-xx)
-- Fast
numeric digits Digits()+3
if xx = 0 then
return 0
if xx = 1 then
return 0.25*Pi()
-- Formula
if xx > 3 then do
xx = 1/xx
return Halfpi()-Arcsin(xx/Sqrt(xx*xx+1))
end
-- Reduce argument towards zero
r = 2-(xx<2)-(xx<0.70)
do r
xx = xx/(1+Sqrt(xx*xx+1))
end
-- Taylor
t = xx; zz = xx; xx = xx*xx
do n = 3 by 2 until zz = v
v = zz
t = -t*xx; zz = zz+t/n
end
-- Undo reduce
return 2**r*zz
Arctan2:
-- Inverse tangent right quadrant
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'n') then say argument
-- Fast
numeric digits Digits()+1
if yy = 0 then
return 0
if xx = 0 then
return 0.5*Sign(yy)*Pi()
-- Formula
if xx > 0 then
return Arctan(yy/xx)
if yy < 0 then
return Arctan(yy/xx)-Pi()
else
return Arctan(yy/xx)+Pi()
Arsech:
-- Inverse secant hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx <= 0 then say argument
if xx > 1 then say argument
-- Fast
if xx = 1 then
return 0
-- Formula
numeric digits Digits()-(xx>0.90)*Xpon(1-xx)+4
zz = 1/xx
return Ln(zz+Sqrt(zz*zz-1))
Arsinh:
-- Inverse sine hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Reflect
if xx < 0 then
return -Arsinh(-xx)
-- Fast
if xx = 0 then
return 0
-- Formula
numeric digits Digits()-(xx<1)*Xpon(xx)+2
if xx > 0.50 then
return Ln(xx+Sqrt(xx*xx+1))
-- Taylor
t = xx; zz = t; x2 = xx*xx
numeric fuzz 2
do n = 2 by 2 until zz = v
v = zz
t = -t*x2*(n-1)*(n-1)/(n*n+n); zz = zz+t
end
return zz
Artanh:
-- Inverse tangent hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if Abs(xx) >= 1 then say argument
-- Fast
if xx = 0 then
return 0
-- Reflect
if xx < 0 then
return -Artanh(-xx)
-- Formula
numeric digits Digits()-Xpon(1-xx)+1
if xx > 0.40 then
return 0.5*Ln((1+xx)/(1-xx))
-- Taylor
numeric digits Digits()+2
t = xx; zz = t; x2 = xx*xx
numeric fuzz 2
do n = 3 by 2 until zz = v
v = zz
t = t*x2*(n-2)/n; zz = zz+t
end
return zz
Beta:
-- Beta
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'n') then say argument
-- Formula
numeric digits Digits()+1
return Gamma(xx)*Gamma(yy)/Gamma(xx+yy)
Base10:
-- Base n -> 10
procedure
arg xx,yy
if xx = '' then say argument
xx = Upper(xx)
if \ Datatype(yy,'w') then say argument
if yy < 2 then say argument
if yy > 36 then say argument
-- Fast
if yy = 2 then
return X2d(B2x(xx))+0
if yy = 16 then
return X2d(xx)+0
-- Convert
b = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'; l = Length(xx)
zz = 0
do i = 1 to l
d = Pos(Substr(xx,i,1),b)-1
if d = -1 | d >= yy then say argument
zz = zz*yy+d
end
return zz
Basenn:
-- Base 10 -> n
procedure
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 2 then say argument
if yy > 36 then say argument
-- Fast
if xx = 1 then
return xx
if yy = 2 then
return X2b(D2x(xx))+0
if yy = 16 then
return D2x(xx)
-- Convert
b = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'
zz = ''
do while xx > 0
r = xx//yy; xx = xx%yy; zz = Substr(b,r+1,1)zz
end
return zz
Cbrt:
-- Cubic root
-- x^(1/3)
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 0
if xx = 1 then
return 1
-- Sign
sx = Sign(xx); xx = Abs(xx)
-- Reduce argument to 1...1000
e = Xpon(xx); e = (e-(2*(e<0)))%3; f = -3*e; xx = Std(xx)'E'f
-- First guess 2 digits accurate
dd = Digits()+2
numeric digits 2
select
when xx < 8 then
zz = 1+(xx-1)/7
when xx < 27 then
zz = 2+(xx-8)/19
when xx < 64 then
zz = 3+(xx-27)/37
when xx < 125 then
zz = 4+(xx-64)/61
when xx < 216 then
zz = 5+(xx-125)/91
when xx < 343 then
zz = 6+(xx-216)/127
when xx < 512 then
zz = 7+(xx-343)/169
when xx < 729 then
zz = 8+(xx-512)/217
otherwise
zz = 9+(xx-729)/271
end
-- Dynamic precision
do n = 1 while dd > 4
d.n = dd; dd = Trunc(0.42*dd)
end
d.n = 4
-- Halley
do i = n to 1 by -1
numeric digits d.i
z3 = zz*zz*zz; zz = zz-zz*(z3-xx)/(2*z3+xx)
end
-- Undo reduce
return sx*zz'E'e
Ceil:
-- Ceiling
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
if Integer(xx) then
return xx
else
return Trunc(xx)+(xx>=0)
Comb:
-- Combination
-- Binomial coefficient
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
if xx-yy < yy then
yy = xx-yy
-- Fast
if yy = 0 then
return 1
if glob.comb.xx.yy <> '' then
return glob.comb.xx.yy
-- Recurring loop
numeric digits Digits()+2
zz = xx
do i = 1 to yy-1
zz = zz*(xx-i)/(i+1)
end
glob.comb.xx.yy = zz
return zz
Cos:
-- Cosine
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 1
-- Reflect
if xx < 0 then
return Cos(-xx)
-- Reduce argument to 0...2pi
numeric digits Digits()+Max(Xpon(xx),0)+2
xx = xx//Tau()
-- Reduce argument to 0...pi/2
select
when xx < Halfpi() then
s = 1
when xx < Pi() then do
s = -1; xx = Pi()-xx
end
when xx < Pi()+Halfpi() then do
s = -1; xx = xx-Pi()
end
otherwise do
s = 1; xx = tau()-xx
end
end
-- Taylor
t = 1; zz = 1; xx = xx*xx
do n = 2 by 2 until zz = v
v = zz
t = -t*xx/(n*(n-1)); zz = zz+t
end
return s*zz
Cosh:
-- Cosine hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Reflect
if xx < 0 then
return Cosh(-xx)
-- Fast
if xx = 0 then
return 1
-- Taylor
numeric digits Digits()+1
if \ Integer(xx) & xx < 5 then do
xx = xx*xx; zz = 1; f = 1; v = 1
do n = 2 by 2 until zz = v
v = zz
f = f*xx/((n-1)*n); zz = zz+f
end
return zz
end
-- Formula
zz = Exp(xx)
return (zz*zz+1)/(2*zz)
Cot:
-- Cotangent
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
zz = Sin(xx)
if zz = 0 then say argument
return Cos(xx)/zz
Cotan:
-- Cotangent
procedure expose glob.
arg xx
-- Alias
return Cot(xx)
Coth:
-- Cotangent hyperbolic
procedure expose glob.
arg xx
if xx = 0 then say argument
-- Formula
numeric digits Digits()+1
return 1/Tanh(xx)
Csc:
-- Cosecant
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
zz = Sin(xx)
if zz = 0 then say argument
return 1/zz
Csch:
-- Cosecant hyperbolic
procedure expose glob.
arg xx
if xx = 0 then say argument
-- Formula
numeric digits Digits()+1
return 1/Sinh(xx)
Deg:
-- Radians -> Degrees
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
numeric digits Digits()+1
return xx*180/Pi()
Dfact:
-- Double factorial
-- n!!
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if glob.dfact.xx <> '' then
return glob.dfact.xx
if xx < 2 then
return 1
-- Loops cf definition
if xx//2 = 1 then do
zz = 1
do n = 3 to xx by 2
zz = zz*n
end
end
else do
zz = 2
do n = 4 to xx by 2
zz = zz*n
end
end
glob.dfact.xx = zz
return zz
Digamma:
-- Digamma
-- Logarithmic derivative of Gamma
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx <= 0 & Integer(xx) then say argument
-- Formula
numeric digits Digits()+2
if xx = 0.125 then
return -Halfpi()-4*Ln(2)-(Pi()+Ln(Sqrt(2)+1)-Ln(Sqrt(2)-1))/Sqrt(2)-Euler()
if xx = 0.25 then
return -Halfpi()-3*Ln(2)-Euler()
if xx = 0.5 then
return -2*Ln(2)-Euler()
if xx = 1 then
return -Euler()
if xx > 0 & xx < 500 & Integer(xx) then
return Harmonic(xx-1)-Euler()
-- Fast series
if xx > 0 & Half(xx) then do
xx = Trunc(xx)
zz = 0
do n = 1 to xx
zz = zz+2/(2*n-1)
end
zz = zz-2*Ln(2)-Euler()
return zz
end
-- Reduce argument to a larger value
numeric digits Digits()+1
dd = Digits(); a = 0
do while xx < dd
a = a+1/xx; xx = xx+1
end
-- Series
zz = -0.5/xx; v = 0
do n = 2 by 2 until zz = v
v = zz
zz = zz+Zeta(1-n)/(xx**n)
end
-- Undo reduce
return Ln(xx)+zz-a
Digit:
-- Round to specified digits
procedure
arg xx,yy
if yy = '' then
yy = Digits()
-- Formula
return Round(Mant(xx),yy-1)'E'Xpon(xx)
Digitprod:
-- Digital product
-- Product of all digits
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Multiply digits
zz = 1
do n = 1 to Length(xx) while zz <> 0
zz = zz*Substr(xx,n,1)
end
return zz
Digitroot:
-- Digital root
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Formula
return 1+(xx-1)//9
Digitsum:
-- Digital sum
-- Sum of all digits
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Sum digits
zz = 0
do n = 1 to Length(xx)
zz = zz+Substr(xx,n,1)
end
return zz
Divisor:
-- Divisor
-- Som of all divisors^y
procedure expose glob. divi.
arg xx,yy
if yy = '' then
yy = 0
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
if \ Datatype(yy,'n') then say argument
numeric digits Digits()+2
select
-- Count divisors
when yy = 0 then do
n = 0
do while xx//2 = 0
n = n+1; xx = xx/2
end
zz = n+1
do i = 3 by 2 while i*i <= xx
n = 0
do while xx//i = 0
n = n+1; xx = xx/i
end
zz = zz*(n+1)
end
if xx > 2 then
zz = zz+zz
return zz
end
-- Get divisors and sum
when yy = 1 then do
d = Divisors(xx); zz = 0
do n = 1 to d
zz = zz+divi.n
end
end
-- Get divisors and sum integer power
when Integer(yy) then do
d = Divisors(xx); zz = 0
do n = 1 to d
zz = zz+divi.n**yy
end
end
-- Get divisors and sum real power
otherwise do
d = Divisors(xx); zz = 0
do n = 1 to d
zz = zz+Power(divi.n,yy)
end
end
end
return zz
Dxgamma:
-- First derivative of Gamma
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx <= 0 & Integer(xx) then say argument
-- Formula
dd = -Digits()
numeric digits -dd*2
a = 1'E'dd
return (Gamma(xx+a)-Gamma(xx))/a
Dxzeta:
-- First derivative of Zeta
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx = 1 then say argument
-- Formula
dd = Digits()
numeric digits dd+1
if xx = -2 then
return -Zeta(3)/(4*Pi()**2)
if xx = -1 then
return 1/12-Ln(Glaisher())
if xx < 0 & Even(xx) then
return (-1)**(-0.5*xx)*Fact(-xx)*Zeta(1-xx)/(2*(Tau())**-xx)
if xx = 0 then
return -0.5*Ln(Tau())
if xx = 0.5 then
return 0.25*(Pi()+2*Euler()+6*Ln(2)+2*Ln(Pi()))*Zeta(0.5)
if xx = 2 then
return Pi()*Pi()*(Euler()+Ln(2)-12*Ln(Glaisher())+Ln(Pi()))/6
if xx > 5 then do
-- Taylor
zz = 0
a = Power(2,xx); b = Ln(4); c = a-2
do n = 1 to 10000 until zz = v
v = zz
if n//2 = 0 then
t1 = 1
else
t1 = -1
zz = zz+(t1*a*Power(n,-xx))*(b+c*Ln(n))/(c*c)
end
end
else do
-- From derivative definition
numeric digits dd*4
e = -dd+dd%5; a = 1'E'e; zz = (Zeta(xx+a)-Zeta(xx))/a
end
return zz
Erf:
-- Error
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Asymptotic series
numeric digits Digits()+2
a = (2*xx*Exp(-xx*xx))/Sqrt(Pi()); b = 2*xx*xx
zz = 0; f = 1; g = 1
do n = 0 until zz = v
v = zz
zz = zz+f/g
f = f*b; g = g*(2*n+3)
end
return a*zz
Erfc:
-- Error complementary
procedure expose glob.
arg xx
-- Formula
return 1-Erf(xx)
Eval:
-- Evaluation of an expression in x,y,z
procedure expose glob.
arg xx,x,y,z
if xx = '' then say argument
-- Evaluate
interpret 'return' xx
Even:
-- Is a number even?
procedure
arg xx
-- Formula
if Integer(xx) then
return \ Abs(xx//2)
else
return 0
Exp:
-- Exponential
-- e^x
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 1
if xx = 1 then
return E()
dd = Digits()
if glob.exp.xx.dd <> '' then
return glob.exp.xx.dd
-- Formula
numeric digits dd+3
if Integer(xx) then do
zz = E()**xx; glob.exp.xx.dd = zz
return zz
end
-- Reduce argument to -1/2...1/2
if Abs(xx) > 1 then do
a = Trunc(xx); xx = Frac(xx)
end
else
a = 0
if Abs(xx) > 0.5 then do
b = 1; xx = 0.5*xx
end
else
b = 0
-- Taylor
zz = 1; t = 1
do n = 1 until zz = v
v = zz
t = (t*xx)/n; zz = zz+t
end
-- Undo reduce
if b <> 0 then
zz = zz*zz
if a <> 0 then
zz = zz*E()**a
glob.exp.xx.dd = zz
return zz
Fact:
-- Integer factorial
-- n!
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return 1
if glob.fact.xx <> '' then
return glob.fact.xx
w = xx-1
-- Loop cf definition
if glob.fact.w = '' then do
zz = 1
do n = 2 to xx
zz = zz*n
end
glob.fact.xx = zz
end
-- Multiply
else
glob.fact.xx = glob.fact.w*xx
return glob.fact.xx
Factorial:
-- Real factorial
-- x!
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 0 & Integer(xx) then say argument
-- Formula
numeric digits Digits()+1
if Integer(xx) then
return Fact(xx)
else
return Gamma(xx+1)
Ffact:
-- Falling factorial
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if \ Datatype(yy,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return 1
if yy = 0 then
return 1
if glob.ffact.xx.yy <> '' then
return glob.ffact.xx.yy
-- Calculate cf definition
zz = 1
do k = 0 to yy-1
zz = zz*(xx-k)
end
glob.ffact.xx.yy = zz
return zz
Floor:
-- Floor
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
if Integer(xx) then
return xx
else
return Trunc(xx)-(xx<0)
Frac:
-- Fractional part
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
return xx-xx%1
Gamma:
-- Gamma
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 1 & Integer(xx) then say argument
-- Formula
if xx < 0 then do
if Half(xx) then do
i = Abs(Floor(xx)); zz = (-1)**i*2**(2*i)*Fact(i)*Sqrt(Pi())/Fact(2*i)
return zz
end
end
if xx > 0 then do
if Integer(xx) then
return Fact(xx-1)
if Half(xx) then do
i = Floor(xx); zz = Fact(2*i)*Sqrt(Pi())/(2**(2*i)*Fact(i))
return zz
end
end
d1 = Digits()
if d1 < 61 then do
-- Lanczos with predefined coefficients
numeric digits d1+2
-- Reflect
if xx < 0 then
return Pi()/(Gamma(1-xx)*Sin(Pi()*xx))
-- Reduce argument to 0.5...1.5
c = Trunc(xx); xx = xx-c
if xx < 0.5 then do
xx = xx+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
xx = xx-1; s = 0
do k = 52 by -1 to 1
s = s*xx+c.k
end
zz = 1/s
-- Undo reduce
if c = -1 then
zz = zz/xx
else do
do i = 1 to c
zz = (xx+i)*zz
end
end
end
else do
xx = xx-1
-- Spouge
-- Estimate digits and iterations
d2 = Floor(1.5*d1); a = Floor(1.1*d1)
numeric digits d2
-- 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)*(xx+k))
end
s = s+Sqrt(Tau()); zz = Power(xx+a,xx+0.5)*Exp(-a-xx)*s
end
return zz
Gcd:
-- Greatest common divisor
procedure
arg xx,yy
if \ Datatype(xx,'n') then say argument
xx = Abs(xx); yy = Abs(yy)
-- Fast
if yy = 0 then
return xx
-- Euclidian division
do until yy = 0
zz = xx//yy; xx = yy; yy = zz
end
return xx
Gmean:
-- Geometric mean
-- GM
procedure
nn = Arg()
if nn = 0 then say argument
-- Cf definition
numeric digits Digits()+1
zz = 1
do n = 1 to nn
xx = Arg(n)
if \ Datatype(xx,'n') then say argument n
if xx <= 0 then say argument n
zz = zz*xx
end
return Nroot(zz,nn)
Half:
-- Is a number half integer?
procedure
arg xx
-- Formula
return (Frac(Abs(xx))=0.5)
Hmean:
-- Harmonic mean
-- HM
procedure
nn = Arg()
if nn = 0 then say argument
-- Cf definition
numeric digits Digits()+1
zz = 0
do n = 1 to nn
xx = Arg(n)
if \ Datatype(xx,'n') then say argument n
if xx = 0 then say argument n
zz = zz+1/xx
end
if zz = 0 then say argument
return nn/zz
Hypot:
-- Hypothenusa
procedure
arg xx,yy
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'n') then say argument
-- Formula
numeric digits Digits()+1
return Sqrt(xx*xx+yy*yy)
Integer:
-- Is a number integer?
procedure
arg xx
-- Formula
return Datatype(xx,'w')
Isqrt:
-- Integer square root
procedure
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 0 then say argument
-- Fast
if xx < 1 then
return 0
if xx < 4 then
return 1
if xx < 9 then
return 2
-- Loop uses only integers
xx = xx%1; q = 1
do until q > xx
q = q*4
end
zz = 0
do while q > 1
q = q%4; t = xx-zz-q; zz = zz%2
if t >= 0 then do
xx = t; zz = zz+q
end
end
return zz
Lambda:
-- Dirichlet lambda
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
numeric digits Digits()+1
return (1-Power(2,-xx))*Zeta(xx)
Lambertw0:
-- Lambert W0
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 0
-- First guess
dd = Digits()+1
numeric digits 2
select
when xx > E() then do
a = Ln(xx); zz = a-Ln(a)
end
when xx > 0 then
zz = xx/E()
when xx >= -1/E() then do
a = E()*xx; b = a+1; c = Sqrt(b); zz = (a*Ln(1+c))/(b+c)
end
otherwise
say xx argument
end
-- Dynamic precision
do n = 1 while dd > 3
d.n = dd; dd = Trunc(0.50*dd)
end
d.n = 3
-- Halley
do i = n to 1 by -1
numeric digits d.i
a = Exp(zz); b = a*zz; zz = zz-(b-xx)/(a+b-((zz+2)*(b-xx))/(2*zz+2))
end
return zz
Lambertwm1:
-- Lambert W-1
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx < -1/E() then say argument
if xx >= 0 then say argument
-- First guess
dd = Digits()+1
numeric digits 2
select
when xx > -0.25 then do
a = Ln(-xx); zz = a-Ln(-a)
end
when xx > -1/E() then do
zz = -1-Sqrt(2)*Sqrt(E()*xx+1)
end
end
-- Dynamic precision
do n = 1 while dd > 3
d.n = dd; dd = Trunc(0.50*dd)
end
d.n = 3
-- Halley
do i = n to 1 by -1
numeric digits d.i
a = Exp(zz); b = a*zz; zz = zz-(b-xx)/(a+b-((zz+2)*(b-xx))/(2*zz+2))
end
return zz
Lcm:
-- Least common multiple function
procedure
arg xx,yy
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'n') then say argument
xx = Abs(xx); yy = Abs(yy)
-- Fast
if yy = 0 then
return 0
-- Euclidian division
a = xx*yy
do until yy = 0
zz = xx//yy; xx = yy; yy = zz
end
return a%xx
Ln:
-- Natural logarithm base e
-- Log
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx <= 0 then say argument
-- Fast
if xx = 1 then
return 0
if xx = 2 then
return Ln2()
if xx = 4 then
return Ln4()
if xx = 10 then
return Ln10()
dd = Digits()
if glob.ln.xx.dd <> '' then
return glob.ln.xx.dd
-- Reduce argument to 0.5...1.5
numeric digits dd+2
do i = 0 by -1 while xx < 0.5
xx = xx+xx
end
do i = i while xx > 1.5
xx = xx*0.5
end
-- Taylor
q = (xx-1)/(xx+1); f = q; zz = q; v = q; q = q*q
do n = 3 by 2 until zz = v
v = zz
f = f*q; zz = zz+f/n
end
-- Undo reduce
zz = 2*zz+i*Ln2(); glob.ln.xx.dd = zz
return zz
Lngamma:
-- Logarithmic gamma
procedure expose glob.
arg xx
-- Formula
return Ln(Gamma(xx))
Log:
-- Natural logarithm base e
-- Ln
procedure expose glob.
arg xx
-- Alias
return Ln(xx)
Log10:
-- Common logarithm base 10
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return Ln(xx)/Ln(10)
Logxy:
-- Logarithm x base y
procedure expose glob.
arg xx,yy
if yy = 1 then say argument
-- Formula
numeric digits Digits()+1
return Ln(xx)/Ln(yy)
Mant:
-- Mantissa
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 0
-- Formula
zz = xx*1E+99999999
return Left(zz,Pos('E',zz)-1)
Natural:
-- Are all arguments in N?
procedure
do i = 1 to Arg()
if \ Datatype(Arg(i),'w') then say argument
if Arg(i) < 1 then say argument
end
return 1
Nroot:
-- Nth root
-- x^(1/n)
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'w') then say argument
if yy < 1 then say argument
if xx < 0 & Even(yy) then say argument
-- Fast
if xx = 0 then
return 0
if xx = 1 then
return 1
if yy = 2 then
return Sqrt(xx)
if yy = 3 then
return Cbrt(xx)
-- Sign
sx = Sign(xx); xx = Abs(xx)
-- First guess 3 digits accurate
d1 = Digits()
numeric digits 3
zz = Exp(Ln(xx)/yy)/1
-- Dynamic precision
d2 = d1+2
numeric digits d2
d = d2
do k = 1 while d > 5
d.k = d; d = Trunc(0.60*d)
end
d.k = 5
-- Newton
yi = 1/yy; y1 = yy-1
do i = k to 1 by -1
numeric digits d.i
zz = zz - yi*(zz-xx/zz**y1)
end
return sx*zz
Odd:
-- Is a number odd?
procedure
arg xx
-- Formula
if Integer(xx) then
return Abs(xx//2)
else
return 0
Perm:
-- Permutations
procedure
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
if xx < yy then say argument
-- Cf definition
zz = 1
do n = (xx-yy)+1 to xx
zz = zz*n
end
return zz
Power:
-- Real power
-- x^y
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'n') then say argument
if xx < 0 & \ Integer(yy) then say argument
-- Fast
numeric digits Digits()+1
if xx = 0 then
return 0
if xx = 1 then
return 1
if yy = -1 then
return 1/xx
if yy = 0 then
return 1
if yy = 1 then
return xx
-- Formula
if Integer(yy) then
return xx**yy
if yy = 0.5 then
return Sqrt(xx)
if Half(yy) then
return Sqrt(xx)**Sign(yy)*xx**(yy%1)
else
return Exp(yy*Ln(xx))
Powmod:
-- Power modulus
procedure
arg xx,yy,zz
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
if \ Datatype(zz,'w') then say argument
if zz < 1 then say argument
-- Fast
if zz = 1 then
return 0
-- Binary method
numeric digits Max(Length(Format(xx,,,0)),Length(Format(yy,,,0)),2*Length(Format(zz,,,0)))
b = xx//zz; rr = 1
do while yy > 0
if yy//2 then
rr = rr*xx//zz
xx = xx*xx//zz; yy = yy%2
end
return rr
Primecount:
-- Prime count
-- Pi(x) estimated number of primes below x
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx = 1 then
return 0
-- Formula
return Round(xx/(Ln(xx)-1),0)
Primeest:
-- Nth prime number estimate
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Estimate
select
when xx <= 5 then do
zz = xx+xx+1
end
when xx <= 12 then do
a = Ln(xx); b = Ln(a); zz = xx*(a+b)
end
when xx <= 200 then do
a = Ln(xx); b = Ln(a); zz = xx*(a+b-1+1.8*b/a)
end
otherwise do
a = Ln(xx); b = Ln(a); zz = xx*(a+b-0.9385)+100
end
end
return Trunc(zz)
Primorial:
-- Primorial
-- p# product of primes
procedure expose prim. glob.
arg xx
if \ Datatype(xx,'w') then say argument
-- Fast
if xx = 0 then
return 1
if Abs(xx) = 1 then
return 2
-- Get primes
if xx < 0 then do
call Primes(xx); p = -xx
end
else
p = Primes(xx)
-- Calculate product
numeric digits Digits()+2
zz = 1
do n = 1 to p
zz = zz*prim.n
end
return zz
Qmean:
-- Quadratic mean
-- QM, rms
procedure
nn = Arg()
if nn = 0 then say argument
-- Cf definition
numeric digits Digits()+1
zz = 0
do n = 1 to nn
xx = Arg(n)
if \ Datatype(xx,'n') then say argument n
zz = zz+xx*xx
end
return Sqrt(zz/nn)
Qtrt:
-- Quartic root
-- x^(1/4)
procedure
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 0
if xx = 1 then
return 1
-- Reduce argument to 1...10000
e = Xpon(xx); e = (e-(3*(e<0)))%4; f = -4*e; xx = Std(xx)'E'f
-- First guess 2 digits accurate
dd = Digits()+2
numeric digits 2
select
when xx < 16 then
zz = 1+(xx-1)/15
when xx < 81 then
zz = 2+(xx-16)/65
when xx < 256 then
zz = 3+(xx-81)/175
when xx < 625 then
zz = 4+(xx-256)/360
when xx < 1296 then
zz = 5+(xx-625)/671
when xx < 2401 then
zz = 6+(xx-1296)/1105
when xx < 4096 then
zz = 7+(xx-2401)/1695
when xx < 6561 then
zz = 8+(xx-4096)/2465
otherwise
zz = 9+(xx-6561)/3439
end
-- Dynamic precision
do n = 1 while dd > 4
d.n = dd; dd = Trunc(0.48*dd)
end
d.n = 4
-- Halley
do i = n to 1 by -1
numeric digits d.i
z4 = zz*zz*zz*zz; z5 = z4*zz; zz = zz-(2*z5-2*zz*xx)/(3*xx+5*z4)
end
-- Undo reduce
return zz'E'e
Rad:
-- Degrees -> Radians
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
numeric digits Digits()+1
return xx*Pi()/180
Radical:
-- Radical
-- Product of unique factors
procedure expose ufac.
arg xx
-- Get unique factors
n = Ufactors(xx)
-- Calculate product
zz = 1
do i = 1 to n
zz = zz*ufac.factor.i
end
return zz
Rand:
-- Random number
procedure expose glob.
arg xx,yy
if xx <> '' then do
if yy = '' then do
yy = xx; xx = 1
end
if \ Datatype(xx,'w') then say argument
if \ Datatype(yy,'w') then say argument
if xx >= yy then say argument
end
-- Get and save first seed from system Date and Time
numeric digits 30
if glob.rand = '' then do
d = Date('s')
a = Substr(d,3,2); b = Substr(d,5,2); c = Substr(d,7,2)
t = Time('l')
d = Substr(t,1,2); e = Substr(t,4,2); f = Substr(t,7,2)
g = Substr(t,10,1); h = Substr(t,11,1); i = Substr(t,12,1)
glob.rand = h||f||d||b||a||c||e||g||i
end
-- Calculate next random number cf HP calculators
glob.rand = Right((glob.rand*2851130928467),15)
-- Uniform deviate (0,1)
zz = '0.'Left(glob.rand,12)
-- or specified range (x,y)
if xx <> '' then
zz = Floor(zz*(yy+1-xx)+xx)
numeric digits 12
return zz+0
Real:
-- Are all arguments in R?
procedure
do i = 1 to Arg()
if \ Datatype(Arg(i),'n') then say argument
end
return 1
Repunit:
-- Repeated 1's
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Formula
return Copies('1',xx)
Rfact:
-- Rising factorial
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
-- Fast
if xx = 1 then
return 1
if yy = 0 then
return 1
if glob.Rfact.xx.yy <> '' then
return glob.Rfact.xx.yy
-- Cf definition
zz = 1
do k = 0 to yy-1
zz = zz*(xx+k)
end
glob.Rfact.xx.yy = zz
return zz
Round:
-- Round to specified decimals
procedure
arg xx,yy
if yy = '' then
yy = 0
if \ Datatype(xx,'n') then say argument
if \ Datatype(yy,'w') then say argument
-- Formula
if yy < 0 then do
yy = Abs(yy); a = 1'E'yy
return Format(xx/a,,0,0)*a
end
else
return Format(xx,,yy,0)
Sci:
-- Convert to scientific format
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
return Mant(xx)'E'Xpon(xx)
Sec:
-- Secant
procedure expose glob.
arg xx
if xx = 0 then
return 1
-- Formula
numeric digits Digits()+1
zz = Cos(xx)
if zz = 0 then say argument
return 1/zz
Sech:
-- Secant hyperbolic
procedure expose glob.
arg xx
-- Formula
numeric digits Digits()+1
return 1/Cosh(xx)
Sigma:
-- Sigma
-- Sum of divisors
procedure expose divi.
arg xx
-- Formula
return Divisor(xx,1)
Sin:
-- Sine
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 0
-- Reflect
if xx < 0 then
return -Sin(-xx)
-- Reduce argument to 0...2pi
numeric digits Digits()+Max(Xpon(xx),0)+1
xx = xx//tau()
-- Reduce argument to 0...pi/2
select
when xx < Halfpi() then
s = 1
when xx < Pi() then do
s = 1; xx = Pi()-xx
end
when xx < Pi()+Halfpi() then do
s = -1; xx = xx-Pi()
end
otherwise do
s = -1; xx = tau()-xx
end
end
-- Taylor
t = xx; zz = xx; xx = xx*xx
do n = 2 by 2 until zz = v
v = zz
t = -t*xx/(n*(n+1)); zz = zz+t
end
return s*zz
Sinh:
-- Sine hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Reflect
if xx < 0 then
return -Sinh(-xx)
-- Fast
if xx = 0 then
return 0
-- Taylor
numeric digits Digits()+1
if \ Integer(xx) & xx < 5 then do
t = xx; zz = xx; xx = xx*xx
do n = 3 by 2 until zz = v
v = zz
t = (t*xx)/(n*n-n); zz = zz+t
end
return zz
end
-- Formula
zz = Exp(xx)
return (zz*zz-1)/(2*zz)
Sqrt:
-- Square root
-- x^(1/2)
procedure
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 0
if xx = 1 then
return 1
dd = Digits()+2
numeric digits dd
if xx = 2 then
return Sqrt2()
if xx = 3 then
return Sqrt3()
if xx = 5 then
return Sqrt5()
-- Reduce argument to 1...100
e = Xpon(xx); e = (e-(e<0))%2; f = -2*e; xx = Std(xx)'E'f
-- First guess 2 digits accurate
numeric digits 2
select
when xx < 4 then
zz = 1+(xx-1)/3
when xx < 9 then
zz = 2+(xx-4)/5
when xx < 16 then
zz = 3+(xx-9)/7
when xx < 25 then
zz = 4+(xx-16)/9
when xx < 36 then
zz = 5+(xx-25)/11
when xx < 49 then
zz = 6+(xx-36)/13
when xx < 64 then
zz = 7+(xx-49)/15
when xx < 81 then
zz = 8+(xx-64)/17
otherwise
zz = 9+(xx-81)/19
end
-- Dynamic precision
do n = 1 while dd > 4
d.n = dd; dd = Trunc(0.36*dd)
end
d.n = 4
-- Halley
do i = n to 1 by -1
numeric digits d.i
z2 = zz*zz; zz = zz-2*zz*(z2-xx)/(3*z2+xx)
end
-- Undo reduce
return zz'E'e
Std:
-- Convert to standard format
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
return Format(xx,,,0)
Tan:
-- Tangent
procedure expose glob.
arg xx
-- Fast
if xx = 0 then
return 0
-- Formula
numeric digits Digits()+1
zz = Cos(xx)
if zz = 0 then say argument
return Sin(xx)/zz
Tanh:
-- Tangent hyperbolic
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
-- Fast
if xx = 0 then
return 0
-- Reflect
if xx < 0 then
return -Tanh(-xx)
-- Formula
numeric digits Digits()+2
if xx < 1 then
numeric digits Digits()-Xpon(xx)
zz = Exp(2*xx)
return (zz-1)/(zz+1)
Totient:
-- Euler's totient function
-- Phi(x) number of relative primes to x
procedure expose fact.
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 1 then say argument
-- Fast
if xx < 3 then
return 1
if xx < 5 then
return 2
-- Multiplicative property using Factors
f = Factors(xx)+1; fact.factor.f = 0
zz = 1; v = fact.factor.1; n = 1
do i = 2 to f
a = fact.factor.i
if a = v then
n = n+1
else do
zz = zz*v**(n-1)*(v-1)
v = a; n = 1
end
end
return zz
Trigamma:
-- Trigamma
-- Derivative of Digamma
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx <= 0 & Integer(xx) then say argument
-- Fast sum
numeric digits Digits()+2
if xx > 0 & xx < 50 then do
if Integer(xx) then do
s = 0
do i = 1 to xx-1
s = s+1/(i*i)
end
return Pi()*Pi()/6-s
end
if Half(xx) then do
s = 0
do i = 1 to xx
a = 2*i-1
s = s+1/(a*a)
end
return Halfpi()*Pi()-4*s
end
end
-- Reduce argument towards a higher value
d = Digits(); a = 0
do while xx < d
a = a+1/(xx*xx); xx = xx+1
end
-- Series
zz = 0; v = 0
do n = 2 by 2
zz = zz+Bernoulli(n)/xx**(n+1)
if zz = v then
leave
v = zz
end
-- Undo reduce
return zz+1/xx+1/(2*xx*xx)+a
Whole:
-- Are all arguments in Z?
procedure
do i = 1 to Arg()
if \ Datatype(Arg(i),'w') then say argument
end
return 1
Xpon:
-- Exponent
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Formula
numeric digits 9
if xx = 0 then
return 0
else
return Right(xx*1E+99999999,9)-99999999
Zero:
-- Are all arguments <> 0?
procedure
do i = 1 to Arg()
if Arg(i) = 0 then say argument
end
return
Zeta:
-- Zeta
procedure expose glob.
arg xx
if \ Datatype(xx,'n') then say argument
if xx = 1 then say argument
-- Fast
dd = Digits()
numeric digits dd+3
if xx = 0 then
return -0.5
if xx = 3 then
return Apery()
if xx < 0 & Even(xx) then
return 0
-- Formula
if xx < 0 & Odd(xx) then
return ((-1)**-xx*Bernoulli(1-xx))/(1-xx)
if xx > 0 & Even(xx) then
return ((-1)**(xx/2+1)*Bernoulli(xx)*(Tau())**xx)/(2*Fact(xx))
-- Reflection
if xx < 0 then
return -2*Factorial(-xx)*Sin(Pi()*-xx/2)*Zeta(1-xx)/Power(Tau(),1-xx)
-- Selected odd integers cf Plouffe
if xx > 4 & xx < 22 & Odd(xx) then do
select
when xx = 5 then do
a = 1470; b = 5; c = 3024; d = 84
end
when xx = 7 then do
a = 56700; b = 19; c = 113400; d = 0
end
when xx = 9 then do
a = 18523890; b = 625; c = 37122624; d = 74844
end
when xx = 11 then do
a = 425675250; b = 1453; c = 851350500; d = 0
end
when xx = 13 then do
a = 257432175; b = 89; c = 514926720; d = 62370
end
when xx = 15 then do
a = 390769879500; b = 13687; c = 781539759000; d = 0
end
when xx = 17 then do
a = 1904417007743250; b = 6758333; c = 3808863131673600; d = 29116187100
end
when xx = 19 then do
a = 21438612514068750; b = 7708537; c = 42877225028137500; d = 0
end
when xx = 21 then do
a = 1881063815762259253125; b = 68529640373; c = 3762129424572110592000; d = 1793047592085750
end
end
s1 = 0; s2 = 0; v1 = 0; v2 = 0
do n = 1 to 1000
m = n**xx; e = Exp(Tau()*n); q = m*e
s1 = s1+1/(q-m); s2 = s2+1/(q+m)
if s1 = v1 & s2 = v2 then
leave
v1 = s1; v2 = s2
end
zz = (b*Pi()**xx-c*s1-d*s2)/a
return zz
end
-- Standard for all other xx
-- Max iterations estimate
select
when dd < 15 then
n = dd*2
when dd < 89 then
n = dd*3
otherwise
n = dd*4
end
-- Reduce argument by factor
a = (-1)/(2**n*(1-Power(2,1-xx)))
-- Series
s = 0; zz = 0
do j = 0 to 2*n-1
b = (-1)**j; c = Power(j+1,xx)
if j < n then
s = 0
else
s = s+Comb(n,j-n)
zz = zz+(s-2**n)*b/c
end
-- Undo reduce
return a*zz
-- Module Helper.inc - 11 Apr 2025
-- Supporting procedures
Sort:
-- Sort 1 or more key stems, syncing 0 or more data stems
arg keys!,data!
kn! = Words(keys!)
do x! = 1 to kn!
key!.x! = Word(keys!,x!)
end
dn! = Words(data!)
do x! = 1 to dn!
data!.x! = Word(data!,x!)
end
n! = Value(key!.1||0); s! = 1; sl!.1 = 1; sr!.1 = n!
do until s! = 0
l! = sl!.s!; r! = sr!.s!; s! = s!-1
do until l! >= r!
-- Insertion sort
if r!-l! < 20 then do
do i! = l!+1 to r!
do x! = 1 to kn!
k!.x! = Value(key!.x!||i!)
end
do x! = 1 to dn!
d!.x! = Value(data!.x!||i!)
end
do j!=i!-1 to l! by -1
do x! = 1 to kn!
a! = Value(key!.x!||j!)
if a! > k!.x! then
leave x!
if a! = k!.x! then
if x! < kn! then
iterate x!
leave j!
end
k! = j!+1
do x! = 1 to kn!
t! = Value(key!.x!||j!)
call Value key!.x!||k!,t!
end
do x! = 1 to dn!
t! = Value(data!.x!||j!)
call Value data!.x!||k!,t!
end
end
k! = j!+1
do x! = 1 to kn!
call Value key!.x!||k!,k!.x!
end
do x! = 1 to dn!
call Value data!.x!||k!,d!.x!
end
end
if s! = 0 then
leave
l! = sl!.s!; r! = sr!.s!; s! = s!-1
end
-- Quicksort
else do
m! = (l!+r!)%2
do x! = 1 to kn!
a! = Value(key!.x!||l!); b! = Value(key!.x!||m!)
if a! > b! then do
do y! = 1 to kn!
t! = Value(key!.y!||l!); u! = Value(key!.y!||m!)
call Value key!.y!||l!,u!; call Value key!.y!||m!,t!
end
do y! = 1 to dn!
t! = Value(data!.y!||l!); u! = Value(data!.y!||m!)
call Value data!.y!||l!,u!; call Value data!.y!||m!,t!
end
leave
end
if a! < b! then
leave
end
do x! = 1 to kn!
a! = Value(key!.x!||l!); b! = Value(key!.x!||r!)
if a! > b! then do
do y! = 1 to kn!
t! = Value(key!.y!||l!); u! = Value(key!.y!||r!)
call Value key!.y!||l!,u!; call Value key!.y!||r!,t!
end
do y! = 1 to dn!
t! = Value(data!.y!||l!); u! = Value(data!.y!||r!)
call Value data!.y!||l!,u!; call Value data!.y!||r!,t!
end
leave
end
if a! < b! then
leave
end
do x! = 1 to kn!
a! = Value(key!.x!||m!); b! = Value(key!.x!||r!)
if a! > b! then do
do y! = 1 to kn!
t! = Value(key!.y!||m!); u! = Value(key!.y!||r!)
call Value key!.y!||m!,u!; call Value key!.y!||r!,t!
end
do y! = 1 to dn!
t! = Value(data!.y!||m!); u! = Value(data!.y!||r!)
call Value data!.y!||m!,u!; call Value data!.y!||r!,t!
end
leave
end
if a! < b! then
leave
end
i! = l!; j! = r!
do x! = 1 to kn!
p!.x! = Value(key!.x!||m!)
end
do until i! > j!
do i! = i!
do x! = 1 to kn!
a! = Value(key!.x!||i!)
if a! < p!.x! then
leave x!
if a! = p!.x! then
if x! < kn! then
iterate x!
leave i!
end
end
do j! = j! by -1
do x! = 1 to kn!
a! = Value(key!.x!||j!)
if a! > p!.x! then
leave x!
if a! = p!.x! then
if x! < kn! then
iterate x!
leave j!
end
end
if i! <= j! then do
do x! = 1 to kn!
t! = Value(key!.x!||i!); u! = Value(key!.x!||j!)
call Value key!.x!||i!,u!; call Value key!.x!||j!,t!
end
do x! = 1 to dn!
t! = Value(data!.x!||i!); u! = Value(data!.x!||j!)
call Value data!.x!||i!,u!; call Value data!.x!||j!,t!
end
i! = i!+1; j! = j!-1
end
end
if j!-l! < r!-i! then do
if i! < r! then do
s! = s!+1; sl!.s! = i!; sr!.s! = r!
end
r! = j!
end
else do
if l! < j! then do
s! = s!+1; sl!.s! = l!; sr!.s! = j!
end
l! = i!
end
end
end
end
return
-- Module Numbers.inc - 11 Apr 2025
-- Useful numbers and functions in number theory
Abundant:
-- Is a number abundant?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 101 then do
a = '12 18 20 24 30 36 40 42 48 54 56 60 66 70 72 78 80 84 88 90 96 100'
if Wordpos(xx,a) = 0 then
return 0
else
return 1
end
-- Cf definition
if Sigma(xx) > 2*xx then
return 1
else
return 0
Alcuin:
-- Alcuin number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 4 then
return 0
-- Formula
if Even(xx) then
return Round(xx*xx/48)
else
return Round((xx+3)**2/48)
Almostprime:
-- Is a number almost k-prime?
procedure expose glob. fact.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 1 then say argument
-- Cf definition
if yy = 1 then
return Prime(xx)
if yy = 2 then
return Semiprime(xx)
if Factors(xx) = yy then
return 1
else
return 0
Arithmetic:
-- Is a number arithmetic?
procedure expose divi.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 101 then do
a = '1 3 5 6 7 11 13 14 15 17 19 20 21 22 23 27 29 30 31 33 ',
|| '35 37 38 39 41 42 43 44 45 46 47 49 51 53 54 55 56 57 59 60 61 62 65 66 ',
|| '67 68 69 70 71 73 77 78 79 83 85 86 87 89 91 92 93 94 95 96 97 99 '
if Wordpos(xx,a) = 0 then
return 0
else
return 1
end
-- Cf definition
s = Sigma(xx)
if Integer(s/divi.0) then
return 1
else
return 0
Bell:
-- Bell number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
if xx = 1 then
return 1
-- Recurring
numeric digits Digits()+2
if xx < 36 then do
zz = 0
do k = 0 to xx
zz = zz+Stirling2(xx,k)
end
end
else do
zz. = 0; zz.0 = 1; zz.1 = 1
do n = 2 to xx
s = 0
do k = 0 to n-1
s = s+Comb(n-1,k)*zz.k
end
zz.n = s
end
zz = zz.xx
end
return zz
Beraha:
-- Beraha number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 1 then
return 4
if xx = 2 then
return 0
if xx = 3 then
return 1
if xx = 4 then
return 2
if xx = 6 then
return 3
-- Formula
numeric digits Digits()+1
if xx = 8 then
return 2+Sqrt(2)
if xx = 10 then
return 0.5*(5+Sqrt(5))
return 2+2*Cos(Tau()/xx)
Bernoulli:
-- Bernoulli number
procedure expose glob. work.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
if xx = 1 then
return -0.5
if xx//2 = 1 then
return 0
p = Digits()
if glob.bernoulli.xx.p <> '' then
return glob.bernoulli.xx.p
-- Recurring method Akiyama-Tanigawa
numeric digits Max(p,xx+xx)
do m = 0 by 1 to xx
work.m = 1/(m+1)
do j = m by -1 to 1
k = j-1; work.k = j*(work.k-work.j)
end
end
-- Result
zz = work.0
glob.bernoulli.xx.p = zz
return zz
Cake:
-- Cake number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Formula
return (xx*xx*xx+5*xx+6)/6
Catalann:
-- Catalan number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Recurring
numeric digits Digits()+Xpon(xx)
zz = 1
do n = 1 to xx
zz = zz*(4*n-2)/(n+1)
end
return zz
Caterer:
-- Lazy caterer's number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Formula
return 0.5*(xx*xx+xx+2)
Central:
-- Central binomial coefficient number
procedure expose glob.
arg xx
-- Formula
return Comb(2*xx,xx)
Chowla:
-- Chowla number
procedure
arg xx
-- Formula
return Aliquot(xx)-1
Composite:
-- Is a number composite?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx = 1 then
return 0
if xx < 101 then do
c = '4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 33 34 35 36 38 39 ',
|| '40 42 44 45 46 48 49 50 51 52 54 55 56 57 58 60 62 63 64 65 66 68 69 70 72 ',
|| '74 75 76 77 78 80 81 82 84 85 86 87 88 90 91 92 93 94 95 96 98 99 100 '
if Wordpos(xx,c) = 0 then
return 0
else
return 1
end
if xx//2 = 0 then
return 1
if xx//3 = 0 then
return 1
if Right(xx,1) = 5 then
return 1
-- Trial division
do n = 5 by 6 while n*n <= xx
if xx//n = 0 then
return 1
if xx//(n+2) = 0 then
return 1
end
return 0
Conwayn:
-- Conway look-and-say number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Series cf definition
zz = '1'
do n = 1 to xx
ns = ''; c = 1; m = Length(zz)
do j = 1 to m
if j = m | Substr(zz,j,1) <> Substr(zz,j+1,1) then do
ns = ns||c||Substr(zz,j,1); c = 1
end
else
c = c+1
end
zz = ns
end
return zz
Coprime:
-- Are 2 numbers coprime?
procedure
arg xx,yy
-- Formula
return (Gcd(xx,yy) = 1)
Cullen:
-- Cullen number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Formula
return xx*2**xx+1
Deficient:
-- Is a number deficient?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 101 then do
d = '1 2 3 4 5 7 8 9 10 11 13 14 15 16 17 19 21 22 23 25 26 27 29 31 32 33 34 35 ',
|| '37 38 39 41 43 44 45 46 47 49 50 51 52 53 55 57 58 59 61 62 63 64 65 67 68 ',
|| '69 71 73 74 75 76 77 79 81 82 83 85 86 87 89 91 92 93 94 95 97 98 99 '
if Wordpos(xx,d) = 0 then
return 0
else
return 1
end
-- Cf definition
if Sigma(xx) < 2*xx then
return 1
else
return 0
Derangement:
-- Derangement number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Recurring cf definition
numeric digits Digits()+2
a = 1; b = 0; zz = 0
do n = 2 to xx
zz = (a+b)*(n-1)
a = b; b = zz
end
return zz
Favard:
-- Favard number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
-- Formula
numeric digits Digits()+1
if xx = 1 then
return Halfpi()
if xx = 2 then
return Pi()**2*0.125
if xx = 3 then
return Pi()**3/24
if xx = 4 then
return 5*Pi()**4/384
if xx = 5 then
return Pi()**5/240
if xx = 6 then
return 61*Pi()**6/46080
if xx = 7 then
return 17*Pi()**7/40320
if xx = 8 then
return 277*Pi()**8/2064384
if xx = 9 then
return 31*Pi()**9/725760
if xx = 10 then
return 50521*Pi()**10/3715891200
if xx = 11 then
return 691*Pi()**11/159667200
if xx = 12 then
return 540553*Pi()**12/392398110720
-- Recurring
numeric digits Digits()+1
zz.0 = 1; zz.1 = Pi()*0.5
do j = 2 to xx
s = 0
do k = 1 to Floor(j*0.5)
a = 2*k; b = a-1; c = j-a
s = s+zz.b*zz.c
end
zz.j = s*Pi()/(2*j)
end
return zz.xx
Fermat:
-- Fermat number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Formula
numeric digits Digits()+1
zz = 2**(2**xx)
return zz
Ffactor:
-- First prime factor by trial division
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 2 then say argument
-- Fast
if Even(xx) then
return 2
if Right(xx,1) = 5 then
return 5
if Prime(xx) then
return xx
-- Check low factors
n = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i)
if xx//p = 0 then
return p
end
-- Check higher factors
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then
if xx//j = 0 then
return j
if p = 3 then
iterate
k = j+2
if xx//k = 0 then
return k
end
if xx > 1 then
return xx
return 0
Fibonacci:
-- Fibonacci number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return xx
-- Recurring
numeric digits Digits()+2
a = 0; b = 1
do n = 2 to xx
zz = a+b; a = b; b = zz
end
return zz
Goldbach:
-- Goldbach number
procedure expose prim. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 4 then say argument
-- Scan primes
p = Primes(xx); zz = 0
do i = 2 to xx%2
if flag.i then do
j = xx-i
if flag.j then do
zz = zz+1
end
end
end
return zz
Golombn:
-- Golomb-Silverman number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx = 1 then
return 1
-- Recurring
numeric digits Digits()+2
zz. = 0; zz.1 = 1
do n = 2 to xx
a = n-1; b = zz.a; c = zz.b; d = n-c
zz.n = 1+zz.d
end
return zz.xx
Gould:
-- Gould number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
-- Recurring
numeric digits Digits()+2
zz. = 0; zz.0 = 1
do n = 1 to xx
if Even(n) then do
a = n%2; zz = zz.a
end
else do
a = n%2; zz = 2*zz.a
end
zz.n = zz
end
return zz.xx
Gregory1:
-- Gregory number 1
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx = 1 then
return 0.5
-- Recurring
numeric digits Digits()+2
zz.1 = 0.5
do n = 2 to xx
s = 0
do k = 1 to n-1
s = s+zz.k/(n+1-k)
end
zz.n = -s+1/(n+1)
end
return (-1)**(xx-1)*zz.xx
Gregory2:
-- Gregory number 2
procedure
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 1 then say argument
-- Fast
if xx = 1 then
return 0.25*Pi()
-- Formula
return Arctan(1/xx)
Harmonic:
-- Harmonic number
procedure
arg xx,yy
if yy = '' then
yy = 1
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
-- Fast
if xx = 0 then
return 0
if xx < 2 then
return 1
-- From definition
numeric digits Digits()+2
zz = 0
if yy = 1 then do
do n = 1 to xx
zz = zz+1/n
end
end
else do
do n = 1 to xx
zz = zz+1/n**yy
end
end
return zz
Highcomposite:
-- Is a number highly composite?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- ... to do
return 0
Home:
-- Home prime number
procedure expose fact. glob. work.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Run thru chain
numeric digits 25
do while Factors(xx) > 1
xx = ''
do i = 1 to fact.0
xx = xx||fact.i
end
end
return xx
Jacobsthal:
-- Jacobsthal number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 3 then
return 1
-- Formula
numeric digits Digits()+2
zz = (2**xx-(-1)**xx)/3
return zz
Kaprekar:
-- Is a number Kaprekar?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 101 then do
k = '1 9 45 55 99'
if Wordpos(xx,k) = 0 then
return 0
else
return 1
end
-- Casting out nines
p = Max(Digits(),2*Length(xx))
numeric digits p
a = xx//9
if a > 2 then
return 0
s = xx*xx
if a = s//9 then do
do i = 1 to Length(s)%2
parse var s l +(i) r
if xx = l+r then
return 1
end
end
return 0
Lah:
-- Lah number
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 1 then say argument
if xx < yy then say argument
-- Fast
if xx = yy then
return 1
-- Formula
s = (-1)**xx
if yy = 1 then
return s*Fact(xx)
else
return s*Comb(xx-1,yy-1)*Fact(xx)/Fact(yy)
Lambdan:
-- Lambda number
procedure expose efac.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 3 then
return 1
if xx < 5 then
return 2
-- Using Efactors and Lcm
f = Efactors(xx)
zz = 1
do i = 1 to f
p = efac.factor.i
e = efac.exponent.i
if p = 2 & e > 2 then
zz = Lcm(zz,2**(e-2))
else
zz = Lcm(zz,(p-1)*p**(e-1))
end
return zz
Lebesguen:
-- Lebesgue number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
if xx = 1 then
return 1/3+2*Sqrt(3)/Pi()
if xx = 2 then
return 0.2+Sqrt(25-2*Sqrt(5))/Pi()
if xx = 3 then
return 1/7+(22*Sin(Pi()/7)-2*Cos(Pi()/14)+10*Cos(3*Pi()/14))/(3*Pi())
if xx = 4 then
return 13/(2*Sqrt(3)*Pi())+1/9+(7*Sin(Tau()/9)-5*Sin(Pi()/9)-Cos(Pi()/18))/Pi()
-- Series expansion
numeric digits Digits()+2
zz = 0
do n = 1 to xx
zz = zz+Tan(n*Pi()/(2*xx+1))/n
end
zz = 1/(2*xx+1)+2*zz/Pi()
return zz
Leonardo:
-- Leonardo number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 3 then
return 1
-- From definition
a = 1; b = 1
do n = 3 to xx
zz = a+b+1
a = b; b = zz
end
return zz
Lucas:
-- Lucas number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 1 then
return 2
if xx = 2 then
return 1
-- Recurring
numeric digits Digits()+2
a = 2; b = 1; zz = 0
do n = 2 to xx
zz = a+b
a = b; b = zz
end
return zz
Magic:
-- Magic number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Formula
return 0.5*(xx*xx*xx+xx)/1
Metallic:
-- Metallic ratio number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
-- Formula
if xx = 1 then
return Golden()
if xx = 2 then
return Silver()
if xx = 3 then
return Bronze()
return 0.5*(xx+Sqrt(xx*xx+4))
Moebius:
-- Moebius number
procedure expose glob. fact. ufac.
arg xx
-- Using number of (unique) prime factors
call Factors(xx)
call Ufactors(xx)
if fact.0 = ufac.0 then
if Even(fact.0) then
return 1
else
return -1
else
return 0
Motzkin:
-- Motzkin number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 3 then
return 1
-- Recurring
numeric digits Digits()+2
a = 1; b = 1; zz = 0
do n = 2 to xx
zz = ((3*n-3)*a+(2*n+1)*b)/(n+2)
a = b; b = zz
end
return zz
Narayana:
-- Narayana cows number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 3 then
return 1
-- Recurring
numeric digits Digits()+2
a = 1; b = 1; c = 1; zz = 0
do n = 3 to xx
zz = a+c
a = b; b = c; c = zz
end
return zz
Nextprime:
-- Next prime number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 2 then say argument
-- Find next prime
xx = xx+1+odd(xx)
do zz = xx by 2 until prime(zz)
end
return zz
Padovan:
-- Padovan number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 3 then
return 1
-- Recurring
numeric digits Digits()+2
a = 1; b = 1; c = 1; zz = 0
do n = 3 to xx
zz = a+b
a = b; b = c; c = zz
end
return zz
Partition:
-- Partition number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
-- Recurring
numeric digits Digits()+2
yy. = 0; yy.0 = 1
do n = 1 to xx
s = 0
do k = 1 to n
a = n-k*(3*k-1)/2; b = n-0.5*k*(3*k+1)
s = s+(-1)**(k+1)*(yy.a+yy.b)
end
yy.n = s
end
return yy.xx
Pell1:
-- Pell number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return xx
-- Recurring
numeric digits Digits()+2
a = 0; b = 1; zz = 0
do n = 2 to xx
zz = a+2*b
a = b; b = zz
end
return zz
Pell2:
-- Pell-Lucas companion number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return 2
-- Recurring
numeric digits Digits()+2
a = 2; b = 2; zz = 0
do n = 2 to xx
zz = a+2*b
a = b; b = zz
end
return zz
Perfect:
-- Is a number perfect?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
p = '6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128 ',
|| '2658455991569831744654692615953842176 ',
|| '191561942608236107294793378084303638130997321548169216 ',
|| '13164036458569648337239753460458722910223472318386943117783728128 ',
|| '14474011154664524427946373126085988481573677491474835889066354349131199152128 '
if Wordpos(xx,p) = 0 then
return 0
else
return 1
Perrin:
-- Perrin number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 3
if xx = 1 then
return 0
if xx = 2 then
return 2
-- Recurring
numeric digits Digits()+2
a = 3; b = 0; c = 2; zz = 0
do n = 3 to xx
zz = a+b
a = b; b = c; c = zz
end
return zz
Persistence:
-- Additive persistence
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 0
-- Cf definition
do zz = 1 until xx < 10
xx = Digitsum(xx)
end
return zz
Pisanoprime:
-- Pisano period
procedure
arg xx,yy
if \ Datatype(xx,'w') then say argument
if \ prime(xx) then say argument
if \ Datatype(yy,'w') then say argument
if yy < 1 then say argument
-- Pisano period
zz = 0
return zz
Pollardrho:
-- Find a prime factor by Pollard-Rho
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 2 then say argument
-- Fast
-- Fast
pr = '2 3 5 7 11 13 17 19 23 29'
do i = 1 to Words(pr)
p = Word(pr,i)
if xx//p = 0 then
return p
end
if Prime(xx) then
return xx
-- Square?
s = Isqrt(xx)
if s*s = xx then
return Pollardrho(s)
-- Pollard Rho
a = Random(1,9); b = a; c = Random(1,9); zz = 1; n = 0
-- a = 2; b = a; c = 1; zz = 1; n = 0
m = Sqrt(xx)
do while zz = 1 & n < m
a = (a*a+c)//xx; b = (b*b+c)//xx; b = (b*b+c)//xx
zz = Gcd(a-b,xx)
if zz = xx then
return 0
n = n+1
end
if Prime(zz) then
return zz
else
return Pollardrho(zz)
Practical:
-- Is a number practical?
procedure expose fact.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 101 then do
p = '1 2 4 6 8 12 16 18 20 24 28 30 32 36 40 42 ',
|| '48 54 56 60 64 66 72 78 80 84 88 90 96 100 '
if Wordpos(xx,p) = 0 then
return 0
else
return 1
end
if Odd(xx) then
return 0
-- Method Srinivasan-Stewart-Sierpinski
n = Factors(xx)
f = 2; s = 2
do i = 2 to n
if fact.i > s+1
then return 0
f = f*fact.i; s = Sigma(f)
end
return 1
Prevprime:
-- Previous prime number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 3 then say argument
-- Fast
if xx = 3 then
return 2
-- Find previous prime
xx = xx-1-odd(xx)
do zz = xx by -2 until prime(zz)
end
return zz
Prime:
-- Is a number prime?
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return 0
xx = xx+0; lp = '2 3 5 7 11 13 17 19 23 29 31 37 41'
if xx < 42 then do
if Wordpos(xx,lp) = 0 then
return 0
else
return 1
end
if Right(xx,1) = 0 then
return 0
if Right(xx,1) = 5 then
return 0
if xx//2 = 0 then
return 0
if xx//3 = 0 then
return 0
-- Miller-Rabin primality test
numeric digits 2*Length(xx)
d = xx-1; e = d
-- Reduce xx-1 by factors of 2
do s = -1 while d//2 = 0
d = d%2
end
-- Thresholds deterministic witnesses
th = '2047 1373653 25326001 3215031751 2152302898747 3474749660383 341550071728321 0 ',
|| '3825123056546413051 0 0 318665857834031151167461 3317044064679887385961981 '
w = Words(th)
-- Up to 13 deterministic trials
if xx < Word(th,w) then do
do k = 1 to w
if xx < Word(th,k) then
leave
end
end
-- or 3 probabilistic trials
else do
lp = ''
do k = 1 to 3
r = Rand(2,e)/1; lp = lp||r||' '
end
k = k-1
end
-- Algorithm using generated witnesses
do i = 1 to k
a = Word(lp,i); zz = Powmod(a,d,xx)
if zz = 1 then
iterate
if zz = e then
iterate
do s
zz = (zz*zz)//xx
if zz = 1 then
return 0
if zz = e then
leave
end
if zz <> e then
return 0
end
return 1
Primell:
-- Is a Mersenne number prime?
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Prime(xx) then say argument
-- Fast
if xx = 2 then
return 1
-- Lucas-Lehmer
m = 2**xx-1; a = 4
do xx-2
a = (a*a-2)//m
end
if a//m = 0 then
return 1
else
return 0
Pronic:
-- Pronic number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Formula
numeric digits Digits()+2
zz = xx*xx+xx
return zz
Recaman:
-- Recaman number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 0
-- Recurring cf definition
zz. = 0; zz.0 = 0; s. = 0; s.0 = 1
do n = 1 to xx
m = n-1; a = zz.m-n
if a > 0 & s.a = 0 then do
zz.n = a; s.a = 1
end
else do
b = zz.m+n
zz.n = b; s.b = 1
end
end
return zz.xx
Semiprime:
-- Is a number semi prime?
procedure expose prim. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
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 '
if xx < 4 then
return 0
if xx < 101 then do
if Wordpos(xx,s) = 0 then
return 0
else
return 1
end
-- Wheeled scan
do i = 2 for 2
if xx//i = 0 then
if Prime(xx%i) then
return 1
else
return 0
end
do i = 5 by 6 until i*i > xx
do j = i by 2 for 2
if xx//j = 0 then
if Prime(xx%j) then
return 1
else
return 0
end
end
return 0
Sorting:
-- Sorting number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 3
if xx = 1 then
return 0
if xx = 2 then
return 2
-- Recurring
zz = 0; i = xx-1; z = 1
do while i >= 0
zz = zz+i; i = i-z; z = z+z
end
return zz
Squarefree:
-- Is a number square free?
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 4 then
return 1
-- Recurring
o = Odd(xx)
do k = 2+o by 1+o to Isqrt(xx)
if xx//(k*k) = 0 then
return 0
end
return 1
Squfof:
-- Find a prime factor by square form factorization
procedure expose mult. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 2 then say argument
-- Fast
pr = '2 3 5 7 11 13 17 19 23 29'
do i = 1 to Words(pr)
p = Word(pr,i)
if xx//p = 0 then
return p
end
if Prime(xx) then
return xx
-- Square?
s = Isqrt(xx)
if s*s = xx then
return Squfof(s)
-- Multipliers
m = '1 3 5 7 11 15 21 33 35 55 77 105 165 231 385 1155'
do i = 1 to Words(m)
mult.i = Word(m,i)
end
mult.0 = i-1
-- Heuristic classic method
bb = Isqrt(s+s)*6
do m = 1 for mult.0
k = mult.m
d = xx*k; po = Isqrt(d); qq = d-po*po
if qq = 0 then
iterate m
p = po; pp = po; qp = 1
do i = 2 while i < bb
b = (po+p)%qq; p = b*qq-p; q = qq; qq = qp+b*(pp-p); r = Isqrt(qq)
if i//2 = 0 then
if r*r = qq then
leave
qp = q; pp = p
end
if i >= bb then
iterate
b = (po-p)%r; p = b*r+p; pp = p; qp = r; qq = (d-pp*pp)%qp
do until p = pp
pp = p; b = (po+p)%qq; q = qq; p = b*qq-p; qq = qp+b*(pp-p); qp = q
end
r = Gcd(xx,qp)
if r <> 1 then
if r <> xx then
if Prime(r) then
return r
end
return 0
Stirling1:
-- Stirling number 1
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
-- Fast
if yy = xx then
return 1
if yy > xx then
return 0
if yy = 0 then
return 0
-- Formula
s = (-1)**(xx-yy)
if yy = 1 then
return s*Fact(xx-1)
if yy = 2 then
return s*Round(Fact(xx-1)*Harmonic(xx-1))
if yy = 3 then
return s*Round(0.5*Fact(xx-1)*(Harmonic(xx-1)**2-Harmonic(xx-1,2)))
if yy = xx-3 then
return s*Comb(xx,2)*Comb(xx,4)
if yy = xx-2 then
return s*Round(0.25*(3*xx-1)*Comb(xx,3))
if yy = xx-1 then
return s*Comb(xx,2)
-- Definition
numeric digits Digits()+2
zz. = 0; zz.1 = 1
do i = 1 to xx-1
a = i
do j = i+1 to 1 by -1
b = j-1; zz.j = zz.b+a*zz.j
end
end
return s*zz.yy
Stirling2:
-- Stirling number 2
procedure expose glob.
arg xx,yy
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
if xx < yy then say argument
-- Fast
if xx = 0 then
return 0
if yy = 0 then
return 0
if xx = yy then
return 1
if yy = 1 then
return 1
-- Formula
if yy = 2 then
return 2**(xx-1)-1
if yy = 3 then
return (3**xx-3*2**xx+3)/6
if yy = xx-1 then
return Comb(xx,2)
-- Definition
numeric digits Digits()+2
zz. = 1; zz = xx-yy
do i = 2 to yy
do j = 1 to zz
a = j-1; zz.j = zz.j+i*zz.a
end
end
return zz.zz
Sylvester:
-- Sylvester number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
if xx > 33 then say argument
-- Definition
numeric digits Digits()+2
zz = 2
do n = 2 to xx
zz = zz*zz-zz+1
end
return zz
Taun:
-- Ramanujan Tau number
procedure expose divi.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Niebur series
numeric digits Digits()+4
s = 0; a = 18*xx*xx
do n = 1 to xx-1
b = n*n
s = s+b*(35*b-52*n*xx+a)*Sigma(n)*Sigma(xx-n)
end
return xx**4*Sigma(xx)-24*s
Tetranaccin:
-- Tetranacci number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 4 then
return 0
-- Recurring
a = 0; b = 0; c = 0; d = 1
do n = 3 to xx
zz = a+b+c; a = b; b = c; c = d; d = zz
end
return zz
Thue:
-- Thue-Morse number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Concatenate 2s complements
zz = 0
do i = 1 to xx
zz = zz||Translate(zz,10,01)
end
return zz
Tribonaccin:
-- Tribonacci number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 3 then
return 0
-- Recurring
a = 0; b = 0; c = 1
do n = 3 to xx
zz = a+b+c; a = b; b = c; c = zz
end
return zz
Ulam:
-- Ulam number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
if xx < 3 then
return xx
-- Series
numeric digits Digits()+2
zz.1 = 1; zz.2 = 2; d = 2; e.= 0; e.1 = 1; e.2 = 1; z = 3
do until d = xx
n = 0
do j = 1 to d
b = z-zz.j
if e.b then do
if zz.j <> b then do
n = n+1
if n > 2 then
leave j
end
end
end
if n = 2 then do
d = d+1; zz.d = z; e.z = 1
end
z = z+1
end
zz = zz.d
return zz
Wedderburn:
-- Wedderburn number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx < 2 then
return xx
-- Recurring
numeric digits Digits()+2
zz. = 0; zz.1 = 1
do n = 2 to xx
if Even(n) then do
a = n%2; s = 0.5*zz.a*(zz.a+1)
end
else do
a = n%2+1; s = 0
end
do i = 1 to a-1
j = n-i; s = s+zz.i*zz.j
end
zz.n = s
end
return zz.xx
Woodall:
-- Woodall number
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Formula
numeric digits Digits()+2
zz = xx*2**xx-1
return zz
ZigZag:
-- Zigzag alternating permutations number
procedure expose glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
if xx = 0 then
return 1
if xx = 1 then
return 1
-- Recurring
numeric digits Digits()+2
zz. = 0; zz.0 = 1; zz.1 = 1
do n = 2 to xx
s = 0
do k = 0 to n-1
a = Comb(n-1,k); b = n-k-1
s = s+a*zz.k*zz.b
end
zz.n = 0.5*s
end
zz = zz.xx
return zz
-- Module Polynomial.inc - 22 Mar 2025
-- Polynomial arithmetic, functions and tools
-- Coefficients are coded as string, highest order first
-- Examples 'a 0 c 0' = ax^3+cx, '1 0 0' = x^2, '2' = 2
Padd:
-- Addition
-- p1+p2
procedure
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
-- Add
wx = Words(xx); wy = Words(yy); wm = Max(wx,wy)
if wx < wy then
xx = Copies('0 ',wm-wx) xx
else
yy = Copies('0 ',wm-wy) yy
zz = ''
do i = 1 to wm
wx = Word(xx,i); wy = Word(yy,i)
if \ Datatype(wx,'n') then say argument i
if \ Datatype(wy,'n') then say argument i
zz = zz wx+wy
end
return zz
Parr2lst:
-- Array -> List
procedure expose poly.
-- Compose list
zz = ''
do i = 1 to poly.0
zz = zz poly.coef.i
end
return zz
Pcube:
-- Cube
-- p^3
procedure
arg xx
-- Formula
return Ppow(xx,3)
Pdif:
-- Differentiation
-- p'
procedure
arg xx
if xx = '' then say argument
-- Differentiate
wx = Words(xx); wm = wx-1
if wm = 0 then
return 0
zz = ''
do i = 1 to wm
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
zz = zz (wm-i+1)*w
end
return zz
Pdiv:
-- Long division
-- p1/p2
procedure
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
wx = Words(xx); wy = Words(yy)
do i = 1 to wx
if \ Datatype(Word(xx,i),'n') then say argument
end
do i = 1 to wy
if \ Datatype(Word(yy,i),'n') then say argument
end
-- Fast
if wx < wy then
return 0','xx
if wy = 1 then
return Pmul(1/yy,xx)','0
-- Divide
numeric digits Digits()+2
nx = xx; dx = yy Copies('0 ',wx-wy); qx = ''
do until Words(nx) < wy
q = Word(nx,1)/Word(dx,1); qx = qx q; z = ''
do i = 1 to Words(nx)
z = z Word(nx,i)-q*Word(dx,i)
end
nx = Subword(z,2); dx = Subword(dx,1,Words(dx)-1)
end
return qx','nx
Peval:
-- Evaluation
-- Value of p1 for y
procedure
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
if \ Datatype(yy,'n') then say argument
-- Evaluate
numeric digits Digits()+2
zz = 0
do i = 1 to Words(xx)
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
zz = zz*yy+w
end
return zz
Plst2form:
-- List -> Formula
procedure
arg xx
-- Fast
if xx = '' then
return 0
-- Compose formula
zz = ''; wm = Words(xx)
do i = 1 to wm
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
if w <> 0 then do
select
when w < 0 then
s = '-'
when i > 1 then
s = '+'
otherwise
s = ''
end
w = Abs(w); e = wm-i
if w = 1 & e > 0 then
w = ''
else
w = Std(w)
select
when e > 1 then
p = 'x^'e
when e = 1 then
p = 'x'
otherwise
p = ''
end
zz = zz||s||w||p
end
end
if zz = '' then
zz = 0
return zz
Pmul:
-- Multiplication
-- p1*p2
procedure expose poly. work.
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
-- Multiply
wx = Words(xx); wy = Words(yy)
do i = 1 to wx
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
work.coefx.i = w
end
do i = 1 to wy
w = Word(yy,i)
if \ Datatype(w,'n') then say argument i
work.coefy.i = w
end
numeric digits Digits()+2
zz = wx+wy-1
poly. = 0; poly.0 = zz
do i = 1 to wx
do j = 1 to wy
k = i+j-1; poly.coef.k = poly.coef.k+work.coefx.i*work.coefy.j
end
end
return Parr2lst()
Pneg:
-- Negation
-- -p
procedure
arg xx
if xx = '' then say argument
-- Negate
zz = ''
do i = 1 to Words(xx)
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
zz = zz (-w)
end
return Strip(zz)
Pnorm:
-- Normalization
procedure
arg xx
if xx = '' then say argument
-- Normalize
zz = ''
do i = 1 to Words(xx)
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
zz = zz w/1
end
return Strip(zz)
Ppow:
-- Integer power
-- p^n
procedure expose poly. comb. work.
arg xx,yy
if xx = '' then say argument
if \ Datatype(yy,'w') then say argument
if yy < 0 then say argument
-- Fast
poly. = 0
if yy = 0 then
return 1
-- Calculate power
numeric digits Digits()+2
wx = Words(xx); zz = wx*yy-yy+1
poly.0 = zz
select
when wx = 1 then do
-- Power of a real number
if \ Datatype(xx,'n') then say argument
poly.coef.1 = xx**yy
end
when wx = 2 then do
-- Newton's binomial
w1 = Word(xx,1); w2 = Word(xx,2)
if \ Datatype(w1,'n') then say argument
if \ Datatype(w2,'n') then say argument
call Combinations -yy
do i = 1 to zz
j = yy-i+1; k = i-1
poly.coef.i = comb.yy.k*w1**j*w2**k
end
end
otherwise do
-- Repeated multiplication
do i = 1 to wx
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
poly.coef1.i = w; poly.coef2.i = w
end
wy = wx
do i = 2 to yy
work. = 0
do j = 1 to wx
do k = 1 to wy
l = j+k-1; work.coef.l = work.coef.l+poly.coef1.j*poly.coef2.k
end
end
if i = yy then
leave i
wx = wx+wy-1
do j = 1 to wx
poly.coef1.j = work.coef.j
end
end
do i = 1 to zz
poly.coef.i = work.coef.i
end
end
end
return Parr2lst()
Pquartic:
-- Quartic
-- p^4
procedure
arg xx
-- Formula
return Psquare(Psquare(xx))
Prev:
-- Coefficient reversion
procedure
arg xx
if xx = '' then say argument
-- Revert
zz = ''
do i = 1 to Words(xx)
zz = Word(xx,i) zz
end
return zz
Psci:
-- Scientific notation
procedure
arg xx
if xx = '' then say argument
-- Convert
zz = ''
do i = 1 to Words(xx)
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
zz = zz Sci(w)
end
return Strip(zz)
Psquare:
-- Square
-- p^2
procedure
arg xx
-- Formula
return Pmul(xx,xx)
Pstd:
-- Standard notation
procedure
arg xx
if xx = '' then say argument
-- Convert
zz = ''
do i = 1 to Words(xx)
w = Word(xx,i)
if \ Datatype(w,'n') then say argument i
zz = zz Std(w)
end
return Strip(zz)
Psub:
-- Subtraction
-- p1-p2
procedure
arg xx,yy
if xx = '' then say argument
if yy = '' then say argument
-- Subtract
wx = Words(xx); wy = Words(yy); wm = Max(wx,wy)
if wx < wy then
xx = Copies('0 ',wm-wx) xx
else
yy = Copies('0 ',wm-wy) yy
zz = ''
do i = 1 to wm
wx = Word(xx,i); wy = Word(yy,i)
if \ Datatype(wx,'n') then say argument i
if \ Datatype(wy,'n') then say argument i
zz = zz wx-wy
end
return zz
-- Module Rational.inc - 22 Mar 2025
-- Rational arithmetic, functions and tools
-- Fractions are coded as string 'numerator denominator'
-- Examples '7 3' = 7/3; '13' = '13 1' = 13; '1 7' = 1/7
Rabs:
-- Absolute value
-- |r|
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Formula
return Rsimp((Abs(xn)) (Abs(xd)))
Radd:
-- Addition
-- r1+r2+r3+...
procedure
if Arg() = 0 then say argument
-- Add
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
if n = 1 then do
zn = xn; zd = xd
end
else do
zn = zn*xd+zd*xn; zd = zd*xd
end
end
return Rsimp((zn) (zd))
Rational:
-- Are all arguments in Q?
procedure
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
end
return 1
Rceil:
-- Ceiling
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Ceiling on floating value
return Ceil(xn/xd)
Rcontfrac:
-- Continued fraction
procedure
arg a0,an,b1,bn,c
if \ Datatype(a0,'w') then say argument
if \ Datatype(b1,'w') then say argument
if \ Datatype(c,'w') then say argument
if an = '' then
say argument
if bn = '' then
say argument
-- Calculate
zz = 0
do n = c by -1 to 1
interpret 'a =' an
if n = 1 then
b = b1
else
interpret 'b =' bn
zz = b/(a+zz)
end
return a0+zz
Rcube:
-- Cube
-- r^3
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Formula
return Rsimp((xn*xn*xn) (xd*xd*xd))
Rden:
-- Denominator
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
return (xd)
Rdiv:
-- Division
-- r1/r2/r3/...
procedure
if Arg() = 0 then say argument
-- Divide
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
if n = 1 then do
zn = xn; zd = xd
end
else do
zn = zn*xd; zd = zd*xn
end
end
return Rsimp((zn) (zd))
Req:
-- Equal?
procedure
arg xx,yy; xx = xx 1; yy = yy 1; parse var xx xn xd .; parse var yy yn yd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
if \ Datatype(yn,'w') then say argument
if \ Datatype(yd,'w') then say argument
if yd = 0 then say argument
-- r1 = r2?
return (xn*yd = xd*yn)
Rfloat:
-- Rational -> Float
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Divide
return xn/xd
Rfloor:
-- Floor
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Floor on floating value
return Floor(xn/xd)
Rge:
-- Greater than or equal?
procedure
arg xx,yy
-- r1 >= r2?
return \ Rlt(xx,yy)
Rgt:
-- Greater than?
procedure
arg xx,yy; xx = xx 1; yy = yy 1; parse var xx xn xd .; parse var yy yn yd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
if \ Datatype(yn,'w') then say argument
if \ Datatype(yd,'w') then say argument
if yd = 0 then say argument
-- r1 > r2?
return (xn*yd > xd*yn)
Rinv:
-- Inversion
-- 1/r
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if xn = 0 then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Reverse
return Rsimp((xd) (xn))
Rle:
-- Less than or equal?
procedure
arg xx,yy
-- r1 <= r2?
return \ Rgt(xx,yy)
Rlst2form:
-- List -> Formula
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Compose formula
return (xn)||'/'||(xd)
Rlt:
-- Less than
procedure
arg xx,yy; xx = xx 1; yy = yy 1; parse var xx xn xd .; parse var yy yn yd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
if \ Datatype(yn,'w') then say argument
if \ Datatype(yd,'w') then say argument
if yd = 0 then say argument
-- r1 < r2?
return (xn*yd < xd*yn)
Rmax:
-- Maximum
-- max(r1,r2,r3,...)
procedure
if Arg() = 0 then say argument
-- Find maximum
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
if n = 1 then do
zn = xn; zd = xd
end
else do
if xn*zd > xd*zn then do
zn = xn; zd = xd
end
end
end
return Rsimp((zn) (zd))
Rmin:
-- Minimum
-- min(r1,r2,r3,...)
procedure
if Arg() = 0 then say argument
-- Find minimum
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
if n = 1 then do
zn = xn; zd = xd
end
else do
if xn*zd < xd*zn then do
zn = xn; zd = xd
end
end
end
return Rsimp((zn) (zd))
Rmod:
-- Modulo division
-- r1-r2*floor(r1/r2)
procedure
arg xx,yy; xx = xx 1; yy = yy 1; parse var xx xn xd .; parse var yy yn yd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
if \ Datatype(yn,'w') then say argument
if yn = 0 then say argument
if \ Datatype(yd,'w') then say argument
if yd = 0 then say argument
-- Formula
return Rsub(xx,Rmul(yy,Floor(Rfloat(xx)/Rfloat(yy))))
Rmul:
-- Multiplication
-- r1*r2*r3*...
procedure
if Arg() = 0 then say argument
-- Multiply
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
if n = 1 then do
zn = xn; zd = xd
end
else do
zn = zn*xn; zd = zd*xd
end
end
return Rsimp((zn) (zd))
Rne:
-- Not equal
procedure
arg xx,yy
-- r1 <> r2?
return \ Req(xx,yy)
Rneg:
-- Negation
-- -r
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Negate
return Rsimp((-xn) (xd))
Rnegative:
-- Less than zero?
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- r1 < 0?
return (xn*xd < 0)
Rnum:
-- Numerator
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Numerator
return (xn)
Rperiod:
-- Length of repeating decimals
procedure
arg xx
zz = Rrepetend(xx)
if zz = 0 then
return 0
else
return Length(zz)
Rperiodinv:
-- Length of repeating decimals 1/x
procedure
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
-- Fast
if xx < 3 then
return 0
-- Recurring
r = 1
do xx
r = 10*r//xx
end
s = r
do zz = 1 until r = s
r = 10*r//xx
end
return zz
Rpositive:
-- Greater than zero?
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- r1 > 0?
return (xn*xd > 0)
Rpow:
-- Integer exponentiation
-- r^n, n in Z
procedure
arg xx,yy; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
if \ Datatype(yy,'w') then say argument
-- Formula
if yy < 0 then
return Rpow(Rinv(xx),-yy)
else
return Rsimp((xn**yy) (xd**yy))
Rquartic:
-- Quartic
-- r^4
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Formula
return Rsimp((xn*xn*xn*xn) (xd*xd*xd*xd))
Rquot:
-- Integer quotient
-- r%n
procedure
arg xx,yy; xx = xx 1; yy = yy 1; parse var xx xn xd .; parse var yy yn yd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
if \ Datatype(yn,'w') then say argument
if yn = 0 then say argument
if \ Datatype(yd,'w') then say argument
if yd = 0 then say argument
-- Formula
return Rfloat(xx)%Rfloat(yy)
Rrat:
-- Float -> Rational
procedure
arg xx
if \ Datatype(xx,'n') then say argument
-- Starting values
in = xx; d = Digits()-1; h = 1'E'd
zn = 0; zd = 1; tn = 1; td = 0
-- Up to precision
do while tn <= h & td <= h
-- Calculate numerator and denominator
x = Trunc(xx)
y = tn; tn = x*tn+zn; zn = y
y = td; td = x*td+zd; zd = y
-- Check precision
if x = xx | tn/td = in then do
if tn > h | td > h then
leave
-- Results
zn = tn; zd = td
leave
end
xx = 1/(xx-x)
end
return Rsimp((zn) (zd))
Rrepetend:
-- Repeating decimals
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Long division
r = xn//xd; r. = ''; p = 1; d = ''
do while r > 0
if r.r <> '' then
return Substr(d,r.r)
r.r = p; r = r*10
a = r%xd; d = d||a; r = r//xd
p = p+1
end
return 0
Rsimp:
-- Simplification
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Simplify
if xd < 0 then do
xn = -xn; xd = -xd
end
g = Gcd(xn,xd)
return (xn/g) (xd/g)
Rsquare:
-- Square
-- r^2
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- Formula
return Rsimp((xn*xn) (xd*xd))
Rsub:
-- Subtraction
-- r1-r2-r3-...
procedure
if Arg() = 0 then say argument
-- Subtract
do n = 1 to Arg()
xx = Arg(n) 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument n
if \ Datatype(xd,'w') then say argument n
if xd = 0 then say argument n
if n = 1 then do
zn = xn; zd = xd
end
else do
zn = zn*xd-zd*xn; zd = zd*xd
end
end
return Rsimp((zn) (zd))
Rzero:
-- Zero?
procedure
arg xx; xx = xx 1; parse var xx xn xd .
if \ Datatype(xn,'w') then say argument
if \ Datatype(xd,'w') then say argument
if xd = 0 then say argument
-- r1 = 0?
return (xn = 0)
-- Module Roots.inc - 22 Mar 2025
-- Solvers for real and complex roots of polynomial equations
-- Real coefficients are coded as string, highest first
-- Example ax^3+cx = 'a 0 c 0'
Ccbeq:
-- Complex roots of cubic equation
-- az^3+bz^2+cz+d = 0 with real coefs
procedure expose cbrt. glob.
arg xx; xx = xx 0 0 0; parse var xx a b c d .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
if \ Datatype(c,'n') then say argument
if \ Datatype(d,'n') then say argument
-- Method Cardano et al
numeric digits Digits()*2
-- Reduce to x^3+ax^2+bx+c
parse value b/a c/a d/a with a1 b1 c1
-- Reduce to x^3+ax+b
h1 = a1*a1; h2 = h1*a1
a2 = b1-h1/3; b2 = -h2/27+h2/9-a1*b1/3+c1
-- Discriminant
s = Round(0.25*b2*b2+a2*a2*a2/27,0.5*Digits())
-- Roots
cbrt. = 0
select
when s < 0 then do
h1 = Abs(a2)/3
f = Arccos(-b2/(2*Sqrt(h1*h1*h1)))
h1 = 2*Sqrt(h1); h2 = f/3; h3 = Pi()/3; h4 = a1/3
cbrt.root.1 = h1*Cos(h2)-h4
cbrt.root.2 = -h1*Cos(h2-h3)-h4
cbrt.root.3 = -h1*Cos(h2+h3)-h4
end
when s = 0 then do
h1 = Cbrt(-0.5*b2); h2 = a1/3
cbrt.root.1 = 2*h1-h2
cbrt.root.2 = -h1-h2
cbrt.root.3 = cbrt.root.2
end
when s > 0 then do
h1 = 0.5*b2
f = Sqrt((h1)**2+(a2*a2*a2/27)); g = Cbrt(-h1+f); h = Cbrt(-h1-f)
h1 = g+h; h2 = -0.5*h1; h3 = a1/3; h4 = 0.5*Sqrt(3)*(g-h)
cbrt.root.1 = h1-h3
cbrt.root.2 = (h2-h3) (h4)
cbrt.root.3 = (h2-h3) (-h4)
end
end
-- Discard roundoff digits
do i = 1 to 3
cbrt.root.i = Cround(cbrt.root.i,Digits()-4)
end
numeric digits Digits()/2
-- Normalize
do i = 1 to 3
cbrt.root.i = Cnorm(cbrt.root.i)
end
-- Number of roots
cbrt.0 = 3
return 3
Cqteq:
-- Complex roots of quartic equation
-- az^4+bz^3+cz^2+dx+e = 0 with real coefs
procedure expose qtrt. glob.
arg xx; xx = xx 0 0 0 0
parse var xx a b c d e .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
if \ Datatype(c,'n') then say argument
if \ Datatype(d,'n') then say argument
if \ Datatype(e,'n') then say argument
-- Method Abramowitz and Stegun, courtesy PlanetCalc
numeric digits Digits()*2
-- Reduce to x^4+ax^3+bx^2+cx+d
parse value b/a c/a d/a e/a with a3 a2 a1 a0
-- Solve resolvent cubic x^3+ax^2+bx+c for real roots
n = Cbeq((1) (-a2) (a1*a3-4*a0) (4*a0*a2-a1**2-a0*a3**2))
-- Use max root
u1 = cbrt.root.1
do i = 2 to n
u1 = Max(cbrt.root.i,u1)
end
-- Compute coefficients square equations
h = Sqrt(Round(0.25*a3**2+u1-a2,Digits()-2)); p1 = 0.5*a3+h; p2 = 0.5*a3-h
h = Sqrt(Round(0.25*u1**2-a0,Digits()-2)); q1 = 0.5*u1+h; q2 = 0.5*u1-h
if Round(p1*q2+p2*q1-a1,0.5*Digits()) <> 0 then
parse value q2 q1 with q1 q2
-- Find roots by solving ax^2+bx+c for complex roots
qtrt. = 0
call Csqeq 1 p1 q1
qtrt.root.1 = sqrt.root.1
qtrt.root.2 = sqrt.root.2
call Csqeq 1 p2 q2
qtrt.root.3 = sqrt.root.1
qtrt.root.4 = sqrt.root.2
-- Discard roundoff digits
do i = 1 to 4
qtrt.root.i = Cround(qtrt.root.i,Digits()-4)
end
-- Normalize
numeric digits Digits()/2
do i = 1 to 4
qtrt.root.i = Cnorm(qtrt.root.i)
end
-- Number of roots
qtrt.0 = 4
return 4
Csqeq:
-- Complex roots of square equation
-- az^2+bz+c = 0 with real coefs
procedure expose sqrt. glob.
arg xx; xx = xx 0 0; parse var xx a b c .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
if \ Datatype(c,'n') then say argument
-- Abc formula
numeric digits Digits()+5
-- Discriminant
s = Round(b*b-4*a*c,Digits()-3)
-- Roots
sqrt. = 0; h1 = -1/(2*a)
select
when s < 0 then do
h2 = Sqrt(-s)
sqrt.root.1 = (h1*b) (h1*h2)
sqrt.root.2 = (h1*b) (-h1*h2)
end
when s = 0 then do
sqrt.root.1 = h1*b
sqrt.root.2 = sqrt.root.1
end
when s > 0 then do
h2 = Sqrt(s)
sqrt.root.1 = h1*(b+h2)
sqrt.root.2 = h1*(b-h2)
end
end
-- Discard roundoff digits
do i = 1 to 2
sqrt.root.i = Cround(sqrt.root.i,Digits()-3)
end
numeric digits Digits()-5
-- Normalize
do i = 1 to 2
sqrt.root.i = Cnorm(sqrt.root.i)
end
-- Number of roots
sqrt.0 = 2
return 2
Cbeq:
-- Real roots of cubic equation
-- ax^3+bx^2+cx+d = 0
procedure expose cbrt. glob.
arg xx; xx = xx 0 0 0 0
parse var xx a b c d .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
if \ Datatype(c,'n') then say argument
if \ Datatype(d,'n') then say argument
-- Method Cardano et al
numeric digits Digits()*2+2
-- Reduce to x^3+ax^2+bx+c
parse value b/a c/a d/a with a1 b1 c1
-- Reduce to x^3+ax+b
h1 = a1*a1; h2 = h1*a1
a2 = b1-h1/3; b2 = -h2/27+h2/9-a1*b1/3+c1
-- Discriminant
s = Round(b2*b2/4+a2*a2*a2/27,Digits()-2)
-- Roots
cbrt. = 0
select
when s < 0 then do
h1 = Abs(a2)/3
f = Arccos(-b2/(2*Sqrt(h1*h1*h1)))
h1 = 2*Sqrt(h1); h2 = f/3; h3 = Pi()/3; h4 = a1/3
x1 = h1*Cos(h2)-h4; x2 = -h1*Cos(h2-h3)-h4; x3 = -h1*Cos(h2+h3)-h4
zz = 1; cbrt.root.zz = x1
if x2 <> x1 then do
zz = zz+1; cbrt.root.zz = x2
end
if x3 <> x1 & x3 <> x2 then do
zz = zz+1; cbrt.root.zz = x3
end
end
when s = 0 then do
h1 = Cbrt(-0.5*b2); h2 = a1/3
x1 = 2*h1-h2; x2 = -h1-h2
zz = 1; cbrt.root.zz = 2*h1-h2
if x2 <> x1 then do
zz = zz+1; cbrt.root.zz = -h1-h2
end
end
when s > 0 then do
h1 = 0.5*b2
f = Sqrt((h1)**2+(a2*a2*a2/27)); g = Cbrt(-h1+f); h = Cbrt(-h1-f)
h1 = g+h; h3 = a1/3
zz = 1; cbrt.root.zz = h1-h3
end
end
-- Normalize
numeric digits Digits()/2
do i = 1 to zz
cbrt.root.i = Round(cbrt.root.i,Digits()-2)
end
-- Number of roots
cbrt.0 = zz
return zz
Lieq:
-- Real roots of linear equation
-- ax+b = 0
procedure expose lirt.
arg xx; xx = xx 0 0
parse var xx a b .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
-- Solution = -b/a
numeric digits Digits()*2
lirt. = 0; lirt.root.1 = -b/a
numeric digits Digits()/2
-- Normalize
lirt.root.1 = lirt.root.1+0
-- Number of roots
lirt.0 = 1
return 1
Qteq:
-- Real roots of quartic equation
-- az^4+bz^3+cz^2+dx+e = 0 with real coefs
procedure expose qtrt. glob.
arg xx; xx = xx 0 0 0 0
parse var xx a b c d e .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
if \ Datatype(c,'n') then say argument
if \ Datatype(d,'n') then say argument
if \ Datatype(e,'n') then say argument
-- Method Abramowitz and Stegun, courtesy PlanetCalc
numeric digits Digits()*2
-- Reduce to x^4+ax^3+bx^2+cx+d
parse value b/a c/a d/a e/a with a3 a2 a1 a0
-- Solve resolvent cubic x^3+ax^2+bx+c for real roots
n = Cbeq((1) (-a2) (a1*a3-4*a0) (4*a0*a2-a1**2-a0*a3**2))
-- Use max root
u1 = cbrt.root.1
do i = 2 to n
u1 = Max(cbrt.root.i,u1)
end
-- Compute coefficients square equations
h = Sqrt(Round(0.25*a3**2+u1-a2,Digits()-2)); p1 = 0.5*a3+h; p2 = 0.5*a3-h
h = Sqrt(Round(0.25*u1**2-a0,Digits()-2)); q1 = 0.5*u1+h; q2 = 0.5*u1-h
if Round(p1*q2+p2*q1-a1,0.5*Digits()) <> 0 then
parse value q2 q1 with q1 q2
-- Find roots by solving ax^2+bx+c for real roots
qtrt. = 0; zz = 0
call Sqeq 1 p1 q1
do i = 1 to sqrt.0
zz = zz+1; qtrt.root.zz = sqrt.root.i
end
call Sqeq 1 p2 q2
do i = 1 to sqrt.0
zz = zz+1; qtrt.root.zz = sqrt.root.i
end
qtrt.0 = zz
-- Discard roundoff digits
do i = 1 to zz
qtrt.root.i = Round(qtrt.root.i,Digits()-4)
end
-- Normalize
numeric digits Digits()/2
do i = 1 to zz
qtrt.root.i = qtrt.root.i/1
end
-- Number of roots
return zz
Sqeq:
-- Real roots of square equation
-- ax^2+bx+c = 0
procedure expose sqrt.
arg xx; xx = xx 0 0 0
parse var xx a b c .
if \ Datatype(a,'n') then say argument
if a = 0 then say argument
if \ Datatype(b,'n') then say argument
if \ Datatype(c,'n') then say argument
-- Discriminant
numeric digits Digits()+2
d = Trunc(b*b-4*a*c,Digits())
-- Abc formula
h1 = -0.5/a
select
when d < 0 then do
zz = 0
end
when d = 0 then do
zz = 1
sqrt.root.1 = h1*b
end
when d > 0 then do
zz = 2; h2 = Sqrt(d)
sqrt.root.1 = h1*(b+h2)
sqrt.root.2 = h1*(b-h2)
end
end
sqrt.0 = zz
numeric digits Digits()-2
-- Normalize
do i = 1 to zz
sqrt.root.i = Trunc(sqrt.root.i,Digits())/1
end
-- Number of roots
return zz
-- Module Sequences.inc - 00 Xxx 9999
-- Generates and saves number sequences
-- Results are stored in 4-letter named stems <ssss.>
-- <ssss>.0 holds the number of entries
Abundants:
-- Abundant numbers
procedure expose abun.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
abun. = 0
if xx < 101 then do
a = '12 18 20 24 30 36 40 42 48 54 56 60 66 70 72 78 80 84 88 90 96 100 999'
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
abun.zz = w
end
zz = zz-1; abun.0 = zz
return zz
end
-- Check if sum(divisors) > 2*number
zz = 0
do i = 1 to xx
s = Sigma(i)
if s > 2*i then do
zz = zz+1; abun.zz = i
end
end
abun.0 = zz
return zz
Additives:
-- Additive primes
procedure expose addi. prim.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
addi. = 0
if xx < 2 then
return 0
if xx < 101 then do
a = '2 3 5 7 11 23 29 41 43 47 61 67 83 89 999'
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
addi.zz = w
end
zz = zz-1; addi.0 = zz
return zz
end
-- Get primes
p = Primes(xx)
-- Collect additive primes
zz = 0
do i = 1 to p
q = prim.i; s = 0
do j = 1 to Length(q)
s = s+Substr(q,j,1)
end
if Prime(s) then do
zz = zz+1; addi.zz = q
end
end
-- Return number of additive primes
addi.0 = zz
return zz
Amicables:
-- Amicable numbers
procedure expose glob. amic.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Scan for amicable pairs
amic. = 0; zz = 0
do i = 1 to xx
s = Sigma(i)-i; glob.Sigma.i = s
if i = glob.Sigma.s then do
if s = i then
iterate
zz = zz+1; amic.1.zz = s; amic.2.zz = i
end
end
amic.0 = zz
return zz
Arithmetics:
-- Arithmetic numbers
procedure expose arit. divi.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
arit. = 0
if xx < 101 then do
a = '1 3 5 6 7 11 13 14 15 17 19 20 21 22 23 27 29 30 31 33 ',
|| '35 37 38 39 41 42 43 44 45 46 47 49 51 53 54 55 56 57 59 60 61 62 65 66 ',
|| '67 68 69 70 71 73 77 78 79 83 85 86 87 89 91 92 93 94 95 96 97 99 999 '
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
arit.zz = w
end
zz = zz-1; arit.0 = zz
return zz
end
-- Check if number is Arithmetic
zz = 0
do i = 1 to xx
if Arithmetic(i) then do
zz = zz+1; arit.zz = i
end
end
arit.0 = zz
return zz
Bernoullis:
-- Bernoulli numbers
procedure expose bern. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- From generating function
numeric digits Max(Digits(),xx+xx)
bern. = 0; bern.0 = 1; bern.1 = -0.5
do n = 2 by 2 to xx
s = 0; n1 = n+1
do k = 0 to n-1
s = s+bern.k*Comb(n1,k)
end
bern.n = -s/n1
end
bern.0 = xx
return xx
Carmichaels:
-- Carmichael 3 strong primes
procedure expose carm.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Janeson
carm. = 0; zz = 0
do p1 = 3 by 2 to xx
if \ Prime(p1) then
iterate p1
pm = p1-1; ps = -p1*p1
do h3 = 1 to pm
t1 = (h3+p1)*pm; t2 = ps//h3
if t2 < 0 then
t2 = t2+h3
do d = 1 to h3+pm
if t1//d <> 0 then
iterate d
if t2 <> d//h3 then
iterate d
p2 = 1+t1%d
if \ Prime(p2) then
iterate d
p3 = 1+(p1*p2%h3)
if \ Prime(p3) then
iterate d
if (p2*p3)//pm <> 1 then
iterate d
zz = zz+1; carm.1.zz = p1; carm.2.zz = p2; carm.3.zz = p3
end
end
end
carm.0 = zz
-- Return count
return zz
Catalans:
-- Catalan numbers
procedure expose cata.
arg xx
if \ Datatype(xx,'n') then say argument
if xx = 0 then say argument
-- Count or threshold?
if xx < 0 then do
xm = Maxnumber(); nm = Abs(xx)
numeric digits nm
end
else do
xm = xx; nm = Maxnumber()
numeric digits Xpon(xm)+2
end
-- Recurring
cata. = 0; c = 1
do zz = 1 to nm until c > xm
c = c*(4*zz-2)/(zz+1); cata.zz = c
end
zz = zz-1; cata.0 = zz
return zz
Chowlas:
-- Chowla numbers
procedure expose chow.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Call aliquot
chow. = 0
do i = 1 to xx
chow.i = Aliquot(i)-1
end
chow.0 = xx
-- Return count
return xx
Combinations:
-- Combinations
procedure expose comb.
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
-- Recurring definition
comb. = 1
-- Comb(xx,n) for n = 1...xx
if xx < 0 then do
xx = -xx; m = xx%2; a = 1
do i = 1 to m
a = a*(xx-i+1)/i; comb.xx.i = a
end
do i = m+1 to xx-1
j = xx-i; comb.xx.i = comb.xx.j
end
zz = xx+1
end
-- Comb(m,n) for m = 1...xx n = 1...xx
else do
do i = 1 to xx
i1 = i-1
do j = 1 to i1
j1 = j-1; comb.i.j = comb.i1.j1+comb.i1.j
end
end
zz = (xx*xx+3*xx+2)/2
end
comb.0 = zz
return zz
Composites:
-- Composite numbers = not prime
procedure expose comp.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
comp. = 0
if xx < 101 then do
c = '4 6 8 9 10 12 14 15 16 18 20 21 22 24 25 26 27 28 30 32 33 34 35 36 38 39 ',
|| '40 42 44 45 46 48 49 50 51 52 54 55 56 57 58 60 62 63 64 65 66 68 69 70 72 ',
|| '74 75 76 77 78 80 81 82 84 85 86 87 88 90 91 92 93 94 95 96 98 99 100 999 '
do zz = 1 to Words(c)
w = Word(c,zz)
if w > xx then
leave
comp.zz = w
end
zz = zz-1; comp.0 = zz
return zz
end
-- Sieve of Eratosthenes
zz. = 0
do i = 2 to xx while i*i <= xx
do j = i+i by i to xx
zz.j = 1
end
end
zz = 0
do i = 4 to xx
if zz.i = 1 then do
zz = zz+1; comp.zz = i
end
end
comp.0 = zz
return zz
Cubans:
-- Cuban primes
procedure expose cuba.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Generate
cuba. = 0; zz = 0
a = 1; b = 8; c = 2
do forever
v = b-a
if v > xx then
leave
if Prime(v) then do
zz = zz+1; cuba.zz = v
end
a = b; c = c+1; b = c*c*c
end
cuba.0 = zz
return zz
Cubics:
-- Cubic primes
procedure expose cubi. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Get primes
cubi. = 0; cubi.1 = 2; zz = 1
call Primes xx
-- Generate
a = 2; b = 3; k = 1
do while b <= xx
b = a + k*k*k
if flag.b = 1 then do
zz = zz+1; cubi.zz = b
a = b; k = 1
end
k = k+1
end
cubi.0 = zz
return zz
Cyclops:
-- Cyclop numbers
procedure expose cycl.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Generate
cycl. = 0; m = 0; zz = 1
do f = 1
h = m+1; m = zz
do i = 1 to 9
do j = h to m
a = cycl.j
do k = 1 to 9
b = i||a||k
if b > xx then
leave f
zz = zz+1; cycl.zz = b
end
end
end
end
cycl.0 = zz
return zz
Deficients:
-- Deficient numbers
procedure expose defi.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
defi. = 0
if xx < 101 then do
a = '1 2 3 4 5 7 8 9 10 11 13 14 15 16 17 19 21 22 23 25 26 27 29 31 32 33 34 35 ',
|| '37 38 39 41 43 44 45 46 47 49 50 51 52 53 55 57 58 59 61 62 63 64 65 67 68 ',
|| '69 71 73 74 75 76 77 79 81 82 83 85 86 87 89 91 92 93 94 95 97 98 99 999'
do zz = 1 to Words(a)
w = Word(a,zz)
if w > xx then
leave
defi.zz = w
end
zz = zz-1; defi.0 = zz
return zz
end
-- Check if sum(divisors) < 2*number
zz = 0
do i = 1 to xx
s = Sigma(i)
if s < 2*i then do
zz = zz+1; defi.zz = i
end
end
defi.0 = zz
return zz
Divisors:
-- Divisors of an integer
procedure expose divi. wrkl. wrkr.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Fast
divi. = 0
if xx = 1 then do
divi.1 = 1; divi.0 = 1
return 1
end
-- Euclid's method
a = 1; wrkl.1 = 1; b = 1; wrkr.1 = xx
m = xx//2
do j = 2+m by 1+m while j*j < xx
if xx//j = 0 then do
a = a+1; wrkl.a = j; b = b+1; wrkr.b = xx%j
end
end
if j*j = xx then do
a = a+1; wrkl.a = j
end
-- Save in table
zz = 0
do i = 1 to a
zz = zz+1; divi.zz = wrkl.i
end
do i = b by -1 to 1
zz = zz+1; divi.zz = wrkr.i
end
divi.0 = zz
-- Return number of divisors
return zz
Efactors:
-- Prime factors with exponent
procedure expose efac.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
efac. = 0
if xx = 1 then
return 0
if Prime(xx) then do
efac.factor.1 = xx; efac.exponent.1 = 1; efac.0 = 1
return 1
end
-- Check low factors
zz = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i); e = 0
do while xx//p = 0
e = e+1; xx = xx%p
end
if e > 0 then do
zz = zz+1; efac.factor.zz = p; efac.exponent.zz = e
end
end
-- Check higher factors
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then do
e = 0
do while xx//j = 0
e = e+1; xx = xx%j
end
end
if e > 0 then do
zz = zz+1; efac.factor.zz = j; efac.exponent.zz = e
end
if p = 3 then
iterate
yy = j+2; e = 0
do while xx//yy = 0
e = e+1; xx = xx%yy
end
if e > 0 then do
zz = zz+1; efac.factor.zz = yy; efac.exponent.zz = e
end
end
if xx > 1 then do
if e = 0 then do
zz = zz+1; efac.factor.zz = xx; efac.exponent.zz = 1
end
end
efac.0 = zz
-- Return number of factor pairs
return zz
Emirps:
-- Reversed primes
procedure expose emir. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Get primes
emir. = 0; zz = 0
call Primes xx
-- Generate
do i = 13 by 2 to xx
a = Reverse(i)
if a = i then
iterate i
if \ flag.i then
iterate i
if \ flag.a then
iterate i
zz = zz+1; emir.zz = i
end
emir.0 = zz
return zz
Erdoss:
-- Erdos primes
procedure expose erdo. prim. flag. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Get primes
erdo. = 0; zz = 0
call Primes xx
-- Generate
do i = 1 to prim.0
a = prim.i
do j = 1
k = Fact(j)
if k >= a then
leave j
p = a-k
if flag.p then
iterate i
end
zz = zz+1; erdo.zz = a
end
erdo.0 = zz
return zz
Euclids:
-- Euclid-Mullin numbers
procedure expose eucl.
arg xx
numeric digits 100
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
if xx > 17 then say argument
-- Using Ffactor = first prime factor
eucl. = 0; eucl.1 = 2; eucl.0 = 1
p = 2
do i = 2 to xx
z = p+1; t = Ffactor(z); eucl.i = t; p = p*t
end
-- Save and return number
eucl.0 = xx
return xx
Extras:
-- Extra primes
procedure expose extr.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Generate
extr. = 0; extr.1 = 2; zz = 1
do i = 3 by 2 to xx
do j = 1 to Length(i)
if Pos(Substr(i,j,1),'014689') > 0 then
iterate i
end
if \ Prime(Digitsum(i)) then
iterate i
if \ Prime(i) then
iterate i
zz = zz+1; extr.zz = i
end
extr.0 = zz
return zz
Factors:
-- Prime factors
procedure expose fact. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 2 then say argument
-- Fast
fact. = 0
if Prime(xx) then do
fact.1 = xx; fact.0 = 1
return 1
end
-- Check low factors
zz = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i)
do while xx//p = 0
zz = zz+1; fact.zz = p
xx = xx%p
end
end
-- Check higher factors
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then do
do while xx//j = 0
zz = zz+1; fact.zz = j
xx = xx%j
end
end
if p = 3 then
iterate
k = j+2
do while xx//k = 0
zz = zz+1; fact.zz = k
xx = xx%k
end
end
-- Last factor
if xx > 1 then do
zz = zz+1; fact.zz = xx
end
fact.0 = zz
-- Return number of factors
return zz
Fareys:
-- Farey fractions
procedure expose fare.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Recurring
fare. = 0; fare.1 = '0/1'; zz = 1
parse value 0 1 1 xx with aa bb cc dd
do while cc <= xx
k = Trunc((xx+bb)/dd)
parse value cc dd k*cc-aa k*dd-bb with aa bb cc dd
zz = zz+1; fare.zz = aa'/'bb
end
fare.0 = zz
-- Number of entries
return zz
Fibonaccis:
-- Fibonacci numbers
procedure expose fibo.
arg xx
if \ Datatype(xx,'w') then say argument
-- Count or threshold?
if xx < 0 then do
xm = Maxnumber(); nm = Abs(xx)
numeric digits nm%2
end
else do
xm = xx; nm = Maxnumber()
numeric digits xpon(xm)+2
end
-- Fast
fibo. = 0
if xm = 0 then
return 0
fibo.1 = 1
if xm = 1 | nm = 1 then do
fibo.0 = 1
return 1
end
-- Recurring
b = 1; c = 1
do zz = 2 until zz = nm | c > xm
fibo.zz = c; a = b; b = c; c = a+b
end
fibo.0 = zz
return zz
Frobeniuss:
-- Frobenius primes
procedure expose frob. prim.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Get primes
frob. = 0; zz = 0
p = Primes(xx%10)
-- Generate
do i = 1 to p-1
j = i+1; f = prim.i * prim.j - prim.i - prim.j
if f > xx then
leave
zz = zz+1; frob.zz = f
end
frob.0 = zz
return zz
Giugas:
-- Giuga numbers
procedure expose giug.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
giug. = 0
if xx < 30 then
return 0
g = '30 858 1722 66198 2214408306 24423128562 432749205173838 14737133470010574 ',
|| '550843391309130318 244197000982499715087866346 554079914617070801288578559178 ',
|| '1910667181420507984555759916338506 999999999999999999999999999999999 '
do zz = 1 to Words(g)
w = Word(g,zz)
if w > xx then do
zz = zz-1
leave
end
giug.zz = w
end
giug.0 = zz
return zz
Hammings:
-- Hamming numbers
procedure expose hamm.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Ensure enough digits
hamm. = 1
numeric digits 2**(Length(Std(xx))+1)
-- Dijkstra
x2 = 2; x3 = 3; x5 = 5
i2 = 1; i3 = 1; i5 = 1
do zz = 2
h = x2
if x3 < h then
h = x3
if x5 < h then
h = x5
hamm.zz = h
if zz = xx then
leave
if x2 = h then do
i2 = i2+1; x2 = hamm.i2+hamm.i2
end
if x3 = h then do
i3 = i3+1; x3 = hamm.i3+hamm.i3+hamm.i3
end
if x5 = h then do
i5 = i5+1; x5 = hamm.i5*5
end
end
hamm.0 = zz
return zz
Highcomposites:
-- Highly composite numbers
procedure expose high.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Thresholds and increments
high. = 0
a = '1 2 6 60 840 55440 720720 61261200 2327925600 321253732800 9999999999999'
b = '1 2 6 60 420 27720 360360 12252240 232792560 80313433200 9999999999999'
c = Words(a)-1; m = 0; zz = 0
-- Collect cf definition
do i = 1 to c
do j = Word(a,i) by Word(b,i) to Min(xx,Word(a,i+1)-1)
d = Divisors(j)
if d > m then do
zz = zz+1; high.zz = j
m = d
end
end
end
high.0 = zz
-- Return count
return zz
Humbles:
-- Humble numbers
procedure expose humb.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Ensure enough digits
humb. = 1
numeric digits 2**(Length(xx)+1)
-- Dijkstra
x2 = 2; x3 = 3; x5 = 5; x7 = 7
i2 = 1; i3 = 1; i5 = 1; i7 = 1
do zz = 2
h = x2
if x3 < h then
h = x3
if x5 < h then
h = x5
if x7 < h then
h = x7
humb.zz = h
if zz = xx then
leave
if x2 = h then do
i2 = i2+1; x2 = humb.i2+humb.i2
end
if x3 = h then do
i3 = i3+1; x3 = humb.i3+humb.i3+humb.i3
end
if x5 = h then do
i5 = i5+1; x5 = humb.i5*5
end
if x7 = h then do
i7 = i7+1; x7 = humb.i7*7
end
end
humb.0 = zz
return zz
Kaprekars:
-- Kaprekar numbers
procedure expose kapr.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
kapr. = 0
if xx < 101 then do
k = '1 9 45 55'
do zz = 1 to Words(k)
w = Word(k,zz)
if w > xx then
leave
kapr.zz = w
end
zz = zz-1; kapr.0 = zz
return zz
end
-- Ensure enough digits
xx = xx+1
p = Max(Digits(),2*Length(xx))
numeric digits p
zz = 1; kape.1 = 1
-- Method casting out nines
do i = 9 for xx-9
a = i//9
if a > 2 then
iterate i
s = i*i
if a = s//9 then do
do j = 1 to Length(s)%2
parse var s l +(j) r
if i <> l+r then
iterate j
zz = zz+1; kapr.zz = i
leave
end
end
end
kapr.0 = zz
return zz
Ludics:
-- Ludic numbers
procedure expose ludi. work.
arg xx
if \ Datatype(xx,'n') then say argument
if xx < 1 then say argument
-- Use a sparse array
work. = 1
do i = 2 to xx%2
if work.i then do
n = 0
do j = i+1 to xx
if work.j then do
n = n+1
if n//i = 0 then do
work.j = 0
end
end
end
end
end
-- Save results
zz = 0
do i = 1 to xx
if work.i then do
zz = zz+1; ludi.zz = i
end
end
-- Count
ludi.0 = zz
return zz
Magics:
-- Magic numbers
procedure expose magi.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Formula
magi. = 0
do zz = 1 until m > xx
m = Magic(zz); magi.zz = m
end
magi.0 = zz
return zz
Magnanimouss:
-- Magnanimous numbers
procedure expose magn. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 0 then say argument
-- Get primes below xx
call Primes xx
magn. = 0; zz = 0
do i = 0 to xx
-- 0...9 always apply
if i < 10 then do
zz = zz+1; magn.zz = i
iterate i
end
-- Numbers having parity do not apply
if i > 1001 then do
parse var i 1 l 2 '' -1 r
if l//2 = r//2 then
iterate i
end
-- Generate candidates and test primality
do j = 1 to Length(i)-1
parse var i 1 l +(j) r
s = l+r
if \ flag.s then
iterate i
end
zz = zz+1; magn.zz = i
end
magn.0 = zz
return zz
Motzkins:
-- Motzkin numbers
procedure expose motz.
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
-- Fast
motz. = 0; motz.0 = 1; motz.1 = 1; motz.0 = 1
if xx = 1 then
return 1
-- Recurring
a = 1; b = 1; c = 0
do n = 2 while c < xx
c = ((3*n-3)*a+(2*n+1)*b)/(n+2)
motz.n = c
a = b; b = c
end
motz.0 = n-1
return motz.0
Perfects:
-- Perfect numbers
procedure expose perf.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
perf. = 0
if xx < 6 then
return 0
p = '6 28 496 8128 33550336 8589869056 137438691328 2305843008139952128 ',
|| '2658455991569831744654692615953842176 ',
|| '191561942608236107294793378084303638130997321548169216 ',
|| '13164036458569648337239753460458722910223472318386943117783728128 ',
|| '14474011154664524427946373126085988481573677491474835889066354349131199152128 ',
|| '99999999999999999999999999999999999999999999999999999999999999999999999999999'
do zz = 1 to Words(p)
w = Word(p,zz)
if w > xx then
leave
perf.zz = w
end
zz = zz-1; perf.0 = zz
return zz
Practicals:
-- Practical numbers
procedure expose prac. fact.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
prac. = 0
if xx < 101 then do
p = '1 2 4 6 8 12 16 18 20 24 28 30 32 36 40 42 ',
|| '48 54 56 60 64 66 72 78 80 84 88 90 96 100 999 '
do zz = 1 to Words(p)
w = Word(p,zz)
if w > xx then
leave
prac.zz = w
end
zz = zz-1; prac.0 = zz
return zz
end
-- Using function practical()
zz = 1; prac.1 = 1
do i = 2 by 2 to xx
if Practical(i) then do
zz = zz+1; prac.zz = i
end
end
prac.0 = zz
return zz
Primes:
-- Primes
procedure expose prim. work. flag. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
-- Fast
prim. = 0; flag. = 0
if xx = 1 then
return 0
-- Nth prime or threshold?
if xx < 0 then
xx = Primeest(Abs(xx))
-- Sieve of Eratosthenes odd integers
work. = 1
do i = 3 by 2 while i*i <= xx
if work.i then do
do j = i*i by i+i to xx
work.j = 0
end
end
end
-- Save and flag results
prim.1 = 2; flag.2 = 1; zz = 1
do i = 3 by 2 to xx
if work.i then do
zz = zz+1; prim.zz = i; flag.i = 1
end
end
prim.0 = zz
return zz
Primorials:
-- Primorials
procedure expose prim. prmo. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
if xx < -50000 then say argument
if xx > 1e250000 then say argument
-- Get primes
numeric digits Max(Digits(),Xpon(xx)+1)
if xx < 0 then
call Primes(xx)
else
call Primes(600000)
-- Multiply primes
prmo. = 0
p = 1; zz = 0
do i = 1 to prim.0
numeric digits xpon(p)+6
p = p*prim.i
if xx > 0 then
if p > xx then
leave
zz = zz+1; prmo.zz = p
end
prmo.0 = zz
return zz
Quadratics:
-- Quadratic primes
procedure expose quad. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Get primes
quad. = 0; quad.1 = 2; zz = 1
call Primes xx
-- Generate
a = 2; b = 3; k = 1
do while b <= xx
b = a + k*k
if flag.b = 1 then do
zz = zz+1; quad.zz = b
a = b; k = 1
end
k = k+1
end
quad.0 = zz
return zz
Reversables:
-- Reversable primes
procedure expose reve. prim. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Get primes
reve. = 0; zz = 0
e = Xpon(xx)+1; p = Primes(1'E'e)
-- Generate
do i = 1 to p
a = prim.i
if a > xx then
leave i
b = Reverse(a)
if \ flag.b then
iterate i
zz = zz+1; reve.zz = a
end
reve.0 = zz
return zz
Smooths:
-- Smooth numbers
procedure expose prim. smoo. numb. indi.
arg xx,yy
if \ Datatype(xx,'n') then say argument
if Abs(xx) < 1 then say argument
if \ Prime(yy) then say argument
-- Count or threshold?
if xx < 0 then do
xm = Maxnumber(); nm = Abs(xx)
end
else do
xm = xx; nm = Maxnumber()
end
-- Get and number primes
call Primes yy
numb. = 0
do i = 1 to prim.0
p = prim.i; numb.p = i
end
-- Collect smooth numbers
smoo. = 0; smoo.1 = 1; smoo.0 = 1; indi. = 1
n = numb.yy
do j = 2 to nm
i = indi.1; z = prim.1*smoo.i
do k = 2 to n
i = indi.k; v = prim.k*smoo.i
if v < z then
z = v
end
if z > xm then
leave j
smoo.j = z
do d = 1 to n
i = indi.d
if prim.d*smoo.i = z then
indi.d = indi.d+1
end
end
zz = j-1; smoo.0 = zz
return zz
Squarefrees:
-- Square free numbers
procedure expose squa. work. prim. glob. flag.
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
-- Fast
squa. = 0
if xx = 1 then do
squa.1 = 1; flag.1 = 1; squa.0 = 1
return 1
end
-- Precision
numeric digits Max(Digits(),Xpon(xx)+1)
-- Get primes
n = Primes(Isqrt(xx))
-- Get rid of multiples of squared primes
work. = 1
do i = 1 to n
p = prim.i; p = p*p
do j = p by p to xx
work.j = 0
end
end
-- Save results
zz = 0
do i = 1 to xx
if work.i then do
zz = zz+1; squa.zz = i; flag.i = 1
end
end
squa.0 = zz
return zz
Squarefrees:
-- Square free numbers
procedure expose squa. work. prim. glob.
arg xx
if \ Datatype(xx,'w') then say argument
if xx = 0 then say argument
-- Fast
squa. = 0
if xx = 1 then do
squa.free.1 = 1; squa.flag.1 = 1; squa.0 = 1
return 1
end
-- Precision
numeric digits Max(Digits(),Xpon(xx)+1)
-- Get primes
n = Primes(Isqrt(xx))
-- Get rid of multiples of squared primes
work. = 1
do i = 1 to n
p = prim.i; p = p*p
do j = p by p to xx
work.j = 0
end
end
-- Save results
zz = 0
do i = 1 to xx
if work.i then do
zz = zz+1; squa.free.zz = i; squa.flag.i = 1
end
end
squa.0 = zz
return zz
Totients:
-- Euler's totient numbers
procedure expose toti.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Recurring sequence
toti. = 0
do i = 1 to xx
toti.i = i
end
do i = 2 to xx
if toti.i < i then
iterate i
do j = i by i to xx
toti.j = toti.j-toti.j/i
end
end
toti.0 = xx
return xx
Ufactors:
-- Unique prime factors
procedure expose ufac.
arg xx
if \ Datatype(xx,'w') then say argument
if xx < 1 then say argument
-- Fast
ufac. = 0
if xx = 1 then
return 0
if xx < 4 then do
ufac.1 = xx; ufac.0 = 1
return 1
end
if Prime(xx) then do
ufac.1 = xx; ufac.0 = 1
return 1
end
-- Check low factors
zz = 0
pr = '2 3 5 7 11 13 17 19 23'
do i = 1 to Words(pr)
p = Word(pr,i); v = xx
do while xx//p = 0
xx = xx%p
end
if xx <> v then do
zz = zz+1; ufac.zz = p
end
end
-- Check higher factors
do j = 29 by 6 while j*j <= xx
p = Right(j,1)
if p <> 5 then do
v = xx
do while xx//j = 0
xx = xx%j
end
if xx <> v then do
zz = zz+1; ufac.zz = j
end
end
if p = 3 then
iterate
yy = j+2; v = xx
do while xx//yy = 0
xx = xx%yy
end
if xx <> v then do
zz = zz+1; ufac.zz = yy
end
end
-- Last factor
if xx > 1 then do
zz = zz+1; ufac.zz = xx
end
ufac.0 = zz
-- Return number of factors
return zz
-- Module Settings.inc - 28 Jan 2025
-- Settings used by many other modules
call Lines 'dummy'
signal on halt name Abend
signal on notready name Abend
signal on novalue name Abend
signal on syntax name Abend
parse version version; parse source . . source
glob. = ''; call Time('r')