diff --git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet
index f8f8084..e764baf 100644
--- a/books/bookvol10.4.pamphlet
+++ b/books/bookvol10.4.pamphlet
@@ -156782,6 +156782,581 @@ UTSodetools(F, UP, L, UTS): Exports == Implementation where
 
 \end{chunk}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{package POLYVEC U32VectorPolynomialOperations}
+\begin{chunk}{TransSolvePackage.input}
+)set break resume
+)sys rm -f TransSolvePackage.output
+)spool TransSolvePackage.output
+)set message test on
+)set message auto off
+)clear all
+)spool
+\end{chunk}
+
+\begin{chunk}{U32VectorPolynomialOperations.help}
+====================================================================
+U32VectorPolynomialOperations examples
+====================================================================
+
+This is a low-level package which implements operations
+on vectors treated as univariate modular polynomials.  Most
+operations takes modulus as parameter.  Modulus is machine
+sized prime which should be small enough to avoid overflow
+in intermediate calculations.
+
+See Also:
+o )show U32VectorPolynomialOperations
+
+\end{chunk}
+
+\pagehead{U32VectorPolynomialOperations}{POLYVEC}
+\pagepic{ps/v104u32vectorpolynomialoperations.eps}{POLYVEC}{1.00}
+
+{\bf Exports:}\\
+\begin{tabular}{lll}
+\cross{POLYVEC}{copyfirst} &
+\cross{POLYVEC}{copyslice} &
+\cross{POLYVEC}{degree} \\
+\cross{POLYVEC}{differentiate} &
+\cross{POLYVEC}{differentiate} &
+\cross{POLYVEC}{divide!} \\
+\cross{POLYVEC}{evalat} &
+\cross{POLYVEC}{extendedgcd} &
+\cross{POLYVEC}{gcd} \\
+\cross{POLYVEC}{gcd} &
+\cross{POLYVEC}{lcm} &
+\cross{POLYVEC}{mul} \\
+\cross{POLYVEC}{mulbybinomial} &
+\cross{POLYVEC}{mulbybinomial} &
+\cross{POLYVEC}{mulbyscalar} \\
+\cross{POLYVEC}{pow} &
+\cross{POLYVEC}{remainder!} &
+\cross{POLYVEC}{resultant} \\
+\cross{POLYVEC}{tomodpa} &
+\cross{POLYVEC}{truncatedmuladd} &
+\cross{POLYVEC}{truncatedmultiplication} \\
+\cross{POLYVEC}{vectoraddmul} &
+\cross{POLYVEC}{vectorcombination} &
+\end{tabular}
+
+\begin{chunk}{package POLYVEC U32VectorPolynomialOperations}
+)abbrev package POLYVEC U32VectorPolynomialOperations
+++ Description: This is a low-level package which implements operations
+++  on vectors treated as univariate modular polynomials.  Most
+++  operations takes modulus as parameter.  Modulus is machine
+++  sized prime which should be small enough to avoid overflow
+++  in intermediate calculations.
+U32VectorPolynomialOperations() : Export == Implementation where
+    PA ==> U32Vector
+    Export ==> with
+        copy_first : (PA, PA, Integer) -> Void
+          ++ copy_first(v1, v2, n) copies first n elements
+          ++ of v2 into n first positions in v1.
+        copy_slice : (PA, PA, Integer, Integer) -> Void
+          ++ copy_first(v1, v2, m, n) copies the slice of v2 starting
+          ++ at m elements and having n elements into corresponding
+          ++ positions in v1.
+        eval_at : (PA, Integer, Integer, Integer) -> Integer
+          ++ eval_at(v, deg, pt, p) treats v as coefficients of
+          ++ polynomial of degree deg and evaluates the
+          ++ polynomial at point pt modulo p
+          ++
+          ++X a:=new(3,1)$U32VEC
+          ++X a.1:=2
+          ++X eval_at(a,2,3,1024)
+          ++X eval_at(a,2,2,8)
+          ++X eval_at(a,2,3,10)
+        vector_add_mul : (PA, PA, Integer, Integer, Integer, Integer) _
+            -> Void
+          ++ vector_add_mul(v1, v2, m, n, c, p) sets v1(m), ...,
+          ++ v1(n) to corresponding extries in v1 + c*v2
+          ++ modulo p.
+        mul_by_binomial : (PA, Integer, Integer) -> Void
+          ++ mul_by_binomial(v,  pt, p) treats v a polynomial
+          ++ and multiplies in place this polynomial by binomial (x + pt).
+          ++ Highest coefficient of product is ignored.
+        mul_by_binomial : (PA, Integer, Integer, Integer) -> Void
+          ++ mul_by_binomial(v, deg, pt, p) treats v as
+          ++ coefficients of polynomial of degree deg and
+          ++ multiplies in place this polynomial by binomial (x + pt).
+          ++ Highest coefficient of product is ignored.
+        mul_by_scalar : (PA, Integer, Integer, Integer) -> Void
+          ++ mul_by_scalar(v, deg, c, p) treats v as
+          ++ coefficients of polynomial of degree deg and
+          ++ multiplies in place this polynomial by scalar c
+        mul : (PA, PA, Integer) -> PA
+          ++ Polynomial multiplication.
+        truncated_multiplication : (PA, PA, Integer, Integer) -> PA
+          ++ truncated_multiplication(x, y, d, p) computes
+          ++ x*y truncated after degree d
+        truncated_mul_add : (PA, PA, PA, Integer, Integer) -> Void
+          ++ truncated_mul_add(x, y, z, d, p) adds to z
+          ++ the produce x*y truncated after degree d
+        pow : (PA, PositiveInteger, NonNegativeInteger, Integer) -> PA
+          ++ pow(u, n, d, p) returns u^n truncated after degree d, except if
+          ++ n=1, in which case u itself is returned
+        differentiate : (PA, Integer) -> PA
+          ++ Polynomial differentiation.
+        differentiate : (PA, NonNegativeInteger, Integer) -> PA
+          ++ Polynomial differentiation.
+        divide! : (PA, PA, PA, Integer) -> Void
+          ++ Polynomial division.
+        remainder! : (PA, PA, Integer) -> Void
+          ++ Polynomial remainder
+        vector_combination : (PA, Integer, PA, Integer, _
+                               Integer, Integer, p : Integer) -> Void
+          ++ vector_combination(v1, c1, v2, c2, n, delta, p) replaces
+          ++ first n + 1 entires of v1 by corresponding entries of
+          ++ c1*v1+c2*x^delta*v2 mod p.
+        to_mod_pa : (SparseUnivariatePolynomial Integer, Integer) -> PA
+          ++ to_mod_pa(s, p) reduces coefficients of polynomial
+          ++ s modulo prime p and converts the result to vector
+        gcd : (PA, PA, Integer) -> PA
+          ++ gcd(v1, v2, p) computes monic gcd of v1 and v2 modulo p.
+        gcd : (PrimitiveArray PA, Integer, Integer, Integer) -> PA
+          ++ gcd(a, lo, hi, p) computes gcd of elements
+          ++ a(lo), a(lo+1), ..., a(hi).
+        lcm : (PrimitiveArray PA, Integer, Integer, Integer) -> PA
+          ++ lcm(a, lo, hi, p) computes lcm of elements
+          ++ a(lo), a(lo+1), ..., a(hi).
+        degree : PA -> Integer
+          ++ degree(v) is degree of v treated as polynomial
+        extended_gcd : (PA, PA, Integer) -> List(PA)
+          ++ extended_gcd(v1, v2, p) gives [g, c1, c2] such
+          ++ that g is \spad{gcd(v1, v2, p)}, \spad{g = c1*v1 + c2*v2}
+          ++ and degree(c1) < max(degree(v2) - degree(g), 0) and
+          ++ degree(c2) < max(degree(v1) - degree(g), 1)
+        resultant : (PA, PA, Integer) -> Integer
+          ++ resultant(v1, v2, p) computes resultant of v1 and v2
+          ++ modulo p.
+
+    Implementation ==> add
+
+        Qmuladdmod ==> QSMULADDMOD6432$Lisp
+        Qmuladd    ==> QSMULADD6432$Lisp
+        Qmul       ==> QSMULMOD32$Lisp
+        Qdot2      ==> QSDOT2MOD6432$Lisp
+        Qrem       ==> QSMOD6432$Lisp
+        modInverse ==> invmod
+
+        copy_first(np : PA, op : PA, n : Integer) : Void ==
+            ns := n pretend SingleInteger
+            for j in 0..(ns - 1) repeat
+                np(j) := op(j)
+
+        copy_slice(np : PA, op : PA, m : Integer, _
+                    n : Integer) : Void ==
+            ms := m pretend SingleInteger
+            ns := n pretend SingleInteger
+            for j in ms..(ms + ns - 1) repeat
+                np(j) := op(j)
+
+        eval_at(v : PA, deg : Integer, pt : Integer, _
+               p : Integer) : Integer ==
+            i : SingleInteger := deg::SingleInteger
+            res : Integer := 0
+            while not(i < 0) repeat
+                res := Qmuladdmod(pt, res, v(i), p)
+                i := i - 1
+            res
+
+        to_mod_pa(s : SparseUnivariatePolynomial Integer, p : Integer) : PA ==
+            zero?(s) => new(1, 0)$PA
+            n0 := degree(s) pretend SingleInteger
+            ncoeffs := new((n0+1) pretend NonNegativeInteger, 0)$PA
+            while not(zero?(s)) repeat
+                n := degree(s)
+                ncoeffs(n) := positiveRemainder(leadingCoefficient(s), p)
+                s := reductum(s)
+            ncoeffs
+
+        vector_add_mul(v1 : PA, v2 : PA, m : Integer, n : Integer, _
+                         c : Integer, p : Integer) : Void ==
+            ms := m pretend SingleInteger
+            ns := n pretend SingleInteger
+            for i in ms..ns repeat
+                v1(i) := Qmuladdmod(c, v2(i), v1(i), p)
+
+        mul_by_binomial(v : PA, n : Integer, pt : Integer, _
+                          p : Integer) : Void ==
+            prev_coeff : Integer := 0
+            ns := n pretend SingleInteger
+            for i in 0..(ns - 1) repeat
+                pp := v(i)
+                v(i) := Qmuladdmod(pt, pp, prev_coeff, p)
+                prev_coeff := pp
+
+        mul_by_binomial(v : PA, pt : Integer, _
+                          p : Integer) : Void ==
+            mul_by_binomial(v, #v, pt, p)
+
+        mul_by_scalar(v : PA, n : Integer, c : Integer, _
+                        p : Integer) : Void ==
+            ns := n pretend SingleInteger
+            for i in 0..ns repeat
+                v(i) := Qmul(c, v(i), p)
+
+        degree(v : PA) : Integer ==
+            n := #v
+            for i in (n - 1)..0 by -1 repeat
+                not(v(i) = 0) => return i
+            -1
+
+        vector_combination(v1 : PA, c1 : Integer, _
+                            v2 : PA, c2 : Integer, _
+                            n : Integer, delta : Integer, _
+                            p : Integer) : Void ==
+            ns := n pretend SingleInteger
+            ds := delta pretend SingleInteger
+            if not(c1 = 1) then
+                ns + 1 < ds =>
+                    for i in 0..ns repeat
+                        v1(i) := Qmul(v1(i), c1, p)
+                for i in 0..(ds - 1) repeat
+                    v1(i) := Qmul(v1(i), c1, p)
+                for i in ds..ns repeat
+                    v1(i) := Qdot2(v1(i), c1, v2(i - ds), c2, p)
+            else
+                for i in ds..ns repeat
+                    v1(i) := Qmuladdmod(c2, v2(i - ds), v1(i), p)
+
+        divide!(r0 : PA, r1 : PA, res : PA, p: Integer) : Void ==
+            dr0 := degree(r0) pretend SingleInteger
+            dr1 := degree(r1) pretend SingleInteger
+            c0 := r1(dr1)
+            c0 := modInverse(c0, p)
+            while not(dr0 < dr1) repeat
+                delta := dr0 - dr1
+                c1 := Qmul(c0, r0(dr0), p)
+                res(delta) := c1
+                c1 := p - c1
+                r0(dr0) := 0
+                dr0 := dr0 - 1
+                if dr0 < 0 then break
+                vector_combination(r0, 1, r1, c1, dr0, delta, p)
+                while r0(dr0) = 0 repeat
+                    dr0 := dr0 - 1
+                    if dr0 < 0 then break
+
+        remainder!(r0 : PA, r1 : PA, p: Integer) : Void ==
+            dr0 := degree(r0) pretend SingleInteger
+            dr1 := degree(r1) pretend SingleInteger
+            c0 := r1(dr1)
+            c0 := modInverse(c0, p)
+            while not(dr0 < dr1) repeat
+                delta := dr0 - dr1
+                c1 := Qmul(c0, r0(dr0), p)
+                c1 := p - c1
+                r0(dr0) := 0
+                dr0 := dr0 - 1
+                if dr0 < 0 then break
+                vector_combination(r0, 1, r1, c1, dr0, delta, p)
+                while r0(dr0) = 0 repeat
+                    dr0 := dr0 - 1
+                    if dr0 < 0 then break
+
+        gcd(x : PA, y : PA, p : Integer) : PA ==
+            dr0 := degree(y) pretend SingleInteger
+            dr1 : SingleInteger
+            if dr0 < 0 then
+                tmpp := x
+                x := y
+                y := tmpp
+                dr1 := dr0
+                dr0 := degree(y) pretend SingleInteger
+            else
+                dr1 := degree(x) pretend SingleInteger
+            dr0 < 0 => return new(1, 0)$PA
+            r0 := new((dr0 + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(r0, y, dr0 + 1)
+            dr1 < 0 =>
+                c := r0(dr0)
+                c := modInverse(c, p)
+                mul_by_scalar(r0, dr0, c, p)
+                return r0
+            r1 := new((dr1 + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(r1, x, dr1 + 1)
+            while 0 < dr1 repeat
+                while not(dr0 < dr1) repeat
+                    delta := dr0 - dr1
+                    c1 := sub_SI(p, r0(dr0))$Lisp
+                    c0 := r1(dr1)
+                    if c0 ~= 1 and delta > 30 then
+                        c0 :=  modInverse(c0, p)
+                        mul_by_scalar(r1, dr1, c0, p)
+                        c0 := 1
+                    r0(dr0) := 0
+                    dr0 := dr0 - 1
+                    vector_combination(r0, c0, r1, c1, dr0, delta, p)
+                    while r0(dr0) = 0 repeat
+                        dr0 := dr0 - 1
+                        if dr0 < 0 then break
+                tmpp := r0
+                tmp := dr0
+                r0 := r1
+                dr0 := dr1
+                r1 := tmpp
+                dr1 := tmp
+            not(dr1 < 0) =>
+                r1(0) := 1
+                return r1
+            c := r0(dr0)
+            c := modInverse(c, p)
+            mul_by_scalar(r0, dr0, c, p)
+            r0
+
+        gcd(a : PrimitiveArray PA, lo : Integer, hi: Integer, p: Integer) _
+              : PA ==
+            res := a(lo)
+            for i in (lo + 1)..hi repeat
+                res := gcd(a(i), res, p)
+            res
+
+        lcm2(v1 : PA, v2 : PA, p : Integer) : PA ==
+            pp := gcd(v1, v2, p)
+            dv2 := degree(v2)
+            dpp := degree(pp)
+            dv2 = dpp =>
+                v1
+            dpp = 0 => mul(v1, v2, p)
+            tmp1 := new((dv2 + 1) pretend NonNegativeInteger, 0)$PA
+            tmp2 := new((dv2 - dpp + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(tmp1, v2, dv2 + 1)
+            divide!(tmp1, pp, tmp2, p)
+            mul(v1, tmp2, p)
+
+        lcm(a : PrimitiveArray PA, lo : Integer, hi: Integer, p: Integer) _
+              : PA ==
+            res := a(lo)
+            for i in (lo + 1)..hi repeat
+                res := lcm2(a(i), res, p)
+            res
+
+        inner_mul : (PA, PA, PA,  SingleInteger, SingleInteger, _
+                      SingleInteger, Integer) -> Void
+
+        mul(x : PA, y : PA, p : Integer) : PA ==
+            xdeg := degree(x) pretend SingleInteger
+            ydeg := degree(y) pretend SingleInteger
+            if xdeg > ydeg then
+                tmpp := x
+                tmp := xdeg
+                x := y
+                xdeg := ydeg
+                y := tmpp
+                ydeg := tmp
+            xcoeffs := x
+            ycoeffs := y
+            xdeg < 0 => x
+            xdeg = 0 and xcoeffs(0) = 1 => copy(y)
+            zdeg : SingleInteger := xdeg + ydeg
+            zdeg0 := ((zdeg + 1)::Integer) pretend NonNegativeInteger
+            zcoeffs := new(zdeg0, 0)$PA
+            inner_mul(xcoeffs, ycoeffs, zcoeffs, xdeg, ydeg, zdeg, p)
+            zcoeffs
+
+        inner_mul(x, y, z, xdeg, ydeg, zdeg, p) ==
+            if ydeg < xdeg then
+                tmpp := x
+                tmp := xdeg
+                x := y
+                xdeg := ydeg
+                y := tmpp
+                ydeg := tmp
+            xdeg :=
+                zdeg < xdeg => zdeg
+                xdeg
+            ydeg :=
+                zdeg < ydeg => zdeg
+                ydeg
+            ss : Integer
+            i : SingleInteger
+            j : SingleInteger
+            for i in 0..xdeg repeat
+                ss := z(i)
+                for j in 0..i repeat
+                    ss := Qmuladd(x(i - j), y(j), ss)
+                z(i) := Qrem(ss, p)
+            for i in (xdeg+1)..ydeg repeat
+                ss := z(i)
+                for j in 0..xdeg repeat
+                    ss := Qmuladd(x(j), y(i-j), ss)
+                z(i) := Qrem(ss, p)
+            for i in (ydeg+1)..zdeg repeat
+                ss := z(i)
+                for j in (i-xdeg)..ydeg repeat
+                    ss := Qmuladd(x(i - j), y(j), ss)
+                z(i) := Qrem(ss, p)
+
+        truncated_mul_add(x, y, z, m, p) ==
+            xdeg := (#x - 1) pretend SingleInteger
+            ydeg := (#y - 1) pretend SingleInteger
+            inner_mul(x, y, z, xdeg, ydeg, m pretend SingleInteger, p)
+
+        truncated_multiplication(x, y, m, p) ==
+            xdeg := (#x - 1) pretend SingleInteger
+            ydeg := (#y - 1) pretend SingleInteger
+            z := new((m pretend SingleInteger + 1)
+                        pretend NonNegativeInteger, 0)$PA
+            inner_mul(x, y, z, xdeg, ydeg, m pretend SingleInteger, p)
+            z
+
+        pow(x : PA, n : PositiveInteger, d: NonNegativeInteger, _
+            p : Integer) : PA ==
+            one? n => x
+            odd?(n)$Integer =>
+                truncated_multiplication(x,
+                    pow(truncated_multiplication(x, x, d, p),
+                        shift(n,-1) pretend PositiveInteger,
+                        d,
+                        p),
+                    d,
+                    p)
+            pow(truncated_multiplication(x, x, d, p),
+                shift(n,-1) pretend PositiveInteger,
+                d,
+                p)
+
+        differentiate(x: PA, p: Integer): PA ==
+            d := #x - 1
+            if zero? d then empty()$PA
+            else
+                r := new(d::NonNegativeInteger, 0)$PA
+                for i in 0..d-1 repeat
+                    i1 := i+1
+                    r.i := Qmul(i1, x.i1, p)
+                r
+
+        differentiate(x: PA, n: NonNegativeInteger, p: Integer): PA ==
+            zero? n => x
+            d := #x - 1
+            if d < n then empty()$PA
+            else
+                r := new((d-n+1) pretend NonNegativeInteger, 0)$PA
+                for i in n..d repeat
+                    j := i-n
+                    f := j+1
+                    for k in j+2..i repeat f := Qmul(f, k, p)
+                    r.(j pretend NonNegativeInteger) := Qmul(f, x.i, p)
+                r
+
+        extended_gcd(x : PA, y : PA, p : Integer) : List(PA) ==
+            dr0 := degree(x) pretend SingleInteger
+            dr1 : SingleInteger
+            swapped : Boolean := false
+            t0 : PA
+            if dr0 < 0 then
+                (x, y) := (y, x)
+                dr1 := dr0
+                dr0 := degree(x) pretend SingleInteger
+                swapped := true
+            else
+                dr1 := degree(y) pretend SingleInteger
+            dr1 < 0 =>
+                dr0 < 0 =>
+                    return [new(1, 0)$PA, new(1, 0)$PA, new(1, 1)$PA]
+                r0 := new((dr0 + 1) pretend NonNegativeInteger, 0)$PA
+                copy_first(r0, x, dr0 + 1)
+                c := r0(dr0)
+                c := modInverse(c, p)
+                mul_by_scalar(r0, dr0, c, p)
+                t0 := new(1, c)$PA
+                if swapped then
+                    return [r0, new(1, 0)$PA, t0]
+                else
+                    return [r0, t0, new(1, 0)$PA]
+            swapped => error "impossible"
+            dt := (dr0 > 0 => dr0 - 1 ; 0)
+            ds := (dr1 > 0 => dr1 - 1 ; 0)
+            -- invariant: r0 = s0*x + t0*y, r1 = s1*x + t1*y
+            r0 := new((dr0 + 1) pretend NonNegativeInteger, 0)$PA
+            t0 := new((dt + 1) pretend NonNegativeInteger, 0)$PA
+            s0 := new((ds + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(r0, x, dr0 + 1)
+            s0(0) := 1
+            r1 := new((dr1 + 1) pretend NonNegativeInteger, 0)$PA
+            t1 := new((dt + 1) pretend NonNegativeInteger, 0)$PA
+            s1 := new((ds + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(r1, y, dr1 + 1)
+            t1(0) := 1
+            while dr1 > 0 repeat
+                while dr0 >= dr1 repeat
+                    delta := dr0 - dr1
+                    c1 := sub_SI(p, r0(dr0))$Lisp
+                    c0 := r1(dr1)
+                    if c0 ~= 1 and delta > 30 then
+                        c0 :=  modInverse(c0, p)
+                        mul_by_scalar(r1, dr1, c0, p)
+                        mul_by_scalar(t1, dt, c0, p)
+                        mul_by_scalar(s1, ds, c0, p)
+                        c0 := 1
+                    r0(dr0) := 0
+                    dr0 := dr0 - 1
+                    vector_combination(r0, c0, r1, c1, dr0, delta, p)
+                    vector_combination(t0, c0, t1, c1, dt, delta, p)
+                    vector_combination(s0, c0, s1, c1, ds, delta, p)
+                    while r0(dr0) = 0 repeat
+                        dr0 := dr0 - 1
+                        if dr0 < 0 then break
+                (r0, r1) := (r1, r0)
+                (dr0, dr1) := (dr1, dr0)
+                (s0, s1) := (s1, s0)
+                (t0, t1) := (t1, t0)
+            dr1 >= 0 =>
+                c := r1(0)
+                c := modInverse(c, p)
+                r1(0) := 1
+                mul_by_scalar(s1, ds, c, p)
+                mul_by_scalar(t1, dt, c, p)
+                return [r1, s1, t1]
+            c := r0(dr0)
+            c := modInverse(c, p)
+            mul_by_scalar(r0, dr0, c, p)
+            mul_by_scalar(s0, ds, c, p)
+            mul_by_scalar(t0, dt, c, p)
+            [r0, s0, t0]
+
+        resultant(x : PA, y : PA, p : Integer) : Integer ==
+            dr0 := degree(x) pretend SingleInteger
+            dr0 < 0 => 0
+            dr1 := degree(y) pretend  SingleInteger
+            dr1 < 0 => 0
+            r0 := new((dr0 + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(r0, x, dr0 + 1)
+            r1 := new((dr1 + 1) pretend NonNegativeInteger, 0)$PA
+            copy_first(r1, y, dr1 + 1)
+            res : SingleInteger := 1
+            repeat
+                dr0 < dr1 =>
+                    (r0, r1) := (r1, r0)
+                    (dr0, dr1) := (dr1, dr0)
+                c0 := r1(dr1)
+                dr1 = 0 =>
+                    while 0 < dr0 repeat
+                        res := Qmul(res, c0, p)
+                        dr0 := dr0 - 1
+                    return res
+                delta := dr0 - dr1
+                c1 := sub_SI(p, r0(dr0))$Lisp
+                if c0 ~= 1 then
+                    c1 :=  Qmul(c1, modInverse(c0, p), p)
+                r0(dr0) := 0
+                dr0 := dr0 - 1
+                vector_combination(r0, 1, r1, c1, dr0, delta, p)
+                res := Qmul(res, c0, p)
+                while r0(dr0) = 0 repeat
+                    dr0 := dr0 - 1
+                    dr0 < 0 => return 0
+                    res := Qmul(res, c0, p)
+
+\end{chunk}
+\begin{chunk}{POLYVEC.dotabb}
+"POLYVEC" [color="#FF4488",href="bookvol10.4.pdf#nameddest=POLYVEC"]
+"A1AGG" [color="#4488FF",href="bookvol10.2.pdf#nameddest=A1AGG"]
+"POLYVEC" -> "A1AGG"
+
+\end{chunk}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \chapter{Chapter V}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \section{package VECTOR2 VectorFunctions2}
diff --git a/books/bookvol5.pamphlet b/books/bookvol5.pamphlet
index b2c79c2..fb3483a 100644
--- a/books/bookvol5.pamphlet
+++ b/books/bookvol5.pamphlet
@@ -24650,6 +24650,7 @@ otherwise the new algebra won't be loaded by the interpreter when needed.
    (|U8Vector| . U8VEC)
    (|U16Vector| . U16VEC)
    (|U32Vector| . U32VEC)
+   (|U32VectorPolynomialOperations| . POLYVEC)
    (|Vector| . VECTOR)
    (|VectorFunctions2| . VECTOR2)
    (|ViewDefaultsPackage| . VIEWDEF)
@@ -39502,6 +39503,66 @@ Given a form, $u$, we try to recover the input line that created it.
 
 \end{chunk}
 
+\section{U32VectorPolynomialOperations}
+
+\defmacro{qsMulAdd6432}
+\begin{chunk}{defmacro qsMulAdd6432}
+(defmacro qsMulAdd6432 (x y z)
+  `(the (unsigned-byte 64)
+     (+ (the (unsigned-byte 64)
+          (* (the (unsigned-byte 32) ,x)
+             (the (unsigned-byte 32) ,y)))
+        (the (unsigned-byte 64) ,z))))
+
+\end{chunk}
+
+\defmacro{qsMulMod32}
+\begin{chunk}{defmacro qsMulMod32}
+(defmacro qsMulMod32 (x y)
+  `(the (unsigned-byte 64)
+     (* (the (unsigned-byte 32) ,x)
+        (the (unsigned-byte 32) ,y))))
+
+\end{chunk}
+
+\defmacro{qsMod6432}
+\begin{chunk}{defmacro qsMod6432}
+(defmacro qsMod6432 (x p)
+  `(the (unsigned-byte 32)
+     (rem (the (unsigned-byte 64) ,x) (the (unsigned-byte 32) ,p))))
+
+\end{chunk}
+
+\defmacro{qsMulAddMod6432}
+\begin{chunk}{defmacro qsMulAddMod6432}
+(defmacro qsMulAddMod6432 (x y z p)
+  `(qsMod6432 (qsMulAdd6432 ,x ,y ,z) ,p))
+
+\end{chunk}
+
+\defmacro{qsMul6432}
+\begin{chunk}{defmacro qsMul6432}
+(defmacro qsMul6432 (x y)
+  `(the (unsigned-byte 64)
+     (* (the (unsigned-byte 32) ,x)
+        (the (unsigned-byte 32) ,y))))
+
+\end{chunk}
+
+\defmacro{qsDot26432}
+\begin{chunk}{defmacro qsDot26432}
+(defmacro qsDot26432 (a1 b1 a2 b2)
+  `(qsMulAdd6432 ,a1 ,b1 (qsMul6432 ,a2 ,b2)))
+
+\end{chunk}
+
+\defmacro{qsDot2Mod6432}
+\begin{chunk}{defmacro qsDot2Mod6432}
+(defmacro qsDot2Mod6432 (a1 b1 a2 b2 p)
+  `(qsMod6432 (qsDot26432 ,a1 ,b1 ,a2 ,b2) ,p))
+
+\end{chunk}
+
 \section{DirectProduct}
 \defun{vec2list}{vec2list}
 \begin{chunk}{defun vec2list}
@@ -44071,6 +44132,13 @@ This needs to work off the internal exposure list, not the file.
 \getchunk{defmacro makeMatrix1U16}
 \getchunk{defmacro makeMatrixU32}
 \getchunk{defmacro makeMatrix1U32}
+\getchunk{defmacro qsDot26432}
+\getchunk{defmacro qsDot2Mod6432}
+\getchunk{defmacro qsMod6432}
+\getchunk{defmacro qsMulAdd6432}
+\getchunk{defmacro qsMulAddMod6432}
+\getchunk{defmacro qsMul6432}
+\getchunk{defmacro qsMulMod32}
 \getchunk{defmacro qvlenU8}
 \getchunk{defmacro qvlenU16}
 \getchunk{defmacro qvlenU32}
diff --git a/books/ps/v104u32vectorpolynomialoperations.eps b/books/ps/v104u32vectorpolynomialoperations.eps
new file mode 100644
index 0000000..2df0fa9
--- /dev/null
+++ b/books/ps/v104u32vectorpolynomialoperations.eps
@@ -0,0 +1,278 @@
+%!PS-Adobe-3.0 EPSF-3.0
+%%Creator: graphviz version 2.26.3 (20100126.1600)
+%%Title: pic
+%%Pages: 1
+%%BoundingBox: 36 36 124 152
+%%EndComments
+save
+%%BeginProlog
+/DotDict 200 dict def
+DotDict begin
+
+/setupLatin1 {
+mark
+/EncodingVector 256 array def
+ EncodingVector 0
+
+ISOLatin1Encoding 0 255 getinterval putinterval
+EncodingVector 45 /hyphen put
+
+% Set up ISO Latin 1 character encoding
+/starnetISO {
+        dup dup findfont dup length dict begin
+        { 1 index /FID ne { def }{ pop pop } ifelse
+        } forall
+        /Encoding EncodingVector def
+        currentdict end definefont
+} def
+/Times-Roman starnetISO def
+/Times-Italic starnetISO def
+/Times-Bold starnetISO def
+/Times-BoldItalic starnetISO def
+/Helvetica starnetISO def
+/Helvetica-Oblique starnetISO def
+/Helvetica-Bold starnetISO def
+/Helvetica-BoldOblique starnetISO def
+/Courier starnetISO def
+/Courier-Oblique starnetISO def
+/Courier-Bold starnetISO def
+/Courier-BoldOblique starnetISO def
+cleartomark
+} bind def
+
+%%BeginResource: procset graphviz 0 0
+/coord-font-family /Times-Roman def
+/default-font-family /Times-Roman def
+/coordfont coord-font-family findfont 8 scalefont def
+
+/InvScaleFactor 1.0 def
+/set_scale {
+       dup 1 exch div /InvScaleFactor exch def
+       scale
+} bind def
+
+% styles
+/solid { [] 0 setdash } bind def
+/dashed { [9 InvScaleFactor mul dup ] 0 setdash } bind def
+/dotted { [1 InvScaleFactor mul 6 InvScaleFactor mul] 0 setdash } bind def
+/invis {/fill {newpath} def /stroke {newpath} def /show {pop newpath} def} bind def
+/bold { 2 setlinewidth } bind def
+/filled { } bind def
+/unfilled { } bind def
+/rounded { } bind def
+/diagonals { } bind def
+
+% hooks for setting color 
+/nodecolor { sethsbcolor } bind def
+/edgecolor { sethsbcolor } bind def
+/graphcolor { sethsbcolor } bind def
+/nopcolor {pop pop pop} bind def
+
+/beginpage {	% i j npages
+	/npages exch def
+	/j exch def
+	/i exch def
+	/str 10 string def
+	npages 1 gt {
+		gsave
+			coordfont setfont
+			0 0 moveto
+			(\() show i str cvs show (,) show j str cvs show (\)) show
+		grestore
+	} if
+} bind def
+
+/set_font {
+	findfont exch
+	scalefont setfont
+} def
+
+% draw text fitted to its expected width
+/alignedtext {			% width text
+	/text exch def
+	/width exch def
+	gsave
+		width 0 gt {
+			[] 0 setdash
+			text stringwidth pop width exch sub text length div 0 text ashow
+		} if
+	grestore
+} def
+
+/boxprim {				% xcorner ycorner xsize ysize
+		4 2 roll
+		moveto
+		2 copy
+		exch 0 rlineto
+		0 exch rlineto
+		pop neg 0 rlineto
+		closepath
+} bind def
+
+/ellipse_path {
+	/ry exch def
+	/rx exch def
+	/y exch def
+	/x exch def
+	matrix currentmatrix
+	newpath
+	x y translate
+	rx ry scale
+	0 0 1 0 360 arc
+	setmatrix
+} bind def
+
+/endpage { showpage } bind def
+/showpage { } def
+
+/layercolorseq
+	[	% layer color sequence - darkest to lightest
+		[0 0 0]
+		[.2 .8 .8]
+		[.4 .8 .8]
+		[.6 .8 .8]
+		[.8 .8 .8]
+	]
+def
+
+/layerlen layercolorseq length def
+
+/setlayer {/maxlayer exch def /curlayer exch def
+	layercolorseq curlayer 1 sub layerlen mod get
+	aload pop sethsbcolor
+	/nodecolor {nopcolor} def
+	/edgecolor {nopcolor} def
+	/graphcolor {nopcolor} def
+} bind def
+
+/onlayer { curlayer ne {invis} if } def
+
+/onlayers {
+	/myupper exch def
+	/mylower exch def
+	curlayer mylower lt
+	curlayer myupper gt
+	or
+	{invis} if
+} def
+
+/curlayer 0 def
+
+%%EndResource
+%%EndProlog
+%%BeginSetup
+14 default-font-family set_font
+1 setmiterlimit
+% /arrowlength 10 def
+% /arrowwidth 5 def
+
+% make sure pdfmark is harmless for PS-interpreters other than Distiller
+/pdfmark where {pop} {userdict /pdfmark /cleartomark load put} ifelse
+% make '<<' and '>>' safe on PS Level 1 devices
+/languagelevel where {pop languagelevel}{1} ifelse
+2 lt {
+    userdict (<<) cvn ([) cvn load put
+    userdict (>>) cvn ([) cvn load put
+} if
+
+%%EndSetup
+setupLatin1
+%%Page: 1 1
+%%PageBoundingBox: 36 36 124 152
+%%PageOrientation: Portrait
+0 0 1 beginpage
+gsave
+36 36 88 116 boxprim clip newpath
+1 1 set_scale 0 rotate 40 41 translate
+0.16355 0.45339 0.92549 graphcolor
+newpath -4 -5 moveto
+-4 112 lineto
+85 112 lineto
+85 -5 lineto
+closepath fill
+1 setlinewidth
+0.16355 0.45339 0.92549 graphcolor
+newpath -4 -5 moveto
+-4 112 lineto
+85 112 lineto
+85 -5 lineto
+closepath stroke
+% POLYVEC
+gsave
+[ /Rect [ 0 72 80 108 ]
+  /Border [ 0 0 0 ]
+  /Action << /Subtype /URI /URI (bookvol10.4.pdf#nameddest=POLYVEC) >>
+  /Subtype /Link
+/ANN pdfmark
+0.93939 0.73333 1 nodecolor
+newpath 80 108 moveto
+0 108 lineto
+0 72 lineto
+80 72 lineto
+closepath fill
+1 setlinewidth
+filled
+0.93939 0.73333 1 nodecolor
+newpath 80 108 moveto
+0 108 lineto
+0 72 lineto
+80 72 lineto
+closepath stroke
+0 0 0 nodecolor
+14 /Times-Roman set_font
+7.5 86.4 moveto 65 (POLYVEC) alignedtext
+grestore
+% A1AGG
+gsave
+[ /Rect [ 6 0 74 36 ]
+  /Border [ 0 0 0 ]
+  /Action << /Subtype /URI /URI (bookvol10.2.pdf#nameddest=A1AGG) >>
+  /Subtype /Link
+/ANN pdfmark
+0.60606 0.73333 1 nodecolor
+newpath 74 36 moveto
+6 36 lineto
+6 0 lineto
+74 0 lineto
+closepath fill
+1 setlinewidth
+filled
+0.60606 0.73333 1 nodecolor
+newpath 74 36 moveto
+6 36 lineto
+6 0 lineto
+74 0 lineto
+closepath stroke
+0 0 0 nodecolor
+14 /Times-Roman set_font
+14 14.4 moveto 52 (A1AGG) alignedtext
+grestore
+% POLYVEC->A1AGG
+gsave
+1 setlinewidth
+0 0 0 edgecolor
+newpath 40 71.83 moveto
+40 64.13 40 54.97 40 46.42 curveto
+stroke
+0 0 0 edgecolor
+newpath 43.5 46.41 moveto
+40 36.41 lineto
+36.5 46.41 lineto
+closepath fill
+1 setlinewidth
+solid
+0 0 0 edgecolor
+newpath 43.5 46.41 moveto
+40 36.41 lineto
+36.5 46.41 lineto
+closepath stroke
+grestore
+endpage
+showpage
+grestore
+%%PageTrailer
+%%EndPage: 1
+%%Trailer
+end
+restore
+%%EOF
diff --git a/changelog b/changelog
index b0bf79c..6e3851d 100644
--- a/changelog
+++ b/changelog
@@ -1,3 +1,8 @@
+20130227 tpd src/axiom-website/patches.html 20130227.02.tpd.patch
+20130227 tpd src/algebra/Makefile add U32VectorPolynomialOperations
+20130227 tpd books/bookvol5 add support for U32VectorPolynomialOperations
+20130227 tpd books/bookvol10.4 add U32VectorPolynomialOperations
+20130227 tpd books/ps/v104u32vectorpolynomialoperations.eps added
 20130227 tpd src/axiom-website/patches.html 20130227.01.tpd.patch
 20130227 tpd src/input/machinearithmetic.input unit test U8Matrix
 20130227 tpd src/algebra/Makefile add U8Matrix
diff --git a/src/algebra/Makefile.pamphlet b/src/algebra/Makefile.pamphlet
index 8a56e3a..c8a0d2c 100644
--- a/src/algebra/Makefile.pamphlet
+++ b/src/algebra/Makefile.pamphlet
@@ -4109,7 +4109,8 @@ LAYER8=\
   ${OUT}/LODOCAT-.o ${OUT}/LWORD.o    ${OUT}/MATCAT.o   ${OUT}/MATCAT-.o  \
   ${OUT}/MATSTOR.o  ${OUT}/ORESUP.o   ${OUT}/OREPCTO.o  ${OUT}/OREUP.o    \
   ${OUT}/PACFFC.o   \
-  ${OUT}/PLOT3D.o   ${OUT}/PR.o       ${OUT}/PREASSOC.o ${OUT}/PRIMARR2.o \
+  ${OUT}/PLOT3D.o   ${OUT}/POLYVEC.o \
+  ${OUT}/PR.o       ${OUT}/PREASSOC.o ${OUT}/PRIMARR2.o \
   ${OUT}/PROJSP.o   \
   ${OUT}/REDORDER.o ${OUT}/SRAGG.o    ${OUT}/SRAGG-.o   ${OUT}/STREAM.o   \
   ${OUT}/SYMPOLY.o  ${OUT}/TS.o       ${OUT}/TUPLE.o    ${OUT}/UPSCAT.o   \
@@ -4544,6 +4545,22 @@ LAYER8=\
 "PLOT3D" -> "TRANFUN"
 /*"PLOT3D" -> {"TRIGCAT"; "ATRIG"; "HYPCAT"; "AHYP"; "ELEMFUN"; "SPFCAT"}*/
 
+"POLYVEC" [color="#FF4488",href="bookvol10.4.pdf#nameddest=POLYVEC"]
+/*"POLYVEC" -> {"SINT" "NNI" "INT"}*/
+"POLYVEC" -> "A1AGG"
+/*"POLYVEC" -> {"FLAGG" "LNAGG" "IXAGG" "HOAGG" "AGG" "TYPE" "SETCAT"}*/
+/*"POLYVEC" -> {"BASTYPE" "KOERCE" "EVALAB" "IEVALAB" "ELTAGG" "ELTAB"}*/
+/*"POLYVEC" -> {"CLAGG" "KONVERT" "ORDSET" "BOOLEAN" "INS" "UFD" "GCDDOM"} */
+/*"POLYVEC" -> {"INTDOM" "COMRING" "RING" "RNG" "ABELGRP" "CABMON"}*/
+/*"POLYVEC" -> {"ABELMON" "ABELSG" "SGROUP" "MONOID" "LMODULE" "BMODULE"}*/
+/*"POLYVEC" -> {"RMODULE" "ALGEBRA" "MODULE" "ENTIRER" "EUCDOM" "PID"}*/
+/*"POLYVEC" -> {"OINTDOM" "ORDRING" "OAGROUP" "OCAMON" "OAMON" "OASGP"}*/
+/*"POLYVEC" -> {"DIFRING" "RETRACT" "LINEXP" "PATMAB" "CFCAT" "REAL"}*/
+/*"POLYVEC" -> {"CHARZ" "STEP" "INS-" "EUCDOM-" "UFD-" "GCDDOM-" "INTDOM-"}*/
+/*"POLYVEC" -> {"ALGEBRA-" "DIFRING-" "ORDRING-" "MODULE-""RING-"}*/
+/*"POLYVEC" -> {"ABELGRP-" "ABELMON-" "MONOID-" "ORDSET-" "ABELSG-"}*/
+/*"POLYVEC" -> {"SGROUP-" "SETCAT-" "RETRACT-" "BASTYPE-" "PI" "PRIMARR"}*/
+
 "PR" [color="#88FF44",href="bookvol10.3.pdf#nameddest=PR"]
 "PR" -> "FAMR"
 /*"PR" -> {"AMR"; "RING"; "RNG"; "ABELGRP"; "CABMON"; "ABELMON"; "ABELSG"}*/
@@ -18399,6 +18416,7 @@ REGRESS= \
  U8Matrix.regress \
  U16Matrix.regress \
  U32Matrix.regress \
+ U32VectorPolynomialOperations.regress \
  U8Vector.regress \
  U16Vector.regress \
  U32Vector.regress \
diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html
index 431ead0..063d677 100644
--- a/src/axiom-website/patches.html
+++ b/src/axiom-website/patches.html
@@ -3993,5 +3993,7 @@ books/bookvol10.4 comment out bad code in GUESS
 books/bookvol10.3 add U16Matrix
 <a href="patches/20130227.01.tpd.patch">20130227.01.tpd.patch</a>
 books/bookvol10.3 add U8Matrix
+<a href="patches/20130227.02.tpd.patch">20130227.02.tpd.patch</a>
+books/bookvol10.4 add U32VectorPolynomialOperations
  </body>
 </html>
