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

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

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.

Apr 2023
Status: Bleeding edge

### Imports

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

### Functions

#### `Function: polynomialdeclare function polynomial(\$coefficients as map(xs:string,item()*)*) as map(xs:string,item()*)`

polynomial()
Construct a representation of a polynomial

##### Params
• coefficients as map(xs:string,item()*)*: 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))
##### Returns
• map(xs:string,item()*)
```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 :)
}
}```

#### ```Function: polynomialdeclare function polynomial(\$coefficients as map(xs:string,item()*)*, \$op as xs:string?) as map(xs:string,item()*)```

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

##### Params
• coefficients as map(xs:string,item()*)*: coefficients of polynomial in decreasing order
• op as xs:string?: operator, one of "sin", "cos", "sinh", "cosh" or empty or "none"
##### Returns
• map(xs:string,item()*)
```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)
}```

#### `Function: as-polynomialdeclare function as-polynomial(\$polynomial as map(xs:string,item()*)) as map(xs:string,item()*)`

as-polynomial()
Construct a polynomial with complex coefficents from one with real coefficients.
If the polynomial already is a complex polynomial, just return it

##### Params
• polynomial as map(xs:string,item()*): polynomial
##### Returns
• map(xs:string,item()*): 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))
}```

#### `Function: as-real-polynomialdeclare function as-real-polynomial(\$polynomial as map(xs:string,item()*)) as map(xs:string,item()*)`

as-real-polynomial()
Convert to a polynomial with real coefficients, if possible. Raises UTIL-CAST
if that is not possible.

##### Params
• polynomial as map(xs:string,item()*): polynomial
##### Returns
• map(xs:string,item()*): 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))
}```

#### `Function: is-realdeclare function is-real(\$polynomial as map(xs:string,item()*)) as xs:boolean`

is-real()
Return true if every coefficient is real-valued.

##### Params
• polynomial as map(xs:string,item()*): polynomial
##### Returns
• xs:boolean: 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)
)
}```

#### `Function: coefficientsdeclare function coefficients(\$poly as map(xs:string,item()*)) as map(xs:string,item()*)*`

coefficients()
Accessor for coefficients of polynomial

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• map(xs:string,item()*)*
```declare function this:coefficients(\$poly as map(xs:string,item()*)) as map(xs:string,item()*)*
{
reverse(\$poly("coefficients"))
}```

#### `Function: orderdeclare function order(\$poly as map(xs:string,item()*)) as xs:integer`

order()
Largest exponent. e.g. order(3x³ - 2) = 3

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• xs:integer
```declare function this:order(\$poly as map(xs:string,item()*)) as xs:integer
{
count(\$poly("coefficients")) - 1
}```

#### `Function: opdeclare function op(\$poly as map(xs:string,item()*)) as xs:string?`

op()
Operator.

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• xs:string?
```declare function this:op(\$poly as map(xs:string,item()*)) as xs:string?
{
\$poly("op")
}```

#### ```Function: valuedeclare function value(\$poly as map(xs:string,item()*), \$x as xs:double) as map(xs:string,item()*)```

value()
Compute f(x) where f is a polynomial in x (a double)

##### Params
• poly as map(xs:string,item()*): the polynomial
• x as xs:double the value to compute the function of
##### Returns
• map(xs:string,item()*)
```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) ))
}```

#### ```Function: zvaluedeclare function zvalue(\$poly as map(xs:string,item()*), \$z as map(xs:string,item()*)) as map(xs:string,item()*)```

zvalue()
Compute f(z) where f is a polynomial in z (a complex number)

##### Params
• poly as map(xs:string,item()*): the polynomial
• z as map(xs:string,item()*): the complex value to compute the function of
##### Returns
• map(xs:string,item()*)
```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))
}```

#### ```Function: zvaluevdeclare function zvaluev(\$poly as map(xs:string,item()*), \$z as xs:double*) as xs:double*```

zvaluev()
Compute f(z) where f is a polynomial in z (a complex number vector)

##### Params
• poly as map(xs:string,item()*): the polynomial
• z as xs:double*: the complex value vector to compute the function of
##### Returns
• xs:double*
```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)))
}```

#### `Function: antiderivativedeclare function antiderivative(\$poly as map(xs:string,item()*)) as map(xs:string,item()*)`

antiderivative()
Anti-derivative of a polynomial, which is another polynomial. We use 0 as the
integration constant.

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• map(xs:string,item()*)
```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
)
}```

#### `Function: derivativedeclare function derivative(\$poly as map(xs:string,item()*)) as map(xs:string,item()*)`

derivative()
Derivative of a polynomial, which is another polynomial.

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• map(xs:string,item()*)
```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
)
)
}```

#### ```Function: definite-integraldeclare function definite-integral(\$poly as map(xs:string,item()*), \$a as xs:double, \$b as xs:double) as map(xs:string,item()*)```

definite-integral()
Definite integral of a polynomial (exact)

##### Params
• poly as map(xs:string,item()*): the polynomial
• a as xs:double: start of range of integration
• b as xs:double: end of range of integration
##### Returns
• map(xs:string,item()*)
```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) )
}```

#### `Function: describedeclare function describe(\$poly as map(xs:string,item()*)) as xs:string`

describe()
Describe the polynomial

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• xs:string
```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: functiondeclare function function(\$poly as map(xs:string,item()*)) as function(xs:double) as map(xs:string,item()*)`

function()
Return a function over values that takes values of the polynomial, e.g. for curve plotting

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• function(xs:double)asmap(xs:string,item()*)
```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)
}
}```

#### `Function: zfunctiondeclare function zfunction(\$poly as map(xs:string,item()*)) as function(map(xs:string,item()*)) as map(xs:string,item()*)`

zfunction()
Return a function over complex values that takes values of the polynomial, e.g. for curve plotting

##### Params
• poly as map(xs:string,item()*): the polynomial
##### Returns
• function(map(xs:string,item()*))asmap(xs:string,item()*)
```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)
}
}```

#### ```Function: adddeclare function add(\$p1 as map(xs:string,item()*), \$p2 as map(xs:string,item()*)) as map(xs:string,item()*)```

##### Params
• p1 as map(xs:string,item()*): one polynomial
• p2 as map(xs:string,item()*): another polynomial
##### Returns
• map(xs:string,item()*)
```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)
)
)
}```

#### ```Function: subdeclare function sub(\$p1 as map(xs:string,item()*), \$p2 as map(xs:string,item()*)) as map(xs:string,item()*)```

sub()
Subtract two polynomials

##### Params
• p1 as map(xs:string,item()*): one polynomial
• p2 as map(xs:string,item()*): another polynomial
##### Returns
• map(xs:string,item()*)
```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)
)
)
}```

#### ```Function: multiplydeclare function multiply(\$p1 as map(xs:string,item()*), \$p2 as map(xs:string,item()*)) as map(xs:string,item()*)```

multiply()
Multiply two polynomials

##### Params
• p1 as map(xs:string,item()*): one polynomial
• p2 as map(xs:string,item()*): another polynomial
##### Returns
• map(xs:string,item()*)
```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)
)
)
}```

#### ```Function: timesdeclare function times(\$p1 as map(xs:string,item()*), \$k as xs:double) as map(xs:string,item()*)```

times()
Multiply polynomial by double constant

##### Params
• p1 as map(xs:string,item()*): one polynomial
• k as xs:double: constant
##### Returns
• map(xs:string,item()*)
```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)
)
}```

#### ```Function: ztimesdeclare function ztimes(\$p1 as map(xs:string,item()*), \$kz as map(xs:string,item()*)) as map(xs:string,item()*)```

ztimes()
Multiply polynomial by complex constant

##### Params
• p1 as map(xs:string,item()*): one polynomial
• kz as map(xs:string,item()*): complex constant
##### Returns
• map(xs:string,item()*)
```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.
:
: @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 two polynomials
:
: @param \$p1: one polynomial
: @param \$p2: another polynomial
:)
\$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)
)
};```