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

http://mathling.com/core/roots


Roots of linear, quadratic, and cubic equations.

Copyright© Mary Holstege 2021-2023
CC-BY (https://creativecommons.org/licenses/by/4.0/)

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-root
declare 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-complex
declare 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-complex
declare 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-roots
declare function quadratic-real-roots($a as xs:double,
  $b as xs:double,
  $c as xs:double) as xs:double*


quadratic-real-roots()
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-complex
declare 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*


quadratic-real-roots-complex()
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:add($sqrt))=>z:divide($a=>z:times(2)),
      (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-roots
declare function quadratic-roots($a as xs:double,
  $b as xs:double,
  $c as xs:double) as map(xs:string,item()*)*


quadratic-roots()
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:add($sqrt))=>z:times(1 div (2*$a)),
      (z:as-complex(-$b)=>z:sub($sqrt))=>z:times(1 div (2*$a))
    )
  )
}

Function: quadratic-roots-complex
declare 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()*)*


quadratic-roots-complex()
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:add($sqrt))=>z:divide($a=>z:times(2)),
      (z:minus($b)=>z:sub($sqrt))=>z:divide($a=>z:times(2))
    )
  )
}

Function: cubic-real-roots
declare 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 (
    this:quadratic-real-roots($b, $c, $d)
  ) 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-complex
declare 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 (
    this:quadratic-real-roots-complex($b, $c, $d)
  ) 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)
      )=>z:add(
        $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-roots
declare 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 (
    this:quadratic-roots($b, $c, $d)
  ) 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(
              $qz=>z:add(
                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 (
        ($bz=>z:add($rζk)=>z:add( $pz=>z:divide($rζk) ))=>z:times(-1 div (3*$a))
      )
    )
  )
}

Function: cubic-roots-complex
declare 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 (
    this:quadratic-roots-complex($b, $c, $d)
  ) 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)
      )=>z:add(
        $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(
            $q=>z:add(
              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 (
          ($b=>z:add($rζk)=>z:add( $p=>z:divide($rζk) ))=>z:times(-1 div 3)=>z:divide($a)
        )
      )
    )
  )
}

Function: real-roots
declare 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: roots
declare 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.
 :
 : Copyright© Mary Holstege 2021-2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @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
 :)

(:~
 : quadratic-real-roots()
 : 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
 :)
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 :)
  )
};

(:~
 : quadratic-real-roots-complex()
 : Same as quadratic-real-roots() except over complex coefficients. 
 :)
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:add($sqrt))=>z:divide($a=>z:times(2)),
      (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()
  )
};

(:~
 : quadratic-roots()
 : 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
 :)
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:add($sqrt))=>z:times(1 div (2*$a)),
      (z:as-complex(-$b)=>z:sub($sqrt))=>z:times(1 div (2*$a))
    )
  )
};

(:~
 : quadratic-roots-complex()
 : Same as quadratic-roots() but over complex coefficients.
 :)
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:add($sqrt))=>z:divide($a=>z:times(2)),
      (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 (
    this:quadratic-real-roots($b, $c, $d)
  ) 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 (
    this:quadratic-real-roots-complex($b, $c, $d)
  ) 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)
      )=>z:add(
        $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 (
    this:quadratic-roots($b, $c, $d)
  ) 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(
              $qz=>z:add(
                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 (
        ($bz=>z:add($rζk)=>z:add( $pz=>z:divide($rζk) ))=>z:times(-1 div (3*$a))
      )
    )
  )
};

(:~
 : 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 (
    this:quadratic-roots-complex($b, $c, $d)
  ) 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)
      )=>z:add(
        $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(
            $q=>z:add(
              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 (
          ($b=>z:add($rζk)=>z:add( $p=>z:divide($rζk) ))=>z:times(-1 div 3)=>z:divide($a)
        )
      )
    )
  )
};

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

| Generated by XQDoc https://xqdoc.org/ with variants of xquerydoc https://github.com/xquery/xquerydoc scripts. |