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

Imports

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 errors="http://mathling.com/core/errors"
       at "../core/errors.xqy"

Functions

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

declare function this:polynomial(
  $coefficients as xs:double*,
  $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:is-constant(
  $poly as map(xs:string,item()*)
) as xs:boolean
{
  this:order($poly) eq 0
}

declare function this:as-real(
  $poly as map(xs:string,item()*)
) as xs:double
{
  if (this:is-constant($poly)) then (
    $poly("coefficients")
  ) else (
    errors:error("UTIL-CAST", ($poly, "xs:double"))
  )
}

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

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:times($c))
  case "cos" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cos($z), $i - 1)=>z:times($c))
  case "sinh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:sinh($z), $i - 1)=>z:times($c))
  case "cosh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cosh($z), $i - 1)=>z:times($c))
  default return
    z:sum(for $c at $i in $poly("coefficients") return z:pow($z, $i - 1)=>z:times($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:times($c))
  case "cos" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cos($z), $i - 1)=>zv:times($c))
  case "sinh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:sinh($z), $i - 1)=>zv:times($c))
  case "cosh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cosh($z), $i - 1)=>zv:times($c))
  default return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow($z, $i - 1)=>zv:times($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 $k * $c div $i), 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 0
      else reverse(for $c at $i in tail($poly("coefficients")) return $k * $c * $i),
      $op
    )
  )
}

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

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

declare function this:function(
  $poly as map(xs:string,item()*)
) as function(xs:double) as xs:double
{
  function($v as xs:double) as xs:double {
    $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 + ($c2[$i],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")!(-.)
  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 + ($c2[$i],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]*$d),
          function($s as xs:double*, $i as xs:integer) as xs:double* {
            let $sub := (
              for $j in $o1 - $i + 1 to $o1 - 1 return 0,
              for $d in $c1 return $c2[$i]*$d
            )
            return (
              for $j in 1 to count($sub) return $sub[$j] + ($s[$j],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")!(. * $k) ),
    this:op($p1)
  )
}

Original Source Code

xquery version "3.1";
(:~
 : Polynomials with real 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. Polynomials are also allowed to be in sin(x), cos(x),
 : sinh(x), and cosh(x) because derivatives and anti-derivatives of these are
 : still representable as polynomials in this way.
 :
 : Copyright© Mary Holstege 2022-2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @since May 2022
 : @custom:Status Active
 :)
module namespace this="http://mathling.com/type/polynomial"; 

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";

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 = (3, 0, 0, -2)
 :)
declare function this:polynomial($coefficients as xs:double*) as map(xs:string,item()*)
{
  (: Prune leading zeros :)
  let $useful :=
    let $first-nonzero := head(for $c at $i in $coefficients where $c!=0 return $i)
    return $coefficients[position() >= $first-nonzero]
  (: Make sure we have something :)
  let $useful :=
    if (empty($useful)) then (0) else $useful
  return (
    map {
      "kind": "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( (3, 0, 0, -2), "sin" ) = 3sin(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 xs:double*,
  $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)
};

(:~
 : is-constant()
 : @param $poly: the polynomial
 : @return true() if this is a 0 order (i.e constant) polynomial
 :)
declare function this:is-constant(
  $poly as map(xs:string,item()*)
) as xs:boolean
{
  this:order($poly) eq 0
};

(:~
 : as-real()
 : Convert a constant polynomial to the constant value. Raises UTIL-CAST if
 : it is not a constant polynomial.
 : @param $poly: the polynomial
 : @return the constant
 :)
declare function this:as-real(
  $poly as map(xs:string,item()*)
) as xs:double
{
  if (this:is-constant($poly)) then (
    $poly("coefficients")
  ) else (
    errors:error("UTIL-CAST", ($poly, "xs:double"))
  )
};

(:~
 : coefficients()
 : Accessor for coefficients of polynomial
 :
 : @param $poly: the polynomial
 :)
declare function this:coefficients($poly as map(xs:string,item()*)) as xs:double*
{
  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
 :
 : @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 xs:double
{
  switch ($poly("op"))
  case "sin" return
    sum(for $c at $i in $poly("coefficients") return math:pow(math:sin($x), $i - 1) * $c)
  case "cos" return
    sum(for $c at $i in $poly("coefficients") return math:pow(math:cos($x), $i - 1) * $c)
  case "sinh" return
    sum(for $c at $i in $poly("coefficients") return math:pow(util:sinh($x), $i - 1) * $c)
  case "cosh" return
    sum(for $c at $i in $poly("coefficients") return math:pow(util:cosh($x), $i - 1) * $c)
  default return
    sum(for $c at $i in $poly("coefficients") return math:pow($x, $i - 1) * $c)
};

(:~
 : 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:times($c))
  case "cos" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cos($z), $i - 1)=>z:times($c))
  case "sinh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:sinh($z), $i - 1)=>z:times($c))
  case "cosh" return
    z:sum(for $c at $i in $poly("coefficients") return z:pow(z:cosh($z), $i - 1)=>z:times($c))
  default return
    z:sum(for $c at $i in $poly("coefficients") return z:pow($z, $i - 1)=>z:times($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:times($c))
  case "cos" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cos($z), $i - 1)=>zv:times($c))
  case "sinh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:sinh($z), $i - 1)=>zv:times($c))
  case "cosh" return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow(zv:cosh($z), $i - 1)=>zv:times($c))
  default return
    zv:sum(for $c at $i in $poly("coefficients") return zv:pow($z, $i - 1)=>zv:times($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 $k * $c div $i), 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 0
      else reverse(for $c at $i in tail($poly("coefficients")) return $k * $c * $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 xs:double
{
  let $anti := this:antiderivative($poly)
  return $anti=>this:value($b) - $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 xs:string($poly("coefficients"))
  else (
    let $x := if (exists($poly("op"))) then $poly("op")||"(x)" else "x"
    return (
      string-join(
        reverse(
          for $c at $i in $poly("coefficients") return (
            if ($c = 0) then ()
            else (
              let $xton :=
                if ($i - 1 = 0) then ""
                else if ($i - 1 = 1) then $x
                else $x||"^"||($i - 1)
              let $k :=
                if ($c = 1 and $i > 1) then ""
                else if ($c = -1 and $i > 1) then "-"
                else xs:string($c)
              return $k||$xton (: c*x^(i-1) :)
            )
          )
        ),
        " + "
      )
    )
  )
};

(:~
 : function()
 : Return a function 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 xs:double
{
  function($v as xs:double) as xs:double {
    $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 + ($c2[$i],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")!(-.)
  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 + ($c2[$i],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]*$d),
          function($s as xs:double*, $i as xs:integer) as xs:double* {
            let $sub := (
              for $j in $o1 - $i + 1 to $o1 - 1 return 0,
              for $d in $c1 return $c2[$i]*$d
            )
            return (
              for $j in 1 to count($sub) return $sub[$j] + ($s[$j],0)[1]
            )
          }
        )
      ),
      this:op($p1)
    )
  )
};

(:~
 : times()
 : Multiply polynomial by 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")!(. * $k) ),
    this:op($p1)
  )
};