cyclotomicMat.at

{ update one entry in a list of polynomials, for instance in a row of a matrix }
set update_row ([i_poly] R, int j,i_poly v) = [i_poly]: R[j]:=v; R
{ update one row of a matrix (replacing it by a given new row) }
set update_matrix_row (CyclotomicMat M, int i, [i_poly] row)= CyclotomicMat:
   M[i]:=row; M
{ update one entry in a matrix }
set update_matrix_entry(CyclotomicMat M, int i, int j, i_poly v)= CyclotomicMat:
  let row=M[i] in row[j]:=v; M[i]:=row; M

set zero_poly_row (int n) = [i_poly]: for :n do poly_0 od
set zero_poly_matrix (int rows, int columns) = CyclotomicMat:
  let !row=zero_poly_row(columns) in for :rows do row od
set zero_poly_matrix(int n) = CyclotomicMat:zero_poly_matrix(n,n)
set scalar_poly_matrix (int n, int c) = CyclotomicMat:
  let !row=zero_poly_row(n) { model row, almost equal to each wanted row } in
  for i:n do update_row(row,i,constant_poly(c)) od
set identity_poly_matrix (int n) = CyclotomicMat: scalar_poly_matrix(n,1)
set mat_as_poly_mat (mat M) = CyclotomicMat:
 let (r,c)=shape(M) in for i:r do for j:c do [M[i,j]] od od

set + (CyclotomicMat M, i_poly p) = CyclotomicMat: { add polynomial scalar matrix }
  for row@i in M  do for e@j in row do if i=j then e+p else e fi od od
set - (CyclotomicMat M, i_poly p) = CyclotomicMat: M+-p

set ^ = (i_poly,int->i_poly): { power of a polynomial }
  let one=i_poly:[1] { use local constant for efficiency } in
  (i_poly p,int n) i_poly: { the actual powering function }
    if n=0 then [1]
    else { no point in trying to do this by repeated squaring: is not faster }
      let f=p in while (n-:=1)>0 do f *:= p od; f
    fi

{ equality of polynomial matrices }
set = (CyclotomicMat A,CyclotomicMat B) = bool:
  all(#A, (int i): let r=A[i],s=B[i] in #r=#s and all(#r, (int j): r[j]=s[j] ))

{ polynomial matrix is 0 if all its entries are zero-length polynomials }
set is_zero (CyclotomicMat M) = bool:
  all(#M, (int i): let r=M[i] in all(#r, (int j): =#r[j]) )

set upper_unitriangular_inverse (CyclotomicMat P) = CyclotomicMat:
  let n=#P { number of rows in matrix }, dummy_row=zero_poly_row(0)
  then Q= for :n do dummy_row od { prepare |n| rows, but leave them empty now }
  in { solve equation $P*Q=I$ for $Q$, proceeding from bottom to top }
  for i:n { reverse traversal loop with reverse accumulation -> straight value }
  ~do let Pi=P[i] { we use only this part of P during this iteration } in
    Q[i] := { computed row is stored as |Q[i]|, and contributed to result }
    (zero_poly_row(i) # poly_1) { beginning of row |Q[i]| is somewhat boring }
    ## { concatenate with more interesting newly computed part of |Q[i]| }
      for j: n-i-1 from i+1 { left to right, but computations are independent }
      do let sum=poly_0 in for k: j-i from i+1 do sum -:= Pi[k] * Q[k][j] od
      ; sum    { contribute $Q_{i,j} = - \sum_{k=i+1}^j P_{i,k}Q_{k,j}$ }
      od
  ~od

set lower_unitriangular_inverse(CyclotomicMat M)=transpose(upper_unitriangular_inverse(transpose(M)))

{ compute PMP^{-1} where P is a permutation of [1,,,n]
  and M is a polynomial matrix
}
set poly_permute_basis([int] P, CyclotomicMat A) = CyclotomicMat: { P is permutation }
  for i:#A do let row=A[P[i]] in for j:#row do row[P[j]] od od

{ some printing functions }

{ make a printable string out of a polynomial and name of indeterminate,
  for use in other print functions }
set poly_format (i_poly v, string q)= string:
   if =v then "0"
   else
    ##for c@k in v
      do
         if =c then ""
         else
            let coef = if c<0 then "-" else "+" fi +
                if abs(c)=1 and k>0 then "" else int_format(abs(c)) fi
	    ,  term = if k>1 then q+"^"+int_format(k) elif k=0 then "" else q fi
	    in string: coef+term
         fi
      od
   fi

{ for |int_poly| and |rat_poly| that value names the indeterminate; use it }
set to_string (int_poly(X,cs,)) = string:
   if =#cs then "0"
   else
      let pow(int k) = string: { print X^k, possibly shortened }
          if k>1 then X+"^"+k.to_string elif k=1 then X else "" fi
      in
    ##for c@k in cs
      do case c
         then "-" + if c!=minus_1 or =k then (-c).to_string else "" fi + pow(k)
         in "" { nothing for terms with a zero coefficient }
         else "+" + if c>1 or =k then c.to_string else "" fi + pow(k)
         esac
      od
   fi

set to_string (rat_poly(X,cs,)) = string:
   if =#cs then "0"
   else
      let pow(int k) = string: { print X^k, possibly shortened }
          if k>1 then X+"^"+k.to_string elif k=1 then X else "" fi
      in
    ##for c@k in cs
      do if c.= then ""
         else let (n,d)=%c in
         ##[if n.< then "-" next n:=-n else "+" fi
           ,if d>1 then c.to_string elif n>1 or k=0 then n.to_string else "" fi
           ,if k>1 then X+"^"+k.to_string elif k=1 then X else "" fi
           ]
         fi
      od
   fi


{ nice printing of a polynomial }
set printPoly (i_poly v) = void: prints(poly_format(v,"q"))

{ M=poly_mat is a polynomial matrix; e.g.
  expr> set M= poly_mat: [[ [0,1], [0,0,1]], [[0], [1]]]
  expr> printPolyMatrix (M)
   +q +q^2
   0 +1
}

set printPolyMatrix(CyclotomicMat M,int space_size,string char)= void:
  let space= " "*space_size, wd = max(for v in ##M do #poly_format(v,char) od)
  in
    for row in M
    do prints(##for v in row do space+l_adjust(wd,poly_format(v,char)) od)
    od
set printPolyMatrix(CyclotomicMat M,int space_size)= void:printPolyMatrix(M,space_size,"q")
set printPolyMatrix (CyclotomicMat M,string char) = void: printPolyMatrix(M,3,char)
set printPolyMatrix (CyclotomicMat M) = void: printPolyMatrix(M,3)

{ this provides maybe more meaningful names for |flex_add| and friends }
set poly_plus    = +@(i_poly,i_poly) { this binding will remain exported }
set poly_minus   = -@(i_poly,i_poly) { this binding will remain exported }
set poly_product = *@(i_poly,i_poly) { this binding will remain exported }

{ restore normal vector arithmetic }
forget +@(vec,vec)
forget -@(vec,vec)
forget *@(vec,vec)
set + = fix_add@(vec,vec)
set - = fix_sub@(vec,vec)
set * = dot_product@(vec,vec)

set sgn_poly (int k) = i_poly:
  assert(!=k,"taking sign of 0"); [if k.< then -1 else 1 fi]

set divide_by (int k,i_poly v) = i_poly: ratvec_as_vec(v/k)

{i_polynomial has only non-negative coefficients}
set is_non_negative = (i_poly->bool): >=@vec { same condition as "dominance" }

set top_left_square_block (CyclotomicMat M,int size) = CyclotomicMat:
  for row in M[:size] do row[:size] od