http://mathling.com/type/polynomial/complex  library module

Imports

import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy"

import module namespace poly="http://mathling.com/type/polynomial"
       at "../types/polynomial.xqy"

import module namespace z="http://mathling.com/core/complex"
       at "../core/complex.xqy"

import module namespace zv="http://mathling.com/core/complex/vector"
       at "../core/vcomplex.xqy"

import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy"

Functions

declare function this:polynomial($coefficients as map(xs:string,item()*)*) as map(xs:string,item()*)
{
  (: Prune leading zeros :)
  let $useful :=
    let $first-nonzero := head(for $c at $i in $coefficients where not(z:eq($c,0)) return $i)
    return $coefficients[position() >= $first-nonzero]
  (: Make sure we have something :)
  let $useful :=
    if (empty($useful)) then (z:as-complex(0)) else $useful
  return
    map {
      "kind": "complex-polynomial",
      "coefficients": reverse($useful) (: Most operations need to use reverse :)
    }
}

declare function this:polynomial(
  $coefficients as map(xs:string,item()*)*,
  $op as xs:string?
) as map(xs:string,item()*)
{
  if ($op = ("sin","cos","sinh","cosh")) then (
    this:polynomial($coefficients)=>map:put("op", $op)
  ) else if (empty($op) or $op="none") then (
    this:polynomial($coefficients)
  ) else errors:error("POLY-BADOP", $op)
}

declare function this:as-polynomial(
  $polynomial as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if ($polynomial("kind")="complex-polynomial") then $polynomial
  else this:polynomial(poly:coefficients($polynomial)!z:as-complex(.), poly:op($polynomial))
}

declare function this:as-real-polynomial(
  $polynomial as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if ($polynomial("kind")="polynomial") then $polynomial
  else poly:polynomial(this:coefficients($polynomial)!z:as-real(.), this:op($polynomial))
}

declare function this:is-real(
  $polynomial as map(xs:string,item()*)
) as xs:boolean
{
  ($polynomial("kind")="polynomial") or
  (
    every $c in this:coefficients($polynomial) satisfies z:is-real($c)
  )
}

declare function this:coefficients($poly as map(xs:string,item()*)) as map(xs:string,item()*)*
{
  reverse($poly("coefficients"))
}

declare function this:order($poly as map(xs:string,item()*)) as xs:integer
{
  count($poly("coefficients")) - 1
}

declare function this:op($poly as map(xs:string,item()*)) as xs:string?
{
  $poly("op")
}

declare function this:value(
  $poly as map(xs:string,item()*),
  $x as xs:double
) as map(xs:string,item()*)
{
  switch ($poly("op"))
  case "sin" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(math:sin($x), $i - 1) ))
  case "cos" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(math:cos($x), $i - 1) ))
  case "sinh" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(util:sinh($x), $i - 1) ))
  case "cosh" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(util:cosh($x), $i - 1) ))
  default return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow($x, $i - 1) ))
}

declare function this:zvalue(
  $poly as map(xs:string,item()*),
  $z as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  switch ($poly("op"))
  case "sin" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:sin($z), $i - 1)=>z:multiply($c))
  case "cos" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cos($z), $i - 1)=>z:multiply($c))
  case "sinh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:sinh($z), $i - 1)=>z:multiply($c))
  case "cosh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cosh($z), $i - 1)=>z:multiply($c))
  default return
    z:sum(for $c at $i in $poly("coefficients") return z:pow($z, $i - 1)=>z:multiply($c))
}

declare function this:zvaluev(
  $poly as map(xs:string,item()*),
  $z as xs:double*
) as xs:double*
{
  switch ($poly("op"))
  case "sin" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:sin($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  case "cos" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cos($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  case "sinh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:sinh($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  case "cosh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cosh($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  default return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow($z, $i - 1)=>zv:multiply(z:as-vector($c)))
}

declare function this:antiderivative($poly as map(xs:string,item()*)) as map(xs:string,item()*)
{
  let $k := if ($poly("op") = "sin") then -1 else 1
  let $op :=
    switch ($poly("op"))
    case "sin" return "cos"
    case "cos" return "sin"
    case "sinh" return "cosh"
    case "cosh" return "sinh"
    default return ()
  return 
    this:polynomial(
      (reverse(for $c at $i in $poly("coefficients") return $c=>z:times($k div $i)), z:as-complex(0)),
      $op
    )
}

declare function this:derivative($poly as map(xs:string,item()*)) as map(xs:string,item()*)
{
  let $k := if ($poly("op") = "cos") then -1 else 1
  let $op :=
    switch ($poly("op"))
    case "sin" return "cos"
    case "cos" return "sin"
    case "sinh" return "cosh"
    case "cosh" return "sinh"
    default return ()
  return (
    this:polynomial(
      if (count($poly("coefficients")) = 1) then z:as-complex(0)
      else reverse(for $c at $i in tail($poly("coefficients")) return $c=>z:times($k*$i)),
      $op
    )
  )
}

declare function this:definite-integral(
  $poly as map(xs:string,item()*),
  $a as xs:double,
  $b as xs:double
) as map(xs:string,item()*)
{
  let $anti := this:antiderivative($poly)
  return $anti=>this:value($b)=>z:sub( $anti=>this:value($a) )
}

declare function this:describe(
  $poly as map(xs:string,item()*)
) as xs:string
{
  if (count($poly("coefficients"))=1) then z:quote($poly("coefficients"))
  else (
    let $z := if (exists($poly("op"))) then $poly("op")||"(z)" else "z"
    return (
      string-join(
        reverse(
          for $c at $i in $poly("coefficients") return (
            if (z:eq($c,0)) then ()
            else (
              let $xton :=
                if ($i - 1 = 0) then ""
                else if ($i - 1 = 1) then $z
                else $z||"^"||($i - 1)
              let $k :=
                if (z:eq($c,1) and $i > 1) then ""
                else if (z:eq($c,-1) and $i > 1) then "-"
                else if (z:is-real($c)) then xs:string(z:pr($c))
                else "("||z:quote($c)||")"
              return $k||$xton (: c*x^(i-1) :)
            )
          )
        ),
        " + "
      )
    )
  )
}

declare function this:function(
  $poly as map(xs:string,item()*)
) as function(xs:double) as map(xs:string,item()*)
{
  function($v as xs:double) as map(xs:string,item()*) {
    $poly=>this:value($v)
  }
}

declare function this:zfunction(
  $poly as map(xs:string,item()*)
) as function(map(xs:string,item()*)) as map(xs:string,item()*)
{
  function($v as map(xs:string,item()*)) as map(xs:string,item()*) {
    $poly=>this:zvalue($v)
  }
}

declare function this:add(
  $p1 as map(xs:string,item()*),
  $p2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if (this:op($p1) != this:op($p2))
  then errors:error("POLY-INCOMPOP", (this:op($p1), this:op($p2), "add"))
  else (),
  let $c1 :=
    if (this:order($p1) > this:order($p2))
    then $p1("coefficients")
    else $p2("coefficients")
  let $c2 := 
    if (this:order($p1) > this:order($p2))
    then $p2("coefficients")
    else $p1("coefficients")
  return (
    this:polynomial(
      reverse(
        for $c at $i in $c1 return $c=>z:add(($c2[$i],z:as-complex(0))[1])
      ),
      this:op($p1)
    )
  )
}

declare function this:sub(
  $p1 as map(xs:string,item()*),
  $p2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if (this:op($p1) != this:op($p2))
  then errors:error("POLY-INCOMPOP", (this:op($p1), this:op($p2), "add"))
  else (),
  let $c1 :=
    if (this:order($p1) > this:order($p2))
    then $p1("coefficients")
    else $p2("coefficients")!z:minus(.)
  let $c2 := 
    if (this:order($p1) > this:order($p2))
    then $p2("coefficients")!z:minus(.)
    else $p1("coefficients")
  return (
    this:polynomial(
      reverse(
        for $c at $i in $c1 return $c=>z:add( ($c2[$i],z:as-complex(0))[1] )
      ),
      this:op($p1)
    )
  )
}

declare function this:multiply(
  $p1 as map(xs:string,item()*),
  $p2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if (this:op($p1) != this:op($p2))
  then errors:error("POLY-INCOMPOP", (this:op($p1), this:op($p2), "add"))
  else (),
  let $o1 :=
    if (this:order($p1) > this:order($p2))
    then this:order($p1)
    else this:order($p2)
  let $c1 := (
    if (this:order($p1) > this:order($p2))
    then $p1("coefficients")
    else $p2("coefficients")
  )
  let $c2 := (
    if (this:order($p1) > this:order($p2))
    then $p2("coefficients")
    else $p1("coefficients")
  )
  return (
    this:polynomial(
      reverse(
        fold-left(2 to count($c2), (for $d in $c1 return $c2[1]=>z:multiply($d)),
          function($s as map(xs:string,item()*)*, $i as xs:integer) as map(xs:string,item()*)* {
            let $sub := (
              for $j in $o1 - $i + 1 to $o1 - 1 return z:as-complex(0),
              for $d in $c1 return $c2[$i]=>z:multiply($d)
            )
            return (
              for $j in 1 to count($sub) return $sub[$j]=>z:add( ($s[$j],z:as-complex(0))[1] )
            )
          }
        )
      ),
      this:op($p1)
    )
  )
}

declare function this:times(
  $p1 as map(xs:string,item()*),
  $k as xs:double
) as map(xs:string,item()*)
{
  this:polynomial(
    reverse( $p1("coefficients")!(.=>z:times($k)) ),
    this:op($p1)
  )
}

declare function this:ztimes(
  $p1 as map(xs:string,item()*),
  $kz as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:polynomial(
    reverse( $p1("coefficients")!(.=>z:multiply($kz)) ),
    this:op($p1)
  )
}

Original Source Code

xquery version "3.1";
(:~
 : Polynomials with complex coefficients; represented as vector of coefficients
 : Most operations require us to process them in reverse order, so we 
 : store them that way. The accessor flips the order back, so internally we 
 : generally do not use it.
 :
 : Copyright© Mary Holstege 2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @since Apr 2023
 : @custom:Status Bleeding edge
 :)
module namespace this="http://mathling.com/type/polynomial/complex"; 

import module namespace errors="http://mathling.com/core/errors"
       at "../core/errors.xqy";
import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy";
import module namespace z="http://mathling.com/core/complex"
       at "../core/complex.xqy";
import module namespace zv="http://mathling.com/core/complex/vector"
       at "../core/vcomplex.xqy";
import module namespace poly="http://mathling.com/type/polynomial"
       at "../types/polynomial.xqy";

declare namespace map="http://www.w3.org/2005/xpath-functions/map";
declare namespace math="http://www.w3.org/2005/xpath-functions/math";

(:~
 : polynomial()
 : Construct a representation of a polynomial
 : @param $coefficients: coefficients of polynomial in decreasing order, e.g.
 : 3x³ - 2 = (z:as-complex(3), z:as-complex(0), z:as-complex(0), z:as-complex(-1))
 :)
declare function this:polynomial($coefficients as map(xs:string,item()*)*) as map(xs:string,item()*)
{
  (: Prune leading zeros :)
  let $useful :=
    let $first-nonzero := head(for $c at $i in $coefficients where not(z:eq($c,0)) return $i)
    return $coefficients[position() >= $first-nonzero]
  (: Make sure we have something :)
  let $useful :=
    if (empty($useful)) then (z:as-complex(0)) else $useful
  return
    map {
      "kind": "complex-polynomial",
      "coefficients": reverse($useful) (: Most operations need to use reverse :)
    }
};

(:~
 : polynomial()
 : Construct a representation of a polynomial in op(x) rather than x.
 : e.g. , polynomial( (z:complex(3,1), 0, 0, z:complex(-1,0)), "sin" ) =
 : (3+i)sin(x)³ - 2
 : @param $coefficients: coefficients of polynomial in decreasing order
 : @param $op: operator, one of "sin", "cos", "sinh", "cosh" or empty or "none"
 :)
declare function this:polynomial(
  $coefficients as map(xs:string,item()*)*,
  $op as xs:string?
) as map(xs:string,item()*)
{
  if ($op = ("sin","cos","sinh","cosh")) then (
    this:polynomial($coefficients)=>map:put("op", $op)
  ) else if (empty($op) or $op="none") then (
    this:polynomial($coefficients)
  ) else errors:error("POLY-BADOP", $op)
};

(:~
 : as-polynomial()
 : Construct a polynomial with complex coefficents from one with real coefficients.
 : If the polynomial already is a complex polynomial, just return it
 :
 : @param $polynomial: polynomial 
 : @return polynomial with complex coefficients
 :)
declare function this:as-polynomial(
  $polynomial as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if ($polynomial("kind")="complex-polynomial") then $polynomial
  else this:polynomial(poly:coefficients($polynomial)!z:as-complex(.), poly:op($polynomial))
};

(:~
 : as-real-polynomial()
 : Convert to a polynomial with real coefficients, if possible. Raises UTIL-CAST
 : if that is not possible.
 :
 : @param $polynomial: polynomial 
 : @return polynomial with real coefficients
 :)
declare function this:as-real-polynomial(
  $polynomial as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if ($polynomial("kind")="polynomial") then $polynomial
  else poly:polynomial(this:coefficients($polynomial)!z:as-real(.), this:op($polynomial))
};

(:~
 : is-real()
 : Return true if every coefficient is real-valued.
 :
 : @param $polynomial: polynomial 
 : @return true() if polynomial has real coefficients
 :)
declare function this:is-real(
  $polynomial as map(xs:string,item()*)
) as xs:boolean
{
  ($polynomial("kind")="polynomial") or
  (
    every $c in this:coefficients($polynomial) satisfies z:is-real($c)
  )
};

(:~
 : coefficients()
 : Accessor for coefficients of polynomial
 :
 : @param $poly: the polynomial
 :)
declare function this:coefficients($poly as map(xs:string,item()*)) as map(xs:string,item()*)*
{
  reverse($poly("coefficients"))
};

(:~
 : order()
 : Largest exponent. e.g. order(3x³ - 2) = 3
 :
 : @param $poly: the polynomial
 :)
declare function this:order($poly as map(xs:string,item()*)) as xs:integer
{
  count($poly("coefficients")) - 1
};

(:~
 : op()
 : Operator.
 :
 : @param $poly: the polynomial
 :)
declare function this:op($poly as map(xs:string,item()*)) as xs:string?
{
  $poly("op")
};

(:~
 : value()
 : Compute f(x) where f is a polynomial in x (a double)
 :
 : @param $poly: the polynomial
 : @param $x the value to compute the function of
 :)
declare function this:value(
  $poly as map(xs:string,item()*),
  $x as xs:double
) as map(xs:string,item()*)
{
  switch ($poly("op"))
  case "sin" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(math:sin($x), $i - 1) ))
  case "cos" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(math:cos($x), $i - 1) ))
  case "sinh" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(util:sinh($x), $i - 1) ))
  case "cosh" return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow(util:cosh($x), $i - 1) ))
  default return
    z:sum(for $c at $i in $poly("coefficients") return $c=>z:times( math:pow($x, $i - 1) ))
};

(:~
 : zvalue()
 : Compute f(z) where f is a polynomial in z (a complex number)
 :
 : @param $poly: the polynomial
 : @param $z: the complex value to compute the function of
 :)
declare function this:zvalue(
  $poly as map(xs:string,item()*),
  $z as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  switch ($poly("op"))
  case "sin" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:sin($z), $i - 1)=>z:multiply($c))
  case "cos" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cos($z), $i - 1)=>z:multiply($c))
  case "sinh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:sinh($z), $i - 1)=>z:multiply($c))
  case "cosh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cosh($z), $i - 1)=>z:multiply($c))
  default return
    z:sum(for $c at $i in $poly("coefficients") return z:pow($z, $i - 1)=>z:multiply($c))
};

(:~
 : zvaluev()
 : Compute f(z) where f is a polynomial in z (a complex number vector)
 :
 : @param $poly: the polynomial
 : @param $z: the complex value vector to compute the function of
 :)
declare function this:zvaluev(
  $poly as map(xs:string,item()*),
  $z as xs:double*
) as xs:double*
{
  switch ($poly("op"))
  case "sin" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:sin($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  case "cos" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cos($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  case "sinh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:sinh($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  case "cosh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cosh($z), $i - 1)=>zv:multiply(z:as-vector($c)))
  default return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow($z, $i - 1)=>zv:multiply(z:as-vector($c)))
};

(:~
 : antiderivative()
 : Anti-derivative of a polynomial, which is another polynomial. We use 0 as the
 : integration constant.
 : 
 : @param $poly: the polynomial
 :)
declare function this:antiderivative($poly as map(xs:string,item()*)) as map(xs:string,item()*)
{
  let $k := if ($poly("op") = "sin") then -1 else 1
  let $op :=
    switch ($poly("op"))
    case "sin" return "cos"
    case "cos" return "sin"
    case "sinh" return "cosh"
    case "cosh" return "sinh"
    default return ()
  return 
    this:polynomial(
      (reverse(for $c at $i in $poly("coefficients") return $c=>z:times($k div $i)), z:as-complex(0)),
      $op
    )
};

(:~
 : derivative()
 : Derivative of a polynomial, which is another polynomial.
 : 
 : @param $poly: the polynomial
 :)
declare function this:derivative($poly as map(xs:string,item()*)) as map(xs:string,item()*)
{
  let $k := if ($poly("op") = "cos") then -1 else 1
  let $op :=
    switch ($poly("op"))
    case "sin" return "cos"
    case "cos" return "sin"
    case "sinh" return "cosh"
    case "cosh" return "sinh"
    default return ()
  return (
    this:polynomial(
      if (count($poly("coefficients")) = 1) then z:as-complex(0)
      else reverse(for $c at $i in tail($poly("coefficients")) return $c=>z:times($k*$i)),
      $op
    )
  )
};

(:~
 : definite-integral()
 : Definite integral of a polynomial (exact)
 : 
 : @param $poly: the polynomial
 : @param $a: start of range of integration
 : @param $b: end of range of integration
 :)
declare function this:definite-integral(
  $poly as map(xs:string,item()*),
  $a as xs:double,
  $b as xs:double
) as map(xs:string,item()*)
{
  let $anti := this:antiderivative($poly)
  return $anti=>this:value($b)=>z:sub( $anti=>this:value($a) )
};

(:~
 : describe()
 : Describe the polynomial
 :
 : @param $poly: the polynomial
 :)
declare function this:describe(
  $poly as map(xs:string,item()*)
) as xs:string
{
  if (count($poly("coefficients"))=1) then z:quote($poly("coefficients"))
  else (
    let $z := if (exists($poly("op"))) then $poly("op")||"(z)" else "z"
    return (
      string-join(
        reverse(
          for $c at $i in $poly("coefficients") return (
            if (z:eq($c,0)) then ()
            else (
              let $xton :=
                if ($i - 1 = 0) then ""
                else if ($i - 1 = 1) then $z
                else $z||"^"||($i - 1)
              let $k :=
                if (z:eq($c,1) and $i > 1) then ""
                else if (z:eq($c,-1) and $i > 1) then "-"
                else if (z:is-real($c)) then xs:string(z:pr($c))
                else "("||z:quote($c)||")"
              return $k||$xton (: c*x^(i-1) :)
            )
          )
        ),
        " + "
      )
    )
  )
};

(:~
 : function()
 : Return a function over values that takes values of the polynomial, e.g. for curve plotting
 :
 : @param $poly: the polynomial
 :)
declare function this:function(
  $poly as map(xs:string,item()*)
) as function(xs:double) as map(xs:string,item()*)
{
  function($v as xs:double) as map(xs:string,item()*) {
    $poly=>this:value($v)
  }
};

(:~
 : zfunction()
 : Return a function over complex values that takes values of the polynomial, e.g. for curve plotting
 :
 : @param $poly: the polynomial
 :)
declare function this:zfunction(
  $poly as map(xs:string,item()*)
) as function(map(xs:string,item()*)) as map(xs:string,item()*)
{
  function($v as map(xs:string,item()*)) as map(xs:string,item()*) {
    $poly=>this:zvalue($v)
  }
};

(:======================================================================
 : Operations between two polynomials
 :======================================================================:)

(:~
 : add()
 : Add two polynomials
 :
 : @param $p1: one polynomial
 : @param $p2: another polynomial
 :)
declare function this:add(
  $p1 as map(xs:string,item()*),
  $p2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if (this:op($p1) != this:op($p2))
  then errors:error("POLY-INCOMPOP", (this:op($p1), this:op($p2), "add"))
  else (),
  let $c1 :=
    if (this:order($p1) > this:order($p2))
    then $p1("coefficients")
    else $p2("coefficients")
  let $c2 := 
    if (this:order($p1) > this:order($p2))
    then $p2("coefficients")
    else $p1("coefficients")
  return (
    this:polynomial(
      reverse(
        for $c at $i in $c1 return $c=>z:add(($c2[$i],z:as-complex(0))[1])
      ),
      this:op($p1)
    )
  )
};

(:~
 : sub()
 : Subtract two polynomials
 :
 : @param $p1: one polynomial
 : @param $p2: another polynomial
 :)
declare function this:sub(
  $p1 as map(xs:string,item()*),
  $p2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if (this:op($p1) != this:op($p2))
  then errors:error("POLY-INCOMPOP", (this:op($p1), this:op($p2), "add"))
  else (),
  let $c1 :=
    if (this:order($p1) > this:order($p2))
    then $p1("coefficients")
    else $p2("coefficients")!z:minus(.)
  let $c2 := 
    if (this:order($p1) > this:order($p2))
    then $p2("coefficients")!z:minus(.)
    else $p1("coefficients")
  return (
    this:polynomial(
      reverse(
        for $c at $i in $c1 return $c=>z:add( ($c2[$i],z:as-complex(0))[1] )
      ),
      this:op($p1)
    )
  )
};

(:~
 : multiply()
 : Multiply two polynomials
 :
 : @param $p1: one polynomial
 : @param $p2: another polynomial
 :)
declare function this:multiply(
  $p1 as map(xs:string,item()*),
  $p2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  if (this:op($p1) != this:op($p2))
  then errors:error("POLY-INCOMPOP", (this:op($p1), this:op($p2), "add"))
  else (),
  let $o1 :=
    if (this:order($p1) > this:order($p2))
    then this:order($p1)
    else this:order($p2)
  let $c1 := (
    if (this:order($p1) > this:order($p2))
    then $p1("coefficients")
    else $p2("coefficients")
  )
  let $c2 := (
    if (this:order($p1) > this:order($p2))
    then $p2("coefficients")
    else $p1("coefficients")
  )
  return (
    this:polynomial(
      reverse(
        fold-left(2 to count($c2), (for $d in $c1 return $c2[1]=>z:multiply($d)),
          function($s as map(xs:string,item()*)*, $i as xs:integer) as map(xs:string,item()*)* {
            let $sub := (
              for $j in $o1 - $i + 1 to $o1 - 1 return z:as-complex(0),
              for $d in $c1 return $c2[$i]=>z:multiply($d)
            )
            return (
              for $j in 1 to count($sub) return $sub[$j]=>z:add( ($s[$j],z:as-complex(0))[1] )
            )
          }
        )
      ),
      this:op($p1)
    )
  )
};


(:~
 : times()
 : Multiply polynomial by double constant
 :
 : @param $p1: one polynomial
 : @param $k: constant
 :)
declare function this:times(
  $p1 as map(xs:string,item()*),
  $k as xs:double
) as map(xs:string,item()*)
{
  this:polynomial(
    reverse( $p1("coefficients")!(.=>z:times($k)) ),
    this:op($p1)
  )
};

(:~
 : ztimes()
 : Multiply polynomial by complex constant 
 :
 : @param $p1: one polynomial
 : @param $kz: complex constant
 :)
declare function this:ztimes(
  $p1 as map(xs:string,item()*),
  $kz as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  this:polynomial(
    reverse( $p1("coefficients")!(.=>z:multiply($kz)) ),
    this:op($p1)
  )
};