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=
     if dup radical * to radical
     then to prime 
  loop radical ;

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* ;