## Wednesday, June 30, 2021

### 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
\ 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

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
zst>> cs 0
do ad i cells + +!
loop rrzs ;

: pnegate \ vect1 -- vect2
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
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
n1 1
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
\ 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

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
zst>> cs 0
do ad i cells + +!
loop ;

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

: 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
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
k 0= if 1 exit then 0 k 0
?do i cells ad + @
k i choose * +
loop ;

: polyacoeff \ ad1 n1 -- vect
n1 1
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 ;