############################################################################
#
#	File:     polynom.icn
#
#	Subject:  Procedures to manipulate multi-variate polynomials
#
#	Author:   Ralph E. Griswold
#
#	Date:     October 1, 2001
#
############################################################################
#
#   This file is in the public domain.
#
############################################################################
#
#  The format for strings omits symbols for multiplication and
#  exponentiation.  For example, 3*a^2 is entered as 3a2.
#
#  A polynomial is represented by a table in which each term, such as 3xy,
#  the xy is #  a key and the corresponding value is the coefficient, 3 in
#  this case.  If a variable is raised to a power, such as x^3, the key
#  is the product of the individual variables, xxx in this case.
#
############################################################################
#
#  Links:  strings, tables
#
############################################################################

link strings
link tables

procedure str2poly(str)			#: convert string to polynomial
   local poly, var, vars, term, factor, power

   poly := table(0)

   str ? {
      while term := (move(1) || tab(upto('-+') | 0)) do {	# possible sign
         term ? {
            factor := 1			# default
            factor := tab(many(&digits ++ '+.-'))
            tab(0) ? {
               vars := ""
               while var := move(1) do {
                  power := 1		# default
                  power := integer(tab(many(&digits)))
                  vars ||:= repl(var, power)
                  }
               }
            poly[csort(vars)] +:= numeric(factor) | fail
            }
         }
      }

   return poly

end

procedure polyadd(poly1, poly2)		#: add polynomials
   local poly, keys, k

   keys := sort(set(keylist(poly1)) ++ set(keylist(poly2)))

   poly := table(0)

   every k := !keys do
      poly[k] := poly1[k] + poly2[k]

   return poly

end

procedure polymod(poly, i)		#: polynomial modular reduction
   local poly1, keys, k

   keys := keylist(poly)

   poly1 := table(0)

   every k := !keys do
      poly1[k] := poly[k] % i

   return poly1

end

procedure polysub(poly1, poly2)		#: subtract polynomials
   local poly, keys, k

   keys := sort(set(keylist(poly1)) ++ set(keylist(poly2)))

   poly := table(0)

   every k := !keys do
      poly[k] := poly1[k] - poly2[k]

   return poly

end

procedure polymul(poly1, poly2)		#: multiply polynomials
   local poly, keys1, keys2, k1, k2

   keys1 := keylist(poly1)
   keys2 := keylist(poly2)

   poly := table(0)

   every k1 := !keys1 do
      every k2 := !keys2 do
         poly[csort(k1 || k2)] +:= poly1[k1] * poly2[k2]

   return poly

end

procedure polyexp(poly1, i)		#: exponentiate polynomial
   local poly

   poly := copy(poly1)

   every 1 to i - 1 do 
      poly := polymul(poly, poly1)

   return poly

end

procedure poly2str(poly)		#: polynomial to string
   local str, keys, k, count, var

   keys := keylist(poly)

   str := ""

   every k := !keys do {
      if poly[k] = 0 then next		# skip term
      else if poly[k] > 0 then str ||:= "+" || ((poly[k] ~= 1) | "")
      else if poly[k] < 0 then str ||:= ((poly[k] ~= 1) | "")
      k ? {
         while var := move(1) do {
            count := 1
            count +:= *tab(many(var))
            if count = 1 then str ||:= var
            else str ||:= var || count
            }
         }
      }

   return str[2:0] | "0"

end

procedure polydiff(poly, var)		#: polynomial differentiation
   local poly_new, keys, k, nvars, newk

   poly_new := table()

   keys := keylist(poly)

   every k := !keys do {
      k ? {
         if newk := tab(upto(var)) then {
            nvars := *tab(many(var))
            newk ||:= repl(var, nvars - 1) || tab(0)
            poly_new[newk] := nvars * poly[k]
            }
         }
      }

   return poly_new

end

procedure polyintg(poly, var)		#: polynomial integration
   local poly_new, keys, k, nvars, newk

   poly_new := table()

   keys := keylist(poly)

   every k := !keys do {
      k ? {
         if newk := tab(upto(var)) then {
            nvars := *tab(many(var))
            newk ||:= repl(var, nvars + 1) || tab(0)
            poly_new[newk] := poly[k] / real(nvars + 1) 
            }
         }
      }

   return poly_new

end

procedure peval(str)			#: string polynomial simplification

   while str ?:= 2(="(", tab(bal(')')), =")", pos(0))

   return poper(str) | str2poly(str)

end

procedure poper(str)			#: find polynomial operation

   return str ? {
      pform(tab(bal('-+*^:|%')), move(1), tab(0))
      }

end

procedure pform(str1, op, str2)		#: polynomial formation

   return case op of {
      "+"  :  polyadd(peval(str1), peval(str2))
      "-"  :  polysub(peval(str1), peval(str2))
      "*"  :  polymul(peval(str1), peval(str2))
      "^"  :  polyexp(peval(str1), str2)
      ":"  :  polydiff(peval(str1), str2)
      "|"  :  polyintg(peval(str1), str2)
      "%"  :  polymod(peval(str1), str2)
      }

end

procedure poly2profile(poly)		#: polynomial to profile sequence
   local str, keys, k, count, vara, i, seg

   keys := keylist(poly)

   str := ""

   every k := !keys do {
      i := poly[k]
      if i < 0 then {			# if negative, reverse sequence
         i := abs(i)
         k := reverse(k)
         }
      str ||:= left(repl(k, i + 1), i * *k)
      }

   return str

end

procedure poly2profilelen(poly)		#: polynomial to profile sequence
   local i, keys, k, count, var

   keys := keylist(poly)

   i := 0

   every k := !keys do
       i +:= *repl(k, abs(poly[k]))	# treat negative as if positive

   return i

end

procedure basepolystr(clist, plist)	# base polynomial string

   return "(" || poly2str(basepoly(clist, plist)) || ")"

end

procedure basepoly(clist, plist)	# base polynomial
   local poly, i, c, p
   static vlist

   initial vlist := string(&lcase)

   poly := table()

   i := 1

   while c := get(clist) & p := get(plist) do {
      poly[repl(vlist[i], (0 <= p))] := (0 ~= c)
      i +:= 1
      }

   return poly

end