# http://mathling.com/core/roots  library module

http://mathling.com/core/roots

Roots of linear, quadratic, and cubic equations.

December 2021
Status: Stable, with bleeding edge handling of complex polynomials

### 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/config
```import module namespace config="http://mathling.com/core/config"
at "../core/config.xqy"```
http://mathling.com/type/polynomial/complex
```import module namespace zpoly="http://mathling.com/type/polynomial/complex"
at "../types/cpolynomial.xqy"```
http://mathling.com/core/errors
```import module namespace errors="http://mathling.com/core/errors"
at "../core/errors.xqy"```

### Functions

#### ```Function: linear-rootdeclare function linear-root(\$a as xs:double, \$b as xs:double) as xs:double?```

linear-root()
Return the root of the linear equation of the form ax + b = 0
If a is 0 there is no solution unless b = 0 in which case there are
an infinite number of solutions, so callers should check on what to
do in that case.

##### Params
• a as xs:double: coefficient of linear equation ax + b = 0
• b as xs:double: coefficient of linear equation ax + b = 0
##### Returns
• xs:double?: root
```declare function this:linear-root(
\$a as xs:double,
\$b as xs:double
) as xs:double?
{
if (\$a = 0) then ()
else -\$b div \$a
}```

#### ```Function: linear-root-complexdeclare function linear-root-complex(\$a as map(xs:string,item()*), \$b as map(xs:string,item()*)) as map(xs:string,item()*)?```

linear-root-complex()
Same as linear-root() but over complex coefficients.

##### Params
• a as map(xs:string,item()*)
• b as map(xs:string,item()*)
##### Returns
• map(xs:string,item()*)?
```declare function this:linear-root-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*)
) as map(xs:string,item()*)?
{
if (z:eq(\$a,0)) then ()
else z:minus(\$b)=>z:divide(\$a)
}```

#### ```Function: linear-real-root-complexdeclare function linear-real-root-complex(\$a as map(xs:string,item()*), \$b as map(xs:string,item()*)) as xs:double?```

linear-real-root-complex()
Same as linear-root() but over complex coefficients and only returns a value
if it has no imaginary part (within epsilon).

##### Params
• a as map(xs:string,item()*)
• b as map(xs:string,item()*)
##### Returns
• xs:double?
```declare function this:linear-real-root-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*)
) as xs:double?
{
if (z:eq(\$a,0)) then ()
else (
let \$val := z:minus(\$b)=>z:divide(\$a)
return if (z:is-real(\$val)) then \$val=>z:as-real() else ()
)
}```

#### ```Function: quadratic-real-rootsdeclare function quadratic-real-roots(\$a as xs:double, \$b as xs:double, \$c as xs:double) as xs:double*```

Return all the real-valued (within ε) roots of the quadratic equation
of the form ax² + bx + c = 0

##### Params
• a as xs:double: coefficient of quadratic equation ax² + bx + c = 0
• b as xs:double: coefficient of quadratic equation ax² + bx + c = 0
• c as xs:double: coefficient of quadratic equation ax² + bx + c = 0
##### Returns
• xs:double*: real roots, if any
```declare function this:quadratic-real-roots(
\$a as xs:double,
\$b as xs:double,
\$c as xs:double
) as xs:double*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (\$a = 0) then (
this:linear-root(\$b, \$c)
) else if (\$b*\$b >= 4*\$a*\$c) then (
(-\$b + math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a),
(-\$b - math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a)
) else (
(: √ negative => not a real root :)
)
}```

#### ```Function: quadratic-real-roots-complexdeclare function quadratic-real-roots-complex(\$a as map(xs:string,item()*), \$b as map(xs:string,item()*), \$c as map(xs:string,item()*)) as xs:double*```

Same as quadratic-real-roots() except over complex coefficients.

##### Params
• a as map(xs:string,item()*)
• b as map(xs:string,item()*)
• c as map(xs:string,item()*)
##### Returns
• xs:double*
```declare function this:quadratic-real-roots-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*)
) as xs:double*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (z:eq(\$a,0)) then (
this:linear-real-root-complex(\$b, \$c)
) else (
let \$sqrt := (: √(b²-4ac) :)
z:sqrt(
z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(4) )
)
let \$vals := (
(z:minus(\$b)=>z:sub(\$sqrt))=>z:divide(\$a=>z:times(2))
)
for \$val in \$vals where z:is-real(\$val) return \$val=>z:as-real()
)
}```

#### ```Function: quadratic-rootsdeclare function quadratic-roots(\$a as xs:double, \$b as xs:double, \$c as xs:double) as map(xs:string,item()*)*```

Return all the roots of quadratic equation of the form ax² + bx + c = 0
Values returned as complex numbers.

##### Params
• a as xs:double: coefficient of quadratic equation ax² + bx + c = 0
• b as xs:double: coefficient of quadratic equation ax² + bx + c = 0
• c as xs:double: coefficient of quadratic equation ax² + bx + c = 0
##### Returns
• map(xs:string,item()*)*: roots
```declare function this:quadratic-roots(
\$a as xs:double,
\$b as xs:double,
\$c as xs:double
) as map(xs:string,item()*)*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (\$a = 0) then (
this:linear-root(\$b, \$c)!z:as-complex(.)
) else if (\$b*\$b >= 4*\$a*\$c) then (
z:as-complex((-\$b + math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a)),
z:as-complex((-\$b - math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a))
) else (
let \$sqrt := z:sqrt(z:as-complex(\$b*\$b - 4*\$a*\$c))
return (
(z:as-complex(-\$b)=>z:sub(\$sqrt))=>z:times(1 div (2*\$a))
)
)
}```

#### ```Function: quadratic-roots-complexdeclare function quadratic-roots-complex(\$a as map(xs:string,item()*), \$b as map(xs:string,item()*), \$c as map(xs:string,item()*)) as map(xs:string,item()*)*```

Same as quadratic-roots() but over complex coefficients.

##### Params
• a as map(xs:string,item()*)
• b as map(xs:string,item()*)
• c as map(xs:string,item()*)
##### Returns
• map(xs:string,item()*)*
```declare function this:quadratic-roots-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*)
) as map(xs:string,item()*)*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (z:eq(\$a,0)) then (
this:linear-root-complex(\$b, \$c)
) else (
let \$sqrt := (: √(b²-4ac) :)
z:sqrt(
z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(4) )
)
return (
(z:minus(\$b)=>z:sub(\$sqrt))=>z:divide(\$a=>z:times(2))
)
)
}```

#### ```Function: cubic-real-rootsdeclare function cubic-real-roots(\$a as xs:double, \$b as xs:double, \$c as xs:double, \$d as xs:double) as xs:double*```

cubic-real-roots()
Return all the real-valued (within ε) roots of cubic equation
of the form ax³ + bx² + cx + d = 0

##### Params
• a as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
• b as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
• c as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
• d as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
##### Returns
• xs:double*: real roots, if any
```declare function this:cubic-real-roots(
\$a as xs:double,
\$b as xs:double,
\$c as xs:double,
\$d as xs:double
) as xs:double*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (\$a = 0) then (
) else (
let \$p := \$b*\$b - 2*\$a*\$c
let \$q := 2*\$b*\$b*\$b - 9*\$a*\$b*\$c + 27*\$a*\$a*\$d
return (
if (\$p = \$q and \$p = 0) then (
-\$b div (3*\$a)
) else (
for \$root in this:cubic-roots(\$a, \$b, \$c, \$d)
where z:is-real(\$root)
return \$root=>z:as-real()
)
)
)
}```

#### ```Function: cubic-real-roots-complexdeclare function cubic-real-roots-complex(\$a as map(xs:string,item()*), \$b as map(xs:string,item()*), \$c as map(xs:string,item()*), \$d as map(xs:string,item()*)) as xs:double*```

cubic-real-roots-complex()
Same as cubic-real-roots() but over complex coefficients.

##### Params
• a as map(xs:string,item()*)
• b as map(xs:string,item()*)
• c as map(xs:string,item()*)
• d as map(xs:string,item()*)
##### Returns
• xs:double*
```declare function this:cubic-real-roots-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*),
\$d as map(xs:string,item()*)
) as xs:double*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (z:eq(\$a,0)) then (
) else (
let \$p := z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(2) ) (: √(b²-4ac) :)
let \$q := (: 2b³-9abc+27a²d :)
(z:pow(\$b,3)=>z:times(2))=>z:sub(
\$a=>z:multiply(\$b)=>z:multiply(\$c)=>z:times(9)
\$a=>z:multiply(\$a)=>z:multiply(\$d)=>z:times(27)
)
return (
if (z:eq(\$p,0) and z:eq(\$q,0)) then (
let \$root := z:minus(\$b)=>z:divide(\$a=>z:times(3)) (: -b/3a :)
where z:is-real(\$root)
return \$root=>z:as-real()
) else (
for \$root in this:cubic-roots-complex(\$a, \$b, \$c, \$d)
where z:is-real(\$root)
return \$root=>z:as-real()
)
)
)
}```

#### ```Function: cubic-rootsdeclare function cubic-roots(\$a as xs:double, \$b as xs:double, \$c as xs:double, \$d as xs:double) as map(xs:string,item()*)*```

cubic-roots()
Return all the roots of cubic equation of the form ax³ + bx² + cx + d = 0
Values returned as complex numbers.

##### Params
• a as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
• b as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
• c as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
• d as xs:double: coefficient of cubic equation ax³ + bx² + cx + d = 0
##### Returns
• map(xs:string,item()*)*: roots
```declare function this:cubic-roots(
\$a as xs:double,
\$b as xs:double,
\$c as xs:double,
\$d as xs:double
) as map(xs:string,item()*)*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (\$a = 0) then (
) else (
let \$p := \$b*\$b - 3*\$a*\$c
let \$q := 2*\$b*\$b*\$b - 9*\$a*\$b*\$c + 27*\$a*\$a*\$d
let \$qz := z:as-complex(\$q)
let \$pz := z:as-complex(\$p)
return if (\$p = \$q and \$p = 0) then z:as-complex(-\$b div (3*\$a)) else (
let \$ζ0 := z:as-complex(1)
let \$ζ1 := z:complex(-1 div 2, math:sqrt(3) div 2)
let \$ζ2 := z:multiply(\$ζ1, \$ζ1)
(: r = cbrt((p ± sqrt(q² - 4p³))/2) :)
let \$r :=
let \$val :=
try {
util:cbrt(
(\$q + math:sqrt(\$q*\$q - 4*\$p*\$p*\$p)) div 2
)
} catch * {xs:double("NaN")}
return (
if (\$val=\$val) then z:as-complex(\$val) (: i.e. not NaN :)
else (
z:cbrt(
z:sqrt(
z:multiply(\$qz,\$qz)=>
z:sub(
z:multiply(z:multiply(z:times(\$pz, 4), \$pz), \$pz)
)
)
)=>z:times(0.5)
)
)
)
let \$r :=
if (z:pr(\$r)=0 and z:pi(\$r)=0) then (
let \$val :=
try {
util:cbrt(
(\$q - math:sqrt(\$q*\$q - 4*\$p*\$p*\$p)) div 2
)
} catch * {xs:double("NaN")}
return (
if (\$val=\$val) then z:as-complex(\$val) (: i.e. not NaN :)
else (
z:cbrt(
\$qz=>z:sub(
z:sqrt(
z:multiply(\$qz,\$qz)=>
z:sub(
z:multiply(z:multiply(z:times(\$pz, 4), \$pz), \$pz)
)
)
)=>z:times(0.5)
)
)
)
) else \$r
let \$rζ0 := \$r=>z:multiply(\$ζ0)
let \$rζ1 := \$r=>z:multiply(\$ζ1)
let \$rζ2 := \$r=>z:multiply(\$ζ2)
let \$bz := z:as-complex(\$b)
(: x = (-1/(3a))(b + (ζ^k)r + p/((ζ^k)r) for k=0 to 2 :)
for \$rζk in (\$rζ0, \$rζ1, \$rζ2) return (
)
)
)
}```

#### ```Function: cubic-roots-complexdeclare function cubic-roots-complex(\$a as map(xs:string,item()*), \$b as map(xs:string,item()*), \$c as map(xs:string,item()*), \$d as map(xs:string,item()*)) as map(xs:string,item()*)*```

cubic-roots-complex()
Same as cubic-roots() but over complex coefficients

##### Params
• a as map(xs:string,item()*)
• b as map(xs:string,item()*)
• c as map(xs:string,item()*)
• d as map(xs:string,item()*)
##### Returns
• map(xs:string,item()*)*
```declare function this:cubic-roots-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*),
\$d as map(xs:string,item()*)
) as map(xs:string,item()*)*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (z:eq(\$a,0)) then (
) else (
let \$p := z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(3) ) (: b²-3ac :)
let \$q := (: 2b³-9abc+27a²d :)
(z:pow(\$b,3)=>z:times(2))=>z:sub(
\$a=>z:multiply(\$b)=>z:multiply(\$c)=>z:times(9)
\$a=>z:multiply(\$a)=>z:multiply(\$d)=>z:times(27)
)
return (
if (z:eq(\$p,0) and z:eq(\$q,0)) then (
z:minus(\$b)=>z:divide(\$a=>z:times(3)) (: -b/3a :)
) else (
let \$ζ0 := z:as-complex(1)
let \$ζ1 := z:complex(-1 div 2, math:sqrt(3) div 2)
let \$ζ2 := z:multiply(\$ζ1, \$ζ1)
(: r = cbrt((p ± sqrt(q² - 4p³))/2) :)
let \$r := (
z:cbrt(
z:sqrt(
z:multiply(\$q,\$q)=>
z:sub(
z:multiply(z:multiply(z:times(\$p, 4), \$p), \$p)
)
)
)=>z:times(0.5)
)
)
let \$r := (
if (z:eq(\$r, 0)) then (
z:cbrt(
\$q=>z:sub(
z:sqrt(
z:multiply(\$q,\$q)=>
z:sub(
z:multiply(z:multiply(z:times(\$p, 4), \$p), \$p)
)
)
)=>z:times(0.5)
)
) else \$r
)
let \$rζ0 := \$r=>z:multiply(\$ζ0)
let \$rζ1 := \$r=>z:multiply(\$ζ1)
let \$rζ2 := \$r=>z:multiply(\$ζ2)
(: x = (-1/(3a))(b + (ζ^k)r + p/((ζ^k)r) for k=0 to 2 :)
for \$rζk in (\$rζ0, \$rζ1, \$rζ2) return (
)
)
)
)
}```

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

real-roots()
Return all the real-valued (within ε) roots of the polynomial.

##### Params
• polynomial as map(xs:string,item()*): the polynomial
##### Returns
• xs:double*: real roots, if any
```declare function this:real-roots(
\$polynomial as map(xs:string,item()*)
) as xs:double*
{
let \$op := poly:op(\$polynomial)
let \$raw-roots := (
if (util:kind(\$polynomial)="polynomial") then (
let \$cfs as xs:double* := poly:coefficients(\$polynomial)
return switch (poly:order(\$polynomial))
case 1 return this:linear-root(\$cfs[1], \$cfs[2])
case 2 return this:quadratic-real-roots(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-real-roots(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
) else if (util:kind(\$polynomial)="complex-polynomial") then (
let \$cfs as map(xs:string,item()*)* := zpoly:coefficients(\$polynomial)
return switch (zpoly:order(\$polynomial))
case 1 return this:linear-real-root-complex(\$cfs[1], \$cfs[2])
case 2 return this:quadratic-real-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-real-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
) else (
errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
)
)
return (
if (empty(\$op)) then (
\$raw-roots
) else if (\$op="sin") then (
(:
: We know values of sin(x) for which poly(sin(x)) = 0
: Values of x are therefor asin of those
: e.g. sin²(x) - 1 = 0 => sin(x) = {1, -1} => x = {asin(1), asin(-1)}
:)
for \$v in \$raw-roots return math:asin(\$v)
) else if (\$op="cos") then (
(:
: We know values of cos(x) for which poly(cos(x)) = 0
: Values of x are therefor acos of those
: e.g. cos²(x) - 1 = 0 => cos(x) = {1, -1} => x = {acos(1), acos(-1)}
:)
for \$v in \$raw-roots return math:acos(\$v)
) else if (\$op="sinh") then (
for \$v in \$raw-roots return util:asinh(\$v)
) else if (\$op="cosh") then (
for \$v in \$raw-roots return util:acosh(\$v)
) else (
errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
)
)
}```
##### Errors

UTIL-CANNOT if the polynomial is order 0 or more than 3

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

roots()
Return all the roots of the polynomial.

##### Params
• polynomial as map(xs:string,item()*): the polynomial
##### Returns
• map(xs:string,item()*)*: roots
```declare function this:roots(
\$polynomial as map(xs:string,item()*)
) as map(xs:string,item()*)*
{
let \$op := poly:op(\$polynomial)
let \$raw-roots := (
if (util:kind(\$polynomial)="polynomial") then (
let \$cfs as xs:double* := poly:coefficients(\$polynomial)
return switch (poly:order(\$polynomial))
case 1 return this:linear-root(\$cfs[1], \$cfs[2])!z:as-complex(.)
case 2 return this:quadratic-roots(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-roots(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("roots", \$polynomial))
) else if (util:kind(\$polynomial)="complex-polynomial") then (
let \$cfs as map(xs:string,item()*)* := zpoly:coefficients(\$polynomial)
return switch (zpoly:order(\$polynomial))
case 1 return this:linear-root-complex(\$cfs[1], \$cfs[2])
case 2 return this:quadratic-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("roots", \$polynomial))
) else (
errors:error("UTIL-CANNOT", ("roots", \$polynomial))
)
)
return (
if (empty(\$op)) then (
\$raw-roots
) else if (\$op="sin") then (
(:
: We know values of sin(x) for which poly(sin(x)) = 0
: Values of x are therefor asin of those
: e.g. sin²(x) - 1 = 0 => sin(x) = {1, -1} => x = {asin(1), asin(-1)}
:)
for \$v in \$raw-roots return z:asin(\$v)
) else if (\$op="cos") then (
(:
: We know values of cos(x) for which poly(cos(x)) = 0
: Values of x are therefor acos of those
: e.g. cos²(x) - 1 = 0 => cos(x) = {1, -1} => x = {acos(1), acos(-1)}
:)
for \$v in \$raw-roots return z:acos(\$v)
) else if (\$op="sinh") then (
for \$v in \$raw-roots return z:asinh(\$v)
) else if (\$op="cosh") then (
for \$v in \$raw-roots return z:acosh(\$v)
) else (
errors:error("UTIL-CANNOT", ("roots", \$polynomial))
)
)
}```
##### Errors

UTIL-CANNOT if the polynomial is order 0 or more than 3

### Original Source Code

```xquery version "3.1";
(:~
: Roots of linear, quadratic, and cubic equations.
:
: @since December 2021
: @custom:Status Stable, with bleeding edge handling of complex polynomials
:)
module namespace this="http://mathling.com/core/roots";

import module namespace config="http://mathling.com/core/config"
at "../core/config.xqy";
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 poly="http://mathling.com/type/polynomial"
at "../types/polynomial.xqy";
import module namespace zpoly="http://mathling.com/type/polynomial/complex"
at "../types/cpolynomial.xqy";

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

(:
: Generic solution to linear equation ax + b = 0
: x = -b/a
:)

(:~
: linear-root()
: Return the root of the linear equation of the form ax + b = 0
: If a is 0 there is no solution unless b = 0 in which case there are
: an infinite number of solutions, so callers should check on what to
: do in that case.
:
: @param \$a: coefficient of linear equation ax + b = 0
: @param \$b: coefficient of linear equation ax + b = 0
: @return root
:)
declare function this:linear-root(
\$a as xs:double,
\$b as xs:double
) as xs:double?
{
if (\$a = 0) then ()
else -\$b div \$a
};

(:~
: linear-root-complex()
: Same as linear-root() but over complex coefficients.
:)
declare function this:linear-root-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*)
) as map(xs:string,item()*)?
{
if (z:eq(\$a,0)) then ()
else z:minus(\$b)=>z:divide(\$a)
};

(:~
: linear-real-root-complex()
: Same as linear-root() but over complex coefficients and only returns a value
: if it has no imaginary part (within epsilon).
:)
declare function this:linear-real-root-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*)
) as xs:double?
{
if (z:eq(\$a,0)) then ()
else (
let \$val := z:minus(\$b)=>z:divide(\$a)
return if (z:is-real(\$val)) then \$val=>z:as-real() else ()
)
};

(:
: Generic solution to quadratic equation ax² + bx + c = 0
: x = (-b ± √(b² - 4ac))/2a
:)

(:~
: Return all the real-valued (within ε) roots of the quadratic equation
: of the form ax² + bx + c = 0
:
: @param \$a: coefficient of quadratic equation ax² + bx + c = 0
: @param \$b: coefficient of quadratic equation ax² + bx + c = 0
: @param \$c: coefficient of quadratic equation ax² + bx + c = 0
: @return real roots, if any
:)
\$a as xs:double,
\$b as xs:double,
\$c as xs:double
) as xs:double*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (\$a = 0) then (
this:linear-root(\$b, \$c)
) else if (\$b*\$b >= 4*\$a*\$c) then (
(-\$b + math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a),
(-\$b - math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a)
) else (
(: √ negative => not a real root :)
)
};

(:~
: Same as quadratic-real-roots() except over complex coefficients.
:)
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*)
) as xs:double*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (z:eq(\$a,0)) then (
this:linear-real-root-complex(\$b, \$c)
) else (
let \$sqrt := (: √(b²-4ac) :)
z:sqrt(
z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(4) )
)
let \$vals := (
(z:minus(\$b)=>z:sub(\$sqrt))=>z:divide(\$a=>z:times(2))
)
for \$val in \$vals where z:is-real(\$val) return \$val=>z:as-real()
)
};

(:~
: Return all the roots of quadratic equation of the form ax² + bx + c = 0
: Values returned as complex numbers.
:
: @param \$a: coefficient of quadratic equation ax² + bx + c = 0
: @param \$b: coefficient of quadratic equation ax² + bx + c = 0
: @param \$c: coefficient of quadratic equation ax² + bx + c = 0
: @return roots
:)
\$a as xs:double,
\$b as xs:double,
\$c as xs:double
) as map(xs:string,item()*)*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (\$a = 0) then (
this:linear-root(\$b, \$c)!z:as-complex(.)
) else if (\$b*\$b >= 4*\$a*\$c) then (
z:as-complex((-\$b + math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a)),
z:as-complex((-\$b - math:sqrt(\$b*\$b - 4*\$a*\$c)) div (2*\$a))
) else (
let \$sqrt := z:sqrt(z:as-complex(\$b*\$b - 4*\$a*\$c))
return (
(z:as-complex(-\$b)=>z:sub(\$sqrt))=>z:times(1 div (2*\$a))
)
)
};

(:~
: Same as quadratic-roots() but over complex coefficients.
:)
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*)
) as map(xs:string,item()*)*
{
(: If a = 0 this is a linear equation, use the simpler math :)
if (z:eq(\$a,0)) then (
this:linear-root-complex(\$b, \$c)
) else (
let \$sqrt := (: √(b²-4ac) :)
z:sqrt(
z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(4) )
)
return (
(z:minus(\$b)=>z:sub(\$sqrt))=>z:divide(\$a=>z:times(2))
)
)
};

(:
: Per: https://en.wikipedia.org/wiki/Cubic_equation
: Generic solution to cubic equation ax³ + bx² + cx + d = 0
:
: Δ0 = b² - 3ac
: Δ1 = 2b³ - 9abc + 27a²d
: C = ((Δ1 ± √(Δ1² - 4Δ0³))/2)^(1/3)
:
: x[k] = -(1/(3a))(b + (ζ^k)C + Δ0/((ζ^k)C)) for k in 0 to 2
:
: ζ = (-1 + √(-3))/2 = (-1 + i√3)/2 = -1/2 + i√3/2
:)

(:~
: cubic-real-roots()
: Return all the real-valued (within ε) roots of cubic equation
: of the form ax³ + bx² + cx + d = 0
:
: @param \$a: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @param \$b: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @param \$c: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @param \$d: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @return real roots, if any
:)
declare function this:cubic-real-roots(
\$a as xs:double,
\$b as xs:double,
\$c as xs:double,
\$d as xs:double
) as xs:double*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (\$a = 0) then (
) else (
let \$p := \$b*\$b - 2*\$a*\$c
let \$q := 2*\$b*\$b*\$b - 9*\$a*\$b*\$c + 27*\$a*\$a*\$d
return (
if (\$p = \$q and \$p = 0) then (
-\$b div (3*\$a)
) else (
for \$root in this:cubic-roots(\$a, \$b, \$c, \$d)
where z:is-real(\$root)
return \$root=>z:as-real()
)
)
)
};

(:~
: cubic-real-roots-complex()
: Same as cubic-real-roots() but over complex coefficients.
:)
declare function this:cubic-real-roots-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*),
\$d as map(xs:string,item()*)
) as xs:double*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (z:eq(\$a,0)) then (
) else (
let \$p := z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(2) ) (: √(b²-4ac) :)
let \$q := (: 2b³-9abc+27a²d :)
(z:pow(\$b,3)=>z:times(2))=>z:sub(
\$a=>z:multiply(\$b)=>z:multiply(\$c)=>z:times(9)
\$a=>z:multiply(\$a)=>z:multiply(\$d)=>z:times(27)
)
return (
if (z:eq(\$p,0) and z:eq(\$q,0)) then (
let \$root := z:minus(\$b)=>z:divide(\$a=>z:times(3)) (: -b/3a :)
where z:is-real(\$root)
return \$root=>z:as-real()
) else (
for \$root in this:cubic-roots-complex(\$a, \$b, \$c, \$d)
where z:is-real(\$root)
return \$root=>z:as-real()
)
)
)
};

(:~
: cubic-roots()
: Return all the roots of cubic equation of the form ax³ + bx² + cx + d = 0
: Values returned as complex numbers.
:
: @param \$a: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @param \$b: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @param \$c: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @param \$d: coefficient of cubic equation ax³ + bx² + cx + d = 0
: @return roots
:)
declare function this:cubic-roots(
\$a as xs:double,
\$b as xs:double,
\$c as xs:double,
\$d as xs:double
) as map(xs:string,item()*)*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (\$a = 0) then (
) else (
let \$p := \$b*\$b - 3*\$a*\$c
let \$q := 2*\$b*\$b*\$b - 9*\$a*\$b*\$c + 27*\$a*\$a*\$d
let \$qz := z:as-complex(\$q)
let \$pz := z:as-complex(\$p)
return if (\$p = \$q and \$p = 0) then z:as-complex(-\$b div (3*\$a)) else (
let \$ζ0 := z:as-complex(1)
let \$ζ1 := z:complex(-1 div 2, math:sqrt(3) div 2)
let \$ζ2 := z:multiply(\$ζ1, \$ζ1)
(: r = cbrt((p ± sqrt(q² - 4p³))/2) :)
let \$r :=
let \$val :=
try {
util:cbrt(
(\$q + math:sqrt(\$q*\$q - 4*\$p*\$p*\$p)) div 2
)
} catch * {xs:double("NaN")}
return (
if (\$val=\$val) then z:as-complex(\$val) (: i.e. not NaN :)
else (
z:cbrt(
z:sqrt(
z:multiply(\$qz,\$qz)=>
z:sub(
z:multiply(z:multiply(z:times(\$pz, 4), \$pz), \$pz)
)
)
)=>z:times(0.5)
)
)
)
let \$r :=
if (z:pr(\$r)=0 and z:pi(\$r)=0) then (
let \$val :=
try {
util:cbrt(
(\$q - math:sqrt(\$q*\$q - 4*\$p*\$p*\$p)) div 2
)
} catch * {xs:double("NaN")}
return (
if (\$val=\$val) then z:as-complex(\$val) (: i.e. not NaN :)
else (
z:cbrt(
\$qz=>z:sub(
z:sqrt(
z:multiply(\$qz,\$qz)=>
z:sub(
z:multiply(z:multiply(z:times(\$pz, 4), \$pz), \$pz)
)
)
)=>z:times(0.5)
)
)
)
) else \$r
let \$rζ0 := \$r=>z:multiply(\$ζ0)
let \$rζ1 := \$r=>z:multiply(\$ζ1)
let \$rζ2 := \$r=>z:multiply(\$ζ2)
let \$bz := z:as-complex(\$b)
(: x = (-1/(3a))(b + (ζ^k)r + p/((ζ^k)r) for k=0 to 2 :)
for \$rζk in (\$rζ0, \$rζ1, \$rζ2) return (
)
)
)
};

(:~
: cubic-roots-complex()
: Same as cubic-roots() but over complex coefficients
:)
declare function this:cubic-roots-complex(
\$a as map(xs:string,item()*),
\$b as map(xs:string,item()*),
\$c as map(xs:string,item()*),
\$d as map(xs:string,item()*)
) as map(xs:string,item()*)*
{
(: If a = 0 this is a quadratic equation, use the simpler math :)
if (z:eq(\$a,0)) then (
) else (
let \$p := z:multiply(\$b,\$b)=>z:sub( z:multiply(\$a,\$c)=>z:times(3) ) (: b²-3ac :)
let \$q := (: 2b³-9abc+27a²d :)
(z:pow(\$b,3)=>z:times(2))=>z:sub(
\$a=>z:multiply(\$b)=>z:multiply(\$c)=>z:times(9)
\$a=>z:multiply(\$a)=>z:multiply(\$d)=>z:times(27)
)
return (
if (z:eq(\$p,0) and z:eq(\$q,0)) then (
z:minus(\$b)=>z:divide(\$a=>z:times(3)) (: -b/3a :)
) else (
let \$ζ0 := z:as-complex(1)
let \$ζ1 := z:complex(-1 div 2, math:sqrt(3) div 2)
let \$ζ2 := z:multiply(\$ζ1, \$ζ1)
(: r = cbrt((p ± sqrt(q² - 4p³))/2) :)
let \$r := (
z:cbrt(
z:sqrt(
z:multiply(\$q,\$q)=>
z:sub(
z:multiply(z:multiply(z:times(\$p, 4), \$p), \$p)
)
)
)=>z:times(0.5)
)
)
let \$r := (
if (z:eq(\$r, 0)) then (
z:cbrt(
\$q=>z:sub(
z:sqrt(
z:multiply(\$q,\$q)=>
z:sub(
z:multiply(z:multiply(z:times(\$p, 4), \$p), \$p)
)
)
)=>z:times(0.5)
)
) else \$r
)
let \$rζ0 := \$r=>z:multiply(\$ζ0)
let \$rζ1 := \$r=>z:multiply(\$ζ1)
let \$rζ2 := \$r=>z:multiply(\$ζ2)
(: x = (-1/(3a))(b + (ζ^k)r + p/((ζ^k)r) for k=0 to 2 :)
for \$rζk in (\$rζ0, \$rζ1, \$rζ2) return (
)
)
)
)
};

(:~
: real-roots()
: Return all the real-valued (within ε) roots of the polynomial.
:
: @param \$polynomial: the polynomial
: @return real roots, if any
: @error UTIL-CANNOT if the polynomial is order 0 or more than 3
:)
declare function this:real-roots(
\$polynomial as map(xs:string,item()*)
) as xs:double*
{
let \$op := poly:op(\$polynomial)
let \$raw-roots := (
if (util:kind(\$polynomial)="polynomial") then (
let \$cfs as xs:double* := poly:coefficients(\$polynomial)
return switch (poly:order(\$polynomial))
case 1 return this:linear-root(\$cfs[1], \$cfs[2])
case 2 return this:quadratic-real-roots(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-real-roots(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
) else if (util:kind(\$polynomial)="complex-polynomial") then (
let \$cfs as map(xs:string,item()*)* := zpoly:coefficients(\$polynomial)
return switch (zpoly:order(\$polynomial))
case 1 return this:linear-real-root-complex(\$cfs[1], \$cfs[2])
case 2 return this:quadratic-real-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-real-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
) else (
errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
)
)
return (
if (empty(\$op)) then (
\$raw-roots
) else if (\$op="sin") then (
(:
: We know values of sin(x) for which poly(sin(x)) = 0
: Values of x are therefor asin of those
: e.g. sin²(x) - 1 = 0 => sin(x) = {1, -1} => x = {asin(1), asin(-1)}
:)
for \$v in \$raw-roots return math:asin(\$v)
) else if (\$op="cos") then (
(:
: We know values of cos(x) for which poly(cos(x)) = 0
: Values of x are therefor acos of those
: e.g. cos²(x) - 1 = 0 => cos(x) = {1, -1} => x = {acos(1), acos(-1)}
:)
for \$v in \$raw-roots return math:acos(\$v)
) else if (\$op="sinh") then (
for \$v in \$raw-roots return util:asinh(\$v)
) else if (\$op="cosh") then (
for \$v in \$raw-roots return util:acosh(\$v)
) else (
errors:error("UTIL-CANNOT", ("real-roots", \$polynomial))
)
)
};

(:~
: roots()
: Return all the roots of the polynomial.
:
: @param \$polynomial: the polynomial
: @return roots
: @error UTIL-CANNOT if the polynomial is order 0 or more than 3
:)
declare function this:roots(
\$polynomial as map(xs:string,item()*)
) as map(xs:string,item()*)*
{
let \$op := poly:op(\$polynomial)
let \$raw-roots := (
if (util:kind(\$polynomial)="polynomial") then (
let \$cfs as xs:double* := poly:coefficients(\$polynomial)
return switch (poly:order(\$polynomial))
case 1 return this:linear-root(\$cfs[1], \$cfs[2])!z:as-complex(.)
case 2 return this:quadratic-roots(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-roots(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("roots", \$polynomial))
) else if (util:kind(\$polynomial)="complex-polynomial") then (
let \$cfs as map(xs:string,item()*)* := zpoly:coefficients(\$polynomial)
return switch (zpoly:order(\$polynomial))
case 1 return this:linear-root-complex(\$cfs[1], \$cfs[2])
case 2 return this:quadratic-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3])
case 3 return this:cubic-roots-complex(\$cfs[1], \$cfs[2], \$cfs[3], \$cfs[4])
default return errors:error("UTIL-CANNOT", ("roots", \$polynomial))
) else (
errors:error("UTIL-CANNOT", ("roots", \$polynomial))
)
)
return (
if (empty(\$op)) then (
\$raw-roots
) else if (\$op="sin") then (
(:
: We know values of sin(x) for which poly(sin(x)) = 0
: Values of x are therefor asin of those
: e.g. sin²(x) - 1 = 0 => sin(x) = {1, -1} => x = {asin(1), asin(-1)}
:)
for \$v in \$raw-roots return z:asin(\$v)
) else if (\$op="cos") then (
(:
: We know values of cos(x) for which poly(cos(x)) = 0
: Values of x are therefor acos of those
: e.g. cos²(x) - 1 = 0 => cos(x) = {1, -1} => x = {acos(1), acos(-1)}
:)
for \$v in \$raw-roots return z:acos(\$v)
) else if (\$op="sinh") then (
for \$v in \$raw-roots return z:asinh(\$v)
) else if (\$op="cosh") then (
for \$v in \$raw-roots return z:acosh(\$v)
) else (
errors:error("UTIL-CANNOT", ("roots", \$polynomial))
)
)
};```