## Monday, November 23, 2015

### Some arithmetical functions

The radical of a natural number is the product of all it's distinct prime factors. An easy way to calculate the radical is to factorize and sort:

: radical1 ( n -- r )
1 dup loc{ n radical prime }
n pollard# sort 0
do dup prime = 0=
then to prime

But instead of using the local variable radical to build up the product one can use the stack:

: radical2 ( n -- r )
1 loc{ n prime }
n pollard# sort 1 swap 0 \ 1 will be at tos in loop
do over prime =
if nip
else over to prime *
then
loop ;

And if the code don't look for the first a occurrence of a prime factor but for the last occurrence of that factor, then there is no need for any local variables.

: undequ ( a b c -- a b c f )
>r 2dup = r> swap ;

The word undequ tests if a=b and put a flag on top of stack.

: radical ( n -- r )
1 swap                \ a non prime dummy at bottom for undequ
pollard# sort
1 swap 0              \ 1 will be at tos in loop
do undequ if nip else * then
loop nip ;            \ nip drops the dummy

The Euler totient function gives the number of positive coprimes less than the number. Euler's product formula implies that
: totients ( n -- t )
1 swap                \ a non prime dummy at bottom for undequ
pollard# sort
1 swap 0              \ 1 will be at tos in loop
do undequ if * else swap 1- * then
loop nip ;            \ nip drops the dummy

Some words that are needed:

: drops ( a1...an n -- )
0 ?do drop loop ;

: ** ( b e -- b^e )
dup 0= if 2drop 1 exit then
1 swap 0 do over * loop nip ;

: m** ( b e -- d )
dup 0= if 2drop 1. exit then
1 swap 1. rot 0
do 2over m*/ loop 2nip ;

: -1** ( n -- i )     \ exponentiation of (-1)
1 and if -1 else 1 then ;

: unit* ( i j -- k )  \ multiplication of units ±1
xor 1+ ;

: ufaculty ( u -- u! )
dup 2 u< if drop 1 exit then
dup 1- recurse * ;

: umfaculty ( u -- ud )
dup 2 u< if drop 1. exit then
dup 1- recurse rot 1 m*/ ;

: oddpart ( a -- i m ) \ m is odd and a=m*2^i
0 swap
begin 2/mod swap 0=
while 1 under+
repeat 2* 1+ ;

The Möbius function is defined for positive integers and is zero if it contains a squared prime factor; otherwise 1 if the number of prime factors is even and else it is -1.

: mobius ( n -- my )
dup 2 u< if drop 1 exit then
1 swap
pollard# sort
dup true loc{ facts sqrfree } 0
do 2dup =
if false to sqrfree
facts i - drops leave
then drop
loop sqrfree
if facts 1 and
if -1
else 1
then
else 0
then nip ;

The word bigomega counts the number of prime factors of a number, and smallomega counts the number of distinct prime factors.

: bigomega ( n -- b )
dup 2 u< if drop 0 exit then
pollard# dup >r drops r> ;

: smallomega ( n -- s )
dup 2 u< if drop 0 exit then
1 swap
pollard# sort 0 swap 0
do undequ 0= if 1+ then nip
loop nip ;

The Legendre symbol (a/p), where p is an odd prime and a is positive, is defined as zero if p divides a; otherwise 1 if there exist a x with a=x^2(mod p) else -1.

The Euler criterion implies that (a/p)=a^((p-1)/2)(mod p):

: legendre ( a p -- i )
tuck dup 1- 2/ swap u**mod dup 1 u>
if swap - else nip then ;

A generalization is the Jacobi symbol (a/n) defined for all odd numbers n and positive integers a:

: jacobi ( a n -- j )
1 loc{ a n unit }
a n ugcd 1 = 0= if 0 exit then
begin n 1 = a 1 = or if unit exit then
a n mod ?dup 0= if 0 exit then
oddpart to a 1 and
if n dup * 1- 3 rshift -1** unit unit* to unit
then a n ugcd 1 = 0= if 0 exit then
n a to n to a
a 1- n 1- * 2 rshift -1** unit unit* to unit
again ;

And the Kronecker symbol (a/n) is a further generalization defined for all integers a and n:

: kronecker ( a n -- j )
0 loc{ a n unit }
n 0= if a abs 1 = if 1 else 0 then exit then
n dup abs to n 0<
if a 0< if -1 else 1 then
else 1
then to unit
a dup abs to a                         \ old_a
n oddpart to n dup                     \ old_a i i
if a 1 and 0=                          \ old_a i f
if 2drop 0 exit                     \ old_a
else 1 and                          \ old_a f
if a 7 and 1+ 4 and              \ old_a f
if unit negate to unit then   \ old_a
then
then
else drop
then a n jacobi dup 0=                 \ old_a i f
if nip exit then                       \ old_a
unit unit* to unit 0<                  \ f
if n 1- 2/ -1** else 1 then            \ ±1
unit unit* ;