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

http://mathling.com/type/polynomial

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.

May 2022
Status: Active

### Imports

http://mathling.com/core/utilities
```import module namespace util="http://mathling.com/core/utilities"
at "../core/utilities.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 xs:double*) as map(xs:string,item()*)`

polynomial()
Construct a representation of a polynomial

##### Params
• coefficients as xs:double*: coefficients of polynomial in decreasing order, e.g. 3x³ - 2 = (3, 0, 0, -2)
##### Returns
• map(xs:string,item()*)
```declare function this:polynomial(\$coefficients as xs:double*) as map(xs:string,item()*)
{
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 :)
}
)
}```

#### ```Function: polynomialdeclare function polynomial(\$coefficients as xs:double*, \$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( (3, 0, 0, -2), "sin" ) = 3sin(x)³ - 2 =

##### Params
• coefficients as xs:double*: 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 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)
}```

#### `Function: is-constantdeclare function is-constant(\$poly as map(xs:string,item()*)) as xs:boolean`

is-constant()

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

#### `Function: as-realdeclare function as-real(\$poly as map(xs:string,item()*)) as xs:double`

as-real()
Convert a constant polynomial to the constant value. Raises UTIL-CAST if
it is not a constant polynomial.

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

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

coefficients()
Accessor for coefficients of polynomial

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

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

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

#### ```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: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))
}```

#### ```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: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))
}```

#### `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 \$k * \$c div \$i), 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 0
else reverse(for \$c at \$i in tail(\$poly("coefficients")) return \$k * \$c * \$i),
\$op
)
)
}```

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

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
• xs:double
```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)
}```

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

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

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

#### `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))
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)
),
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))
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)
),
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))
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*\$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)
)
}
)
),
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 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")!(. * \$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.
:
: @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()*)
{
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)
};

(:~
: 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
:======================================================================:)

(:~
:
: @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))
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)
),
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))
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)
),
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))
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*\$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)
)
}
)
),
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)
)
};```