text:
println("Loading required packages...")
using LinearAlgebra, Random, Printf, Statistics
using BenchmarkTools, JSON

using AccurateArithmetic
using AccurateArithmetic.Summation: accumulate, two_sum, fast_two_sum
using AccurateArithmetic.Summation: sumAcc, dotAcc, compSumAcc, compDotAcc
using AccurateArithmetic.Test

output(x) = @printf "%.2e " x
err(val::T, ref::T) where {T} = min(1, max(eps(T), abs((val-ref)/ref)))

FAST_TESTS = false

# * Accuracy

function run_accuracy(gen, funs, labels, title, filename)
    println("-> $title")

    n = 100     # Vector size
    c1 = 2.     # Condition number (min)
    c2 = 1e45   # -                (max)
    logstep = 2 # Step for condition number increase (log scale)

    data = [Float64[] for _ in 1:(1+length(funs))]

    c = c1
    while c < c2
        i = 1
        (x, d, C) = gen(rand(n:n+10), c)
        output(C)
        push!(data[i], C)

        for fun in funs
            i += 1
            r = fun(x...)
            ε = err(r, d)
            output(ε)
            push!(data[i], ε)
        end

        println()
        c *= logstep
    end

    open(filename*".json", "w") do f
        JSON.print(f, Dict(
            :type   => :accuracy,
            :title  => title,
            :labels => labels,
            :data   => data))
    end
end

function run_accuracy()
    println("Running accuracy tests...")

    function gen_sum(n, c)
        (x, d, c) = generate_sum(n, c)
        ((x,), d, c)
    end
    run_accuracy(gen_sum,
                 (sum,        sum_naive, sum_oro, sum_kbn),
                 ("pairwise", "naive",   "oro",   "kbn"),
                 "Error of summation algorithms",
                 "sum_accuracy")

    function gen_dot(n, c)
        (x, y, d, c) = generate_dot(n, c)
        ((x, y), d, c)
    end
    run_accuracy(gen_dot,
                 (dot,    dot_naive, dot_oro),
                 ("blas", "naive",   "oro"),
                 "Error of dot product algorithms",
                 "dot_accuracy")
end

# * Optimal u_shift

function run_ushift(gen, acc, title, filename)
    println("-> $title")

    if FAST_TESTS
        title *= " [FAST]"
        sizes = [2^(3*i) for i in 3:4]
        ushifts = (0,2)
        BenchmarkTools.DEFAULT_PARAMETERS.evals = 1
        BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
    else
        sizes = [2^(3*i) for i in 2:6]
        ushifts = 0:4
        BenchmarkTools.DEFAULT_PARAMETERS.evals = 2
        BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5.0
    end

    data = [[] for _ in 1:(1+length(sizes))]
    for ushift in ushifts
        i = 1
        print(ushift, " ")
        push!(data[i], ushift)

        for n in sizes
            i += 1
            x = gen(n)

            b = @benchmark accumulate($x, $acc, $(Val(:scalar)), $(Val(ushift)))
            t = minimum(b.times) / n
            output(t)
            push!(data[i], t)
        end
        println()
    end

    open(filename*".json", "w") do f
        JSON.print(f, Dict(
            :type   => :ushift,
            :title  => title,
            :labels => sizes,
            :data   => data))
    end
end

function run_ushift()
    BenchmarkTools.DEFAULT_PARAMETERS.evals = 2
    println("Finding optimal ushift...")

    run_ushift(n->(rand(n),), sumAcc,
               "Performance of naive summation",
               "sum_naive_ushift")

    run_ushift(n->(rand(n),), compSumAcc(two_sum),
               "Performance of ORO summation",
               "sum_oro_ushift")

    run_ushift(n->(rand(n),), compSumAcc(fast_two_sum),
               "Performance of KBN summation",
               "sum_kbn_ushift")

    run_ushift(n->(rand(n), rand(n)), dotAcc,
               "Performance of naive dot product",
               "dot_naive_ushift")

    run_ushift(n->(rand(n), rand(n)), compDotAcc,
               "Performance of compensated dot product",
               "dot_oro_ushift")
end

# * Performance comparisons

function run_performance(n2, gen, funs, labels, title, filename)
    println("-> $title")

    if FAST_TESTS
        title *= " [FAST]"
        logstep = 100.
        BenchmarkTools.DEFAULT_PARAMETERS.evals = 1
        BenchmarkTools.DEFAULT_PARAMETERS.seconds = 0.5
    else
        logstep = 1.1
        BenchmarkTools.DEFAULT_PARAMETERS.evals = 2
        BenchmarkTools.DEFAULT_PARAMETERS.seconds = 5.0
    end

    data = [Float64[] for _ in 1:(1+length(funs))]
    n = 32
    while n < n2
        i = 1
        x = gen(n)
        output(n)
        push!(data[i], n)

        for fun in funs
            i += 1
            b = @benchmark $fun($(x)...)
            t = minimum(b.times) / n
            output(t)
            push!(data[i], t)
        end

        println()
        N = Int(round(n*logstep))
        N = 32*div(N, 32)
        n = max(N, n+32)
    end

    open(filename*".json", "w") do f
        JSON.print(f, Dict(
            :type      => :performance,
            :title     => title,
            :labels    => labels,
            :elem_size => sizeof(gen(1)),
            :data      => data))
    end
end

function run_performance()
    println("Running performance tests...")

    run_performance(1e8, n->(rand(n),),
                    (sum,        sum_naive, sum_oro, sum_kbn),
                    ("pairwise", "naive",   "oro",   "kbn"),
                    "Performance of summation implementations",
                    "sum_performance")

    BLAS.set_num_threads(1)
    run_performance(3e7, n->(rand(n), rand(n)),
                    (dot,    dot_naive, dot_oro),
                    ("blas", "naive",   "oro"),
                    "Performance of dot product implementations",
                    "dot_performance")
end

# * All tests

function run_tests(fast=false)
    global FAST_TESTS = fast
    print("Running tests")
    FAST_TESTS && print(" in FAST mode")
    println("...\n")

    run_accuracy()
    sleep(5)
    run_ushift()
    sleep(5)
    run_performance()
    println("Normal end of the performance tests")
end


text:
#**************************************************************************************
# make.jl
# =============== part of the GeoEfficiency.jl package.
#
# script for building documentaion of the GeoEfficiency.jl package.
#
#**************************************************************************************

using  Documenter,  DocumenterMarkdown #, DocumenterLaTeX
using GeoEfficiency

const PAGES = Any[
    #"Home" => "index.md",
    "Home" => "index.md",
    "Manual" => [
		"manual/GeoEfficiency.md",
        "manual/Error.md",
        "manual/Input_Console.md",
        "manual/Physics_Model.md",
        "manual/Input_Batch.md",
        "manual/Calculations.md",
        "manual/Output_Interface.md",
     ],
    "Index" => "list.md",
]

const formats = Any[
    Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        canonical = "https://DrKrar.github.io/GeoEfficiency.jl/docs/build/dev/",

    ),

]

isdefined(@__MODULE__,:DocumenterMarkdown) &&  push!(formats, Markdown())

if "pdf" in ARGS
    Sys.iswindows() ?   push!(formats, LaTeX(platform = "native")) :
                        push!(formats, LaTeX(platform = "docker"))
end

makedocs(
    format  = formats,
    modules = [GeoEfficiency],
    clean   = "clean" in ARGS,
    doctest = "doctest" in ARGS,
    sitename= "GeoEfficiency.jl",
    authors = "Mohamed E. Krar",
    pages   = PAGES,
    assets  = ["assets/custom.css"],
)

mktempdir() do tmp
    # Hide the PDF from html-deploydocs
    build = joinpath(@__DIR__, "build")
    files = readdir(build)
    idx = findfirst(f -> startswith(f, "GeoEfficiency.jl") && endswith(f, ".pdf"), files)
    pdf = idx === nothing ? nothing : joinpath(build, files[idx])
    if pdf !== nothing
        pdf = mv(pdf, joinpath(tmp, basename(pdf)))
    end
    # Deploy HTML pages
    @info "Deploying HTML pages"
    deploydocs(
        repo = "github.com/DrKrar/GeoEfficiency.jl.git",
        versions = ["stable" => "v^", "v#.#", "dev" => "dev"],
    )
    if isdefined(@__MODULE__,:DocumenterMarkdown)
        # Deploy Markup pages
        @info "Deploying MarkUp pages"
        deploydocs(
            repo = "github.com/DrKrar/GeoEfficiency.jl.git",
            target = "build/Mrk",
            versions = ["stable" => "v^", "v#.#", "dev" => "dev"],
        )
    end #if
    # Put back PDF into docs/build/pdf
    mkpath(joinpath(build, "pdf"))
    if pdf !== nothing
        pdf = mv(pdf, joinpath(build, "pdf", basename(pdf)))
    end
    # Deploy PDF
    @info "Deploying PDF"
    deploydocs(
        repo = "github.com/DrKrar/GeoEfficiency.jl.git",
        target = "build/pdf",
        branch = "gh-pages-pdf",
        forcepush = true,
    )
end

ext: jl


text:
 #=================================================================================
                                     Utilities
 ==================================================================================#

function subplot_gridsize(num)
    poss = [[x,Int(ceil(num/x))] for x in range(1,stop=Int(floor(sqrt(num)))+1)]
    choice = findmin([sum(x) for x in poss])[2]  # argmin
    return sort(poss[choice])
end

 #=================================================================================
                                   Plotting zoo
 ==================================================================================#

function pairplot(X::AbstractArray{T, 2}; figsize=(10,10), alpha=0.5, bins=50,
    axs=nothing, hvlines0=false, kwargs...) where T <: AbstractFloat
    n, d = size(X)
    (d > 20) && throw("will not plot for d > 20")

    (axs == nothing) && ((fig, axs) = PyPlot.subplots(d, d, figsize=figsize))
    for ix = 1:d, iy = 1:d
        if ix != iy
            if alpha isa Number
                axs[iy, ix].scatter(X[:, ix], X[:, iy], alpha=alpha; kwargs...)
            elseif alpha isa AbstractArray
                scatter_alpha(X[:,ix], X[:,iy], alpha, ax=axs[iy,ix]; kwargs...)
            end
            hvlines0 && [x(0, linestyle=":", color="grey") for x in [axs[iy,ix].axhline, axs[iy,ix].axvline]];
        else
            axs[ix, iy].hist(X[:, ix], bins=bins)
        end
    end
    return axs
end

function scatter_arrays(xs...)
    n = length(xs)
    sz = subplot_gridsize(n)
    f, axs = PyPlot.subplots(sz..., figsize=(5 + (sz[2]>1), sz[1]*3))
    if n == 1   # axs is not an array!
        axs.scatter(unpack_arr(xs[1])...)
        return
    else
        for i in eachindex(xs)
            ax = axs[i]; x = xs[i]
            ax.scatter(unpack_arr(x)...)
        end
    end
end

function scatter_alpha(x1::Vector{T}, x2::Vector{T}, alpha::Vector{T2}; cmap_ix::Union{Int, Vector{T3}}=0, cmap::String="tab10",
        rescale_alpha::Bool=true, ax=nothing) where T <: Real where T2 <: AbstractFloat where T3 <: Signed
    n = length(alpha)
    ax = something(ax, gca())
    if isa(cmap_ix, Signed)
        cols = repeat(collect(ColorMap(cmap)(cmap_ix))', n, 1)
    else
        cols = ColorMap(cmap)(cmap_ix)
    end
    rescale_alpha ? (alpha /= maximum(alpha)) : nothing
    cols[:,4] = alpha
    ax.scatter(x1, x2, color=cols)
end

function hinton(matrix; max_weight=nothing, ax=gca())
    """
    Draw Hinton diagram for visualizing a weight matrix.
    https://matplotlib.org/gallery/specialty_plots/hinton_demo.html
    """

    max_weight = something(max_weight, 2^ceil(log(maximum(abs.(matrix[:]))) / log(2)))

    ax.patch.set_facecolor("gray")
    ax.set_aspect("auto")  #("equal", "box")
    ax.xaxis.set_major_locator(plt.NullLocator())
    ax.yaxis.set_major_locator(plt.NullLocator())

    for ix in CartesianIndices(matrix)
        w = matrix[ix]
        color = w > 0 ? "white" : "black"
        sz = √(abs(w) / max_weight)
        rect = plt.Rectangle([ix[1] - sz / 2, ix[2] - sz / 2], sz, sz,
                             facecolor=color, edgecolor=color)
        ax.add_patch(rect)
    end

    ax.autoscale_view()
    ax.invert_yaxis()
end

 #=================================================================================
                                  Axis manipulation
 ==================================================================================#

function ax_lim_one_side(ax, xy; limstart=nothing, limend=nothing, type="constant")
    #=
    Ported from pyalexutil. Manipulates one current axis limit without changing others.
    The interfaces `x_lim_one_side`, `y_lim_one_side` use this functionality and may
    be preferred to this function as they read better.
    :param ax           - axis object to manipulate
    :param xy           - either "x", "y", the dimension of the axis to manipulate
    :param limstart     - the argument for changing the start number (nothing = no change)
    :param limend       - the     "     "     "      "   end    "
    :param type         - what to do with the start/end axis: "constant" specifies
                          overriding current value with limstart/limend, "multiply"/"*"
                          and "add"/"+" also accepted which multiply/add curr. number.
    =#
    lims = (xy == "x") ? ax.get_xlim() : ax.get_ylim()
    lims = collect(lims)  # make mutable (is tuple typed)

    if type == "m" || type == "multiply" || type == "*"
        f = *
    elseif type == "a" || type == "add" || type == "+"
        f = +
    elseif type == "c" || type == "constant"
        f = (x, y) -> y
    else
        throw("Unexpected limtype (expecting 'constant', 'add', 'multiply')")
    end

    if limstart != nothing; lims[1] = f(lims[1], limstart); end
    if limend != nothing; lims[2] = f(lims[2], limend); end

    (xy == "x") ? ax.set_xlim(lims) : ax.set_ylim(lims);

end

# convenience wrappers for clean code
function x_lim_one_side(ax; s=nothing, e=nothing, type="constant")
    #=
    Change either the start or end of the specified x-axis. See `ax_lim_one_side`.
    =#
    ax_lim_one_side(ax, "x", limstart=s, limend=e, type=type)
end

function y_lim_one_side(ax; s=nothing, e=nothing, type="constant")
    #=
    Change either the start or end of the specified y-axis. See `ax_lim_one_side`.
    =#
    ax_lim_one_side(ax, "y", limstart=s, limend=e, type=type)
end

x_lim_start_zero(ax) = x_lim_one_side(ax; s=0.)
y_lim_start_zero(ax) = y_lim_one_side(ax; s=0.)

rmaxislabel_y(ax::PyCall.PyObject) = ax.tick_params(axis="y",which="both",left=false,labelleft=false)
rmaxislabel_x(ax::PyCall.PyObject) = ax.tick_params(axis="x",which="both",bottom=false,labelbottom=false)

@deprecate shiftedColorMap(cmap; start=0, midpoint=0.5, stop=1.0, name="shiftedcmap") AxPlot.utils.shiftedColorMap(cmap; start=0, midpoint=0.5, stop=1.0, name="shiftedcmap")

max_stars_repo_name: ornithos/AxPlot

max_issues_repo_name: ornithos/AxPlot

max_forks_repo_name: ornithos/AxPlot


text:
module Package

export hello, domath

"""
    hello(who::String)

Return "Hello, `who`".
"""
hello(who::String) = "Hello, $who"

"""
    domath(x::Number)

Return `x + 5`.
"""
domath(x::Number) = x + 5

"""
    compute_square(x)

Compute the square of something.
"""
compute_square(x) = x^2

end # module


text:
global const ETCD_VER = "v3.2.12"
global const ETCD_DEPS_PATH = joinpath(dirname(dirname(@__FILE__)), "deps")
global const ETCD_SRC_PATH = if is_apple()
    joinpath(ETCD_DEPS_PATH, "etcd-$(ETCD_VER)-darwin-amd64.zip")
elseif is_linux()
    joinpath(ETCD_DEPS_PATH, "etcd-$(ETCD_VER)-linux-amd64.tar.gz")
end

global const ETCD_DEST_PATH = if is_apple()
    joinpath(ETCD_DEPS_PATH, "etcd-$(ETCD_VER)-darwin-amd64")
elseif is_linux()
    joinpath(ETCD_DEPS_PATH, "etcd-$(ETCD_VER)-linux-amd64")
end

global const ETCD_BIN = joinpath(ETCD_DEST_PATH, "etcd")

global const ETCD_DOWNLOAD_URI = if is_apple()
    "https://github.com/coreos/etcd/releases/download/$(ETCD_VER)/etcd-$(ETCD_VER)-darwin-amd64.zip"
elseif is_linux()
    "https://github.com/coreos/etcd/releases/download/$(ETCD_VER)/etcd-$(ETCD_VER)-linux-amd64.tar.gz"
else
    ""
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct AVRFixed <: AVR
        Emf::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a AVR that returns a fixed voltage to the rotor winding

# Arguments
- `Emf::Float64`: Fixed voltage to the rotor winding, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: Fixed AVR has no states
- `n_states::Int64`: Fixed AVR has no states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AVRFixed <: AVR
    "Fixed voltage to the rotor winding"
    Emf::Float64
    ext::Dict{String, Any}
    "Fixed AVR has no states"
    states::Vector{Symbol}
    "Fixed AVR has no states"
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AVRFixed(Emf, ext=Dict{String, Any}(), )
    AVRFixed(Emf, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function AVRFixed(; Emf, ext=Dict{String, Any}(), )
    AVRFixed(Emf, ext, )
end

# Constructor for demo purposes; non-functional.
function AVRFixed(::Nothing)
    AVRFixed(;
        Emf=0,
        ext=Dict{String, Any}(),
    )
end

"""Get AVRFixed Emf."""
get_Emf(value::AVRFixed) = value.Emf
"""Get AVRFixed ext."""
get_ext(value::AVRFixed) = value.ext
"""Get AVRFixed states."""
get_states(value::AVRFixed) = value.states
"""Get AVRFixed n_states."""
get_n_states(value::AVRFixed) = value.n_states
"""Get AVRFixed internal."""
get_internal(value::AVRFixed) = value.internal


text:
"""
    MeasureSNR(signal, noisy; db=false)

Measure the signal-to-noise ratio between the clean input `signal` and the
contaminated input `noisy`.

# Arguments
* `signal::Array{Real, N}`: N-dimensional clean signal. `N` must be <= 5.
* `noisy::Array{Real, N}`: N-dimensional noisy signal of same size as `signal`.

# Keyword arguments
* `db::Bool=false`: `db=false` if the signal-to-noise ratio is measured by
amplitude, or `db=true` if the signal-to-noise ratio is measure in dB.

# Example
```julia
julia> d, extent = SeisHypEvents(); dnoisy = SeisAddNoise(d, 2);
MeasureSNR(d, dnoisy)
```
"""

function MeasureSNR{Ts<:Real, Tn<:Real, N}(signal::Array{Ts, N},
                                           noisy::Array{Tn, N}; db::Bool=false)

    if db==false
        snr = norm(signal)/norm(signal-noisy)
    elseif db==true
        snr = 20*log10(norm(signal)/norm(signal-noisy))
    end
    return snr

end


text:
using LinearAlgebra
import NFFT3

mutable struct NFFT3Plan{T,D} <: AbstractNFFTPlan{T,D,1}
  parent::NFFT3.NFFT{D}
end

mutable struct NFCT3Plan{T,D} <: AbstractNFCTPlan{T,D,1}
  parent::NFFT3.NFCT{D}
end

mutable struct NFST3Plan{T,D} <: AbstractNFSTPlan{T,D,1}
  parent::NFFT3.NFST{D}
end

################
# constructors
################

function NFFT3Plan(x::Matrix{T}, N::NTuple{D,Int};
              dims::Union{Integer,UnitRange{Int64}}=1:D,
              precompute::PrecomputeFlags=TENSOR, sortNodes=false,
              fftflags=UInt32(NFFT3.FFTW_ESTIMATE),
              kwargs...) where {D,T}

  if dims != 1:D
    error("NFFT3 does not support directional plans!")
  end

  m, σ, reltol = accuracyParams(; kwargs...)

  if precompute == AbstractNFFTs.LINEAR
    prePsi =  NFFT3.PRE_LIN_PSI
  elseif precompute == AbstractNFFTs.FULL
    prePsi = NFFT3.PRE_FULL_PSI
  elseif precompute == AbstractNFFTs.TENSOR
    prePsi = NFFT3.PRE_PSI
  else
    error("Precompute $(precompute) not supported by NFFT3!")
  end
  sortN = sortNodes ? NFFT3.NFCT_SORT_NODES : UInt32(0)

  f1 = UInt32(
    NFFT3.PRE_PHI_HUT |
    prePsi |
    # NFFT3.MALLOC_X |
    # NFFT3.MALLOC_F_HAT |
    # NFFT3.MALLOC_F |
    NFFT3.FFTW_INIT |
    NFFT3.FFT_OUT_OF_PLACE |
    sortN |
    NFFT3.NFCT_OMP_BLOCKWISE_ADJOINT
     )

  f2 = UInt32(fftflags | NFFT3.FFTW_DESTROY_INPUT)

  n = ntuple(d -> (ceil(Int,σ*N[d])÷2)*2, D) # ensure that n is an even integer

  parent = NFFT3.NFFT(Int32.(reverse(N)), size(x,2), Int32.(reverse(n)), Int32(m), f1, f2)

  xx = Float64.(reverse(x, dims=1))
  parent.x = D == 1 ? vec(xx) : xx

  return NFFT3Plan{T,D}(parent)
end

function NFCT3Plan(x::Matrix{T}, N::NTuple{D,Int};
              dims::Union{Integer,UnitRange{Int64}}=1:D,
              precompute::PrecomputeFlags=TENSOR, sortNodes=false,
              fftflags=UInt32(NFFT3.FFTW_ESTIMATE),
              kwargs...) where {D,T}

  if dims != 1:D
    error("NFFT3 does not support directional plans!")
  end

  m, σ, reltol = accuracyParams(; kwargs...)

  if precompute == AbstractNFFTs.LINEAR
    prePsi =  NFFT3.PRE_LIN_PSI
  elseif precompute == AbstractNFFTs.FULL
    prePsi = NFFT3.PRE_FULL_PSI
  elseif precompute == AbstractNFFTs.TENSOR
    prePsi = NFFT3.PRE_PSI
  else
    error("Precompute $(precompute) not supported by NFFT3!")
  end
  sortN = sortNodes ? NFFT3.NFCT_SORT_NODES : UInt32(0)

  f1 = UInt32(
    NFFT3.PRE_PHI_HUT |
    prePsi |
    # NFFT3.MALLOC_X |
    # NFFT3.MALLOC_F_HAT |
    # NFFT3.MALLOC_F |
    NFFT3.FFTW_INIT |
    NFFT3.FFT_OUT_OF_PLACE |
    sortN |
    NFFT3.NFCT_OMP_BLOCKWISE_ADJOINT
     )

  f2 = UInt32(fftflags | NFFT3.FFTW_DESTROY_INPUT)

  n = ntuple(d -> (ceil(Int,σ*N[d])÷2)*2, D) # ensure that n is an even integer

  parent = NFFT3.NFCT(Int32.(reverse(N)), size(x,2), Int32.(reverse(n)), Int32(m), f1, f2)

  xx = Float64.(reverse(x, dims=1))
  parent.x = D == 1 ? vec(xx) : xx

  return NFCT3Plan{T,D}(parent)
end

function NFST3Plan(x::Matrix{T}, N::NTuple{D,Int};
              dims::Union{Integer,UnitRange{Int64}}=1:D,
              precompute::PrecomputeFlags=TENSOR, sortNodes=false,
              fftflags=UInt32(NFFT3.FFTW_ESTIMATE),
              kwargs...) where {D,T}

  if dims != 1:D
    error("NFFT3 does not support directional plans!")
  end

  m, σ, reltol = accuracyParams(; kwargs...)

  if precompute == AbstractNFFTs.LINEAR
    prePsi =  NFFT3.PRE_LIN_PSI
  elseif precompute == AbstractNFFTs.FULL
    prePsi = NFFT3.PRE_FULL_PSI
  elseif precompute == AbstractNFFTs.TENSOR
    prePsi = NFFT3.PRE_PSI
  else
    error("Precompute $(precompute) not supported by NFFT3!")
  end
  sortN = sortNodes ? NFFT3.NFCT_SORT_NODES : UInt32(0)

  f1 = UInt32(
    NFFT3.PRE_PHI_HUT |
    prePsi |
    # NFFT3.MALLOC_X |
    # NFFT3.MALLOC_F_HAT |
    # NFFT3.MALLOC_F |
    NFFT3.FFTW_INIT |
    NFFT3.FFT_OUT_OF_PLACE |
    sortN |
    NFFT3.NFCT_OMP_BLOCKWISE_ADJOINT
     )

  f2 = UInt32(fftflags | NFFT3.FFTW_DESTROY_INPUT)

  n = ntuple(d -> (ceil(Int,σ*N[d])÷2)*2, D) # ensure that n is an even integer

  parent = NFFT3.NFST(Int32.(reverse(N)), size(x,2), Int32.(reverse(n)), Int32(m), f1, f2)

  xx = Float64.(reverse(x, dims=1))

  parent.x = D == 1 ? vec(xx) : xx

  return NFST3Plan{T,D}(parent)
end

function Base.show(io::IO, p::NFFT3Plan{D}) where {D}
  print(io, "NFFT3Plan")
end

function Base.show(io::IO, p::NFCT3Plan{D}) where {D}
  print(io, "NFCT3Plan")
end

function Base.show(io::IO, p::NFST3Plan{D}) where {D}
  print(io, "NFST3Plan")
end

AbstractNFFTs.size_in(p::NFFT3Plan) = reverse(Int.(p.parent.N))
AbstractNFFTs.size_out(p::NFFT3Plan) = (Int(p.parent.M),)

AbstractNFFTs.size_in(p::NFCT3Plan) = reverse(Int.(p.parent.N))
AbstractNFFTs.size_out(p::NFCT3Plan) = (Int(p.parent.M),)

AbstractNFFTs.size_in(p::NFST3Plan) = reverse(Int.(p.parent.N .- 1))
AbstractNFFTs.size_out(p::NFST3Plan) = (Int(p.parent.M),)

function AbstractNFFTs.plan_nfft(::Type{<:Array}, x::Matrix{T}, N::NTuple{D,Int}, rest...;
                   timing::Union{Nothing,TimingStats} = nothing, kargs...) where {T,D}
  t = @elapsed begin
    p = NFFT3Plan(x, N, rest...; kargs...)
  end
  if timing != nothing
    timing.pre = t
  end
  return p
end

function AbstractNFFTs.plan_nfct(::Type{<:Array}, x::Matrix{T}, N::NTuple{D,Int}, rest...;
                   timing::Union{Nothing,TimingStats} = nothing, kargs...) where {T,D}
  t = @elapsed begin
    p = NFCT3Plan(x, N, rest...; kargs...)
  end
  if timing != nothing
    timing.pre = t
  end
  return p
end

function AbstractNFFTs.plan_nfst(::Type{<:Array}, x::Matrix{T}, N::NTuple{D,Int}, rest...;
                   timing::Union{Nothing,TimingStats} = nothing, kargs...) where {T,D}
  t = @elapsed begin
    p = NFST3Plan(x, N, rest...; kargs...)
  end
  if timing != nothing
    timing.pre = t
  end
  return p
end

function LinearAlgebra.mul!(fHat::StridedArray, p::NFFT3Plan{D}, f::AbstractArray;
             verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {D}
  #consistencyCheck(p, f, fHat)

  p.parent.fhat = vec(f)
  p.parent.f = vec(fHat)
  NFFT3.nfft_trafo(p.parent)
  fHat[:] .= p.parent.f

  return fHat
end

function LinearAlgebra.mul!(f::StridedArray, pl::Adjoint{Complex{T},<:NFFT3Plan{T}}, fHat::AbstractArray;
                     verbose=false, timing::Union{Nothing,TimingStats} = nothing) where T
  #consistencyCheck(p, f, fHat)
  p = pl.parent

  p.parent.f = vec(fHat)
  p.parent.fhat = vec(f)
  tadjoint = fApprox = NFFT3.nfft_adjoint(p.parent)
  f[:] .= p.parent.fhat

  return f
end

function LinearAlgebra.mul!(fHat::StridedArray, p::NFCT3Plan{D}, f::AbstractArray;
             verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {D}

  p.parent.fhat = vec(f)
  p.parent.f = vec(fHat)
  NFFT3.nfct_trafo(p.parent)
  fHat[:] .= p.parent.f

  return fHat
end

function LinearAlgebra.mul!(f::StridedArray, pl::Transpose{T,<:NFCT3Plan{T}}, fHat::AbstractArray;
                     verbose=false, timing::Union{Nothing,TimingStats} = nothing) where T
  p = pl.parent

  p.parent.f = vec(fHat)
  p.parent.fhat = vec(f)
  tadjoint = fApprox = NFFT3.nfct_adjoint(p.parent)
  f[:] .= p.parent.fhat

  return f
end

function LinearAlgebra.mul!(fHat::StridedArray, p::NFST3Plan{D}, f::AbstractArray;
             verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {D}

  p.parent.fhat = vec(f)
  p.parent.f = vec(fHat)
  NFFT3.nfst_trafo(p.parent)
  fHat[:] .= p.parent.f

  return fHat
end

function LinearAlgebra.mul!(f::StridedArray, pl::Transpose{T,<:NFST3Plan{T}}, fHat::AbstractArray;
                     verbose=false, timing::Union{Nothing,TimingStats} = nothing) where T
  p = pl.parent

  p.parent.f = vec(fHat)
  p.parent.fhat = vec(f)
  tadjoint = fApprox = NFFT3.nfst_adjoint(p.parent)
  f[:] .= p.parent.fhat

  return f
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((b*coth(c+d*x))^7//2, x) == :(b^7//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+b^7//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-2*b^3*d^-1*(b*coth(c+d*x))^(1/2)+-2//5*b*d^-1*(b*coth(c+d*x))^5//2)
@test integrate((b*coth(c+d*x))^5//2, x) == :(b^5//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-1*b^5//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+-2//3*b*d^-1*(b*coth(c+d*x))^3//2)
@test integrate((b*coth(c+d*x))^3//2, x) == :(b^3//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+b^3//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-2*b*d^-1*(b*coth(c+d*x))^(1/2))
@test integrate((b*coth(c+d*x))^(1/2), x) == :(b^(1/2)*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-1*b^(1/2)*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2)))
@test integrate((b*coth(c+d*x))^-1//2, x) == :(b^-1//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+b^-1//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2)))
@test integrate((b*coth(c+d*x))^-3//2, x) == :(b^-3//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-1*b^-3//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+-2*b^-1*d^-1*(b*coth(c+d*x))^-1//2)
@test integrate((b*coth(c+d*x))^-5//2, x) == :(b^-5//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+b^-5//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-2//3*b^-1*d^-1*(b*coth(c+d*x))^-3//2)
@test integrate((b*coth(c+d*x))^-7//2, x) == :(b^-7//2*d^-1*arctanh(b^-1//2*(b*coth(c+d*x))^(1/2))+-1*b^-7//2*d^-1*arctan(b^-1//2*(b*coth(c+d*x))^(1/2))+-2*b^-3*d^-1*(b*coth(c+d*x))^-1//2+-2//5*b^-1*d^-1*(b*coth(c+d*x))^-5//2)
@test integrate((b*coth(c+d*x))^4//3, x) == :(b^4//3*d^-1*arctanh(b^-1//3*(b*coth(c+d*x))^1//3)+-3*b*d^-1*(b*coth(c+d*x))^1//3+-1//4*b^4//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+-1*b^1//3*(b*coth(c+d*x))^1//3)+1//4*b^4//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+b^1//3*(b*coth(c+d*x))^1//3)+(1/2)*3^(1/2)*b^4//3*d^-1*arctan(1//3*3^(1/2)*(1+2*b^-1//3*(b*coth(c+d*x))^1//3))+-1//2*3^(1/2)*b^4//3*d^-1*arctan(1//3*3^(1/2)*(1+-2*b^-1//3*(b*coth(c+d*x))^1//3)))
@test integrate((b*coth(c+d*x))^2//3, x) == :(b^2//3*d^-1*arctanh(b^-1//3*(b*coth(c+d*x))^1//3)+-1//4*b^2//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+-1*b^1//3*(b*coth(c+d*x))^1//3)+1//4*b^2//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+b^1//3*(b*coth(c+d*x))^1//3)+(1/2)*3^(1/2)*b^2//3*d^-1*arctan(1//3*3^(1/2)*(1+-2*b^-1//3*(b*coth(c+d*x))^1//3))+-1//2*3^(1/2)*b^2//3*d^-1*arctan(1//3*3^(1/2)*(1+2*b^-1//3*(b*coth(c+d*x))^1//3)))
@test integrate((b*coth(c+d*x))^1//3, x) == :(-1//2*b^1//3*d^-1*log(b^2//3+-1*(b*coth(c+d*x))^2//3)+1//4*b^1//3*d^-1*log(b^4//3+(b*coth(c+d*x))^4//3+b^2//3*(b*coth(c+d*x))^2//3)+-1//2*3^(1/2)*b^1//3*d^-1*arctan(1//3*3^(1/2)*b^-2//3*(b^2//3+2*(b*coth(c+d*x))^2//3)))
@test integrate((b*coth(c+d*x))^-1//3, x) == :(-1//2*b^-1//3*d^-1*log(b^2//3+-1*(b*coth(c+d*x))^2//3)+1//4*b^-1//3*d^-1*log(b^4//3+(b*coth(c+d*x))^4//3+b^2//3*(b*coth(c+d*x))^2//3)+(1/2)*3^(1/2)*b^-1//3*d^-1*arctan(1//3*3^(1/2)*b^-2//3*(b^2//3+2*(b*coth(c+d*x))^2//3)))
@test integrate((b*coth(c+d*x))^-2//3, x) == :(b^-2//3*d^-1*arctanh(b^-1//3*(b*coth(c+d*x))^1//3)+-1//4*b^-2//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+-1*b^1//3*(b*coth(c+d*x))^1//3)+1//4*b^-2//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+b^1//3*(b*coth(c+d*x))^1//3)+(1/2)*3^(1/2)*b^-2//3*d^-1*arctan(1//3*3^(1/2)*(1+2*b^-1//3*(b*coth(c+d*x))^1//3))+-1//2*3^(1/2)*b^-2//3*d^-1*arctan(1//3*3^(1/2)*(1+-2*b^-1//3*(b*coth(c+d*x))^1//3)))
@test integrate((b*coth(c+d*x))^-4//3, x) == :(b^-4//3*d^-1*arctanh(b^-1//3*(b*coth(c+d*x))^1//3)+-3*b^-1*d^-1*(b*coth(c+d*x))^-1//3+-1//4*b^-4//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+-1*b^1//3*(b*coth(c+d*x))^1//3)+1//4*b^-4//3*d^-1*log(b^2//3+(b*coth(c+d*x))^2//3+b^1//3*(b*coth(c+d*x))^1//3)+(1/2)*3^(1/2)*b^-4//3*d^-1*arctan(1//3*3^(1/2)*(1+-2*b^-1//3*(b*coth(c+d*x))^1//3))+-1//2*3^(1/2)*b^-4//3*d^-1*arctan(1//3*3^(1/2)*(1+2*b^-1//3*(b*coth(c+d*x))^1//3)))
@test integrate((b*coth(c+d*x)^2)^3//2, x) == :(-1//2*b*d^-1*(b*coth(c+d*x)^2)^(1/2)*coth(c+d*x)+b*d^-1*(b*coth(c+d*x)^2)^(1/2)*log(sinh(c+d*x))*tanh(c+d*x))
@test integrate((b*coth(c+d*x)^2)^(1/2), x) == :(d^-1*(b*coth(c+d*x)^2)^(1/2)*log(sinh(c+d*x))*tanh(c+d*x))
@test integrate((b*coth(c+d*x)^2)^-1//2, x) == :(d^-1*(b*coth(c+d*x)^2)^-1//2*coth(c+d*x)*log(cosh(c+d*x)))
@test integrate((b*coth(c+d*x)^2)^-3//2, x) == :(-1//2*b^-1*d^-1*(b*coth(c+d*x)^2)^-1//2*tanh(c+d*x)+b^-1*d^-1*(b*coth(c+d*x)^2)^-1//2*coth(c+d*x)*log(cosh(c+d*x)))
@test integrate((b*coth(c+d*x)^2)^4//3, x) == :(-3//5*b*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)+b*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*arctanh(coth(c+d*x)^1//3)+-1//4*b*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*b*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*b*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3))+-1//2*b*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^2)^2//3, x) == :(-3*d^-1*(b*coth(c+d*x)^2)^2//3*tanh(c+d*x)+d^-1*(b*coth(c+d*x)^2)^2//3*coth(c+d*x)^-4//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^2)^2//3*coth(c+d*x)^-4//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^2)^2//3*coth(c+d*x)^-4//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^2//3*coth(c+d*x)^-4//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^2//3*coth(c+d*x)^-4//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^2)^1//3, x) == :(d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^1//3*coth(c+d*x)^-2//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^2)^-1//3, x) == :(d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^2)^-2//3, x) == :(-3*d^-1*(b*coth(c+d*x)^2)^-2//3*coth(c+d*x)+d^-1*(b*coth(c+d*x)^2)^-2//3*coth(c+d*x)^4//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^2)^-2//3*coth(c+d*x)^4//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^2)^-2//3*coth(c+d*x)^4//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^-2//3*coth(c+d*x)^4//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^2)^-2//3*coth(c+d*x)^4//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^2)^-4//3, x) == :(-3//5*b^-1*d^-1*(b*coth(c+d*x)^2)^-1//3*tanh(c+d*x)+b^-1*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*arctanh(coth(c+d*x)^1//3)+-1//4*b^-1*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*b^-1*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*b^-1*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*b^-1*d^-1*(b*coth(c+d*x)^2)^-1//3*coth(c+d*x)^2//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^3)^3//2, x) == :(-2//3*b*d^-1*(b*coth(c+d*x)^3)^(1/2)+-2//7*b*d^-1*(b*coth(c+d*x)^3)^(1/2)*coth(c+d*x)^2+b*d^-1*(b*coth(c+d*x)^3)^(1/2)*coth(c+d*x)^-3//2*arctanh(coth(c+d*x)^(1/2))+-1*b*d^-1*(b*coth(c+d*x)^3)^(1/2)*coth(c+d*x)^-3//2*arctan(coth(c+d*x)^(1/2)))
@test integrate((b*coth(c+d*x)^3)^(1/2), x) == :(-2*d^-1*(b*coth(c+d*x)^3)^(1/2)*tanh(c+d*x)+d^-1*(b*coth(c+d*x)^3)^(1/2)*coth(c+d*x)^-3//2*arctan(coth(c+d*x)^(1/2))+d^-1*(b*coth(c+d*x)^3)^(1/2)*coth(c+d*x)^-3//2*arctanh(coth(c+d*x)^(1/2)))
@test integrate((b*coth(c+d*x)^3)^-1//2, x) == :(-2*d^-1*(b*coth(c+d*x)^3)^-1//2*coth(c+d*x)+d^-1*(b*coth(c+d*x)^3)^-1//2*coth(c+d*x)^3//2*arctanh(coth(c+d*x)^(1/2))+-1*d^-1*(b*coth(c+d*x)^3)^-1//2*coth(c+d*x)^3//2*arctan(coth(c+d*x)^(1/2)))
@test integrate((b*coth(c+d*x)^3)^-3//2, x) == :(-2//3*b^-1*d^-1*(b*coth(c+d*x)^3)^-1//2+-2//7*b^-1*d^-1*(b*coth(c+d*x)^3)^-1//2*tanh(c+d*x)^2+b^-1*d^-1*(b*coth(c+d*x)^3)^-1//2*coth(c+d*x)^3//2*arctan(coth(c+d*x)^(1/2))+b^-1*d^-1*(b*coth(c+d*x)^3)^-1//2*coth(c+d*x)^3//2*arctanh(coth(c+d*x)^(1/2)))
@test integrate((b*coth(c+d*x)^3)^4//3, x) == :(-1*b*d^-1*(b*coth(c+d*x)^3)^1//3+b*x*(b*coth(c+d*x)^3)^1//3*tanh(c+d*x)+-1//3*b*d^-1*(b*coth(c+d*x)^3)^1//3*coth(c+d*x)^2)
@test integrate((b*coth(c+d*x)^3)^2//3, x) == :(x*(b*coth(c+d*x)^3)^2//3*tanh(c+d*x)^2+-1*d^-1*(b*coth(c+d*x)^3)^2//3*tanh(c+d*x))
@test integrate((b*coth(c+d*x)^3)^1//3, x) == :(d^-1*(b*coth(c+d*x)^3)^1//3*log(sinh(c+d*x))*tanh(c+d*x))
@test integrate((b*coth(c+d*x)^3)^-1//3, x) == :(d^-1*(b*coth(c+d*x)^3)^-1//3*coth(c+d*x)*log(cosh(c+d*x)))
@test integrate((b*coth(c+d*x)^3)^-2//3, x) == :(x*(b*coth(c+d*x)^3)^-2//3*coth(c+d*x)^2+-1*d^-1*(b*coth(c+d*x)^3)^-2//3*coth(c+d*x))
@test integrate((b*coth(c+d*x)^3)^-4//3, x) == :(-1*b^-1*d^-1*(b*coth(c+d*x)^3)^-1//3+x*b^-1*(b*coth(c+d*x)^3)^-1//3*coth(c+d*x)+-1//3*b^-1*d^-1*(b*coth(c+d*x)^3)^-1//3*tanh(c+d*x)^2)
@test integrate((b*coth(c+d*x)^4)^3//2, x) == :(b*x*(b*coth(c+d*x)^4)^(1/2)*tanh(c+d*x)^2+-1*b*d^-1*(b*coth(c+d*x)^4)^(1/2)*tanh(c+d*x)+-1//3*b*d^-1*(b*coth(c+d*x)^4)^(1/2)*coth(c+d*x)+-1//5*b*d^-1*(b*coth(c+d*x)^4)^(1/2)*coth(c+d*x)^3)
@test integrate((b*coth(c+d*x)^4)^(1/2), x) == :(x*(b*coth(c+d*x)^4)^(1/2)*tanh(c+d*x)^2+-1*d^-1*(b*coth(c+d*x)^4)^(1/2)*tanh(c+d*x))
@test integrate((b*coth(c+d*x)^4)^-1//2, x) == :(x*(b*coth(c+d*x)^4)^-1//2*coth(c+d*x)^2+-1*d^-1*(b*coth(c+d*x)^4)^-1//2*coth(c+d*x))
@test integrate((b*coth(c+d*x)^4)^-3//2, x) == :(x*b^-1*(b*coth(c+d*x)^4)^-1//2*coth(c+d*x)^2+-1*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//2*coth(c+d*x)+-1//3*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//2*tanh(c+d*x)+-1//5*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//2*tanh(c+d*x)^3)
@test integrate((b*coth(c+d*x)^4)^4//3, x) == :(-3*b*d^-1*(b*coth(c+d*x)^4)^1//3*tanh(c+d*x)+-3//7*b*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)+-3//13*b*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^3+b*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*arctanh(coth(c+d*x)^1//3)+-1//4*b*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*b*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*b*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3))+-1//2*b*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^4)^2//3, x) == :(-3//5*d^-1*(b*coth(c+d*x)^4)^2//3*tanh(c+d*x)+d^-1*(b*coth(c+d*x)^4)^2//3*coth(c+d*x)^-8//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^4)^2//3*coth(c+d*x)^-8//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^4)^2//3*coth(c+d*x)^-8//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^2//3*coth(c+d*x)^-8//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^2//3*coth(c+d*x)^-8//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^4)^1//3, x) == :(-3*d^-1*(b*coth(c+d*x)^4)^1//3*tanh(c+d*x)+d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^1//3*coth(c+d*x)^-4//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^4)^-1//3, x) == :(-3*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)+d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^4)^-2//3, x) == :(-3//5*d^-1*(b*coth(c+d*x)^4)^-2//3*coth(c+d*x)+d^-1*(b*coth(c+d*x)^4)^-2//3*coth(c+d*x)^8//3*arctanh(coth(c+d*x)^1//3)+-1//4*d^-1*(b*coth(c+d*x)^4)^-2//3*coth(c+d*x)^8//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*d^-1*(b*coth(c+d*x)^4)^-2//3*coth(c+d*x)^8//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^-2//3*coth(c+d*x)^8//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*d^-1*(b*coth(c+d*x)^4)^-2//3*coth(c+d*x)^8//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3)))
@test integrate((b*coth(c+d*x)^4)^-4//3, x) == :(-3*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)+-3//7*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*tanh(c+d*x)+-3//13*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*tanh(c+d*x)^3+b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*arctanh(coth(c+d*x)^1//3)+-1//4*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*log(1+coth(c+d*x)^2//3+-1*coth(c+d*x)^1//3)+1//4*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*log(1+coth(c+d*x)^1//3+coth(c+d*x)^2//3)+(1/2)*3^(1/2)*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*arctan(1//3*3^(1/2)*(1+-2*coth(c+d*x)^1//3))+-1//2*3^(1/2)*b^-1*d^-1*(b*coth(c+d*x)^4)^-1//3*coth(c+d*x)^4//3*arctan(1//3*3^(1/2)*(1+2*coth(c+d*x)^1//3)))
@test integrate((1+coth(x))^5, x) == :(-8*coth(x)+-2*(1+coth(x))^2+16x+16*log(sinh(x))+-2//3*(1+coth(x))^3+-1//4*(1+coth(x))^4)
@test integrate((1+coth(x))^4, x) == :(-1*(1+coth(x))^2+-4*coth(x)+8x+8*log(sinh(x))+-1//3*(1+coth(x))^3)
@test integrate((1+coth(x))^3, x) == :(-2*coth(x)+4x+4*log(sinh(x))+-1//2*(1+coth(x))^2)
@test integrate((1+coth(x))^2, x) == :(-1*coth(x)+2x+2*log(sinh(x)))
@test integrate((1+coth(x))^-1, x) == :((1/2)*x+-1*(2+2*coth(x))^-1)
@test integrate((1+coth(x))^-2, x) == :(-1*(4+4*coth(x))^-1+-1//4*(1+coth(x))^-2+1//4*x)
@test integrate((1+coth(x))^-3, x) == :(-1*(8+8*coth(x))^-1+-1//6*(1+coth(x))^-3+-1//8*(1+coth(x))^-2+1//8*x)
@test integrate((1+coth(x))^-4, x) == :(-1*(16+16*coth(x))^-1+-1//8*(1+coth(x))^-4+-1//12*(1+coth(x))^-3+-1//16*(1+coth(x))^-2+1//16*x)
@test integrate((1+coth(x))^-5, x) == :(-1*(32+32*coth(x))^-1+-1//10*(1+coth(x))^-5+-1//16*(1+coth(x))^-4+-1//24*(1+coth(x))^-3+-1//32*(1+coth(x))^-2+1//32*x)
@test integrate((1+coth(x))^7//2, x) == :(-8*(1+coth(x))^(1/2)+-4//3*(1+coth(x))^3//2+-2//5*(1+coth(x))^5//2+8*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^5//2, x) == :(-4*(1+coth(x))^(1/2)+-2//3*(1+coth(x))^3//2+4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^3//2, x) == :(-2*(1+coth(x))^(1/2)+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^(1/2), x) == :(2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-1//2, x) == :(-1*(1+coth(x))^-1//2+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-3//2, x) == :(-1//2*(1+coth(x))^-1//2+-1//3*(1+coth(x))^-3//2+1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-5//2, x) == :(-1//4*(1+coth(x))^-1//2+-1//5*(1+coth(x))^-5//2+-1//6*(1+coth(x))^-3//2+1//8*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((a+b*coth(c+d*x))^5, x) == :(a*x*(a^4+5*b^4+10*a^2*b^2)+-1//4*b*d^-1*(a+b*coth(c+d*x))^4+b*d^-1*(b^4+5*a^4+10*a^2*b^2)*log(sinh(c+d*x))+-2//3*a*b*d^-1*(a+b*coth(c+d*x))^3+-1//2*b*d^-1*(a+b*coth(c+d*x))^2*(b^2+3*a^2)+-4*a*b^2*d^-1*(a^2+b^2)*coth(c+d*x))
@test integrate((a+b*coth(c+d*x))^4, x) == :(x*(a^4+b^4+6*a^2*b^2)+-1//3*b*d^-1*(a+b*coth(c+d*x))^3+-1*a*b*d^-1*(a+b*coth(c+d*x))^2+-1*b^2*d^-1*(b^2+3*a^2)*coth(c+d*x)+4*a*b*d^-1*(a^2+b^2)*log(sinh(c+d*x)))
@test integrate((a+b*coth(c+d*x))^3, x) == :(a*x*(a^2+3*b^2)+-1//2*b*d^-1*(a+b*coth(c+d*x))^2+b*d^-1*(b^2+3*a^2)*log(sinh(c+d*x))+-2*a*b^2*d^-1*coth(c+d*x))
@test integrate((a+b*coth(c+d*x))^2, x) == :(x*(a^2+b^2)+-1*b^2*d^-1*coth(c+d*x)+2*a*b*d^-1*log(sinh(c+d*x)))
@test integrate((a+b*coth(c+d*x))^-1, x) == :(a*x*(a^2+-1*b^2)^-1+-1*b*d^-1*(a^2+-1*b^2)^-1*log(a*sinh(c+d*x)+b*cosh(c+d*x)))
@test integrate((a+b*coth(c+d*x))^-2, x) == :(x*(a^2+-1*b^2)^-2*(a^2+b^2)+b*d^-1*(a+b*coth(c+d*x))^-1*(a^2+-1*b^2)^-1+-2*a*b*d^-1*(a^2+-1*b^2)^-2*log(a*sinh(c+d*x)+b*cosh(c+d*x)))
@test integrate((a+b*coth(c+d*x))^-3, x) == :(a*x*(a^2+-1*b^2)^-3*(a^2+3*b^2)+b*d^-1*(a+b*coth(c+d*x))^-2*(-2*b^2+2*a^2)^-1+-1*b*d^-1*(a^2+-1*b^2)^-3*(b^2+3*a^2)*log(a*sinh(c+d*x)+b*cosh(c+d*x))+2*a*b*d^-1*(a+b*coth(c+d*x))^-1*(a^2+-1*b^2)^-2)
@test integrate((a+b*coth(c+d*x))^-4, x) == :(x*(a^2+-1*b^2)^-4*(a^4+b^4+6*a^2*b^2)+b*d^-1*(a+b*coth(c+d*x))^-3*(-3*b^2+3*a^2)^-1+a*b*d^-1*(a+b*coth(c+d*x))^-2*(a^2+-1*b^2)^-2+b*d^-1*(a+b*coth(c+d*x))^-1*(a^2+-1*b^2)^-3*(b^2+3*a^2)+-4*a*b*d^-1*(a^2+-1*b^2)^-4*(a^2+b^2)*log(a*sinh(c+d*x)+b*cosh(c+d*x)))
@test integrate((4+6*coth(c+d*x))^-1, x) == :(-1//5*x+3//10*d^-1*log(2*sinh(c+d*x)+3*cosh(c+d*x)))
@test integrate((4+-6*coth(c+d*x))^-1, x) == :(-1//5*x+-3//10*d^-1*log(-2*sinh(c+d*x)+3*cosh(c+d*x)))
@test integrate((a+b*coth(c+d*x))^(1/2), x) == :(d^-1*(a+b)^(1/2)*arctanh((a+b)^-1//2*(a+b*coth(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^(1/2)*arctanh((a+-1b)^-1//2*(a+b*coth(c+d*x))^(1/2)))
@test integrate((a+b*coth(c+d*x))^-1//2, x) == :(d^-1*(a+b)^-1//2*arctanh((a+b)^-1//2*(a+b*coth(c+d*x))^(1/2))+-1*d^-1*(a+-1b)^-1//2*arctanh((a+-1b)^-1//2*(a+b*coth(c+d*x))^(1/2)))
@test integrate((1+coth(x))^-1*sinh(x)^4, x) == :((8+-8*coth(x))^-1+-3*(16+16*coth(x))^-1+-3//32*(1+coth(x))^-2+-1//24*(1+coth(x))^-3+1//32*(1+-1*coth(x))^-2+5//16*x)
@test integrate((1+coth(x))^-1*sinh(x)^3, x) == :(-4//5*cosh(x)+4//15*cosh(x)^3+-1*(5+5*coth(x))^-1*sinh(x)^3)
@test integrate((1+coth(x))^-1*sinh(x)^2, x) == :((4+4*coth(x))^-1+-1*(8+-8*coth(x))^-1+-3//8*x+1//8*(1+coth(x))^-2)
@test integrate((1+coth(x))^-1*sinh(x), x) == :(2//3*cosh(x)+-1*(3+3*coth(x))^-1*sinh(x))
@test integrate((1+coth(x))^-1*csch(x), x) == :(-1*(1+coth(x))^-1*csch(x))
@test integrate((1+coth(x))^-1*csch(x)^2, x) == :(-1*log(1+coth(x)))
@test integrate((1+coth(x))^-1*csch(x)^3, x) == :(-1*csch(x)+arctanh(cosh(x)))
@test integrate((1+coth(x))^-1*csch(x)^4, x) == :(-1//2*coth(x)^2+coth(x))
@test integrate((a+b*coth(x))^-1*sinh(x)^4, x) == :(-1*b^5*(a^2+-1*b^2)^-3*log(a+b*coth(x))+-1*(-4*b^2+4*a^2)^-1*sinh(x)^4*(b+-1*a*coth(x))+-1//8*(a^2+-1*b^2)^-2*sinh(x)^2*(4*b^3+-1*a*b^2*(7+-3*a^2*b^-2)*coth(x))+-1//16*(a+b)^-3*(3*a^2+8*b^2+9*a*b)*log(1+-1*coth(x))+1//16*(a+-1b)^-3*(3*a^2+8*b^2+-9*a*b)*log(1+coth(x)))
@test integrate((a+b*coth(x))^-1*sinh(x)^3, x) == :(a*(-3*b^2+3*a^2)^-1*cosh(x)^3+-1*a*(a^2+-1*b^2)^-1*cosh(x)+-1*b*(-3*b^2+3*a^2)^-1*sinh(x)^3+-1*b^3*(a^2+-1*b^2)^-2*sinh(x)+-1*b^4*(a^2+-1*b^2)^-5//2*arctanh((a^2+-1*b^2)^-1//2*(b+a*coth(x))*sinh(x))+a*b^2*(a^2+-1*b^2)^-2*cosh(x))
@test integrate((a+b*coth(x))^-1*sinh(x)^2, x) == :(-1*b^3*(a^2+-1*b^2)^-2*log(a+b*coth(x))+-1*(-2*b^2+2*a^2)^-1*sinh(x)^2*(b+-1*a*coth(x))+-1//4*(a+-1b)^-2*(a+-2b)*log(1+coth(x))+1//4*(a+b)^-2*(a+2b)*log(1+-1*coth(x)))
@test integrate((a+b*coth(x))^-1*sinh(x), x) == :(a*(a^2+-1*b^2)^-1*cosh(x)+-1*b*(a^2+-1*b^2)^-1*sinh(x)+-1*b^2*(a^2+-1*b^2)^-3//2*arctanh((a^2+-1*b^2)^-1//2*(b+a*coth(x))*sinh(x)))
@test integrate((a+b*coth(x))^-1*csch(x), x) == :(-1*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(b+a*coth(x))*sinh(x)))
@test integrate((a+b*coth(x))^-1*csch(x)^2, x) == :(-1*b^-1*log(a+b*coth(x)))
@test integrate((a+b*coth(x))^-1*csch(x)^3, x) == :(-1*b^-1*csch(x)+a*b^-2*arctanh(cosh(x))+-1*b^-2*(a^2+-1*b^2)^(1/2)*arctanh((a^2+-1*b^2)^-1//2*(b+a*coth(x))*sinh(x)))
@test integrate((a+b*coth(x))^-1*csch(x)^4, x) == :(-1//2*b^-1*coth(x)^2+a*b^-2*coth(x)+-1*b^-3*(a^2+-1*b^2)*log(a+b*coth(x)))
@test integrate((1+coth(x))^-1*cosh(x)^4, x) == :(-1*(8+-8*coth(x))^-1+-3*(16+16*coth(x))^-1+-1//24*(1+coth(x))^-3+1//16*x+1//32*(1+-1*coth(x))^-2+5//32*(1+coth(x))^-2)
@test integrate((1+coth(x))^-1*cosh(x)^3, x) == :(-1//3*sinh(x)^3+-1//5*sinh(x)^5+1//5*cosh(x)^5)
@test integrate((1+coth(x))^-1*cosh(x)^2, x) == :(-1*(4+4*coth(x))^-1+-1*(8+-8*coth(x))^-1+1//8*x+1//8*(1+coth(x))^-2)
@test integrate((1+coth(x))^-1*cosh(x), x) == :(-1//3*sinh(x)^3+1//3*cosh(x)^3)
@test integrate((1+coth(x))^-1*sech(x), x) == :(-1*sinh(x)+arctan(sinh(x))+cosh(x))
@test integrate((1+coth(x))^-1*sech(x)^2, x) == :(-1*log(1+coth(x))+-1*log(tanh(x))+tanh(x))
@test integrate((1+coth(x))^-1*sech(x)^3, x) == :(-1*sech(x)+-1//2*arctan(sinh(x))+(1/2)*sech(x)*tanh(x))
@test integrate((1+coth(x))^-1*sech(x)^4, x) == :((1/2)*tanh(x)^2+-1//3*tanh(x)^3)
@test integrate((1+coth(x))^(1/2)*sech(x)^2, x) == :((1+coth(x))^(1/2)*tanh(x)+arctanh((1+coth(x))^(1/2)))
@test integrate((a+b*coth(x))^-1*cosh(x)^4, x) == :(-1*(-4*b^2+4*a^2)^-1*sinh(x)^4*(b+-1*a*coth(x))+-1//8*(a^2+-1*b^2)^-2*sinh(x)^2*(4*b*(-1*b^2+2*a^2)+-1*a*(-1*b^2+5*a^2)*coth(x))+-1*b*a^4*(a^2+-1*b^2)^-3*log(a+b*coth(x))+-1//16*a*(a+b)^-3*(b+3a)*log(1+-1*coth(x))+1//16*a*(a+-1b)^-3*(-1b+3a)*log(1+coth(x)))
@test integrate((a+b*coth(x))^-1*cosh(x)^3, x) == :(a*(a^2+-1*b^2)^-1*sinh(x)+a*(-3*b^2+3*a^2)^-1*sinh(x)^3+-1*b*(-3*b^2+3*a^2)^-1*cosh(x)^3+a*b^2*(a^2+-1*b^2)^-2*sinh(x)+b*a^3*(a^2+-1*b^2)^-5//2*arctanh((a^2+-1*b^2)^-1//2*(a*cosh(x)+b*sinh(x)))+-1*b*a^2*(a^2+-1*b^2)^-2*cosh(x))
@test integrate((a+b*coth(x))^-1*cosh(x)^2, x) == :(-1*(-2*b^2+2*a^2)^-1*sinh(x)^2*(b+-1*a*coth(x))+-1//4*a*(a+b)^-2*log(1+-1*coth(x))+1//4*a*(a+-1b)^-2*log(1+coth(x))+-1*b*a^2*(a^2+-1*b^2)^-2*log(a+b*coth(x)))
@test integrate((a+b*coth(x))^-1*cosh(x), x) == :(a*(a^2+-1*b^2)^-1*sinh(x)+-1*b*(a^2+-1*b^2)^-1*cosh(x)+a*b*(a^2+-1*b^2)^-3//2*arctanh((a^2+-1*b^2)^-1//2*(a*cosh(x)+b*sinh(x))))
@test integrate((a+b*coth(x))^-1*sech(x), x) == :(a^-1*arctan(sinh(x))+b*a^-1*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a*cosh(x)+b*sinh(x))))
@test integrate((a+b*coth(x))^-1*sech(x)^2, x) == :(a^-1*tanh(x)+-1*b*a^-2*log(a+b*coth(x))+-1*b*a^-2*log(tanh(x)))
@test integrate((a+b*coth(x))^-1*sech(x)^3, x) == :((1/2)*a^-1*arctan(sinh(x))+(1/2)*a^-1*sech(x)*tanh(x)+-1*b*a^-2*sech(x)+-1*a^-3*b^2*arctan(sinh(x))+b*a^-3*(a^2+-1*b^2)^(1/2)*arctanh((a^2+-1*b^2)^-1//2*(a*cosh(x)+b*sinh(x))))
@test integrate((a+b*coth(x))^-1*sech(x)^4, x) == :(-1//3*a^-1*tanh(x)^3+a^-3*(a^2+-1*b^2)*tanh(x)+(1/2)*b*a^-2*tanh(x)^2+-1*b*a^-4*(a^2+-1*b^2)*log(a+b*coth(x))+-1*b*a^-4*(a^2+-1*b^2)*log(tanh(x)))
@test integrate((im+2*coth(x))^-1*sech(x), x) == :(-1*im*arctan(sinh(x))+-2//5*5^(1/2)*arctanh(1//5*5^(1/2)*(-2*im*sinh(x)+cosh(x))))
@test integrate((1+coth(x))^-1*tanh(x)^4, x) == :(tanh(x)^2+-2*log(cosh(x))+-5//2*tanh(x)+-5//6*tanh(x)^3+5//2*x+(2+2*coth(x))^-1*tanh(x)^3)
@test integrate((1+coth(x))^-1*tanh(x)^3, x) == :(-1*tanh(x)^2+2*log(cosh(x))+-3//2*x+3//2*tanh(x)+(2+2*coth(x))^-1*tanh(x)^2)
@test integrate((1+coth(x))^-1*tanh(x)^2, x) == :(-1*log(cosh(x))+-3//2*tanh(x)+3//2*x+(2+2*coth(x))^-1*tanh(x))
@test integrate((1+coth(x))^-1*tanh(x), x) == :((2+2*coth(x))^-1+-1//2*x+log(cosh(x)))
@test integrate((1+coth(x))^-1, x) == :((1/2)*x+-1*(2+2*coth(x))^-1)
@test integrate((1+coth(x))^-1*coth(x), x) == :((2+2*coth(x))^-1+(1/2)*x)
@test integrate((1+coth(x))^-1*coth(x)^2, x) == :(-1*(2+2*coth(x))^-1+-1//2*x+log(sinh(x)))
@test integrate((1+coth(x))^-1*coth(x)^3, x) == :(-1*log(sinh(x))+-3//2*coth(x)+3//2*x+(2+2*coth(x))^-1*coth(x)^2)
@test integrate((1+coth(x))^-1*coth(x)^4, x) == :(-1*coth(x)^2+2*log(sinh(x))+-3//2*x+3//2*coth(x)+(2+2*coth(x))^-1*coth(x)^3)
@test integrate((1+coth(x))^3//2*coth(x), x) == :(-2*(1+coth(x))^(1/2)+-2//3*(1+coth(x))^3//2+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^(1/2)*coth(x), x) == :(-2*(1+coth(x))^(1/2)+2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-1//2*coth(x), x) == :((1+coth(x))^-1//2+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-3//2*coth(x), x) == :(-1//2*(1+coth(x))^-1//2+1//3*(1+coth(x))^-3//2+1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^3//2*coth(x)^2, x) == :(-2*(1+coth(x))^(1/2)+-2//5*(1+coth(x))^5//2+2*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^(1/2)*coth(x)^2, x) == :(-2//3*(1+coth(x))^3//2+2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-1//2*coth(x)^2, x) == :(-1*(1+coth(x))^-1//2+-2*(1+coth(x))^(1/2)+(1/2)*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((1+coth(x))^-3//2*coth(x)^2, x) == :(-1//3*(1+coth(x))^-3//2+3//2*(1+coth(x))^-1//2+1//4*2^(1/2)*arctanh((1/2)*2^(1/2)*(1+coth(x))^(1/2)))
@test integrate((a+b*coth(x))^-1*tanh(x)^4, x) == :(-1//3*a^-1*tanh(x)^3+a*x*(a^2+-1*b^2)^-1+(1/2)*b*a^-2*tanh(x)^2+-1*a^-3*(a^2+b^2)*tanh(x)+-1*b*a^-4*(a^2+b^2)*log(cosh(x))+-1*a^-4*b^5*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x)))
@test integrate((a+b*coth(x))^-1*tanh(x)^3, x) == :(-1//2*a^-1*tanh(x)^2+b*a^-2*tanh(x)+a^-3*(a^2+b^2)*log(cosh(x))+-1*b*x*(a^2+-1*b^2)^-1+a^-3*b^4*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x)))
@test integrate((a+b*coth(x))^-1*tanh(x)^2, x) == :(-1*a^-1*tanh(x)+a*x*(a^2+-1*b^2)^-1+-1*b*a^-2*log(cosh(x))+-1*a^-2*b^3*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x)))
@test integrate((a+b*coth(x))^-1*tanh(x), x) == :(a^-1*log(cosh(x))+-1*b*x*(a^2+-1*b^2)^-1+a^-1*b^2*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x)))
@test integrate((a+b*coth(x))^-1, x) == :(a*x*(a^2+-1*b^2)^-1+-1*b*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x)))
@test integrate((a+b*coth(x))^-1*coth(x), x) == :(a*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x))+-1*b*x*(a^2+-1*b^2)^-1)
@test integrate((a+b*coth(x))^-1*coth(x)^2, x) == :(b^-1*log(sinh(x))+-1*a*x*b^-2+x*a^3*b^-2*(a^2+-1*b^2)^-1+-1*a^2*b^-1*(a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x)))
@test integrate((a+b*coth(x))^-1*coth(x)^3, x) == :(-1*b^-1*coth(x)+a*(a^2+-1*b^2)^-1*log(sinh(x))+-1*b*x*(a^2+-1*b^2)^-1+a^3*b^-2*(a^2+-1*b^2)^-1*log(a+b*coth(x)))
@test integrate((a+b*coth(x))^-1*coth(x)^4, x) == :(-1//2*b^-1*coth(x)^2+a*x*(a^2+-1*b^2)^-1+a*b^-2*coth(x)+-1*b*(a^2+-1*b^2)^-1*log(sinh(x))+-1*a^4*b^-3*(a^2+-1*b^2)^-1*log(a+b*coth(x)))
@test integrate((a+b*coth(x))^-1*coth(x)^5, x) == :(-1//3*b^-1*coth(x)^3+a*(a^2+-1*b^2)^-1*log(sinh(x))+(1/2)*a*b^-2*coth(x)^2+-1*b*x*(a^2+-1*b^2)^-1+-1*b^-3*(a^2+b^2)*coth(x)+a^5*b^-4*(a^2+-1*b^2)^-1*log(a+b*coth(x)))
@test integrate(x*(a+b*coth(x))^-2*csch(x)^2, x) == :((a^2+-1*b^2)^-1*log(a*sinh(x)+b*cosh(x))+x*b^-1*(a+b*coth(x))^-1+-1*a*x*b^-1*(a^2+-1*b^2)^-1)
@test integrate(x^-1*coth(a+b*log(c*x^n)), x) == :(b^-1*n^-1*log(sinh(a+b*log(c*x^n))))
@test integrate(x^-1*coth(a+b*log(c*x^n))^2, x) == :(-1*b^-1*n^-1*coth(a+b*log(c*x^n))+log(x))
@test integrate(x^-1*coth(a+b*log(c*x^n))^3, x) == :(b^-1*n^-1*log(sinh(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*coth(a+b*log(c*x^n))^2)
@test integrate(x^-1*coth(a+b*log(c*x^n))^4, x) == :(-1*b^-1*n^-1*coth(a+b*log(c*x^n))+-1//3*b^-1*n^-1*coth(a+b*log(c*x^n))^3+log(x))
@test integrate(x^-1*coth(a+b*log(c*x^n))^5, x) == :(b^-1*n^-1*log(sinh(a+b*log(c*x^n)))+-1//2*b^-1*n^-1*coth(a+b*log(c*x^n))^2+-1//4*b^-1*n^-1*coth(a+b*log(c*x^n))^4)
@test integrate(x^-1*coth(a+b*log(c*x^n))^5//2, x) == :(b^-1*n^-1*arctanh(coth(a+b*log(c*x^n))^(1/2))+-1*b^-1*n^-1*arctan(coth(a+b*log(c*x^n))^(1/2))+-2//3*b^-1*n^-1*coth(a+b*log(c*x^n))^3//2)
@test integrate(x^-1*coth(a+b*log(c*x^n))^3//2, x) == :(b^-1*n^-1*arctan(coth(a+b*log(c*x^n))^(1/2))+b^-1*n^-1*arctanh(coth(a+b*log(c*x^n))^(1/2))+-2*b^-1*n^-1*coth(a+b*log(c*x^n))^(1/2))
@test integrate(x^-1*coth(a+b*log(c*x^n))^(1/2), x) == :(b^-1*n^-1*arctanh(coth(a+b*log(c*x^n))^(1/2))+-1*b^-1*n^-1*arctan(coth(a+b*log(c*x^n))^(1/2)))
@test integrate(x^-1*coth(a+b*log(c*x^n))^-1//2, x) == :(b^-1*n^-1*arctan(coth(a+b*log(c*x^n))^(1/2))+b^-1*n^-1*arctanh(coth(a+b*log(c*x^n))^(1/2)))
@test integrate(x^-1*coth(a+b*log(c*x^n))^-3//2, x) == :(b^-1*n^-1*arctanh(coth(a+b*log(c*x^n))^(1/2))+-1*b^-1*n^-1*arctan(coth(a+b*log(c*x^n))^(1/2))+-2*b^-1*n^-1*coth(a+b*log(c*x^n))^-1//2)
@test integrate(x^-1*coth(a+b*log(c*x^n))^-5//2, x) == :(b^-1*n^-1*arctan(coth(a+b*log(c*x^n))^(1/2))+b^-1*n^-1*arctanh(coth(a+b*log(c*x^n))^(1/2))+-2//3*b^-1*n^-1*coth(a+b*log(c*x^n))^-3//2)
@test integrate((a+b*coth(x)^2+c*coth(x)^4)^-1//2*coth(x)^5, x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2a+coth(x)^2*(b+2c)))+-1//2*c^-1*(a+b*coth(x)^2+c*coth(x)^4)^(1/2)+1//4*c^-3//2*(b+-2c)*arctanh((1/2)*c^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2*c*coth(x)^2)))
@test integrate((a+b*coth(x)^2+c*coth(x)^4)^-1//2*coth(x)^3, x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2a+coth(x)^2*(b+2c)))+-1//2*c^-1//2*arctanh((1/2)*c^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2*c*coth(x)^2)))
@test integrate((a+b*coth(x)^2+c*coth(x)^4)^-1//2*coth(x), x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2a+coth(x)^2*(b+2c))))
@test integrate((a+b*coth(x)^2+c*coth(x)^4)^-1//2*tanh(x), x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2a+coth(x)^2*(b+2c)))+-1//2*a^-1//2*arctanh((1/2)*a^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(2a+b*coth(x)^2)))
@test integrate((a+b*coth(x)^2+c*coth(x)^4)^-1//2*tanh(x)^3, x) == :((1/2)*(a+b+c)^-1//2*arctanh((1/2)*(a+b+c)^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2a+coth(x)^2*(b+2c)))+-1//2*a^-1//2*arctanh((1/2)*a^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(2a+b*coth(x)^2))+-1//2*a^-1*(a+b*coth(x)^2+c*coth(x)^4)^(1/2)*tanh(x)^2+1//4*b*a^-3//2*arctanh((1/2)*a^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(2a+b*coth(x)^2)))
@test integrate((a+b*coth(x)^2+c*coth(x)^4)^(1/2)*coth(x), x) == :(-1//2*(a+b*coth(x)^2+c*coth(x)^4)^(1/2)+(1/2)*(a+b+c)^(1/2)*arctanh((1/2)*(a+b+c)^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2a+coth(x)^2*(b+2c)))+-1//4*c^-1//2*(b+2c)*arctanh((1/2)*c^-1//2*(a+b*coth(x)^2+c*coth(x)^4)^-1//2*(b+2*c*coth(x)^2)))
@test integrate((coth(a*c+b*c*x)^2)^5//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-15//4*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*arctanh(exp(c*(a+b*x)))*tanh(a*c+b*c*x)+-4*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-4*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-55//6*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+25//4*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-1*exp(c*(a+b*x))*tanh(a*c+b*c*x)+26//3*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-3*exp(c*(a+b*x))*tanh(a*c+b*c*x))
@test integrate((coth(a*c+b*c*x)^2)^3//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-3*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*arctanh(exp(c*(a+b*x)))*tanh(a*c+b*c*x)+-2*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-2*exp(c*(a+b*x))*tanh(a*c+b*c*x)+3*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*(1+-1*exp(2*c*(a+b*x)))^-1*exp(c*(a+b*x))*tanh(a*c+b*c*x))
@test integrate((coth(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*exp(c*(a+b*x))*tanh(a*c+b*c*x)+-2*b^-1*c^-1*(coth(a*c+b*c*x)^2)^(1/2)*arctanh(exp(c*(a+b*x)))*tanh(a*c+b*c*x))
@test integrate((coth(a*c+b*c*x)^2)^-1//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*coth(a*c+b*c*x)*exp(c*(a+b*x))+-2*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*arctan(exp(c*(a+b*x)))*coth(a*c+b*c*x))
@test integrate((coth(a*c+b*c*x)^2)^-3//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*coth(a*c+b*c*x)*exp(c*(a+b*x))+-3*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*arctan(exp(c*(a+b*x)))*coth(a*c+b*c*x)+-2*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*(1+exp(2*c*(a+b*x)))^-2*coth(a*c+b*c*x)*exp(c*(a+b*x))+3*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*(1+exp(2*c*(a+b*x)))^-1*coth(a*c+b*c*x)*exp(c*(a+b*x)))
@test integrate((coth(a*c+b*c*x)^2)^-5//2*exp(c*(a+b*x)), x) == :(b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*coth(a*c+b*c*x)*exp(c*(a+b*x))+-15//4*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*arctan(exp(c*(a+b*x)))*coth(a*c+b*c*x)+-4*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*(1+exp(2*c*(a+b*x)))^-4*coth(a*c+b*c*x)*exp(c*(a+b*x))+-55//6*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*(1+exp(2*c*(a+b*x)))^-2*coth(a*c+b*c*x)*exp(c*(a+b*x))+25//4*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*(1+exp(2*c*(a+b*x)))^-1*coth(a*c+b*c*x)*exp(c*(a+b*x))+26//3*b^-1*c^-1*(coth(a*c+b*c*x)^2)^-1//2*(1+exp(2*c*(a+b*x)))^-3*coth(a*c+b*c*x)*exp(c*(a+b*x)))
@test integrate(sin(coth(a+b*x))^3, x) == :(-3//8*b^-1*cosintegral(1+-1*coth(a+b*x))*sin(1)+-3//8*b^-1*cosintegral(1+coth(a+b*x))*sin(1)+-1//8*b^-1*cos(3)*sinintegral(3+-3*coth(a+b*x))+-1//8*b^-1*cos(3)*sinintegral(3+3*coth(a+b*x))+1//8*b^-1*cosintegral(3+-3*coth(a+b*x))*sin(3)+1//8*b^-1*cosintegral(3+3*coth(a+b*x))*sin(3)+3//8*b^-1*cos(1)*sinintegral(1+-1*coth(a+b*x))+3//8*b^-1*cos(1)*sinintegral(1+coth(a+b*x)))
@test integrate(sin(coth(a+b*x))^2, x) == :(-1//4*b^-1*log(1+-1*coth(a+b*x))+1//4*b^-1*log(1+coth(a+b*x))+-1//4*b^-1*cos(2)*cosintegral(2+2*coth(a+b*x))+-1//4*b^-1*sin(2)*sinintegral(2+2*coth(a+b*x))+1//4*b^-1*cos(2)*cosintegral(2+-2*coth(a+b*x))+1//4*b^-1*sin(2)*sinintegral(2+-2*coth(a+b*x)))
@test integrate(sin(coth(a+b*x)), x) == :((1/2)*b^-1*cos(1)*sinintegral(1+-1*coth(a+b*x))+(1/2)*b^-1*cos(1)*sinintegral(1+coth(a+b*x))+-1//2*b^-1*cosintegral(1+-1*coth(a+b*x))*sin(1)+-1//2*b^-1*cosintegral(1+coth(a+b*x))*sin(1))
@test integrate(cos(coth(a+b*x))^3, x) == :(-3//8*b^-1*cos(1)*cosintegral(1+-1*coth(a+b*x))+-3//8*b^-1*sin(1)*sinintegral(1+-1*coth(a+b*x))+-1//8*b^-1*cos(3)*cosintegral(3+-3*coth(a+b*x))+-1//8*b^-1*sin(3)*sinintegral(3+-3*coth(a+b*x))+1//8*b^-1*cos(3)*cosintegral(3+3*coth(a+b*x))+1//8*b^-1*sin(3)*sinintegral(3+3*coth(a+b*x))+3//8*b^-1*cos(1)*cosintegral(1+coth(a+b*x))+3//8*b^-1*sin(1)*sinintegral(1+coth(a+b*x)))
@test integrate(cos(coth(a+b*x))^2, x) == :(-1//4*b^-1*log(1+-1*coth(a+b*x))+1//4*b^-1*log(1+coth(a+b*x))+-1//4*b^-1*cos(2)*cosintegral(2+-2*coth(a+b*x))+-1//4*b^-1*sin(2)*sinintegral(2+-2*coth(a+b*x))+1//4*b^-1*cos(2)*cosintegral(2+2*coth(a+b*x))+1//4*b^-1*sin(2)*sinintegral(2+2*coth(a+b*x)))
@test integrate(cos(coth(a+b*x)), x) == :((1/2)*b^-1*cos(1)*cosintegral(1+coth(a+b*x))+(1/2)*b^-1*sin(1)*sinintegral(1+coth(a+b*x))+-1//2*b^-1*cos(1)*cosintegral(1+-1*coth(a+b*x))+-1//2*b^-1*sin(1)*sinintegral(1+-1*coth(a+b*x)))

max_stars_repo_path: test/6 Hyperbolic functions/6.4 Hyperbolic cotangent/6.4.2 Hyperbolic cotangent functions.jl


text:
using DSP
using PyPlot

# Generate a signal u of length m = 50, with each entry a random value that is either −1 or +1
u = rand([1, -1], 50)

c = [1, 0.7, -0.3]
y = conv(c, u)

# Plot u and y
plot(u) # Blue Graph
plot(y) # Orange Graph

h = [0.9, -0.5, 0.5, -0.4, 0.3, -0.3, 0.2, -0.1]
z = conv(h , y)

# Plot for z = h * y
plot(z)


text:
@testset "Polar Decomposition" begin

Random.seed!(1103)

# For each variant, decomposing A to U, H, verify:
# U is unitary
# H is hermitian and positive semi-definite
# backward error is small

for T in [Float64, Complex{Float64}]
    @testset "Newton $T" begin
        for n in [1,3,10]
            A = rand(T,n,n)
            U,H = polar(A, alg = :newton)
            @test U'*U ≈ Matrix(I,n,n) atol=1e-7
            @test ishermitian(H)
            for i in eigvals(H)
                @test i >= 0.
            end
            @test A ≈ U*H

            # also check the properties here (once should suffice)
            r = polar(A, alg = :newton)
            @test U == r.U
            @test H == r.H

            r = polar(A, alg = :newton, verbose = true)
            @test U == r.U
            @test H == r.H

            m = n + 2
            A = rand(T,m,n)
            U,H = polar(A, alg = :newton)
            @test U'*U ≈ Matrix(I,n,n) atol=1e-7
            @test ishermitian(H)
            for i in eigvals(H)
                @test i >= 0.
            end
            @test A ≈ U*H
        end
        A = rand(T,5,7)
        @test_throws ArgumentError U,H = polar(A, alg = :newton)
    end
end

for T in [Float64, Complex{Float64}]
    @testset "Halley $T" begin
        for n in [1,3,10]
            A = rand(T, n, n)
            r = polar(A, alg = :halley)
            U = r.U
            H = r.H
            @test U'*U ≈ Matrix(I,n,n) atol=1e-7
            @test ishermitian(H)
            for i in eigvals(H)
                @test i >= 0.
            end
            @test A ≈ U*H
        end
    end

    @testset "QDWH $T" begin
      for m in [1,3,10]
        B = rand(T, m, m)
        U2,H2 = polar(B, alg = :qdwh)
        @test U2'*U2 ≈ Matrix(I,m,m) atol=1e-7
        @test ishermitian(H2)
        for i in eigvals(H2)
            @test i >= 0.
        end
        @test B ≈ U2*H2
      end
    end
end

for T in [Float64, Complex{Float64}]
    @testset "SVD $T" begin
        for n in [1,3,10]
            m = n + rand(0:5)
            A = rand(T, m, n)
            U,H = polar(A, alg = :svd)
            @test U'*U ≈ Matrix(I,n,n) atol=1e-7
            @test ishermitian(H)
            for i in eigvals(H)
                @test i >= 0.
            end
            @test A ≈ U*H
        end
    end
end

for T in [Float64, Complex{Float64}]
    @testset "Newton-Schulz $T" begin
        for n in [1,3,10]
            m = n
            A = rand(T, m, n)
            nA = opnorm(A)
            A = (1.5 / nA) * A
            U,H = polar(A, alg = :schulz)
            @test U'*U ≈ Matrix(I,n,n) atol=1e-7
            @test ishermitian(H)
            for i in eigvals(H)
                @test i >= 0.
            end
            @test A ≈ U*H
        end
    end
end

for T in [Float64, Complex{Float64}]
    @testset "Hybrid $T" begin
        for n in [1,3,10]
            A = rand(T,n,n)
            U,H = polar(A, alg =:hybrid)
            @test U'*U ≈ Matrix(I,n,n) atol=1e-7
            @test ishermitian(H)
            for i in eigvals(H)
                @test i >= 0.
            end
            @test A ≈ U*H atol=1e-7
        end
    end
end

@testset "Float32" begin
    for n in [1,3,10]
        A = Array{Float32}(undef, n, n)
        copyto!(A, rand(n,n))

        r1 = polar(A, alg =:newton)
        r2 = polar(A, alg =:qdwh)
        r3 = polar(A, alg =:halley)
        r4 = polar(A, alg =:svd)
        r5 = polar(A, alg =:hybrid)

        @test r1.U ≈ r2.U
        @test r2.U ≈ r3.U
        @test r3.U ≈ r4.U
        @test r4.U ≈ r5.U
    end
end

end


text:
using CSSUtil
using WebIO
using Test

@testset "hbox" begin
    el1 = node(:div, "Hello world!")
    el2 = node(:div, "Goodbye world!")
    box = hbox(el1, el2)
    @test props(box)[:style]["display"] == "flex"
    @test props(box)[:style]["flex-direction"] == "row"
    @test el1 ∈ children(box)
    @test el2 ∈ children(box)
end

@testset "vbox" begin
    el1 = node(:div, "Hello world!")
    el2 = node(:div, "Goodbye world!")
    box = vbox(el1, el2)
    @test props(box)[:style]["display"] == "flex"
    @test props(box)[:style]["flex-direction"] == "column"
    @test el1 ∈ children(box)
    @test el2 ∈ children(box)
end

@testset "hline" begin
    line = hline()
    @test haskey(props(line)[:style], "borderBottom")
    @test props(line)[:style]["align-self"] == "stretch"
end

@testset "vline" begin
    line = vline()
    @test haskey(props(line)[:style], "borderLeft")
    @test props(line)[:style]["align-self"] == "stretch"
end


text:
module StrangSplitting
export make_time_step_2d!

function make_time_step_2d!(f, vΔt, advect!)
    for (fi, α) in zip(eachcol(f[1]), vΔt[1](f[2]))
        advect![1](fi, α/2)
    end
    f[2][:] = (f[1])'
    for (fi, α) in zip(eachcol(f[2]), vΔt[2](f[1]))
        advect![2](fi, α)
    end
    f[1][:] = (f[2])'
    for (fi, α) in zip(eachcol(f[1]), vΔt[1](f[2]))
        advect![1](fi, α/2)
    end
end

end


text:
using .Samsara, Conception
using UUIDs, Random

""" LinkedGate is a subtype of Conception.AbstractConcept
LinkedGate is a subtype of Conception.AbstractConcept to mark that it's not supposed to become part of hte history: It is not a temporal element, since a transition throught a gate happens instantaneosly.
CTOR arguments:
    - nodeA : node that can be activated, on one side of gate.
    - nodeB : node that can be activated, on one side of gate.
    - condition : node that can be activated; gate keeper.
"""
mutable struct LinkedGate <: Conception.AbstractConcept
    _id
    pre_node::AbstractConcept
    post_node::Union{AbstractConcept, Nothing}
    _gate_is_open::Bool
    _unlocked_by_trait::Union{AbstractConcept, Nothing}
    function LinkedGate(nodeA::T, nodeB::Union{T, Nothing}=nothing;
                        conditional::Union{L, Nothing}=nothing) where
                            {T <: Conception.AbstractConcept, L <: Conception.AbstractConcept}
        new("gate", nodeA, nodeB, false, conditional)
    end
end

Conception.get_id(it::LinkedGate) = string(it)

""" the_node(gate::LinkedGate)
return the node that is activated when the gate is relevant (pre_node)
"""
function the_node(the_gate::LinkedGate)
    the_gate.pre_node
end

""" show(IO, LinkeGate)
Print LinkeGate as { < > , <|> } for {opened, closed}
"""
function Base.show(io::IO, arg::LinkedGate)
    text_state = gate_is_open(arg) ? ' ' : '|'
    print(io, "<"*text_state*">")
end

""" is_open(LinkedGate)
getter for whether gate is open or not """
function gate_is_open(the_gate::LinkedGate)
    isnothing(the_gate._unlocked_by_trait) && return the_gate._gate_is_open

    return is_active(the_gate._unlocked_by_trait)
end
" convenience: Make for readable code "
gate_is_closed(the_gate::LinkedGate) = !gate_is_open(the_gate)

""" gate_set_connected_node!*(gate, node)
Set connected node in a gate WITH NOTHING on other side // requires an uninitialized gate.
ARG:
    - gate : An uninitialized gate of type LinkedGate
    - node : Any node of type Conception.AbstractConcept ~~ struct that support activate! function.
"""
function gate_set_connected_node!(the_gate::LinkedGate, node::Conception.AbstractConcept)
    !isnothing(the_gate.post_node) && throw(ArgumentError("can't reset an allready set gate"))
    the_gate.post_node = node
end
function gate_set_connected_node!(the_gate::LinkedGate, Nothing)
    the_gate.post_node = nothing
end

""" open_gate(the_gate)
Open up the gate, such that activation of gate causes transmission from nodeA to nodeB or vice versa
"""
function open_gate!(the_gate::LinkedGate, new_state=true)
    the_gate._gate_is_open = new_state
end

""" activate!(LinkedGate)   =>  do the gate : try to enter
activate port, i.e. try to enter
(Extends Conception.activate!(..) -- but not TemporalType, since AbstractConcept is higher.
Therefore, I have to explicitly write Conception.activate! in the next sentence)
"""
function Conception.activate!(the_gate::LinkedGate)
    !gate_is_open(the_gate) && return
    # Letting MuExS handle eactivation, and activate!(nothing) just return,
    #   a closed gate does not cause anything // works by default.
    # ~> isnothing(the_gate.post_node) && return # Ikkje naudsynt!

    if is_active(the_gate.pre_node)
        the_node = the_gate.pre_node
        the_other_node = the_gate.post_node
    else
        the_node = the_gate.post_node
        the_other_node = the_gate.pre_node
    end
    #println("activate gate: "*string(the_node)*"→"*string(the_other_node))
    activate!(the_other_node)
end

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
# cudnnConvolutionMode_t
const CUDNN_CONVOLUTION = Cint(0)
const CUDNN_CROSS_CORRELATION = Cint(1)

# cudnnConvolutionFwdPreference_t
const CUDNN_CONVOLUTION_FWD_NO_WORKSPACE = Cint(0)
const CUDNN_CONVOLUTION_FWD_PREFER_FASTEST = Cint(1)
const CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT = Cint(2)

# cudnnConvolutionFwdAlgo_t
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM = Cint(0)
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM = Cint(1)
const CUDNN_CONVOLUTION_FWD_ALGO_GEMM = Cint(2)
const CUDNN_CONVOLUTION_FWD_ALGO_DIRECT = Cint(3)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT = Cint(4)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING = Cint(5)
const CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD = Cint(6)
const CUDNN_CONVOLUTION_FWD_ALGO_WINOGRAD_NONFUSED = Cint(7)
const CUDNN_CONVOLUTION_FWD_ALGO_COUNT = Cint(8)

# cudnnConvolutionBwdFilterPreference_t
const CUDNN_CONVOLUTION_BWD_FILTER_NO_WORKSPACE = Cint(0)
const CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST = Cint(1)
const CUDNN_CONVOLUTION_BWD_FILTER_SPECIFY_WORKSPACE_LIMIT = Cint(2)

# cudnnConvolutionBwdFilterAlgo_t
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0 = Cint(0)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1 = Cint(1)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT = Cint(2)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3 = Cint(3)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD = Cint(4)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_WINOGRAD_NONFUSED = Cint(5)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT_TILING = Cint(6)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_COUNT = Cint(7)

# cudnnConvolutionBwdDataPreference_t
const CUDNN_CONVOLUTION_BWD_DATA_NO_WORKSPACE = Cint(0)
const CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST = Cint(1)
const CUDNN_CONVOLUTION_BWD_DATA_SPECIFY_WORKSPACE_LIMIT = Cint(2)

# cudnnConvolutionBwdDataAlgo_t
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_0 = Cint(0)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_1 = Cint(1)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT = Cint(2)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING = Cint(3)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD = Cint(4)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_WINOGRAD_NONFUSED = Cint(5)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_COUNT = Cint(6)

mutable struct ConvolutionDesc
    ptr::Cptr
    work

    function ConvolutionDesc(::Type{T}, pads, strides, dilations) where T
        ref = Ref{Cptr}()
        @cudnn :cudnnCreateConvolutionDescriptor (Ptr{Cptr},) ref
        desc = new(ref[], nothing)
        finalizer(desc, x -> @cudnn :cudnnDestroyConvolutionDescriptor (Cptr,) x.ptr)

        @assert length(pads) == length(strides) == length(dilations)
        N = length(pads)
        cpads = Cint[pads[i] for i=N:-1:1]
        cstrides = Cint[strides[i] for i=N:-1:1]
        cdilations = Cint[dilations[i] for i=N:-1:1]
        mode = CUDNN_CROSS_CORRELATION
        @cudnn(:cudnnSetConvolutionNdDescriptor,
            (Cptr,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cint},Cint,Cint),
            desc, N, cpads, cstrides, cdilations, mode, datatype(T))
        desc
    end
end

Base.unsafe_convert(::Type{Cptr}, desc::ConvolutionDesc) = desc.ptr

function convolution(w::CuArray{T,N}, x::CuArray{T,N}, pads, strides, dilations) where {T,N}
    @assert size(w,N-1) == size(x,N-1)
    convdesc = ConvolutionDesc(T, pads, strides, dilations)
    wdesc = FilterDesc(w)
    xdesc = TensorDesc(x)

    ydims = ntuple(N-2) do d
        k = (size(w,d)-1) * dilations[d] + 1
        1 + (size(x,d) + 2pads[d] - k) ÷ strides[d]
    end
    y = similar(x, ydims..., size(w,N), size(x,N))
    ydesc = TensorDesc(y)

    h = gethandle()
    ref = Ref{Cint}()
    preference = CUDNN_CONVOLUTION_FWD_PREFER_FASTEST
    @cudnn(:cudnnGetConvolutionForwardAlgorithm,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Csize_t,Ptr{Cint}),
        h, xdesc, wdesc, convdesc, ydesc, preference, 0, ref)
    algo = ref[]

    ref = Ref{Csize_t}()
    @cudnn(:cudnnGetConvolutionForwardWorkspaceSize,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Csize_t}),
        h, xdesc, wdesc, convdesc, ydesc, algo, ref)
    workspace = CuArray{UInt8}(Int(ref[]))

    @cudnn(:cudnnConvolutionForward,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Cptr,Csize_t,Cptr,Cptr,Cptr),
        h, T[1], xdesc, x, wdesc, w, convdesc, algo, workspace, length(workspace), T[0], ydesc, y)

    convdesc.work = (wdesc,xdesc,ydesc)
    y, convdesc
end

function ∇convolution_bias!(dy::CuArray{T}, db::CuArray{T}) where T
    wdesc,xdesc,ydesc = convdesc.work
    @cudnn(:cudnnConvolutionBackwardBias,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr),
        h, T[1], ydesc, dy, T[1], bdesc, db)
end

function ∇convolution_filter!(convdesc, x::CuArray{T}, dy::CuArray{T}, dw::CuArray{T}) where T
    wdesc,xdesc,ydesc = convdesc.work
    h = gethandle()
    ref = Ref{Cint}()
    preference = CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST
    @cudnn(:cudnnGetConvolutionBackwardFilterAlgorithm,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Csize_t,Ptr{Cint}),
        h, xdesc, ydesc, convdesc, wdesc, preference, 0, ref)
    algo = ref[]

    ref = Ref{Csize_t}()
    @cudnn(:cudnnGetConvolutionBackwardFilterWorkspaceSize,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Csize_t}),
        h, xdesc, ydesc, convdesc, wdesc, algo, ref)
    workspace = CuArray{UInt8}(Int(ref[]))

    @cudnn(:cudnnConvolutionBackwardFilter,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Cptr,Csize_t,Cptr,Cptr,Cptr),
        h, T[1], xdesc, x, ydesc, dy, convdesc, algo, workspace, length(workspace), T[1], wdesc, dw)
end

function ∇convolution_data!(convdesc, w::CuArray{T}, dy::CuArray{T}, dx::CuArray{T}) where T
    wdesc,xdesc,ydesc = convdesc.work
    h = gethandle()
    ref = Ref{Cint}()
    preference = CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST
    @cudnn(:cudnnGetConvolutionBackwardDataAlgorithm,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Csize_t,Ptr{Cint}),
        h, wdesc, ydesc, convdesc, xdesc, preference, 0, ref)
    algo = ref[]

    ref = Ref{Csize_t}()
    @cudnn(:cudnnGetConvolutionBackwardDataWorkspaceSize,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Ptr{Csize_t}),
        h, wdesc, ydesc, convdesc, xdesc, algo, ref)
    workspace = CuArray{UInt8}(Int(ref[]))

    @cudnn(:cudnnConvolutionBackwardData,
        (Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cptr,Cint,Cptr,Csize_t,Cptr,Cptr,Cptr),
        h, T[1], wdesc, w, ydesc, dy, convdesc, algo, workspace, length(workspace), T[1], xdesc, dx)
end