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

http://mathling.com/core/functions


Functions that we need to perform analytical operations on, e.g.
arithmetic ops, derivative(), etc.
A unifying API for things like polynomials, complex polynomials, polynomial
quotients, etc.

Provides for polymorphism over a certain class of functions in the same way
that geo/euclidean provides polymorphism over a class of geometric objects.

Kinds of functions supported:
constants (core/complex)
real polynomials over x, sin(x), and cos(x) (types/polynomial)
complex polynomials over x, sin(x), and cos(x) (types/cpolynomial)

Note: limitations on some quotient operators right now

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

April 2023
Status: Bleeding edge

Imports

http://mathling.com/type/polynomial/quotient
import module namespace quotient="http://mathling.com/type/polynomial/quotient"
       at "../types/quotient.xqy"
http://mathling.com/core/utilities
import module namespace util="http://mathling.com/core/utilities"
       at "../core/utilities.xqy"
http://mathling.com/type/polynomial
import module namespace poly="http://mathling.com/type/polynomial"
       at "../types/polynomial.xqy"
http://mathling.com/core/complex
import module namespace z="http://mathling.com/core/complex"
       at "../core/complex.xqy"
http://mathling.com/core/complex/vector
import module namespace zv="http://mathling.com/core/complex/vector"
       at "../core/vcomplex.xqy"
http://mathling.com/core/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: value
declare function value($f as map(xs:string,item()*), $x as xs:double) as map(xs:string,item()*)


value()
The value of f(x) for a real-valued x. Returns a complex value:
use z:as-real() if you need to.

Params
  • f as map(xs:string,item()*): an object representing a certain kind of function
  • x as xs:double: double value
Returns
  • map(xs:string,item()*): : complex value f(x)
declare function this:value(
  $f as map(xs:string,item()*),
  $x as xs:double
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then $f
    else if ($kind="polynomial") then poly:value($f, $x)=>z:as-complex()
    else if ($kind="complex-polynomial") then zpoly:value($f, $x)
    else if ($kind="polynomial-quotient") then (
      this:value(quotient:u($f), $x)=>z:divide( this:value(quotient:v($f), $x) )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: zvalue
declare function zvalue($f as map(xs:string,item()*), $z as map(xs:string,item()*)) as map(xs:string,item()*)


zvalue()
The value of f(z) for a complex value z

Params
  • f as map(xs:string,item()*): an object representing a certain kind of function
  • z as map(xs:string,item()*): complex value
Returns
  • map(xs:string,item()*): : a complex value f(z)
declare function this:zvalue(
  $f as map(xs:string,item()*),
  $z as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then $f
    else if ($kind="polynomial") then poly:zvalue($f, $z)
    else if ($kind="complex-polynomial") then zpoly:zvalue($f, $z)
    else if ($kind="polynomial-quotient") then (
      this:zvalue(quotient:u($f), $z)=>z:divide( this:zvalue(quotient:v($f), $z) )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: zvaluev
declare function zvaluev($f as map(xs:string,item()*), $z as xs:double*) as xs:double*


zvaluev()
The value of f(z) for a complex value z vector

Params
  • f as map(xs:string,item()*): an object representing a certain kind of function
  • z as xs:double*: complex value vector
Returns
  • xs:double*: : a complex value f(z) vector
declare function this:zvaluev(
  $f as map(xs:string,item()*),
  $z as xs:double*
) as xs:double*
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then $f=>z:as-vector()
    else if ($kind="polynomial") then poly:zvaluev($f, $z)
    else if ($kind="complex-polynomial") then zpoly:zvaluev($f, $z)
    else if ($kind="polynomial-quotient") then (
      this:zvaluev(quotient:u($f), $z)=>zv:divide( this:zvaluev(quotient:v($f), $z) )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: derivative
declare function derivative($f as map(xs:string,item()*)) as map(xs:string,item()*)


derivative()
Compute the derivative of the function.

Params
  • f as map(xs:string,item()*): an object representing a certain kind of function
Returns
  • map(xs:string,item()*): : an object representing the function that is the derivative of f
declare function this:derivative(
  $f as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then (
      z:as-complex(0)
    ) else if ($kind="polynomial") then (
      poly:derivative($f)
    ) else if ($kind="complex-polynomial") then (
      zpoly:derivative($f)
    ) else if ($kind="polynomial-quotient") then (
      (: d(u/v) = (v du - u dv) / v² :)
      let $u := quotient:u($f)
      let $v := quotient:v($f)
      let $du := this:derivative($u)
      let $dv := this:derivative($v)
      return (
        quotient:quotient(
          (this:multiply($v, $du)=>this:sub( this:multiply($u, $dv) )),
          this:multiply($v, $v)
        )
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: antiderivative
declare function antiderivative($f as map(xs:string,item()*)) as map(xs:string,item()*)


anti-derivative()
Compute the anti-derivative of the function.

Params
  • f as map(xs:string,item()*): an object representing a certain kind of function
Returns
  • map(xs:string,item()*): : an object representing the function that is the anti-derivative of f
declare function this:antiderivative(
  $f as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then (
      if (z:eq($f, 0)) then $f else zpoly:polynomial( ($f, z:as-complex(0)) )
    ) else if ($kind="polynomial") then (
      poly:antiderivative($f)
    ) else if ($kind="complex-polynomial") then (
      zpoly:antiderivative($f)
    ) else if ($kind="polynomial-quotient") then (
      (: Not easy:
       : ∫u'/v = u/v + ∫(uv'/v²) so
       : ∫u/v = (∫u)/v + ∫((∫u)v'/v²)
       : Which in theory we could keep going on until the right hand bit gets
       : down to a v' that is 0, but skip this for now
       :)
      errors:error("UTIL-NOTIMPLEMENTED", ("anti-derivative", $f))
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: definite-integral
declare function definite-integral($f as map(xs:string,item()*), $a as xs:double, $b as xs:double) as xs:double


definite-integral()
Definite integral of a function (exact)
Note: if the result includes an imaginary part, this will raise
an error.

Params
  • f as map(xs:string,item()*): the function
  • 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(
  $f as map(xs:string,item()*),
  $a as xs:double,
  $b as xs:double
) as xs:double
{
  let $anti := this:antiderivative($f)
  return ($anti=>this:value($b)=>z:sub( $anti=>this:value($a) ) )=>z:as-real()
}

Function: describe
declare function describe($f as map(xs:string,item()*)) as xs:string


describe()
Describe the function

Params
  • f as map(xs:string,item()*): the function
Returns
  • xs:string
declare function this:describe(
  $f as map(xs:string,item()*)
) as xs:string
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then z:quote($f)
    else if ($kind="polynomial") then poly:describe($f)
    else if ($kind="complex-polynomial") then zpoly:describe($f)
    else if ($kind="polynomial-quotient") then (
      "("||this:describe(quotient:u($f))||")/("||this:describe(quotient:v($f))||")"
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: function
declare function function($f as map(xs:string,item()*)) as function(xs:double) as xs:double


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

Params
  • f as map(xs:string,item()*): the function
Returns
  • function(xs:double)asxs:double
declare function this:function(
  $f as map(xs:string,item()*)
) as function(xs:double) as xs:double
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then
      function($v as xs:double) as xs:double {
        z:as-real($f)
      }
    else if ($kind="polynomial") then poly:function($f)
    else if ($kind="complex-polynomial") then zpoly:function($f)
    else if ($kind="polynomial-quotient") then (
      let $uf := this:function(quotient:u($f))
      let $vf := this:function(quotient:v($f))
      return (
        function($v as xs:double) as xs:double {
          $uf($v) div $vf($v)
        }
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: zfunction
declare function zfunction($f 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 function, e.g. for curve plotting

Params
  • f as map(xs:string,item()*): the function
Returns
  • function(map(xs:string,item()*))asmap(xs:string,item()*)
declare function this:zfunction(
  $f as map(xs:string,item()*)
) as function(map(xs:string,item()*)) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then
      function($v as map(xs:string,item()*)) as map(xs:string,item()*) {
        $f
      }
    else if ($kind="polynomial") then poly:function($f)
    else if ($kind="complex-polynomial") then zpoly:function($f)
    else if ($kind="polynomial-quotient") then (
      let $uf := this:zfunction(quotient:u($f))
      let $vf := this:zfunction(quotient:v($f))
      return (
        function($v as map(xs:string,item()*)) as map(xs:string,item()*) {
          $uf($v)=>z:divide( $vf($v) )
        }
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: add
declare function add($f1 as map(xs:string,item()*), $f2 as map(xs:string,item()*)) as map(xs:string,item()*)


add()
Add two functions

Params
  • f1 as map(xs:string,item()*): one function
  • f2 as map(xs:string,item()*): another function
Returns
  • map(xs:string,item()*)
declare function this:add(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex") then $f1=>z:add($f2) else this:add($f2, $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        zpoly:as-polynomial($f1)=>zpoly:add( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>poly:add($f2)
      ) else if ($kind2="complex-polynomial") then (
        zpoly:as-polynomial($f1)=>zpoly:add( $f2 )
      ) else if ($kind2="polynomial-quotient") then (
        this:add($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        $f1=>zpoly:add( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>zpoly:add( zpoly:as-polynomial($f2) )
      ) else if ($kind2="complex-polynomial") then (
        $f1=>zpoly:add($f2)
      ) else if ($kind2="polynomial-quotient") then (
        this:add($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: b/c + a = (b + ac) / c :)
        quotient:quotient(
          this:add( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient(
          this:add( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1)
        let $v1 := quotient:v($f1)
        let $u2 := quotient:u($f2)
        let $v2 := quotient:v($f2)
        (: a/b + c/d = (ad + cb) / (bd) :)
        return (
          quotient:quotient(
            this:add( this:multiply($u1, $v2), this:multiply($u2, $v1) ),
            this:multiply( $v1, $v2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
}

Function: sub
declare function sub($f1 as map(xs:string,item()*), $f2 as map(xs:string,item()*)) as map(xs:string,item()*)


sub()
Subtract two functions

Params
  • f1 as map(xs:string,item()*): one function
  • f2 as map(xs:string,item()*): another function
Returns
  • map(xs:string,item()*)
declare function this:sub(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex") then $f1=>z:sub($f2) else this:add(this:minus($f2), $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        zpoly:as-polynomial($f1)=>zpoly:sub( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>poly:sub($f2)
      ) else if ($kind2="complex-polynomial") then (
        zpoly:as-polynomial($f1)=>zpoly:sub( $f2 )
      ) else if ($kind2="polynomial-quotient") then (
        this:add(this:minus($f2), $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        $f1=>zpoly:sub( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>zpoly:sub( zpoly:as-polynomial($f2) )
      ) else if ($kind2="complex-polynomial") then (
        $f1=>zpoly:sub($f2)
      ) else if ($kind2="polynomial-quotient") then (
        this:add(this:minus($f2), $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: b/c - a = (b - ac) / c :)
        quotient:quotient(
          this:sub( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="complex-polynomial") then (
        (: b/c - a = (b + ac) / c :)
        quotient:quotient(
          this:sub( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1)
        let $v1 := quotient:v($f1)
        let $u2 := quotient:u($f2)
        let $v2 := quotient:v($f2)
        (: a/b - c/d = (ad - cb) / (bd) :)
        return (
          quotient:quotient(
            this:sub( this:multiply($u1, $v2), this:multiply($u2, $v1) ),
            this:multiply( $v1, $v2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
}

Function: multiply
declare function multiply($f1 as map(xs:string,item()*), $f2 as map(xs:string,item()*)) as map(xs:string,item()*)


multiply()
Multiply two functions

Params
  • f1 as map(xs:string,item()*): one function
  • f2 as map(xs:string,item()*): another function
Returns
  • map(xs:string,item()*)
declare function this:multiply(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex") then $f1=>z:multiply($f2) else this:multiply($f2, $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(z:as-real($f2))
        else zpoly:as-polynomial($f1)=>zpoly:multiply( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>poly:multiply($f2)
      ) else if ($kind2="complex-polynomial") then (
        zpoly:as-polynomial($f1)=>zpoly:multiply( $f2 )
      ) else if ($kind2="polynomial-quotient") then (
        this:multiply($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(z:as-real($f2))
        else $f1=>zpoly:multiply( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>zpoly:multiply( zpoly:as-polynomial($f2) )
      ) else if ($kind2="complex-polynomial") then (
        $f1=>zpoly:multiply($f2)
      ) else if ($kind2="polynomial-quotient") then (
        this:multiply($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: b/c * a = ba / c :)
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1)
        let $v1 := quotient:v($f1)
        let $u2 := quotient:u($f2)
        let $v2 := quotient:v($f2)
        (: a/b * c/d = (ac) / (bd) :)
        return (
          quotient:quotient(
            this:multiply( $u1, $u1 ),
            this:multiply( $v1, $v2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
}

Function: divide
declare function divide($f1 as map(xs:string,item()*), $f2 as map(xs:string,item()*)) as map(xs:string,item()*)


divide()
Divide two functions

Params
  • f1 as map(xs:string,item()*): one function
  • f2 as map(xs:string,item()*): another function
Returns
  • map(xs:string,item()*)
declare function this:divide(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex")
      then $f1=>z:divide($f2)
      else this:multiply(z:reciprocal($f2), $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(1 div z:as-real($f2))
        else zpoly:as-polynomial($f1)=>zpoly:ztimes( z:reciprocal($f2) )
      ) else if ($kind2="polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="polynomial-quotient") then (
        (: a / (b/c) =  ac / b :)
        quotient:quotient(
          this:multiply($f1, quotient:v($f2)),
          quotient:u($f2)
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(1 div z:as-real($f2))
        else $f1=>zpoly:ztimes( z:reciprocal($f2) )
      ) else if ($kind2="polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="polynomial-quotient") then (
        (: a / (b/c) =  ac / b :)
        quotient:quotient(
          this:multiply($f1, quotient:v($f2)),
          quotient:u($f2)
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), z:reciprocal($f2)),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: (b/c) / a = b / ca :)
        quotient:quotient(
          quotient:u($f1),
          this:multiply(quotient:v($f1), $f2)
        )
      ) else if ($kind2="complex-polynomial") then (
        (: (b/c) / a = b / ca :)
        quotient:quotient(
          quotient:u($f1),
          this:multiply(quotient:v($f1), $f2)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1) 
        let $v1 := quotient:v($f1) 
        let $u2 := quotient:u($f2) 
        let $v2 := quotient:v($f2) 
        (: a/b / c/d = a/b * d/c = ad/bc :)
        return (
          quotient:quotient(
            this:multiply( $u1, $v2),
            this:multiply( $v1, $u2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
}

Function: times
declare function times($f as map(xs:string,item()*), $k as xs:double) as map(xs:string,item()*)


times()
Multiply function by constant

Params
  • f as map(xs:string,item()*): the function
  • k as xs:double: constant
Returns
  • map(xs:string,item()*)
declare function this:times(
  $f as map(xs:string,item()*),
  $k as xs:double
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then z:times($f, $k)
    else if ($kind="polynomial") then poly:times($f, $k)
    else if ($kind="complex-polynomial") then zpoly:times($f, $k)
    else if ($kind="polynomial-quotient") then (
      quotient:quotient(
        this:times(quotient:u($f), $k),
        quotient:v($f)
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: minus
declare function minus($f as map(xs:string,item()*)) as map(xs:string,item()*)


minus()
Negate the function

Params
  • f as map(xs:string,item()*): the function
Returns
  • map(xs:string,item()*)
declare function this:minus(
  $f as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then z:minus($f)
    else if ($kind="polynomial") then poly:times($f, -1)
    else if ($kind="complex-polynomial") then zpoly:times($f, -1)
    else if ($kind="polynomial-quotient") then (
      quotient:quotient(
        this:minus(quotient:u($f)),
        quotient:v($f)
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Function: order
declare function order($f as map(xs:string,item()*)) as xs:integer


order()
Order of the function.

Params
  • f as map(xs:string,item()*): the function
Returns
  • xs:integer
declare function this:order(
  $f as map(xs:string,item()*)
) as xs:integer
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then 0
    else if ($kind="polynomial") then poly:order($f)
    else if ($kind="complex-polynomial") then zpoly:order($f)
    else if ($kind="polynomial-quotient") then (
      abs(this:order(quotient:u($f)) - this:order(quotient:v($f)))
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
}

Original Source Code

xquery version "3.1";
(:~
 : Functions that we need to perform analytical operations on, e.g.
 : arithmetic ops, derivative(), etc.
 : A unifying API for things like polynomials, complex polynomials, polynomial
 : quotients, etc.
 :
 : Provides for polymorphism over a certain class of functions in the same way
 : that geo/euclidean provides polymorphism over a class of geometric objects.
 :
 : Kinds of functions supported:
 : constants (core/complex)
 : real polynomials over x, sin(x), and cos(x) (types/polynomial)
 : complex polynomials over x, sin(x), and cos(x) (types/cpolynomial)
 :
 : Note: limitations on some quotient operators right now
 :
 : Copyright© Mary Holstege 2022-2023
 : CC-BY (https://creativecommons.org/licenses/by/4.0/)
 : @since April 2023
 : @custom:Status Bleeding edge
 :)
module namespace this="http://mathling.com/core/functions"; 

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 zv="http://mathling.com/core/complex/vector"
       at "../core/vcomplex.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";
import module namespace quotient="http://mathling.com/type/polynomial/quotient"
       at "../types/quotient.xqy";

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

(:~
 : value()
 : The value of f(x) for a real-valued x. Returns a complex value:
 : use z:as-real() if you need to.
 :
 : @param $f: an object representing a certain kind of function
 : @param $x: double value
 : @return: complex value f(x)
 :)
declare function this:value(
  $f as map(xs:string,item()*),
  $x as xs:double
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then $f
    else if ($kind="polynomial") then poly:value($f, $x)=>z:as-complex()
    else if ($kind="complex-polynomial") then zpoly:value($f, $x)
    else if ($kind="polynomial-quotient") then (
      this:value(quotient:u($f), $x)=>z:divide( this:value(quotient:v($f), $x) )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : zvalue()
 : The value of f(z) for a complex value z
 :
 : @param $f: an object representing a certain kind of function
 : @param $z: complex value
 : @return: a complex value f(z)
 :)
declare function this:zvalue(
  $f as map(xs:string,item()*),
  $z as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then $f
    else if ($kind="polynomial") then poly:zvalue($f, $z)
    else if ($kind="complex-polynomial") then zpoly:zvalue($f, $z)
    else if ($kind="polynomial-quotient") then (
      this:zvalue(quotient:u($f), $z)=>z:divide( this:zvalue(quotient:v($f), $z) )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : zvaluev()
 : The value of f(z) for a complex value z vector
 :
 : @param $f: an object representing a certain kind of function
 : @param $z: complex value vector
 : @return: a complex value f(z) vector
 :)
declare function this:zvaluev(
  $f as map(xs:string,item()*),
  $z as xs:double*
) as xs:double*
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then $f=>z:as-vector()
    else if ($kind="polynomial") then poly:zvaluev($f, $z)
    else if ($kind="complex-polynomial") then zpoly:zvaluev($f, $z)
    else if ($kind="polynomial-quotient") then (
      this:zvaluev(quotient:u($f), $z)=>zv:divide( this:zvaluev(quotient:v($f), $z) )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : derivative()
 : Compute the derivative of the function.
 :
 : @param $f: an object representing a certain kind of function
 : @return: an object representing the function that is the derivative of f
 :)
declare function this:derivative(
  $f as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then (
      z:as-complex(0)
    ) else if ($kind="polynomial") then (
      poly:derivative($f)
    ) else if ($kind="complex-polynomial") then (
      zpoly:derivative($f)
    ) else if ($kind="polynomial-quotient") then (
      (: d(u/v) = (v du - u dv) / v² :)
      let $u := quotient:u($f)
      let $v := quotient:v($f)
      let $du := this:derivative($u)
      let $dv := this:derivative($v)
      return (
        quotient:quotient(
          (this:multiply($v, $du)=>this:sub( this:multiply($u, $dv) )),
          this:multiply($v, $v)
        )
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : anti-derivative()
 : Compute the anti-derivative of the function.
 :
 : @param $f: an object representing a certain kind of function
 : @return: an object representing the function that is the anti-derivative of f
 :)
declare function this:antiderivative(
  $f as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then (
      if (z:eq($f, 0)) then $f else zpoly:polynomial( ($f, z:as-complex(0)) )
    ) else if ($kind="polynomial") then (
      poly:antiderivative($f)
    ) else if ($kind="complex-polynomial") then (
      zpoly:antiderivative($f)
    ) else if ($kind="polynomial-quotient") then (
      (: Not easy:
       : ∫u'/v = u/v + ∫(uv'/v²) so
       : ∫u/v = (∫u)/v + ∫((∫u)v'/v²)
       : Which in theory we could keep going on until the right hand bit gets
       : down to a v' that is 0, but skip this for now
       :)
      errors:error("UTIL-NOTIMPLEMENTED", ("anti-derivative", $f))
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : definite-integral()
 : Definite integral of a function (exact)
 : Note: if the result includes an imaginary part, this will raise
 : an error.
 : 
 : @param $f: the function
 : @param $a: start of range of integration
 : @param $b: end of range of integration
 :)
declare function this:definite-integral(
  $f as map(xs:string,item()*),
  $a as xs:double,
  $b as xs:double
) as xs:double
{
  let $anti := this:antiderivative($f)
  return ($anti=>this:value($b)=>z:sub( $anti=>this:value($a) ) )=>z:as-real()
};

(:~
 : describe()
 : Describe the function
 :
 : @param $f: the function
 :)
declare function this:describe(
  $f as map(xs:string,item()*)
) as xs:string
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then z:quote($f)
    else if ($kind="polynomial") then poly:describe($f)
    else if ($kind="complex-polynomial") then zpoly:describe($f)
    else if ($kind="polynomial-quotient") then (
      "("||this:describe(quotient:u($f))||")/("||this:describe(quotient:v($f))||")"
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : function()
 : Return a function that takes values of the function, e.g. for curve plotting
 :
 : @param $f: the function
 :)
declare function this:function(
  $f as map(xs:string,item()*)
) as function(xs:double) as xs:double
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then
      function($v as xs:double) as xs:double {
        z:as-real($f)
      }
    else if ($kind="polynomial") then poly:function($f)
    else if ($kind="complex-polynomial") then zpoly:function($f)
    else if ($kind="polynomial-quotient") then (
      let $uf := this:function(quotient:u($f))
      let $vf := this:function(quotient:v($f))
      return (
        function($v as xs:double) as xs:double {
          $uf($v) div $vf($v)
        }
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : zfunction()
 : Return a function over complex values that takes values of the function, e.g. for curve plotting
 :
 : @param $f: the function
 :)
declare function this:zfunction(
  $f as map(xs:string,item()*)
) as function(map(xs:string,item()*)) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then
      function($v as map(xs:string,item()*)) as map(xs:string,item()*) {
        $f
      }
    else if ($kind="polynomial") then poly:function($f)
    else if ($kind="complex-polynomial") then zpoly:function($f)
    else if ($kind="polynomial-quotient") then (
      let $uf := this:zfunction(quotient:u($f))
      let $vf := this:zfunction(quotient:v($f))
      return (
        function($v as map(xs:string,item()*)) as map(xs:string,item()*) {
          $uf($v)=>z:divide( $vf($v) )
        }
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:======================================================================
 : Operations 
 :======================================================================:)

(:~
 : add()
 : Add two functions
 :
 : @param $f1: one function
 : @param $f2: another function
 :)
declare function this:add(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex") then $f1=>z:add($f2) else this:add($f2, $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        zpoly:as-polynomial($f1)=>zpoly:add( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>poly:add($f2)
      ) else if ($kind2="complex-polynomial") then (
        zpoly:as-polynomial($f1)=>zpoly:add( $f2 )
      ) else if ($kind2="polynomial-quotient") then (
        this:add($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        $f1=>zpoly:add( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>zpoly:add( zpoly:as-polynomial($f2) )
      ) else if ($kind2="complex-polynomial") then (
        $f1=>zpoly:add($f2)
      ) else if ($kind2="polynomial-quotient") then (
        this:add($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: b/c + a = (b + ac) / c :)
        quotient:quotient(
          this:add( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient(
          this:add( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1)
        let $v1 := quotient:v($f1)
        let $u2 := quotient:u($f2)
        let $v2 := quotient:v($f2)
        (: a/b + c/d = (ad + cb) / (bd) :)
        return (
          quotient:quotient(
            this:add( this:multiply($u1, $v2), this:multiply($u2, $v1) ),
            this:multiply( $v1, $v2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
};

(:~
 : sub()
 : Subtract two functions
 :
 : @param $f1: one function
 : @param $f2: another function
 :)
declare function this:sub(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex") then $f1=>z:sub($f2) else this:add(this:minus($f2), $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        zpoly:as-polynomial($f1)=>zpoly:sub( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>poly:sub($f2)
      ) else if ($kind2="complex-polynomial") then (
        zpoly:as-polynomial($f1)=>zpoly:sub( $f2 )
      ) else if ($kind2="polynomial-quotient") then (
        this:add(this:minus($f2), $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        $f1=>zpoly:sub( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>zpoly:sub( zpoly:as-polynomial($f2) )
      ) else if ($kind2="complex-polynomial") then (
        $f1=>zpoly:sub($f2)
      ) else if ($kind2="polynomial-quotient") then (
        this:add(this:minus($f2), $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: b/c - a = (b - ac) / c :)
        quotient:quotient(
          this:sub( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="complex-polynomial") then (
        (: b/c - a = (b + ac) / c :)
        quotient:quotient(
          this:sub( quotient:u($f1), this:multiply($f2, quotient:v($f1)) ),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1)
        let $v1 := quotient:v($f1)
        let $u2 := quotient:u($f2)
        let $v2 := quotient:v($f2)
        (: a/b - c/d = (ad - cb) / (bd) :)
        return (
          quotient:quotient(
            this:sub( this:multiply($u1, $v2), this:multiply($u2, $v1) ),
            this:multiply( $v1, $v2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
};

(:~
 : multiply()
 : Multiply two functions
 :
 : @param $f1: one function
 : @param $f2: another function
 :)
declare function this:multiply(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex") then $f1=>z:multiply($f2) else this:multiply($f2, $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(z:as-real($f2))
        else zpoly:as-polynomial($f1)=>zpoly:multiply( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>poly:multiply($f2)
      ) else if ($kind2="complex-polynomial") then (
        zpoly:as-polynomial($f1)=>zpoly:multiply( $f2 )
      ) else if ($kind2="polynomial-quotient") then (
        this:multiply($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(z:as-real($f2))
        else $f1=>zpoly:multiply( zpoly:polynomial(($f2)) )
      ) else if ($kind2="polynomial") then (
        $f1=>zpoly:multiply( zpoly:as-polynomial($f2) )
      ) else if ($kind2="complex-polynomial") then (
        $f1=>zpoly:multiply($f2)
      ) else if ($kind2="polynomial-quotient") then (
        this:multiply($f2, $f1)
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: b/c * a = ba / c :)
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), $f2),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1)
        let $v1 := quotient:v($f1)
        let $u2 := quotient:u($f2)
        let $v2 := quotient:v($f2)
        (: a/b * c/d = (ac) / (bd) :)
        return (
          quotient:quotient(
            this:multiply( $u1, $u1 ),
            this:multiply( $v1, $v2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
};

(:~
 : divide()
 : Divide two functions
 :
 : @param $f1: one function
 : @param $f2: another function
 :)
declare function this:divide(
  $f1 as map(xs:string,item()*),
  $f2 as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind1 := util:kind($f1)
  let $kind2 := util:kind($f2)
  return (
    if ($kind1="complex") then (
      if ($kind2="complex")
      then $f1=>z:divide($f2)
      else this:multiply(z:reciprocal($f2), $f1)
    ) 
    else if ($kind1="polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(1 div z:as-real($f2))
        else zpoly:as-polynomial($f1)=>zpoly:ztimes( z:reciprocal($f2) )
      ) else if ($kind2="polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="polynomial-quotient") then (
        (: a / (b/c) =  ac / b :)
        quotient:quotient(
          this:multiply($f1, quotient:v($f2)),
          quotient:u($f2)
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="complex-polynomial") then (
      if ($kind2="complex") then (
        if (z:is-real($f2))
        then $f1=>this:times(1 div z:as-real($f2))
        else $f1=>zpoly:ztimes( z:reciprocal($f2) )
      ) else if ($kind2="polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="complex-polynomial") then (
        quotient:quotient($f1, $f2)
      ) else if ($kind2="polynomial-quotient") then (
        (: a / (b/c) =  ac / b :)
        quotient:quotient(
          this:multiply($f1, quotient:v($f2)),
          quotient:u($f2)
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    )
    else if ($kind1="polynomial-quotient") then (
      if ($kind2="complex") then (
        quotient:quotient(
          this:multiply(quotient:u($f1), z:reciprocal($f2)),
          quotient:v($f1)
        )
      ) else if ($kind2="polynomial") then (
        (: (b/c) / a = b / ca :)
        quotient:quotient(
          quotient:u($f1),
          this:multiply(quotient:v($f1), $f2)
        )
      ) else if ($kind2="complex-polynomial") then (
        (: (b/c) / a = b / ca :)
        quotient:quotient(
          quotient:u($f1),
          this:multiply(quotient:v($f1), $f2)
        )
      ) else if ($kind2="polynomial-quotient") then (
        let $u1 := quotient:u($f1) 
        let $v1 := quotient:v($f1) 
        let $u2 := quotient:u($f2) 
        let $v2 := quotient:v($f2) 
        (: a/b / c/d = a/b * d/c = ad/bc :)
        return (
          quotient:quotient(
            this:multiply( $u1, $v2),
            this:multiply( $v1, $u2 )
          )
        )
      ) else (
        errors:error("UTIL-BADARGS", ("f2", $f2))
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f1", $f1))
    )
  )
};

(:~
 : times()
 : Multiply function by constant
 :
 : @param $f: the function
 : @param $k: constant
 :)
declare function this:times(
  $f as map(xs:string,item()*),
  $k as xs:double
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then z:times($f, $k)
    else if ($kind="polynomial") then poly:times($f, $k)
    else if ($kind="complex-polynomial") then zpoly:times($f, $k)
    else if ($kind="polynomial-quotient") then (
      quotient:quotient(
        this:times(quotient:u($f), $k),
        quotient:v($f)
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : minus()
 : Negate the function
 :
 : @param $f: the function
 :)
declare function this:minus(
  $f as map(xs:string,item()*)
) as map(xs:string,item()*)
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then z:minus($f)
    else if ($kind="polynomial") then poly:times($f, -1)
    else if ($kind="complex-polynomial") then zpoly:times($f, -1)
    else if ($kind="polynomial-quotient") then (
      quotient:quotient(
        this:minus(quotient:u($f)),
        quotient:v($f)
      )
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};

(:~
 : order()
 : Order of the function. 
 :
 : @param $f: the function
 :)
declare function this:order(
  $f as map(xs:string,item()*)
) as xs:integer
{
  let $kind := util:kind($f) return (
    if ($kind="complex") then 0
    else if ($kind="polynomial") then poly:order($f)
    else if ($kind="complex-polynomial") then zpoly:order($f)
    else if ($kind="polynomial-quotient") then (
      abs(this:order(quotient:u($f)) - this:order(quotient:v($f)))
    ) else (
      errors:error("UTIL-BADARGS", ("f", $f))
    )
  )
};