Polynomials in Z[X]

The separate stacks xst, yst and zst used for the implementation of nested sets, can also be used to dynamical allocation of arrays, here arrays for the coefficients of polynomials. Defining a non nested list of integers

( 3 0 -2 7 -3 1 )  ok

works as long as the list not is interpreted as nested, when the negative numbers signalize a count. (All number elements of nested sets or lists must be non negative.)

Interpreted as non nested linear lists all set stack manipulation words works, so the administration of these polynomials works as  usual in Forth. The list above correspond to the polynomial

( 3 0 -2 7 -3 1 ) poly. 3-2x²+7x³-3x⁴+x⁵ ok

\ Polynomials

: >da \ vect -- vect ad n 
  zst @ zst@ cs tuck cells - swap ; 
\ Gives the address to the first coefficient plus the count
\ of the polynomial at top of stack

: >da2 \ vect2 vect1 -- vect2 vect1 ad2 n2  
  adn2 cell- cell / ;
\ Gives address and count to the second polynomial of stack

: >zst+ \ vect1 m -- vect2 
  zst> swap >zst 2 - >zst ; 
\ Add item to the list

: da \ -- vect ad  
  -1 >zst zst @ ; 
\ Initiate an empty list 

: da. \ vect -- 
  >da 0
  do dup i cells + @ .
  loop zdrop drop ;

\ print the coefficients

The word Z. cannot be used since it interpret the list as nested.

( 3 0 -2 7 -3 1 ) da. 3 0 -2 7 -3 1  ok

Some systems can write exponents:

\ Printing polynomials

create potence 

s" "   s, s" x"  s, s" x²" s, s" x³" s, s" x⁴" s, 

s" x⁵" s, s" x⁶" s, s" x⁷" s, s" x⁸" s, s" x⁹" s, 

s" x¹⁰" s, s" x¹¹" s, s" x¹²" s, s" x¹³" s, s" x¹⁴" s, 

true value lowterm 

: .term \ i n -- 

  ?dup 0= if drop exit then

  dup 0<

  if ." -"

  else lowterm 0= if ." +" then

  then abs dup 1 > 2 pick 0= or

  if 0 <# #s #> type 

  else drop 

  then false to lowterm 

  cells potence + count type ; 

: poly. \ vect --

  true to lowterm

  >da 0 

  do i over i cells + @ .term

  loop zdrop drop ;

Since BigZ is limited to non negative integers greatest common divisor is defined for unsigned integers UGCD and a word GCD for all integers have to be defined:

: gcd \ n1 n2 -- n      \ Greatest common divisor

  2dup or 0= 

  if 2drop 1 exit then 

  abs swap abs 

  2dup u< if swap then  \ smallest of a b on top of stack

  ?dup 0= if exit then  \ return second value if tos is zero

  begin tuck            \ y x y first time b a b

     0 swap um/mod      \ y x 0 y --> y r q

     drop dup 0=        \ y r [r=0]

  until drop ;          \ y

: multgcd \ k1...kn n -- gcd

  1 do gcd loop ;

\ Gives multiple greatest common device

\ Calculation with polynomials

: polynom \ ad n m -- m' 

  locals| m | cells over + cell- 

  dup @ -rot cell- 

  ?do m * i @ + -cell +loop ; 

\ m' is the evaluation of m with polynomial at ad n

: polyn \ vect m -- vect m'

  >da rot polynom ;

\ m' is the evaluation of m with the polynomial vect

: gcoeff \ vect -- vect n

  zst @ cell - @ ;

\ Gives the coefficient of the greatest power

: rrzs \ vect1 -- vect2  "reduce right zeroes"

  begin gcoeff 0= 

  while zst> zst> drop 2 + >zst

  repeat ;

\ Eliminate leading coefficient equal to zero

: poly* \ ad1 n1 ad2 n2 -- vect 

  locals| n2 ad2 n1 ad1 | da drop 

  n1 n2 + 1- 0 

  do 0 i 1+ 0 

     do j i - 0 n2 within i n1 < and

        if i cells ad1 + @ 

           j i - cells ad2 + @ * + 

        then

     loop >zst+

  loop rrzs ; 

\ Multiply polynomials given by arrays

: p* \ vect1 vect2 -- vect3 

  >da2 >da poly*

  znip znip ;

\ Multiply polynomials

( 0 -1 2 1 ) zdup poly. -x+2x²+x³ ok

( 2 0 2 3 ) zdup poly. 2+2x²+3x³ ok

p* poly. -2x+4x²+x⁴+8x⁵+3x⁶ ok

: p+ \ vect1 vect2 -- vect3 

  adn2 nip adn1 nip < if zswap then 

  adn2 drop locals| ad | 

  zst>> cs 0 

  do ad i cells + +! 

  loop rrzs ;

\ Add polynomials

: pnegate \ vect1 -- vect2

  adn1 cell- 0

  do dup i + dup @ negate swap ! cell

  +loop drop ;

\ Negate a polynomial

: p- \ vect1 vect2 -- vect3

  pnegate p+ ;

\ Subtract polynomials

: ps/ \ vect1 n -- vect2

  locals| n |

  >da cells over + swap

  do i @ n / i ! cell +loop ;

\ Divide polynomial with integer

: makepoly \ vect ad n --  name of polynomial

  cr ." : " type space 

  zst> zst> .

  cs 1- 1

  do ." over * " zst> . ." + "

  loop ." * " zst> . ." + ; " ;

\ Prints the definition of a polynomial to be pasted

( 3 0 -2 7 -3 1 ) s" poly1" makepoly

: poly1 1 over * -3 + over * 7 + over * -2 + over * 0 + * 3 + ;  ok

Copying and pasting the output easy define the polynomial POLY1, eventually after some cleaning.

There is a nice theory of integer-valued polynomials by Pólya. That is, polynomial with rational coefficients that gives integer outputs for integer inputs. The set of these polynomials is a subring of Q[X], the ring of all rational polynomials. This subring is denoted int(Z).

This can be used to calculate the fixed prime divisors of polynomials with integer coefficients. Even if the coefficients of a polynomial in Z[X] has the greatest common divisor 1, all the outputs of the polynomial might be divisible by a fixed number. For example the polynomial 

6-9x-2x²+5x⁴ 

always give an output that is divisible by 6, and therefore have the fixed prime divisors 2 and 3.

\ Integer valued polynomials 

: bin*sum \ ad k -- sum 

  locals| k ad |

  k 0= if 1 exit then 0 k 0 

  ?do i cells ad + @ 

     k i choose * + 

  loop ; 

\ Calculate the sum in the figure above

: polyacoeff \ ad1 n1 -- vect 

  da locals| ad2 n1 ad1 |

  ad1 @ >zst+ 

  n1 1

  ?do ad1 n1 i polynom 

     ad2 i bin*sum - >zst+

  loop ; 

\ Calculate the vector (c0,...,cn) from 

\ integer polynomial at ad1 n1

: polya \ ad n m -- m'

  swap -rot locals| m ad | 0 swap 0

  ?do i cells ad + @ 

     m i choose * +

  loop ;

\ m' is the evaluation of m with the pólya function at ad n

: coeffgcd \ vect -- n

  zst>> cs       \ CS transform set count into stack count
  multgcd ;

: fixdiv \ vect -- vect n

  >da            \ get address and count for polynomial

  polyacoeff     \ calculate Pólya's coefficients

  coeffgcd ;

\ The multiple GCD of c0,...,cn is the fixed divisor of the

\ corresponding original polynomial with integer coefficients

Eisensteins criteria: If there exist a prime number p which not divides an but a0,...,an-1, and p² not divide a0, then
a0+a1*x+a2*x²+...+an*x^n is irreducible over the rational numbers.

\ Eisenstein's criteria

: iseisenstein \ vect -- vect flag   "is an Eisenstein polynomial?"

  zdup zst> 2 + zst> abs false 0 locals| p flag  an |

  >zst zst>> cs multgcd abs primes ?dup

  if 0

     do to p flag 0=

        if an p umod 0= 0=

           >da drop @ abs p ^2 umod 0= 0= and

           to flag 

        then

     loop 

  then flag ;

Most polynomial are irreducible but very few are Eisenstein.

Theorem 1 (Chen & al): If the polynomial f(X) in Z[X] is reducible then the number of positive primes of the form f(a) is less then or equal the degree of f(X). For all degrees n there is a reducible polynomial f(X) with different a1,...,an such that f(ai) is a positive prime.

So finding a polynomial giving primes for n+1 different values is finding an irreducible polynomial. But a lot of polynomials have a fixed divisor greater than one, and those can't be proved irreducible by the theorem above. A more relevant test is therefore to combine theorem 1 with the Pólya fix divisor test.

: ischen \ vect -- vect flag  "false may be and true is irreducible"
  fixdiv 1 > if false exit then
  degree locals| n |
  0 bits n / 2e s>f f** f>s 1000 min 0
  do i polyn dup 0>
     if isprime -
        dup n > if leave then
     else drop
     then
     i negate polyn dup 0>
     if isprime -
        dup n > if leave then
     else drop
     then
  loop n > ;

\ Polynomials

\ Dynamical allocation of arrays

: >da \ vect -- vect ad n
  zst @ zst@ cs tuck cells - swap ;
\ Gives the address to the first coefficient plus the count
\ of the polynomial at top of stack

: >da2 \ vect2 vect1 -- vect2 vect1 ad2 n2 
  adn2 cell- cell / ;
\ Gives address and count to the second polynomial of stack

: >xst+ \ vect1 m -- vect2
  xst> swap >xst 2 - >xst ;
\ Add item to the xst list

: >zst+ \ vect1 m -- vect2
  zst> swap >zst 2 - >zst ;
\ Add item to the zst list

: da \ -- vect ad 
  -1 >zst zst @ ;
\ Initiate an empty list

: da. \ vect --
  >da 0
  do dup i cells + @ .
  loop zdrop drop ;
\ print the coefficients

\ Printing polynomials
\ 64 bits systems only

create potence
s" " s,    s" x" s,   s" x²" s,  s" x³" s,  s" x⁴" s,
s" x⁵" s,  s" x⁶" s,  s" x⁷" s,  s" x⁸" s,  s" x⁹" s,
s" x¹⁰" s, s" x¹¹" s, s" x¹²" s, s" x¹³" s, s" x¹⁴" s,
s" x¹⁵" s, s" x¹⁶" s, s" x¹⁷" s, s" x¹⁸" s, s" x¹⁹" s,

true value lowterm
: .term \ i n --
  ?dup 0= if drop exit then
  dup 0<
  if ." -"
  else lowterm 0= if ." +" then
  then abs dup 1 > 2 pick 0= or
  if 0 <# #s #> type
  else drop
  then false to lowterm
  cells potence + count type ;

: p. \ vect --
  true to lowterm
  >da 0
  do i over i cells + @ .term
  loop zdrop drop ;

\ Greatest common divisors for multiple integers

: multgcd \ k1...kn n -- gcd
  dup 0= if exit then
  swap abs swap 1
  ?do swap abs ugcd loop ;
\ Gives multiple greatest common device

\ Calculation with polynomials

: polynom \ ad n m -- m'
  locals| m | cells over + cell-
  dup @ -rot cell-
  ?do m * i @ + -cell +loop ;
\ m' is the evaluation of m with polynomial at ad n

: sbs* \ sb m -- sb*m
  dup 0< xs> xor >xs abs bs* ;

: s>sb \ n -- sb
  dup abs s>b 0< >xs ;

: sbpolynom \ ad n m -- sb
  locals| m | cells over + cell-
  dup @ s>sb cell-
  ?do m sbs* i @ s>sb sb+ -cell +loop ;
\ single input and big output

: polyn \ vect m -- vect m'
  >da rot polynom ;
\ m' is the evaluation of m with the polynomial vect

: sbpolyn \ vect m -- vect sb
  >da rot sbpolynom ;
\ m' is the evaluation of m with the polynomial vect

: gcoeff \ vect -- vect n
  zst @ cell - @ ;
\ Gives the coefficient of the greatest power

: lcoeff \ vect -- vect n
  >da drop @ ;
\ Gives the coefficient of the constant term

: rrzs \ vect1 -- vect2  "reduce right zeroes"
  begin gcoeff 0=
  while zst> zst> drop 2 + >zst
  repeat ;
\ Eliminate leading coefficient equal to zero

: poly* \ ad1 n1 ad2 n2 -- vect
  locals| n2 ad2 n1 ad1 | da drop
  n1 n2 + 1- 0
  do 0 i 1+ 0
     do j i - 0 n2 within i n1 < and
        if i cells ad1 + @
           j i - cells ad2 + @ * +
        then
     loop >zst+
  loop rrzs ;
\ Multiply polynomials given by arrays

: p* \ vect1 vect2 -- vect3
  >da2 >da poly*
  znip znip ;
\ Multiply polynomials

: v+ \ vect1 vect2 -- vect3
  adn2 nip adn1 nip < if zswap then
  adn2 drop locals| ad |
  zst>> cs 0
  do ad i cells + +!
  loop ;
\ Add vectors

: p+ \ vect1 vect2 -- vect3
  v+ rrzs ;
\ Add polynomials

: ps* \ vect1 n -- vect2
  locals| n |
  >da cells over + swap
  do i @ n * i ! cell +loop ;
\ Multiply polynomial with integer

: pnegate \ vect1 -- vect2
  -1 ps* ;

false [if]
: pnegate \ vect1 -- vect2
  adn1 cell- 0
  do dup i + dup @ negate swap ! cell
  +loop drop ;
\ Negate a polynomial
[then]

: p- \ vect1 vect2 -- vect3
  pnegate p+ ;
\ Subtract polynomials

: v- \ vect1 vect2 -- vect3
  pnegate v+ ;
\ Subtract vectors

: ps/ \ vect1 n -- vect2
  locals| n |
  >da cells over + swap
  do i @ n / i ! cell +loop ;
\ Divide polynomial with integer

: degree \ vect -- vect n
  zst@ cs 1- ;

\ long division

: vcutr \ vect1 n -- vect2
  degree swap 1- -
  >r zst>> cs r@ - >xs
  r> drops xs> 2* 1+ negate >>zst ;
\ vect2 is the n rightmost coefficients of vect1

: vshiftr \ vect1 -- vect2
  zst> zst> drop 2 + >zst ;
\ drop the rightmost coefficient

: getcoeff \ xvect i -- xvect n
  cells xst @ cell - swap - @ ;

: vor \ vect -- flag
  zst>> cs 1 ?do or loop ;

: ldivide \ -- q r
  zst> zst@ swap >zst
  yst> yst@ swap >yst
  /mod swap ;

: lclean \ --
  xst setdrop yst setdrop ;

: lbuild \ v n q -- v' n+1
  dup >xs 1 under+
  yst zst setcopy ps* v- ;

: lnodiv \ --
  drop 0
  ?do xsdrop loop false ;

: p/ \ v1 v2 -- v1/v2 flag
  false locals| flag |
  degree zst yst setmove
  degree zst xst setcopy
  over 1+ vcutr                           \ w
  2 + 2 under+ swap 0 -rot
  do ldivide
     if true to flag leave then
     lbuild
     vshiftr ( i getcoeff ) zswap vmerge  \ w'
  loop flag if lclean lnodiv exit then
  ldivide if lclean lnodiv exit then
  lbuild vor lclean
  ( over 0 ?do xs> loop ) nip 0= ;
\ flag is true if v2 divides v1
\ else result is irrelevant

\ auto definition of polynomial
: makepoly \ vect ad n --  name of polynomial
  cr ." : " type space
  zst> zst> .
  cs 1- 1
  do ." over * " zst> . ." + "
  loop ." * " zst> . ." + ; " ;
\ Prints the definition of a polynomial to be
\ copied and pasted

\ Integer valued polynomials

: bin*sum \ ad k -- sum
  locals| k ad |
  k 0= if 1 exit then 0 k 0
  ?do i cells ad + @
     k i choose * +
  loop ;

: polyacoeff \ ad1 n1 -- vect
  da locals| ad2 n1 ad1 |
  ad1 @ >zst+
  n1 1
  ?do ad1 n1 i polynom
     ad2 i bin*sum - >zst+
  loop ;
\ Calculate the vector (c0,...,cn) from
\ integer polynomial at ad1 n1

: polya \ ad n m -- m'
  swap -rot locals| m ad | 0 swap 0
  ?do i cells ad + @
     m i choose * +
  loop ;
\ m' is the evaluation of m with the pólya function at ad n

: coeffgcd \ vect -- n
  zst>> cs       \ CS transform set count into stack count
  multgcd ;
\ GCD of the coefficients

: fixdiv \ vect -- vect n
  >da            \ get address and count for polynomial
  polyacoeff     \ calculate Pólya's coefficients
  coeffgcd ;
\ The multiple GCD of c0,...,cn is the fixed divisor of the
\ corresponding original polynomial with integer coefficients

: divcofac \ vect -- vect'
  zdup coeffgcd ps/ ;

: iseisenstein \ vect -- vect flag
  zdup zst> 2 + zst> abs false 0 locals| p flag  an |
  >zst coeffgcd dup an ugcd 1 <>
  if zdrop drop false exit then
  primes ?dup
  if 0
     do to p flag 0=
        if an p umod 0= 0=
           >da drop @ abs p ^2 umod 0= 0= and
           to flag
        then
     loop
  then flag ;

2000 value xlim

: isirr \ vect -- vect flag
  iseisenstein if true exit then
  degree 0= if gcoeff isp exit then
  degree 1 = if zdup coeffgcd 1 = exit then
  fixdiv degree 0 0 locals| posp negp n d |
  0 sbpolyn d bs/mod drop bisprime
  if xs@ if negp 1+ to negp else posp 1+ to posp then
  then xsdrop
  xlim 1
  do i sbpolyn d bs/mod drop bisprime
     if xs@ if negp 1+ to negp else posp 1+ to posp then
     then xsdrop
     i negate sbpolyn d bs/mod drop bisprime
     if xs@ if negp 1+ to negp else posp 1+ to posp then
     then xsdrop
     posp n > negp n > or if leave then
  loop posp n > negp n > or ;

: nopsqr \ x p -- x'     p|x
  begin 2dup /mod swap 0=
  while -rot nip
  repeat drop * ;

: negate? \ |n| -- n
  2 random if negate then ;

: pickprime \ n -- p
  primes dup >r 1 max random
  pick r> drops ;

: geteis0 \ u -- vect p
  ( )
  2 - 1 max random 2 +
  dup pickprime
  tuck nopsqr negate? >zst+ ;

: x/p^n \ an p -- an'
  begin 2dup mod 0=
  while tuck / swap
  repeat drop ;

: geteisvar \ n u -- vect
  dup geteis0 locals| p u | 1- 1 max random 1+ 0
  ?do u p / random 1+ p * negate? >zst+
  loop u 1+ random 1+ p x/p^n dup 0= or
  negate? >zst+
  divcofac ;

: dupderiv \ vect -- vect vect'
  ( >da swap locals| ad | 1
  do ad i cells + @ i * loop ) ;

: deriv \ vect -- vect'
  dupderiv znip ;

\ p(x) --> p(x+d)

: mtransl \ k d ak -- vect
  locals| ak d k |
  ( k 1+ 0
  do k i choose d i ** * ak *
  loop ) ;

: zerovect \ n -- vect
  >r ( r> 0
  do 0 loop ) ;

: ptransl \ vect1 d -- vect2
  locals| d |
  >da 0 over zerovect
  do i over i cells + @ d swap
     mtransl p+
  loop drop znip ;

\ Rational roots

: q* \ a b c d -- ac/(ac,bd) bd/(ac,bd)
  rot * >r * r>         \ ac bd
  2dup abs swap abs
  ugcd tuck             \ ac gcd bd gcd
  / >r / r> ;

: q/  2swap q* ;

: q+ \ a b c d -- (ad+bc)/gcd bd/gcd
  dup 3 pick * >r      \ a b c d  r: bd
  -rot * -rot * +      \ a*d+b*c  r: bd
  dup abs r@ abs
  ugcd r> over         \ a*d+b*c gcd bd gcd
  / >r / r> ;

: q-  negate q+ ;

: qpolynom \ ad n a b -- a' b'
  locals| b a | cells over + cell-
  dup @ >r cell- r> 1 2swap
  do a b q* i @ 1 q+ -cell +loop ;

: getpospairs \ vect -- vect set
  lcoeff abs gcoeff abs divz divz
  cartprod ;

: getypair \ yset -- yset' y x
  yst> drop yst> yst> ;

: haverationalroots \ vect -- vect flag
  lcoeff 0= if true exit then
  getpospairs zst yst setmove
  begin yst@
  while ysplit
     getypair 2dup ugcd 1 =
     if >da 2over qpolynom drop 0=
        if yst setdrop 2drop true exit then
        >r negate >r
        >da r> r> qpolynom drop 0=
        if yst setdrop true exit then
     else 2drop
     then
  repeat yst> ;

: setofroots \ vect -- vect set
  lcoeff 0= if true exit then
  getpospairs
  zst yst setmove xst @
  begin yst@
  while ysplit
     getypair 2dup ugcd 1 =
     if >da 2over qpolynom drop 0=
        if ( 2dup ) zst xst setmove then
        swap negate swap
        >da 2over qpolynom drop 0=
        if ( 2dup ) zst xst setmove then
     then 2drop
  repeat yst> drop
  xst @ - cell / 2* >xst
  xst zst setmove ;

: .root \ b a --  "a/b"
  dup 0= if . drop exit then
  over abs 1 = if . drop exit then
  . 8 emit ." /" . ;

: .roots \ set --
  zst> cs 3 / 0
  do zst> drop zst> zst> .root space loop
;

: isirreducible \ vect -- vect flag
  haverationalroots degree 1 > and
  if false else isirr then ;

About Pollard rho

The Pollard rho algorithm factorize a composite number n in a time proportional to √p, where p is the smallest prime that divides n. In worst case p≈√n so the algorithm factorize proportional to ∜n. Compare with trial-and-error that factorize proportional to √n.

The algorithm is easy to implement and is built upon an the greatest common divisor gcd (which can be calculated by the algorithm of Euclide) and a polynomial P(x) of degree >1. Most common polynomial to use is P(x)=x²+1. Let Pn(x)=P(x) mod n, that is the rest when P(x) is divided by n. Pn is a function Pn:ℕ→ℤn, where ℤn={0,1,2,...,n-1}. The function Pn will act like a simple pseudo random generator in the algorithm.

Define Xi+1=Pn(Xi). Since ℤn is a finite set there are smallest i,j such that Xi=Xi+j, when the sequence start to repeat itself. Mostly the sequence Xi isn't cyclic from the start Xo, but later on at some Xi. There is at trick to find out if Xi is in the loop or not. It's like sending away a turtle and a hare on the same track at the same time. When the hare again is comming besides the turtle they must both be in the loop.

Therefore, define a second sequence Yi+1=Pn(Pn(Yi)), where Xo=Yo. When Yi=Xi, i>0, Xi and Yi are in the cycle. But for Pollard rho the real important sequences are the uncalculated sequences Vi+1=Pm(Vi) and Wi+1=Pm(Pm(Wi)) where Vo=Wo=Xo and m|n (which aren't calculated since we don't know any m|n yet). When Wi=Vi, i>0, then m|Xi-Yi and gcd(Xi-Yi,n)=m. If m=1 the result is neglected and the process go on with Xi+1 and Yi+1.

If m=n the start values Xo=Yo is said to fail and may be increamented by 1 for an other try. Now, gcd(Xi-Yi,n)=n when Xi and Yi become equal before any of the uncalculated sequences Vi and Wi become equal for some non trivial divisor m>1 of n. There could still be non trivial divisors to obtain for gcd(Xj-Yj,n) for j>i when Xi=Yi, but there is no guarranty and Xi=Yi is a natural terminal case.

An example with P(x)=x²+1 that is factorized while continuing after failure for Xo=2 is n=4294952621. An example that can't be factorized after failure with Xo=2 but with Xo=3 is n=4294939069.

Increasing Xo after termination when gcd(Xi-Yi,n)=n will always find a factor of n, because if m is a non trivial factor of n, then Xo=m immediately gives m as a factor. That is, for P(x)=x²+1.

BigZ - new top level instructions

I've reorganized and extended some top level words: 

create-set \ m n xt -- set 
filter-set \ set1 xt -- set2 
build-set \ m n xt -- set 
transform-set \ set1 xt -- set2 

This words should be used with the word :| to define and form sets. 

1 10 :| b a | a ^2 b ^2 + isprime ; create-set zdup cr zet. 
{(1,1),(2,1),(4,1),(6,1),(1,2),(3,2),(5,2),(7,2),(2,3),(8,3),(1,4),(5,4),(9,4),(2,5),(4,5),(6,5),(8,5),(1,6),(5,6),(2,7),(8,7),(3,8),(5,8),(7,8),(4,9)} ok 

:| b a | a ^2 b ^2 + ; transform-set zdup cr zet. {2,5,13,17,29,37,41,53,61,73,89,97,113} ok 

:| n | n 4 mod 3 = ; filter-set zet. 0 ok 

The word :| define a nameless word and count the number of parameters.

: bl# \ ad n -- m    count the number of spaces in the string
  over + swap 0 -rot
  do i c@ bl = -
  loop ;

variable loc#    \ the number of parameters
variable sta#    \ the number of outputs on the stack

: locals# \ --
  >in @ >r
  [char] | parse bl# loc# !
  r> >in !
  1 sta# ! ; immediate

: :| \ --
  :noname
  postpone locals#
  postpone locals| ;

The nameless words, represented by xt on stack, could be of two types:

1. Taking parameters and leaving a flag.
2. Taking parameters and leaving a non negative integer.

In the first case a set with the "dimension" stated by the parameters is the result. This works with CREATE-SET and FILTER-SET.

In the second case a set of values is the result. This works with BUILD-SET and TRANSFORM-SET.

1 10 :| x | x ^2 ; build-set cr zet.
{1,4,9,16,25,36,49,64,81} ok

Normally these nameless words leave one integer on the stack, but for TRANSFORM-SET there is an option (using 2; instead of ;) when the word leaves two non negative integers on the stack

{ 0 10 | all }

creates the set {0,1,2,3,4,5,6,7,8,9} which can be transformed to a two dimensional set

:| n | n n ^2 2; transform-set cr zet.
{(0,0),(1,1),(2,4),(3,9),(4,16),(5,25),(6,36),(7,49),(8,64),(9,81)} ok

The purpose with this words is to be able to quickly inspect conjectures and Diophantine equations. 

a and b are coprime if there exists values x and y such that 

ax+by=1

Since BigZ just deals with sets of non negative numbers we can search x and y in the equation

ax-by=1

1 10 :| b a y x | a x * b y * - 1 = ; create-set cr zet.
{(2,1,1,1),(3,2,1,1),(4,3,1,1),(5,4,1,1),(6,5,1,1),(7,6,1,1),(8,7,1,1),(9,8,1,1),(1,1,2,1),(2,3,2,1),(3,5,2,1),(4,7,2,1),(5,9,2,1),(1,2,3,1),(2,5,3,1),(3,8,3,1),(1,3,4,1),(2,7,4,1),(1,4,5,1),(2,9,5,1),(1,5,6,1),(1,6,7,1),(1,7,8,1),(1,8,9,1),(3,1,1,2),(5,2,1,2),(7,3,1,2),(9,4,1,2),(1,1,3,2),(3,4,3,2),(5,7,3,2),(1,2,5,2),(3,7,5,2),(1,3,7,2),(1,4,9,2),(4,1,1,3),(7,2,1,3),(2,1,2,3),(5,3,2,3),(8,5,2,3),(1,1,4,3),(4,5,4,3),(7,9,4,3),(2,3,5,3),(5,8,5,3),(1,2,7,3),(4,9,7,3),(2,5,8,3),(5,1,1,4),(9,2,1,4),(3,2,3,4),(7,5,3,4),(1,1,5,4),(5,6,5,4),(3,5,7,4),(1,2,9,4),(6,1,1,5),(3,1,2,5),(8,3,2,5),(2,1,3,5),(7,4,3,5),(4,3,4,5),(9,7,4,5),(1,1,6,5),(6,7,6,5),(3,4,7,5),(2,3,8,5),(4,7,9,5),(7,1,1,6),(5,4,5,6),(1,1,7,6),(7,8,7,6),(8,1,1,7),(4,1,2,7),(5,2,3,7),(2,1,4,7),(9,5,4,7),(3,2,5,7),(6,5,6,7),(1,1,8,7),(8,9,8,7),(4,5,9,7),(9,1,1,8),(3,1,3,8),(5,3,5,8),(7,6,7,8),(1,1,9,8),(5,1,2,9),(7,3,4,9),(2,1,5,9),(4,3,7,9),(8,7,8,9)} ok

The vectors in the set are of the form (x,y,a,b), that is the opposite of its appearance as parameters. (Which is logical in Forth with postfix notations and values on a stack, but still a bit awkward).

The idea is a one row code and for that purpose there is a need of short words:

: isp  isprime ;         \ n -- flag
: isq  sqr ;             \ is perfect square: n -- flag
: isqf  sqrfree ;        \ is square free: n -- flag
: isem  bigomega 2 = ;   \ is semi prime: n -- flag
: ispp  smallomega 1 = ; \ is prime power: n -- flag

: 2sqs  2sqsum ;         \ square sum: a b -- sum
: 3sqs  3sqsum ;         \ square sum: a b c -- sum
: 4sqs  4sqsum ;         \ square sum: a b c d -- sum

: cop  coprime ;         \ are coprime: m n -- flag
: div  swap mod 0= ;     \ divides: m n -- flag

: <<  \ i j k -- flag  i<j<k
  over > -rot < and ;

: <<=  \ i j k -- flag  i<=j<=k
  over >= -rot <= and ;

: z. zet. ;

: fi postpone else postpone false postpone then ; immediate 
\ a short version of ELSE FALSE THEN.

When inspecting a Diophantine equation there might be symmetries and trivial cases to weed out.

1 100 :| c b a | a b c << a b cop and if a b 2sqs c ^2 = else false then ; create-set

can be shortened to

1 100 :| c b a | a b c << a b cop and if a b 2sqs c ^2 = fi ; create-set 

cr zet.
{(3,4,5),(5,12,13),(8,15,17),(7,24,25),(20,21,29),(12,35,37),(9,40,41),(28,45,53),(11,60,61),(33,56,65),(16,63,65),(48,55,73),(36,77,85),(13,84,85),(39,80,89),(65,72,97)} ok
Experimenting with sets gives a great opportunity to find and test conjectures:

1 50 :| b a | a b 2sqsum isprime ; create-set  ok
:| b a | a b + ; transform-set cr z.
{2,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61,63,65,67,69,71,73,75,77,79,83,85,87,89,91,93,95} ok 

This suggests that all odd numbers larger than 1 can be written as a sum a+b where a²+b² is prime.

Define 

: sqeq \ m n a b -- ma²+nb²
  dup * rot * -rot
  dup * * + ;

and the set {a+b<50|a,b>0 & ma²+nb² is prime} for some different values of m and n

5 value m
7 value n

1 50 :| b a | m n a b sqeq isprime ; create-set
:| b a | a b + ; transform-set cr zet.
{5,7,11,13,17,19,23,25,29,31,35,37,41,43,47,49,53,55,59,61,65,67,71,73,77,79,83,85,89,91} ok

3 to m
4 to n

1 50 :| b a | m n a b sqeq isprime ; create-set
:| b a | a b + ; transform-set cr zet.
{2,3,4,5,6,8,9,10,11,12,13,15,16,17,18,19,20,22,23,24,25,26,27,29,30,31,32,33,34,36,37,38,39,40,41,43,44,45,46,47,48,50,51,52,53,54,55,57,58,59,60,61,62,64,65,66,67,68,69,71,72,73,74,75,76,78,79,80,81,82,83,86,87,88,89,93,94,95,96} ok

This will show a pattern which suggest that 
{a+b|a,b>0 & ma²+nb² is prime}={k>1|gcd(k,m+n)=1}

However, further tests will show that some changes are needed:

Routines to check the conjecture and print the table:

: prime_partition \ m n k -- flag
  false locals| flag k |
  k 1
  ?do 2dup k i - i sqeq isprime
     if true to flag leave then
  loop 2drop flag ;

\ 100000 value nx     Gives the same table as
1000 value nx

: maximal_exception \ m n -- mx  
  2dup + 1 locals| mx m+n n m |
  nx 3
  do m+n i coprime
     if m n i prime_partition 0=
        if i to mx then
     then
  loop mx ;

: table \ m2 n2 --
  locals| n2 | cr
  2 spaces dup 1 do i 3 .r loop 1
  do n2 1 cr i 2 .r
     do j i coprime
        if j i maximal_exception 
        else 0
        then 3 .r
     loop
  loop ;

The source code for BigZ can be loaded here:

https://github.com/Lehs/BigZ/blob/master/bigzet.txt

Karatsuba multiplication

The time for direct multiplication is proportional to n² where n is the number of figures in the multiplicands x,y. When choosing a big base B=2ⁿ and and writing

xy=(x0+B*x1)(y0+B*y1)=x0y0 + (x0*y1+x1*y0)B + x1*y1*B²

there are four multiplications of smaller numbers, also here the calculation time is proportional to n². However

(x0*y1+x1*y0)=(x0+x1)(y0+y1)-x0*y0-x1*y1

why it's enough to calculate three smaller multiplications 

x0*y0, x1*y1 and (x0+x1)(y0+y1). 

The multiplication with B are fast left shifting and if the shifting and the addition where cost free, the recursive Karatsuba multiplication would be very efficient

but unfortunately the extra math (and in Forth also some stack juggling) takes a lot of time and the method is efficient only for rather big numbers. For very big numbers, however, it's very efficient.

Here is the way I implemented it in ANS Forth:

: bcells* \ big m -- big*C^m
  cells top$ locals| n ad mb |
  ad ad mb + n move
  ad mb erase
  mb bvp @ +! ;
\ C is the number of digits in a cell

: bcells/ \ big m -- big/C^m
  cells top$ locals| n ad mb |
  ad mb + ad n move
  mb negate bvp @ +! ;

: bsplit \ w ad -- u v 
  dup nextfree < 
  if bvp @ dup @ vp+ bvp @ ! ! 
  else drop bzero
  then ;
\ A big number is split on the big stack at address ad

: btransmul \ x y -- x0 x1 y0 y1 m     B=2^bits 
  len1 len2 max cell + lcell 1+ rshift     \ m
  dup >r cells 
  >bx first over + bsplit 
  bx> first + bsplit r> ; 
\ x=x0+x1*B^m  y=y0+y1*B^m 

0x84 value karalim \ break point byte length for termination.

: b* \ x y -- xy

  len1 len2 max karalim < 

  if b* exit then

  btransmul >r                   \ x0 x1 y0 y1

  3 bpick 2 bpick recurse >bx    \ bx: x0*y0

  2 bpick 1 bpick recurse >bx    \ bx: x0*y0 x1*y1

  b+ >bx b+ bx>   recurse        \ (x0+x1)(y0+y1)

  bx b- by b- r@ bcells*         \ z1*C^m

  bx> r> 2* bcells* bx> b+ b+ <top ;

\ Karatsuba multiplication

How to use BigZ - part 3

The binomial coefficients for big integers

The number of possibilities to choose k objects from n objects soon get to big for a single cell number. The word bschoose gives a big integer result for single cell inputs.

: bschoose \ n k -- b
  bone 0
  ?do dup i - bs*
     i 1+ bs/mod drop
  loop drop ;

 ok
2000 500 bschoose cr b.
5648284895675941420424412140748481039502890353942825357221051675360331984776743417002364625179991976070068866284527555107208940603781511988000970130381311935878493235111594076219803768997324618773852975824828528735285833615310777764160933348372329757027402537319600321600269195597902747298520883357267710485334098751949232380773741897267988881873218260056305793069941805234442045890109611836653468404129012879905442075185208447514284775689056520318572740750419026192611832748925888424320  ok

This word produce big integers with single cell factors that can be analysed by the word 

sfacset \ b -- b' set

2000 1000 bschoose sfacset bdrop cr zet.

{2,5,7,11,13,17,19,23,37,41,43,53,59,67,73,79,101,103,113,127,131,149,151,167,173,179,181,211,251,257,263,269,271,277,281,283,337,347,349,353,359,367,373,379,383,389,397,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,1153,1163,1171,1181,1187,1193,1201,1213,1217,1223,1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,1297,1301,1303,1307,1319,1321,1327,1361,1367,1373,1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,1453,1459,1471,1481,1483,1487,1489,1493,1499,1511,1523,1531,1543,1549,1553,1559,1567,1571,1579,1583,1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,1667,1669,1693,1697,1699,1709,1721,1723,1733,1741,1747,1753,1759,1777,1783,1787,1789,1801,1811,1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,1901,1907,1913,1931,1933,1949,1951,1973,1979,1987,1993,1997,1999} ok

The word bsetprod calculates the big product of the singles in set

: bsetprod \ set -- b
  bone                   \ big one
  foreach                \ make ready for do-loop
  ?do zst> bs* loop ;

and can be used to calculate the radical for big integers with single cell factors:

: bsradical \ b -- b'
  sfacset bdrop 
  bsetprod ;

50 25 bschoose bsradical cr b.
1504888171878  ok

Erdős squarefree conjecture (proved 1996) states that the central binomial coefficient (2n)Cn is not squarefree if n>4. The word sqrfacset calculates the set of all factors that occurs more than once: 

: sqrfacset \ b -- set
  bdup bsradical b/
  sfacset bdrop ;

20000 10000 bschoose sqrfacset cr zet.
{2,3,7,11,23,29,41,47,53,61,71,73,79,109,127,137,139} ok

The word

: maxel \ set -- n   non e
  zst> zst@ swap >zst zdrop ;

gives the maximal element in a set of integers.

: erdprime \ n -- p
  dup 2* swap bschoose
  sqrfacset maxel ;

How to use BigZ - part 2

Gaps in increasing sequences of natural numbers

Strictly increasing sequences as

(1,2,3,...), (2,3,5,7,11,13,...) and (1,2,4,8,16,...)

can partly be represented as sorted sets in BigZ. Differences of sequences is an analogy to differentiation of functions. When defining a function gapz, to be the sorted set of all gaps between consecutive numbers of a set, also gapz become an analogy of the derivative of functions.

Apply gapz on an arithmetic sequence gives a set with a single element. If n times apply gapz on an increasing polynomial series of degree n gives a singleton set. On a sequence of exponential function with base two does nothing to the infinite set.

: gapz \ s -- s'
  0 locals| n |        \ counts the number of gaps in s'
  foreach 1+           \ prepare elements of s for the do-loop
  do zst> zst@ - >xst  \ the gap between the largest consecutive's
     n 1+ to n
  loop zst> drop       \ drop the smallest element of s
  n 2* negate >xst     \ calculate the set-count for s'
  xst zst setmove      \ move the set to zst
  set-sort reduce ;    \ sort and eliminate copies

{ 1 1000 | prime } gapz cr zet.

{1,2,4,6,8,10,12,14,18,20} ok

Partitions of a number n into distinct primes

: collincl \ s n -- s'
  0 >xst
  begin zst@
  while zsplit 
     dup >zst zfence zmerge
     set-sort reduce zfence
     xst zst setmove zmerge
     zst xst setmove
  repeat zdrop drop
  xst zst setmove
  reduce ;
\ include n in all sets in s

: xunion \ set --
  xst zst setmove union 
  zst xst setmove ;
\ Union of the top sets on the xst- and zst-stacks
\ is put on the xst-stack

: primeset \ m -- set
  pi dup 1+ 1
  ?do i pnr@ >zst 
  loop 2* negate >zst ;
\ Create the set of all primes < m+1

: memb \ s n -- flag
  false swap
  adn1 over + swap
  ?do dup i @ =
     if -1 under+ leave then cell
  +loop drop ;
\ Faster test if n is a member in the sorted number set s

For T being the set of primes:

The algorithm can be used with corrections for n=2p.

: termcase \ n -- flag
  case 2 of true endof
       3 of true endof
      11 of true endof
     dup of false endof
  endcase ; 
\ terminal cases: prime numbers without additional partitions

I have no proof that there are additional partitions for all primes greater than 11, but as far as the algorithm will go the terminal cases above are correct.

: z2@ \ set -- set n
  zst> zst@ swap >zst ;
\ read the largest element in the set

: lowlim \ set n -- set p
  0 swap adn1 over + swap
  ?do i @ under+ 2dup < 0= 
     if 2drop i @ leave then cell
  +loop ;
\ p is the smallest prime such that 2+3+5+...+p > n
  
: setsum \ set -- sum
  0 foreach ?do zst> + loop ;
\ The sum of all elements in set

: sumcorr \ s n -- s'
  locals| n |
  0 >xst
  begin zst@ 
  while zsplit zdup setsum n =
     if zfence xunion
     else zdrop
     then
  repeat zst> drop 
  xst zst setmove ;
\ Removes all partitions from s such that the sum < n

: dps \ n -- set
  dup 2 < if drop 0 >zst exit then 
  dup termcase if >zst -2 >zst -4 >zst exit then
  0 >xst
  dup primeset
  dup lowlim locals| low n |
  begin zst@ 
     if z2@ low <
        if false else true then
     else false 
     then
  while zsplit n zst> dup >r - ?dup 
     if recurse
        zst@ 
        if r> collincl n sumcorr xunion 
        else zst> drop r> drop 
        then 
     else { { r> } } xunion  
     then 
  repeat zdrop
  xst zst setmove 

  set-sort reduce ;

\ The set of partitions of n>0 into distinct primes

20 dps cr zet.

{{2,7,11},{2,5,13},{7,13},{3,17}} ok

50 dps cr zet.

{{2,7,11,13,17},{2,5,11,13,19},{7,11,13,19},{2,5,7,17,19},{3,11,17,19},{2,5,7,13,23},{3,11,13,23},{2,3,5,17,23},{3,7,17,23},{3,5,19,23},{2,3,5,11,29},{3,7,11,29},{3,5,13,29},{2,19,29},{3,5,11,31},{2,17,31},{19,31},{2,11,37},{13,37},{2,7,41},{2,5,43},{7,43},{3,47}} ok

: A000586 \ n -- 

  ." 1," 1+ 1 

  ?do i dps cardinality 0

     <# [char] , hold #s #> type 

  loop ;

\ List A000586 

100 a000586 cr
1,0,1,1,0,2,0,2,1,1,2,1,2,2,2,2,3,2,4,3,4,4,4,5,5,5,6,5,6,7,6,9,7,9,9,9,11,11,11,13,12,14,15,15,17,16,18,19,20,21,23,22,25,26,27,30,29,32,32,35,37,39,40,42,44,45,50,50,53,55,57,61,64,67,70,71,76,78,83,87,89,93,96,102,106,111,114,119,122,130,136,140,147,150,156,164,170,178,183,188,198, ok

Partitions of a number n into distinct non composites

A variant of the above.

: termcase1 \ n -- flag

  case 1 of true endof

       2 of true endof

     dup of false endof

  endcase ; 

: dps1 \ n -- set

  dup 0= if >zst exit then 

  dup termcase1 if >zst -2 >zst -4 >zst exit then

  0 >xst

  dup { 1 } primeset zmerge

  dup lowlim locals| low n |

  begin zst@ 

     if z2@ low <

        if false else true then

     else false 

     then

  while zsplit n zst> dup >r - ?dup 

     if recurse

        zst@ 

        if r> collincl n sumcorr xunion 

        else zst> drop r> drop 

        then 

     else { { r> } } xunion  

     then 

  repeat zdrop

  xst zst setmove 

  set-sort reduce ;

50 dps1 cr zet. 

{{1,3,5,11,13,17},{2,7,11,13,17},{1,2,3,5,7,13,19},{2,5,11,13,19},{7,11,13,19},{2,5,7,17,19},{1,2,11,17,19},{3,11,17,19},{1,13,17,19},{1,3,5,7,11,23},{2,5,7,13,23},{1,2,11,13,23},{3,11,13,23},{2,3,5,17,23},{1,2,7,17,23},{3,7,17,23},{1,2,5,19,23},{3,5,19,23},{1,7,19,23},{2,3,5,11,29},{1,2,7,11,29},{3,7,11,29},{1,2,5,13,29},{3,5,13,29},{1,7,13,29},{1,3,17,29},{2,19,29},{1,2,5,11,31},{3,5,11,31},{1,7,11,31},{1,2,3,13,31},{1,5,13,31},{2,17,31},{19,31},{1,2,3,7,37},{1,5,7,37},{2,11,37},{13,37},{1,3,5,41},{2,7,41},{2,5,43},{7,43},{1,2,47},{3,47}} ok

: test \ n --    n>0

  1+ 1 

  ?do i dps1 cardinality 0

     <# [char] , hold #s #> type 

  loop ;

100 cr test
1,1,2,1,2,2,2,3,2,3,3,3,4,4,4,5,5,6,7,7,8,8,9,10,10,11,11,11,13,13,15,16,16,18,18,20,22,22,24,25,26,29,30,32,33,34,37,39,41,44,45,47,51,53,57,59,61,64,67,72,76,79,82,86,89,95,100,103,108,112,118,125,131,137,141,147,154,161,170,176,182,189,198,208,217,225,233,241,252,266,276,287,297,306,320,334,348,361,371,386, ok