# 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

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

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

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

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

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

#### ```Function: definite-integraldeclare 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: describedeclare 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 (
)
)
}```

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

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

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

##### 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 (
)
else if (\$kind1="polynomial") then (
if (\$kind2="complex") then (
) else if (\$kind2="polynomial") then (
) else if (\$kind2="complex-polynomial") then (
) else if (\$kind2="polynomial-quotient") then (
) else (
)
)
else if (\$kind1="complex-polynomial") then (
if (\$kind2="complex") then (
) else if (\$kind2="polynomial") then (
) else if (\$kind2="complex-polynomial") then (
) else if (\$kind2="polynomial-quotient") then (
) else (
)
)
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(
quotient:v(\$f1)
)
) else if (\$kind2="complex-polynomial") then (
quotient:quotient(
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 (
)
) else (
)
)
}```

#### ```Function: subdeclare 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 (
) else (
)
)
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 (
) else (
)
)
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 (
)
) else (
)
)
}```

#### ```Function: multiplydeclare 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 (
)
)
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 (
)
)
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 (
)
) else (
)
)
}```

#### ```Function: dividedeclare 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 (
)
)
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 (
)
)
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 (
)
) else (
)
)
}```

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

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

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

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

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

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

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

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

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

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

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

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

(:~
:
: @param \$f1: one function
: @param \$f2: another function
:)
\$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 (
)
else if (\$kind1="polynomial") then (
if (\$kind2="complex") then (
) else if (\$kind2="polynomial") then (
) else if (\$kind2="complex-polynomial") then (
) else if (\$kind2="polynomial-quotient") then (
) else (
)
)
else if (\$kind1="complex-polynomial") then (
if (\$kind2="complex") then (
) else if (\$kind2="polynomial") then (
) else if (\$kind2="complex-polynomial") then (
) else if (\$kind2="polynomial-quotient") then (
) else (
)
)
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(
quotient:v(\$f1)
)
) else if (\$kind2="complex-polynomial") then (
quotient:quotient(
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 (
)
) else (
)
)
};

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

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

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

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

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

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