text:
mutable struct BatchTD <: Optimizer
    α::Float64
    α_mod::Float64
    avg::Bool
    BatchTD(α, α_mod, avg) = new(α, α_mod, avg)
end

function update!(value::VFunction, opt::BatchTD, s_t, s_tp1, r, γ, ρ, terminal)
    α = opt.α*opt.α_mod
    if opt.avg
        α = α/length(ρ)
    end
    δ = r + γ.*(dot.(s_tp1, [value.weights])) - dot.(s_t, [value.weights])
    Δθ = α*sum(ρ.*(δ.*s_t))
    value.weights .+= Δθ
end

mutable struct BatchTDC <: Optimizer
    α::Float64
    β::Float64
    α_mod::Float64
    β_mod::Float64
    avg::Bool
    h::IdDict
    # TDC(α, β) = new(α, β, 1.0, 1.0, true)
    # TDC(α, β, α_mod, β_mod, avg) = new(α, β, α_mod, β_mod, avg)
end

function update!(value::VFunction, opt::BatchTDC, s_t, s_tp1, r, γ, ρ, terminal)
    # α = TD.α*TD.α_mod
    # β = TD.β*TD.β_mod

    h = get!(opt.h, value, zero(value.weight))

    α = opt.α*opt.α_mod
    β = opt.β*opt.β_mod
    if opt.avg
        α = α/length(ρ)
        β = β/length(ρ)
    end

    δ = r + γ.*(dot.(s_tp1, [value.weights])) - dot.(s_t, [value.weights])
    Δθ = α*sum(ρ.*(δ.*s_t - γ.*dot.(s_t, [value.h]).*s_tp1))
    Δh = β*sum(s_t.*(ρ.*δ - dot.(s_t, [value.h])))

    value.weights .+= Δθ
    value.h .+= Δh
end

mutable struct WISBatchTDC <: Optimizer
    α::Float64
    β::Float64
    α_mod::Float64
    β_mod::Float64
    WISTDC(α, β) = new(α, β, 1.0, 1.0)
end

function update!(value::VFunction, opt::WISBatchTDC, s_t, s_tp1, r, γ, ρ, terminal)
    δ = r + γ.*(dot.(s_tp1, [value.weights])) - dot.(s_t, [value.weights])
    Δθ = α*sum(ρ.*(δ.*s_t - γ.*dot.(s_t, [value.h]).*s_tp1))
    Δh = β*sum(s_t.*(ρ.*δ - dot.(s_t, [value.h])))

    value.weights .+= Δθ./sum(ρ)
    value.h .+= Δh./sum(ρ)
end

mutable struct BatchTDC2 <: Optimizer
    α::Float64
    β::Float64
    α_mod::Float64
    β_mod::Float64
    avg::Bool
    TDC2(α, β) = new(α, β, 1.0, 1.0, true)
end

function update!(value::VFunction, opt::BatchTDC2, s_t, s_tp1, r, γ, ρ, terminal)
    # Python code sample = (prev_phi.copy(), phi.copy(), action, reward, state, prev_state, pi[action]/mu[action])
    # td_error = (sample[3] + gamma*(sample[1].dot(weights_tdc_iwer)) - sample[0].dot(weights_tdc_iwer))
    # weights_tdc_iwer = weights_tdc_iwer + alpha*sample[-1]*(td_error*sample[0] - gamma*sample[1]*(sample[0].dot(h_tdc_iwer)))
    # h_tlndc_iwer = h_tdc_iwer + alpha_h*(sample[-1]*td_error - sample[0].dot(h_tdc_iwer))*sample[0]
    α = opt.α*opt.α_mod
    β = opt.β*opt.β_mod
    if opt.avg
        α = α/length(ρ)
        β = β/length(ρ)
    end

    δ = r + γ.*(dot.(s_tp1, [value.weights])) - dot.(s_t, [value.weights])
    Δθ = α*sum(ρ.*(δ.*s_t - γ.*dot.(s_t, [value.h]).*s_tp1))
    Δh = β*sum(s_t.*(ρ.*(δ - dot.(s_t, [value.h]))))

    value.weights .+= Δθ
    value.h .+= Δh
end

mutable struct BatchGTD2 <: Optimizer
    α::Float64
    β::Float64
    α_mod::Float64
    β_mod::Float64
    avg::Bool
    GTD2(α, β) = new(α, β, 1.0, 1.0, true)
    GTD2(α, β, α_mod, β_mod, avg) = new(α, β, α_mod, β_mod, avg)
end

function update!(value::VFunction, opt::BatchGTD2, s_t, s_tp1, r, γ, ρ, terminal)
    # Python code sample = (prev_phi.copy(), phi.copy(), action, reward, state, prev_state, pi[action]/mu[action])
    # td_error = (sample[3] + gamma*(sample[1].dot(weights_gtd2_iwer)) - sample[0].dot(weights_gtd2_iwer))
    # weights_gtd2_iwer = weights_gtd2_iwer + alpha*sample[-1]*(sample[0] - gamma*sample[1])*(sample[0].dot(h_gtd2_iwer))
    # h_gtd2_iwer = h_gtd2_iwer + alpha_h*(sample[-1]*td_error - sample[0].dot(h_gtd2_iwer))*sample[0]
    α = opt.α*opt.α_mod
    β = opt.β*opt.β_mod
    if opt.avg
        α = α/length(ρ)
        β = β/length(ρ)
    end

    δ = r + γ.*(dot.(s_tp1, [value.weights])) - dot.(s_t, [value.weights])
    Δθ = α*sum(ρ.*(s_t - γ.*s_tp1).*(dot.(s_t, [value.h])))
    Δh = β*sum(s_t.*(ρ.*δ - dot.(s_t, [value.h])))

    value.weights .+= Δθ
    value.h .+= Δh
end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
cd(@__DIR__) # changes the directory to the current directory, the default I guess is the HOME
using Pkg; Pkg.activate("."); Pkg.instantiate()

#Pkg.add("OrdinaryDiffEq")
#Pkg.add("ModelingToolkit")
#Pkg.add("DiffEqOperators")
#Pkg.add("DomainSets")
#Pkg.add("Plots")
# neuralpde
#Pkg.add("NeuralPDE")
#Pkg.add("Flux")
#Pkg.add("GalacticOptim")
#Pkg.add("DiffEqFlux")
#Pkg.add("Optim")

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets
# Method of Manufactured Solutions: exact solution
# dot notation, vectorized, see https://julialang.org/blog/2017/01/moredots/
u_exact = (x,t) -> exp.(-t) * cos.(x)

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ cos(x),
        u(t,0) ~ exp(-t),
        u(t,1) ~ exp(-t) * cos(1)]

# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0)]

# PDE system - ModelingToolkit
@named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

using NeuralPDE, Flux, GalacticOptim, Optim, DiffEqFlux
# Neural Network
chain = FastChain(FastDense(2,16, Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))
initθ = Float64.(DiffEqFlux.initial_params(chain))
discretization_nn = PhysicsInformedNN(chain, GridTraining(dx); init_params = initθ)

prob_nn = discretize(pdesys, discretization_nn)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

# optimizer
opt = Optim.BFGS()
res = GalacticOptim.solve(prob_nn,opt; cb = cb, maxiters=1200)
phi = discretization_nn.phi

# Method of lines discretization
dx = 0.1

discretization = MOLFiniteDifference([x=>dx],t)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
using OrdinaryDiffEq
sol = solve(prob,Tsit5(),saveat=0.2)

# Plot results and compare with exact solution
x = (0:dx:1)[2:end-1]
t = sol.t

using Plots
plt = plot()

for i in 1:length(t)
    plot!(x,sol.u[i],label="Numerical, t=$(t[i])")
    plot!(x,res.u[i],label="NN, t=$(t[i])")
    scatter!(x, u_exact(x, t[i]),label="Exact, t=$(t[i])")

end
display(plt)
savefig("plot.png")

ext: jl


text:
"""
**Jackknife** errors for (non-linear) functions of uncertain data, i.e.
`g(<a>, <b>, ...)` where `<a>`, `<b>` are the means of data sets `a` and `b`.
Use `jackknife_full(g, a, b, ...)` to get the jackknife estimate, bias and
error or `jackknife(g, a, b, ...)` to get just the estimate and error.

See: [`jackknife_full`](@ref), [`jackknife`](@ref)
"""
module Jackknife

import Statistics: mean, var

export jackknife_full, jackknife

"""
    jackknife_full(g::Function, a[, b, ...])

Returns the jackknife estimate, bias and error for a given function
`g(<a>, <b>, ...)` acting on the means of the samples `a`, `b`, etc.

Example:
    # Assuming some sample xs is given
    # variance of sample xs
    g(x2, x) = x2 - x^2
    estimate, bias, error = jackknife_full(g, xs.^2, xs)

See also: [`estimate`](@ref), [`bias`](@ref), [`std_error`](@ref)
"""
function jackknife_full(g::Function, samples::AbstractVector{<:Number}...)
    reduced_results = leaveoneout(g, samples...)
    return estimate(g, samples...; reduced_results = reduced_results),
           bias(g, samples...; reduced_results = reduced_results),
           _std_error(reduced_results)
end

"""
    jackknife(g::Function, a[, b, ...])

Returns the jackknife estimate and error for a given function `g(<a>, <b>, ...)`
acting on the means of the samples `a`, `b`, etc.

Example:
    # Assuming some sample xs is given
    # variance of sample xs
    g(x2, x) = x2 - x^2
    error = jackknife(g, xs.^2, xs)

See also: [`estimate`](@ref), [`std_error`](@ref)
"""
function jackknife(g::Function, samples::AbstractVector{<:Number}...)
    reduced_results = leaveoneout(g, samples...)
    return estimate(g, samples...; reduced_results = reduced_results),
           _std_error(reduced_results)
end

"""
    leaveoneout(g::Function, samples::AbstractVector{<:Number}...)

Generates N sub-samples for each sample in `samples`, where one value is left
out and applies it mean to the function `g`. The result is used for a
1-jackknife.
"""
function leaveoneout(g::Function, samples::AbstractVector{<:Number}...)
    @assert(
        length(samples[1]) > 1,
        "The sample must have multiple values! (length(samples[1]) > 1)"
    )

    full_sample_sums = map(sum, samples)
    invN1 = 1 / (length(samples[1]) - 1)
    red_sample_means = map(samples, full_sample_sums) do sample, full_sample_sum
        [(full_sample_sum - sample[i]) * invN1 for i in eachindex(sample)]
    end
    g.(red_sample_means...)
end

"""
    var(g::Function, a[, b, ...])

Returns the jackknife variance for a given function `g(<a>, <b>, ...)` acting on
the means of the samples `a`, `b`, etc.
"""
function var(g::Function, samples::AbstractVector{<:Number}...)
    _var(leaveoneout(g, samples...))
end
_var(gis::AbstractVector{<:Complex}) = _var(real(gis)) + _var(imag(gis))
function _var(reduced_results::AbstractVector{<:Real})
    # With N values each sample, N subsamples are created and N reduced_results
    # are calculated.
    N = length(reduced_results)
    # Eq. (3.35) in QMC Methods book
    return var(reduced_results) * (N - 1)^2 / N
end

"""
    std_error(g::Function, a[, b, ...])

Returns the jackknife error for a given function `g(<a>, <b>, ...)` acting on
the means of the samples `a`, `b`, etc.
"""
function std_error(g::Function, samples::AbstractVector{<:Number}...)
    _std_error(leaveoneout(g, samples...))
end
function _std_error(gis::AbstractVector{<:Complex})
    sqrt(_std_error(real(gis))^2 + _std_error(imag(gis))^2)
end
_std_error(gis::AbstractVector{<:Real}) = sqrt(_var(gis))

"""
    bias(g::Function, a[, b, ...])

Returns the jackknife bias for a given function `g(<a>, <b>, ...)` acting on
the means of the samples `a`, `b`, etc.
"""
function bias(
        g::Function,
        samples::AbstractVector{<:Number}...;
        reduced_results::AbstractVector{<:Number} = leaveoneout(g, samples...)
    )
    # Basically Eq. (3.33)
    (length(reduced_results) - 1) * (
        mean(reduced_results) - g(map(mean, samples)...)
    )
end

"""
    estimate(g::Function, a[, b, ...])

Returns the (bias corrected) jackknife estimate for a given function
`g(<a>, <b>, ...)` acting on the means of the samples `a`, `b`, etc. If
`leaveoneout` has already been calculated it can be supplied via the keyword
argument `reduced_results`.
"""
function estimate(
        g::Function,
        samples::AbstractVector{<:Number}...;
        reduced_results::AbstractVector{<:Number} = leaveoneout(g, samples...)
    )
    # Eq. (3.34) in QMC Methods book
    n = length(reduced_results)
    return n * g(map(mean, samples)...) - (n - 1) * mean(reduced_results)
end

end # module

# TODO: Prebinning + Jackknife


text:
module Vimes

using Dates
using MacroTools
using Printf
using SourceWalk
using StatsBase

using SourceWalk.CSTParser
using SourceWalk: sourcemap

include("library.jl")
include("patch.jl")
include("run.jl")

end # module


text:
ComputeEnvironment = Fermi.Environments.Environment(
                                       Fermi.Environments.No_IC(),
                                       true,
                                       Fermi.Environments.NoAccelerator(),
                                       Fermi.Environments.NoCommunicator(),
                                       Fermi.Environments.TBLIS()
                                      )


text:
immutable KmerSpace
    alphabet
    dict
    ks
    function KmerSpace( alphabet, ks )
        dict = [ alphabet[i]=>i-1 for i in 1:Base.length(alphabet) ]
        new( alphabet, dict, ks )
    end
end

nkmers( k, alphabet ) = Base.length( alphabet ) .^ k

Base.length( ks::KmerSpace ) = sum( nkmers( ks.ks, ks.alphabet ) )

function make_kmer( k, alphabet, idx )
    n = Base.length( alphabet )
    kmer = Array( eltype(alphabet), k )
    #use 1 as starting index for consistency
    idx -= 1
    for i = k:-1:1
        offset = idx % n
        kmer[i] = alphabet[ offset + 1 ]
        idx = floor( idx / n )
    end
    return kmer
end

function kmer_index( string, dict )
    idx = 0
    n = Base.length( dict )
    for i = 1:Base.length(string)
        idx *= n
        idx += dict[ string[i] ]
    end
    return idx + 1
end

type KmerSpaceState
    i
    ki
    ke
end

Base.start( ks::KmerSpace ) = KmerSpaceState( 1, 1, nkmers( ks.ks[1], ks.alphabet ) )
function Base.next( ks::KmerSpace, s::KmerSpaceState )
    i = s.i
    ki = s.ki
    ke = s.ke
    if i > ke
        i = 1
        ki += 1
        ke = nkmers( ks.ks[ki], ks.alphabet )
    end
    kmer = make_kmer( ks.ks[ki], ks.alphabet, i )
    next_state = KmerSpaceState( i+1, ki, ke )
    return kmer, next_state
end
Base.done( ks::KmerSpace, s::KmerSpaceState ) = s.i > s.ke && s.ki >= Base.length( ks.ks )

function nfeatures( len::Int, kmer_space )
    n = 0
    for k in kmer_space.ks
        n += ( len - k + 1 ) * nkmers( k, kmer_space.alphabet )
    end
    return n
end

nfeatures( string, kmer_space ) = nfeatures( Base.length(string), kmer_space )

function featurize( string, kmer_space )
    fs = falses( nfeatures( string, kmer_space ) )
    offset = 0
    for k = kmer_space.ks
        n = nkmers( k, kmer_space.alphabet )
        for p = 1:Base.length(string)-k+1
            kmer = string[p:p+k-1]
            i = kmer_index( kmer, kmer_space.dict )
            fs[ offset + (p-1)*n + i ] = true
        end
        offset += n*(Base.length(string)-k+1)
    end
    return fs
end

function feature_names( len, kmer_space )
    names = Array( ASCIIString, nfeatures( len, kmer_space ) )
    offset = 0
    for k = kmer_space.ks
        n = nkmers( k, kmer_space.alphabet )
        for p = 1:len-k+1
            for i = 1:n
                kmer = join( make_kmer( k, kmer_space.alphabet, i ) )
                names[ offset + (p-1)*n + i ] = string( "[",p,"]",kmer )
            end
        end
        offset += n*(len-k+1)
    end
    return names
end

function featurize_all( xs, kmer_space )
    m = size( xs, 1 )
    if m > 0
        n = nfeatures( xs[1], kmer_space )
        fs = Array( Bool, m, n )
        @parallel for i = 1:m
            fs[i,:] = featurize( xs[i], kmer_space )
        end
        return fs, feature_names( Base.length( xs[1] ), kmer_space )
    end
    return Array( bool, 0, 0 ), Array( ASCIIString, 0 )
end

max_stars_repo_name: tbepler/ets-analysis-scripts

max_issues_repo_name: tbepler/ets-analysis-scripts

max_forks_repo_name: tbepler/ets-analysis-scripts


text:
#
# Functions and types to interpolate ephemeris data with
# a Cubic Hermite interpolator.
#
# This functionality is ~entirely~ enabled by Interpolations.jl!
#

export ContinuousEphemeris

"""
An abstract supertype for all possible interpolations for ephemeris data.
"""
abstract type AbstractEphemerisInterpolation end

"""
A callable structure which holds ephemeris data for a body's positions and velocities between the timepoints specified by `T`.
Call the object with a Julian day input to interpolate the body's positions and velocities using `Interpolations.CubicHermite`.
"""
struct ContinuousEphemeris{D<:Number, TX<:Interpolations.CubicHermite, TY<:Interpolations.CubicHermite, TZ<:Interpolations.CubicHermite} <:AbstractEphemerisInterpolation
    T::Tuple{D,D}
    X::TX
    Y::TY
    Z::TZ
end

"""
Construct a `ContinuousEphemeris` instance from an appropriately formatted
`DataFrame`. This is mostly an internal package method. See `ephemeris(args...; continuous = Val{true}, kwargs...)`
for a more user friendly constructor!
"""
function ContinuousEphemeris(data::DataFrame)

    T = eltype(data.dⱼ)
    timespan = (
        data.dⱼ[1],
        data.dⱼ[end],
    )

    X = Interpolations.CubicHermite(data.dⱼ, data.x, data.ẋ)
    Y = Interpolations.CubicHermite(data.dⱼ, data.y, data.ẏ)
    Z = Interpolations.CubicHermite(data.dⱼ, data.z, data.ż)

    return ContinuousEphemeris{T, typeof(X), typeof(Y), typeof(Z)}(
        timespan, X, Y, Z
    )
end

"""
Interpolate the body's position and velocity by a Julian day input.
"""
function (ephemeris::ContinuousEphemeris)(t::Number; bias = false, type = nothing)
    if bias
        time = t + ephemeris.T[1]
    else
        time = t
    end
    if isnothing(type)
        return vcat(
            ephemeris.X(time),
            ephemeris.Y(time),
            ephemeris.Z(time),
            Interpolations.gradient(ephemeris.X, time),
            Interpolations.gradient(ephemeris.Y, time),
            Interpolations.gradient(ephemeris.Z, time),
        )
    else
        return (
            ephemeris.X(time),
            ephemeris.Y(time),
            ephemeris.Z(time),
            Interpolations.gradient(ephemeris.X, time),
            Interpolations.gradient(ephemeris.Y, time),
            Interpolations.gradient(ephemeris.Z, time),
        ) |> type
    end
end

function Base.show(io::IO, ephemeris::ContinuousEphemeris)
    print(io, "Cubic Hermite interpolation of ephemeris data")
end


text:
"""
    ARDTransform(v::AbstractVector)
    ARDTransform(s::Real, dims::Int)

Multiply every vector of observation by `v` element-wise
```
    v = rand(3)
    tr = ARDTransform(v)
```
"""
struct ARDTransform{Tv<:AbstractVector{<:Real}} <: Transform
    v::Tv
end

ARDTransform(s::Real, dims::Integer) = ARDTransform(fill(s, dims))

@functor ARDTransform

function set!(t::ARDTransform{<:AbstractVector{T}}, ρ::AbstractVector{T}) where {T<:Real}
    @assert length(ρ) == dim(t) "Trying to set a vector of size $(length(ρ)) to ARDTransform of dimension $(dim(t))"
    t.v .= ρ
end

dim(t::ARDTransform) = length(t.v)

(t::ARDTransform)(x::Real) = first(t.v) * x
(t::ARDTransform)(x) = t.v .* x

_map(t::ARDTransform, x::AbstractVector{<:Real}) = t.v' .* x
_map(t::ARDTransform, x::ColVecs) = ColVecs(t.v .* x.X)
_map(t::ARDTransform, x::RowVecs) = RowVecs(t.v' .* x.X)

Base.isequal(t::ARDTransform, t2::ARDTransform) = isequal(t.v, t2.v)

Base.show(io::IO, t::ARDTransform) =
    print(io, "ARD Transform (dims: ", dim(t),")")


text:
"""
discrete Lehmann representation for imaginary-time/Matsubara-freqeuncy correlator
"""
# include("spectral.jl")
# using ..Spectral
# include("./functional/builder.jl")
# include("./discrete/builder.jl")
# include("operation.jl")

"""
struct DLRGrid

    DLR grids for imaginary-time/Matsubara frequency correlators

#Members:
- `isFermi`: bool is fermionic or bosonic
- `symmetry`: particle-hole symmetric :ph, or particle-hole asymmetric :pha, or :none
- `Euv` : the UV energy scale of the spectral density
- `β` : inverse temeprature
- `Λ`: cutoff = UV Energy scale of the spectral density * inverse temperature
- `rtol`: tolerance absolute error
- `size` : number of DLR basis
- `ω` : selected representative real-frequency grid
- `n` : selected representative Matsubara-frequency grid (integer)
- `ωn` : (2n+1)π/β
- `τ` : selected representative imaginary-time grid
"""
struct DLRGrid
    isFermi::Bool
    symmetry::Symbol
    Euv::Float64
    β::Float64
    Λ::Float64
    rtol::Float64

    # dlr grids
    # size::Int # rank of the dlr representation
    ω::Vector{Float64}
    n::Vector{Int} # integers, (2n+1)π/β gives the Matsubara frequency
    ωn::Vector{Float64} # (2n+1)π/β
    τ::Vector{Float64}

    """
    function DLRGrid(Euv, β, rtol, isFermi::Bool; symmetry::Symbol = :none, rebuild = false, folder = nothing, algorithm = :functional)

    Create DLR grids

    #Arguments:
    - `Euv` : the UV energy scale of the spectral density
    - `β` : inverse temeprature
    - `isFermi`: bool is fermionic or bosonic
    - `symmetry`: particle-hole symmetric :ph, or particle-hole asymmetric :pha, or :none
    - `rtol`: tolerance absolute error
    - `rebuild` : set false to load DLR basis from the file, set true to recalculate the DLR basis on the fly
    - `folder` : the folder to load the DLR file if rebuild = false, or the folder to save the DLR file if rebuild = true
    - `algorithm` : if rebuild = true, then set :functional to use the functional algorithm to generate the DLR basis, or set :discrete to use the matrix algorithm.
    """
    function DLRGrid(Euv, β, rtol, isFermi::Bool, symmetry::Symbol = :none; rebuild = false, folder = nothing, algorithm = :functional)
        Λ = Euv * β # dlr only depends on this dimensionless scale
        # println("Get $Λ")
        @assert rtol > 0.0 "rtol=$rtol is not positive and nonzero!"
        @assert Λ > 0 "Energy scale $Λ must be positive!"
        @assert symmetry == :ph || symmetry == :pha || symmetry == :none "symmetry must be :ph, :pha or nothing"

        if Λ > 1e8
            printstyled("Warning: current implementation may cause ~ 3-4 digits loss for Λ ≥ 1e8!\n", color = :red)
        end

        if rtol >= 1e-6
            printstyled("Warning: current implementation may cause ~ 3-4 digits loss for rtol ≥ 1e-6!\n", color = :red)
        end

        if Λ < 100
            Λ = Int(100)
        else
            Λ = 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n
        end

        rtolpower = Int(floor(log10(rtol))) # get the biggest n so that rtol>1e-n
        if abs(rtolpower) < 4
            rtolpower = -4
        end

        if symmetry == :none
            # if isFermi
            filename = "universal_$(Λ)_1e$(rtolpower).dlr"
            # else
            #     error("Generic bosonic dlr has not yet been implemented!")
            # end
        elseif symmetry == :ph
            filename = "ph_$(Λ)_1e$(rtolpower).dlr"
        elseif symmetry == :pha
            filename = "pha_$(Λ)_1e$(rtolpower).dlr"
        else
            error("$symmetry is not implemented!")
        end

        dlr = new(isFermi, symmetry, Euv, β, Λ, rtol, [], [], [], [])
        if rebuild
            _build!(dlr, folder, filename, algorithm)
        else
            _load!(dlr, folder, filename)
        end
        return dlr
    end
end

"""
Base.size(dlrGrid::DLRGrid) = length(dlrGrid.ω)
Base.length(dlrGrid::DLRGrid) = length(dlrGrid.ω)
rank(dlrGrid::DLRGrid) = length(dlrGrid.ω)

get the rank of the DLR grid, namely the number of the DLR coefficients.
"""
Base.size(dlrGrid::DLRGrid) = length(dlrGrid.ω)
Base.length(dlrGrid::DLRGrid) = length(dlrGrid.ω)
rank(dlrGrid::DLRGrid) = length(dlrGrid.ω)

function _load!(dlrGrid::DLRGrid, folder, filename)
    searchdir(path, key) = filter(x -> occursin(key, x), readdir(path))
    function finddlr(folder, filename)
        for dir in folder
            if length(searchdir(dir, filename)) > 0
                #dlr file found
                return joinpath(dir, filename)
            end
        end
        error("Cann't find the DLR file $filename in the folders $folder. Please generate it first!")
    end

    folder = isnothing(folder) ? [] : collect(folder)
    push!(folder, string(@__DIR__, "/../basis/"))

    grid = readdlm(finddlr(folder, filename), comments = true, comment_char = '#')
    # println("reading $filename")

    β = dlrGrid.β
    ω, τ = grid[:, 2], grid[:, 3]

    if dlrGrid.isFermi
        n = Int.(grid[:, 4])
        ωn = @. (2n + 1.0) * π / β
    else
        n = Int.(grid[:, 5])
        ωn = @. 2n * π / β
    end
    for r = 1:length(ω)
        push!(dlrGrid.ω, ω[r] / β)
        push!(dlrGrid.τ, τ[r] * β)
        push!(dlrGrid.n, n[r])
        push!(dlrGrid.ωn, ωn[r])
    end
end

function _build!(dlrGrid::DLRGrid, folder, filename, algorithm)
    isFermi = dlrGrid.isFermi
    β = dlrGrid.β
    if algorithm == :discrete || dlrGrid.symmetry == :none
        ω, τ, nF, nB = Discrete.build(dlrGrid, true)
    elseif algorithm == :functional && (dlrGrid.symmetry == :ph || dlrGrid.symmetry == :pha)
        ω, τ, nF, nB = Functional.build(dlrGrid, true)
    else
        error("$algorithm has not yet been implemented!")
    end
    rank = length(ω)
    if isnothing(folder) == false
        open(joinpath(folder, filename), "w") do io
            @printf(io, "# %5s  %25s  %25s  %25s  %20s\n", "index", "freq", "tau", "fermi n", "bose n")
            for r = 1:rank
                @printf(io, "%5i  %32.17g  %32.17g  %16i  %16i\n", r, ω[r], τ[r], nF[r], nB[r])
            end
        end
    end
    for r = 1:rank
        push!(dlrGrid.ω, ω[r] / β)
        push!(dlrGrid.τ, τ[r] * β)
        n = isFermi ? nF[r] : nB[r]
        push!(dlrGrid.n, n)
        push!(dlrGrid.ωn, isFermi ? (2n + 1.0) * π / β : 2n * π / β)
    end
end


text:
###############################################################################
#
#   fmpq.jl : Flint rationals
#
###############################################################################

export fmpq, FlintQQ, FractionField, Rational, FlintRationalField, height,
       height_bits, isless, reconstruct, next_minimal, next_signed_minimal,
       next_calkin_wilf, next_signed_calkin_wilf, dedekind_sum, harmonic,
       bernoulli, bernoulli_cache

###############################################################################
#
#   Data type and parent methods
#
###############################################################################

fmpq(a::Rational{BigInt}) = fmpq(fmpz(a.num), fmpz(a.den))

function fmpq(a::Rational{Int})
  r = fmpq()
  ccall((:fmpq_set_si, :libflint), Void, (Ref{fmpq}, Int64, UInt64), r, numerator(a), denominator(a))
  return r
end

fmpq(a::Rational{T}) where {T <: Integer} = fmpq(numerator(a), denominator(a))

fmpq(a::Integer) = fmpq(fmpz(a), fmpz(1))

fmpq(a::Integer, b::Integer) = fmpq(fmpz(a), fmpz(b))

fmpq(a::fmpz, b::Integer) = fmpq(a, fmpz(b))

fmpq(a::Integer, b::fmpz) = fmpq(fmpz(a), b)

parent(a::fmpq) = FlintQQ

parent_type(::Type{fmpq}) = FlintRationalField

elem_type(::Type{FlintRationalField}) = fmpq

base_ring(a::FlintRationalField) = FlintZZ

base_ring(a::fmpq) = FlintZZ

isdomain_type(::Type{fmpq}) = true

###############################################################################
#
#   Hashing
#
###############################################################################

function Base.hash(a::fmpq, h::UInt)
   return _hash_integer(a.num, _hash_integer(a.den, h))
end

###############################################################################
#
#   Constructors
#
###############################################################################

function //(x::fmpz, y::fmpz)
   iszero(y) && throw(DivideError())
   g = gcd(x, y)
   return fmpq(divexact(x, g), divexact(y, g))
end

//(x::fmpz, y::Integer) = x//fmpz(y)

//(x::Integer, y::fmpz) = fmpz(x)//y

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function numerator(a::fmpq)
   z = fmpz()
   ccall((:fmpq_numerator, :libflint), Void, (Ref{fmpz}, Ref{fmpq}), z, a)
   return z
end

function denominator(a::fmpq)
   z = fmpz()
   ccall((:fmpq_denominator, :libflint), Void, (Ref{fmpz}, Ref{fmpq}), z, a)
   return z
end

doc"""
    abs(a::fmpq)
> Return the absolute value of $a$.
"""
function abs(a::fmpq)
   z = fmpq()
   ccall((:fmpq_abs, :libflint), Void, (Ref{fmpq}, Ref{fmpq}), z, a)
   return z
end

zero(a::FlintRationalField) = fmpq(0)

one(a::FlintRationalField) = fmpq(1)

function isone(a::fmpq)
   return Bool(ccall((:fmpq_is_one, :libflint), Cint, (Ref{fmpq}, ), a))
end

function iszero(a::fmpq)
   return Bool(ccall((:fmpq_is_zero, :libflint), Cint, (Ref{fmpq}, ), a))
end

isunit(a::fmpq) = !iszero(a)

doc"""
    height(a::fmpq)
> Return the height of the fraction $a$, namely the largest of the absolute
> values of the numerator and denominator.
"""
function height(a::fmpq)
   temp = fmpz()
   ccall((:fmpq_height, :libflint), Void, (Ref{fmpz}, Ref{fmpq}), temp, a)
   return temp
end

doc"""
    height_bits(a::fmpq)
> Return the number of bits of the height of the fraction $a$.
"""
function height_bits(a::fmpq)
   return ccall((:fmpq_height_bits, :libflint), Int, (Ref{fmpq},), a)
end

function deepcopy_internal(a::fmpq, dict::ObjectIdDict)
   z = fmpq()
   ccall((:fmpq_set, :libflint), Void, (Ref{fmpq}, Ref{fmpq}), z, a)
   return z
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(a::fmpq) = a

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FlintRationalField)
   print(io, "Rational Field")
end

function show(io::IO, a::fmpq)
   print(io, numerator(a))
   if denominator(a) != 1
      print(io, "//", denominator(a))
   end
end

needs_parentheses(x::fmpq) = false

isnegative(x::fmpq) = x < 0

show_minus_one(::Type{fmpq}) = false

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(a::fmpq)
   z = fmpq()
   ccall((:fmpq_neg, :libflint), Void, (Ref{fmpq}, Ref{fmpq}), z, a)
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fmpq, b::fmpq)
   z = fmpq()
   ccall((:fmpq_add, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
   return z
end

function -(a::fmpq, b::fmpq)
   z = fmpq()
   ccall((:fmpq_sub, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
   return z
end

function *(a::fmpq, b::fmpq)
   z = fmpq()
   ccall((:fmpq_mul, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function +(a::fmpq, b::Int)
   z = fmpq()
   ccall((:fmpq_add_si, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Int), z, a, b)
   return z
end

function +(a::fmpq, b::fmpz)
   z = fmpq()
   ccall((:fmpq_add_fmpz, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
   return z
end

+(a::Int, b::fmpq) = b + a

+(a::fmpz, b::fmpq) = b + a

+(a::fmpq, b::Rational{T}) where {T <: Integer} = a + fmpq(b)

+(a::Rational{T}, b::fmpq) where {T <: Integer} = b + a

function -(a::fmpq, b::Int)
   z = fmpq()
   ccall((:fmpq_sub_si, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Int), z, a, b)
   return z
end

function -(a::fmpq, b::fmpz)
   z = fmpq()
   ccall((:fmpq_sub_fmpz, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
   return z
end

-(a::fmpq, b::Rational{T}) where {T <: Integer} = a - fmpq(b)

-(a::Rational{T}, b::fmpq) where {T <: Integer} = fmpq(b) - a

function *(a::fmpq, b::fmpz)
   z = fmpq()
   ccall((:fmpq_mul_fmpz, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
   return z
end

*(a::fmpz, b::fmpq) = b*a

function -(a::fmpz, b::fmpq)
   n = a*denominator(b) - numerator(b)
   d = denominator(b)
   g = gcd(n, d)
   return parent(b)(divexact(n, g), divexact(d, g))
end

*(a::fmpq, b::Rational{T}) where {T <: Integer} = a * fmpq(b)

*(a::Rational{T}, b::fmpq) where {T <: Integer} = b * a

//(a::fmpq, b::Rational{T}) where {T <: Integer} = a//fmpq(b)

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::fmpq, b::fmpq)
   return ccall((:fmpq_equal, :libflint), Bool,
                (Ref{fmpq}, Ref{fmpq}), a, b)
end

doc"""
    isless(a::fmpq, b::fmpq)
> Return `true` if $a < b$, otherwise return `false`.
"""
function isless(a::fmpq, b::fmpq)
   return ccall((:fmpq_cmp, :libflint), Cint,
                (Ref{fmpq}, Ref{fmpq}), a, b) < 0
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

function ==(a::fmpq, b::Int)
   return ccall((:fmpq_equal_si, :libflint), Bool, (Ref{fmpq}, Int), a, b)
end

==(a::Int, b::fmpq) = b == a

function ==(a::fmpq, b::fmpz)
   return ccall((:fmpq_equal_fmpz, :libflint), Bool,
                (Ref{fmpq}, Ref{fmpz}), a, b)
end

==(a::fmpz, b::fmpq) = b == a

==(a::fmpq, b::Rational{T}) where {T <: Integer} = a == fmpq(b)

==(a::Rational{T}, b::fmpq) where {T <: Integer} = b == a

doc"""
    isless(a::fmpq, b::Integer)
> Return `true` if $a < b$, otherwise return `false`.
"""
function isless(a::fmpq, b::Integer)
   z = fmpq(b)
   return ccall((:fmpq_cmp, :libflint), Cint,
                (Ref{fmpq}, Ref{fmpq}), a, z) < 0
end

doc"""
    isless(a::Integer, b::fmpq)
> Return `true` if $a < b$, otherwise return `false`.
"""
function isless(a::Integer, b::fmpq)
   z = fmpq(a)
   return ccall((:fmpq_cmp, :libflint), Cint,
                (Ref{fmpq}, Ref{fmpq}), z, b) < 0
end

doc"""
    isless(a::fmpq, b::fmpz)
> Return `true` if $a < b$, otherwise return `false`.
"""
function isless(a::fmpq, b::fmpz)
   z = fmpq(b)
   return ccall((:fmpq_cmp, :libflint), Cint,
                (Ref{fmpq}, Ref{fmpq}), a, z) < 0
end

doc"""
    isless(a::fmpz, b::fmpq)
> Return `true` if $a < b$, otherwise return `false`.
"""
function isless(a::fmpz, b::fmpq)
   z = fmpq(a)
   return ccall((:fmpq_cmp, :libflint), Cint,
                (Ref{fmpq}, Ref{fmpq}), z, b) < 0
end

isless(a::Rational{T}, b::fmpq) where {T <: Integer} = isless(fmpq(a), b)

isless(a::fmpq, b::Rational{T}) where {T <: Integer} = isless(a, fmpq(b))

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpq, b::Int)
   temp = fmpq()
   ccall((:fmpq_pow_si, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Int), temp, a, b)
   return temp
end

###############################################################################
#
#   Shifting
#
###############################################################################

doc"""
    <<(a::fmpq, b::Int)
> Return $2^b/a$.
"""
function >>(a::fmpq, b::Int)
   z = fmpq()
   ccall((:fmpq_div_2exp, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Int), z, a, b)
   return z
end

doc"""
    <<(a::fmpq, b::Int)
> Return $2^b\times a$.
"""
function <<(a::fmpq, b::Int)
   z = fmpq()
   ccall((:fmpq_mul_2exp, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Int), z, a, b)
   return z
end

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fmpq)
   z = fmpq()
   ccall((:fmpq_inv, :libflint), Void, (Ref{fmpq}, Ref{fmpq}), z, a)
   return z
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(a::fmpq, b::fmpq)
   z = fmpq()
   ccall((:fmpq_div, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
   return z
end

div(a::fmpq, b::fmpq) = divexact(a, b)

function rem(a::fmpq, b::fmpq)
   iszero(b) && throw("Divide by zero in rem")
   return fmpq(0)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(a::fmpq, b::fmpz)
   z = fmpq()
   ccall((:fmpq_div_fmpz, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, a, b)
   return z
end

divexact(a::fmpz, b::fmpq) = inv(b)*a

divexact(a::fmpq, b::Integer) = divexact(a, fmpz(b))

divexact(a::Integer, b::fmpq) = inv(b)*a

divexact(a::fmpq, b::Rational{T}) where {T <: Integer} = divexact(a, fmpq(b))

divexact(a::Rational{T}, b::fmpq) where {T <: Integer} = divexact(fmpq(a), b)

###############################################################################
#
#   Modular arithmetic
#
###############################################################################

doc"""
    mod(a::fmpq, b::fmpz)
> Return $a \pmod{b}$ where $b$ is an integer coprime to the denominator of
> $a$.
"""
function mod(a::fmpq, b::fmpz)
   z = fmpz()
   ccall((:fmpq_mod_fmpz, :libflint), Void,
         (Ref{fmpz}, Ref{fmpq}, Ref{fmpz}), z, a, b)
   return z
end

doc"""
    mod(a::fmpq, b::Integer)
> Return $a \pmod{b}$ where $b$ is an integer coprime to the denominator of
> $a$.
"""
mod(a::fmpq, b::Integer) = mod(a, fmpz(b))

###############################################################################
#
#   GCD
#
###############################################################################

function gcd(a::fmpq, b::fmpq)
   z = fmpq()
   ccall((:fmpq_gcd, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, a, b)
   return z
end

################################################################################
#
#   Ad hoc Remove and valuation
#
################################################################################

remove(a::fmpq, b::Integer) = remove(a, fmpz(b))

valuation(a::fmpq, b::Integer) = valuation(a, fmpz(b))

###############################################################################
#
#   Rational reconstruction
#
###############################################################################

doc"""
    reconstruct(a::fmpz, b::fmpz)
> Attempt to find a rational number $n/d$ such that
> $0 \leq |n| \leq \lfloor\sqrt{m/2}\rfloor$ and
> $0 < d \leq \lfloor\sqrt{m/2}\rfloor$ such that gcd$(n, d) = 1$ and
> $a \equiv nd^{-1} \pmod{m}$. If no solution exists, an exception is thrown.
"""
function reconstruct(a::fmpz, b::fmpz)
   c = fmpq()
   if !Bool(ccall((:fmpq_reconstruct_fmpz, :libflint), Cint,
                  (Ref{fmpq}, Ref{fmpz}, Ref{fmpz}), c, a, b))
      error("Impossible rational reconstruction")
   end
   return c
end

doc"""
    reconstruct(a::fmpz, b::Integer)
> Attempt to find a rational number $n/d$ such that
> $0 \leq |n| \leq \lfloor\sqrt{m/2}\rfloor$ and
> $0 < d \leq \lfloor\sqrt{m/2}\rfloor$ such that gcd$(n, d) = 1$ and
> $a \equiv nd^{-1} \pmod{m}$. If no solution exists, an exception is thrown.
"""
reconstruct(a::fmpz, b::Integer) =  reconstruct(a, fmpz(b))

doc"""
    reconstruct(a::Integer, b::fmpz)
> Attempt to find a rational number $n/d$ such that
> $0 \leq |n| \leq \lfloor\sqrt{m/2}\rfloor$ and
> $0 < d \leq \lfloor\sqrt{m/2}\rfloor$ such that gcd$(n, d) = 1$ and
> $a \equiv nd^{-1} \pmod{m}$. If no solution exists, an exception is thrown.
"""
reconstruct(a::Integer, b::fmpz) =  reconstruct(fmpz(a), b)

doc"""
    reconstruct(a::Integer, b::Integer)
> Attempt to find a rational number $n/d$ such that
> $0 \leq |n| \leq \lfloor\sqrt{m/2}\rfloor$ and
> $0 < d \leq \lfloor\sqrt{m/2}\rfloor$ such that gcd$(n, d) = 1$ and
> $a \equiv nd^{-1} \pmod{m}$. If no solution exists, an exception is thrown.
"""
reconstruct(a::Integer, b::Integer) =  reconstruct(fmpz(a), fmpz(b))

###############################################################################
#
#   Rational enumeration
#
###############################################################################

doc"""
    next_minimal(a::fmpq)
> Given $x$, returns the next rational number in the sequence obtained by
> enumerating all positive denominators $q$, and for each $q$ enumerating
> the numerators $1 \le p < q$ in order and generating both $p/q$ and $q/p$,
> but skipping all gcd$(p,q) \neq 1$. Starting with zero, this generates
> every nonnegative rational number once and only once, with the first
> few entries being $0, 1, 1/2, 2, 1/3, 3, 2/3, 3/2, 1/4, 4, 3/4, 4/3, \ldots$.
> This enumeration produces the rational numbers in order of minimal height.
> It has the disadvantage of being somewhat slower to compute than the
> Calkin-Wilf enumeration. If $x < 0$ we throw a `DomainError()`.
"""
function next_minimal(a::fmpq)
   a < 0 && throw(DomainError())
   c = fmpq()
   ccall((:fmpq_next_minimal, :libflint), Void, (Ref{fmpq}, Ref{fmpq}), c, a)
   return c
end

doc"""
    next_signed_minimal(a::fmpq)
> Given a signed rational number $x$ assumed to be in canonical form,
> returns the next element in the minimal-height sequence generated by
> `next_minimal` but with negative numbers interleaved. The sequence begins
> $0, 1, -1, 1/2, -1/2, 2, -2, 1/3, -1/3, \ldots$. Starting with zero, this
> generates every rational number once and only once, in order of minimal
> height.
"""
function next_signed_minimal(a::fmpq)
   c = fmpq()
   ccall((:fmpq_next_signed_minimal, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}), c, a)
   return c
end

doc"""
    next_calkin_wilf(a::fmpq)
> Given $x$ return the next number in the breadth-first traversal of the
> Calkin-Wilf tree. Starting with zero, this generates every nonnegative
> rational number once and only once, with the first few entries being
> $0, 1, 1/2, 2, 1/3, 3/2, 2/3, 3, 1/4, 4/3, 3/5, 5/2, 2/5, \ldots$.
> Despite the appearance of the initial entries, the Calkin-Wilf enumeration
> does not produce the rational numbers in order of height: some small
> fractions will appear late in the sequence. This order has the advantage of
> being faster to produce than the minimal-height order.
"""
function next_calkin_wilf(a::fmpq)
   a < 0 && throw(DomainError())
   c = fmpq()
   ccall((:fmpq_next_calkin_wilf, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}), c, a)
   return c
end

doc"""
    next_signed_calkin_wilf(a::fmpq)
> Given a signed rational number $x$ returns the next element in the
> Calkin-Wilf sequence with negative numbers interleaved. The sequence begins
> $0, 1, -1, 1/2, -1/2, 2, -2, 1/3, -1/3, \ldots$. Starting with zero, this
> generates every rational number once and only once, but not in order of
> minimal height.
"""
function next_signed_calkin_wilf(a::fmpq)
   c = fmpq()
   ccall((:fmpq_next_signed_calkin_wilf, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}), c, a)
   return c
end

###############################################################################
#
#   Special functions
#
###############################################################################

doc"""
    harmonic(n::Int)
> Computes the harmonic number $H_n = 1 + 1/2 + 1/3 + \cdots + 1/n$.
> Table lookup is used for $H_n$ whose numerator and denominator
> fit in a single limb. For larger $n$, a divide and conquer strategy is used.
"""
function harmonic(n::Int)
   n < 0 && throw(DomainError())
   c = fmpq()
   ccall((:fmpq_harmonic_ui, :libflint), Void, (Ref{fmpq}, Int), c, n)
   return c
end

doc"""
    bernoulli(n::Int)
> Computes the Bernoulli number $B_n$ for nonnegative $n$.
"""
function bernoulli(n::Int)
   n < 0 && throw(DomainError())
   c = fmpq()
   ccall((:bernoulli_fmpq_ui, :libarb), Void, (Ref{fmpq}, Int), c, n)
   return c
end

doc"""
    bernoulli_cache(n::Int)
> Precomputes and caches all the Bernoulli numbers up to $B_n$.
> This is much faster than repeatedly calling `bernoulli(k)`.
> Once cached, subsequent calls to `bernoulli(k)` for any $k \le n$
> will read from the cache, making them virtually free.
"""
function bernoulli_cache(n::Int)
   n = n + 1
   n < 0 && throw(DomainError())
   ccall((:bernoulli_cache_compute, :libarb), Void, (Int,), n)
end

doc"""
    dedekind_sum(h::fmpz, k::fmpz)
"""
function dedekind_sum(h::fmpz, k::fmpz)
   c = fmpq()
   ccall((:fmpq_dedekind_sum, :libflint), Void,
         (Ref{fmpq}, Ref{fmpz}, Ref{fmpz}), c, h, k)
   return c
end

doc"""
    dedekind_sum(h::fmpz, k::Integer)
> Computes the Dedekind sum $s(h,k)$ for arbitrary $h$ and $k$.
"""
dedekind_sum(h::fmpz, k::Integer) = dedekind_sum(h, fmpz(k))

doc"""
    dedekind_sum(h::Integer, k::fmpz)
> Computes the Dedekind sum $s(h,k)$ for arbitrary $h$ and $k$.
"""
dedekind_sum(h::Integer, k::fmpz) = dedekind_sum(fmpz(h), k)

doc"""
    dedekind_sum(h::Integer, k::Integer)
> Computes the Dedekind sum $s(h,k)$ for arbitrary $h$ and $k$.
"""
dedekind_sum(h::Integer, k::Integer) = dedekind_sum(fmpz(h), fmpz(k))

###############################################################################
#
#   Unsafe operators and functions
#
###############################################################################

function zero!(c::fmpq)
   ccall((:fmpq_zero, :libflint), Void,
         (Ref{fmpq},), c)
   return c
end

function mul!(c::fmpq, a::fmpq, b::fmpq)
   ccall((:fmpq_mul, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), c, a, b)
   return c
end

function addeq!(c::fmpq, a::fmpq)
   ccall((:fmpq_add, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), c, c, a)
   return c
end

function add!(c::fmpq, a::fmpq, b::fmpq)
   ccall((:fmpq_add, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), c, a, b)
   return c
end

###############################################################################
#
#   Conversions to/from flint Julia rationals
#
###############################################################################

function Rational(z::fmpq)
   r = Rational{BigInt}(0)
   ccall((:fmpq_get_mpz_frac, :libflint), Void,
         (Ref{BigInt}, Ref{BigInt}, Ref{fmpq}), r.num, r.den, z)
   return r
end

function Rational(z::fmpz)
   return Rational{BigInt}(BigInt(z))
end

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

(a::FlintRationalField)() = fmpq(fmpz(0), fmpz(1))

(a::FlintRationalField)(b::Rational) = fmpq(numerator(b), denominator(b))

(a::FlintRationalField)(b::Integer) = fmpq(b)

(a::FlintRationalField)(b::Int, c::Int) = fmpq(b, c)

(a::FlintRationalField)(b::fmpz) = fmpq(b)

(a::FlintRationalField)(b::Integer, c::Integer) = fmpq(b, c)

(a::FlintRationalField)(b::fmpz, c::Integer) = fmpq(b, c)

(a::FlintRationalField)(b::Integer, c::fmpz) = fmpq(b, c)

(a::FlintRationalField)(b::fmpz, c::fmpz) = fmpq(b, c)

(a::FlintRationalField)(b::fmpq) = b

###############################################################################
#
#   Conversions and promotions
#
###############################################################################

convert(::Type{fmpq}, a::Integer) = fmpq(a)

convert(::Type{fmpq}, a::fmpz) = fmpq(a)

Base.promote_rule(::Type{fmpq}, ::Type{T}) where {T <: Integer} = fmpq

Base.promote_rule(::Type{fmpq}, ::Type{fmpz}) = fmpq

promote_rule(::Type{fmpq}, ::Type{fmpz}) = fmpq

Base.promote_rule(::Type{fmpq}, ::Type{Rational{T}}) where {T <: Integer} = fmpq

convert(::Type{Rational{BigInt}}, a::fmpq) = Rational(a)

###############################################################################
#
#   FractionField constructor
#
###############################################################################

FractionField(base::FlintIntegerRing) = FlintQQ


text:
#=
test_collection:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-09
=#

module TestCollection
    # packages
    using Test

    # external modules
    include("tuple.jl")
    include("list.jl")
    include("set.jl")
    include("dictionary.jl")

    # methods
    function test()
        @testset "Collection" begin
            @testset "Tuple" begin
                @test isempty(Tuple.t) == false
                @test isempty(Tuple.t_n) == false
                @test_throws BoundsError Tuple.get_val(Tuple.t, 4)
                @test_throws ErrorException Tuple.t_n[:d]
                @test_nowarn Tuple.main()
            end
            @testset "List" begin
                @test isempty(List.list_empty) == true
                @test isempty(List.list) == false
                @test_nowarn List.main()
            end
            @testset "Dictionary" begin
                @test_nowarn Dictionary.main()
            end
            @testset "Set" begin
                @test_nowarn SetModule.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestCollection
    TestCollection.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
my_tests = ["test_hmc.jl",
            "test_empmctuner.jl",
            "test_syntax.jl",
            "parsers/test_diff.jl",
            "parsers/unit_tests.jl",
            "test_dists.jl",
            "test_ARS1.jl",
            "test_ARS2.jl"]

println("Running tests:")

for my_test in my_tests
    println("  * $(my_test) *")
    include(my_test)
end


text:
using Quaternions
using LinearAlgebra
using StaticArrays

abstract type AbstractCamera end

"""
    PerspectiveCamera
"""
mutable struct PerspectiveCamera <: AbstractCamera
    near::GLfloat
    far::GLfloat
    fov::GLfloat
    movingSpeed::GLfloat
    headingSpeed::GLfloat
    position::SVector{3,GLfloat}
    rotationMatrix::SMatrix{3,3,GLfloat,9}
    viewMatrix::SMatrix{4,4,GLfloat,16}
    quaternion::Quaternion{GLfloat}
    forward::SVector{3,GLfloat}
    right::SVector{3,GLfloat}
    up::SVector{3,GLfloat}
end
PerspectiveCamera(; near=0.1, far=1000.0, fov=deg2rad(67),
                    movingSpeed=5.0, headingSpeed=50, position=[0,0,0],
                    rotationMatrix=Matrix(I,3,3), viewMatrix=Matrix(I,4,4),
                    quaternion=qrotation([0,1,0], 0), forward=[0,0,-1],
                    right=[1,0,0], up=[0,1,0]) = PerspectiveCamera(near, far, fov, movingSpeed, headingSpeed, position,
                                                                   rotationMatrix, viewMatrix, quaternion,
                                                                   forward, right, up)

"""
    OrthographicCamera
"""
mutable struct OrthographicCamera <: AbstractCamera
    near::GLfloat
    far::GLfloat
    movingSpeed::GLfloat
    headingSpeed::GLfloat
    position::SVector{3,GLfloat}
    rotationMatrix::SMatrix{3,3,GLfloat,9}
    viewMatrix::SMatrix{4,4,GLfloat,16}
    quaternion::Quaternion{GLfloat}
    forward::SVector{3,GLfloat}
    right::SVector{3,GLfloat}
    up::SVector{3,GLfloat}
end
OrthographicCamera(; near=0.1, far=1000.0,
                     movingSpeed=5.0, headingSpeed=50, position=[0,0,0],
                     rotationMatrix=Matrix(I,3,3), viewMatrix=Matrix(I,4,4),
                     quaternion=qrotation([0,1,0], 0), forward=[0,0,-1],
                     right=[1,0,0], up=[0,1,0]) = OrthographicCamera(near, far, movingSpeed, headingSpeed, position,
                                                                     rotationMatrix, viewMatrix, quaternion,
                                                                     forward, right, up)

"""
    rotate!(camera, axis, angle; fwd=[0,0,-1], rgt=[1,0,0], up=[0,1,0])
Incrementally rotate `camera` about `axis` by `angle` in degree.
"""
function rotate!(camera::AbstractCamera, axis, angle::Real; fwd=[0,0,-1], rgt=[1,0,0], up=[0,1,0])
    camera.quaternion = qrotation(Vector(axis), deg2rad(angle)) * camera.quaternion    # incrementally update quaternion
    camera.rotationMatrix = rotationmatrix(camera.quaternion)    # update new rotation matrix
    camera.forward = camera.rotationMatrix * fwd    # update new forward direction vector
    camera.right = camera.rotationMatrix * rgt    # update new right direction vector
    camera.up = camera.rotationMatrix * up    # update new up direction vector
    return camera
end

"""
    get_view_matrix(camera)
Retrieve camera's latest view matrix. Generally speaking, there is no need to directly
set view matrix to some hard-coded values. You should update view matrix via [`update_view_matrix!`](@ref).
"""
get_view_matrix(camera::AbstractCamera) = camera.viewMatrix

"""
    update_view_matrix!(camera::AbstractCamera, move)
Translate `camera` about `right` axis by `move[1]`; `up` axis by `move[2]`; `forward` axis by `move[3]`.
"""
function update_view_matrix!(camera::AbstractCamera, move::AbstractVector)
    camera.position += camera.forward * move[3]
    camera.position += camera.up * move[2]
    camera.position += camera.right * move[1]
    camera.rotationMatrix = rotationmatrix(camera.quaternion)
    homogenousMatrix = vcat([camera.rotationMatrix camera.position], SMatrix{1,4,GLfloat,4}(0,0,0,1))
    # note that, view matrix is for transforming objects to eye space, but our goal is
    # to transform camera(eye space), so we need an inverse here.
    camera.viewMatrix = inv(homogenousMatrix)
end

function get_current_aspectratio(window::GLFW.Window)
    width, height = GLFW.GetFramebufferSize(window)
    return width / height
end

function get_projective_matrix(camera::PerspectiveCamera, aspectRatio::Real)
    near = camera.near
    far = camera.far
    fov = camera.fov
    range = tan(0.5*fov) * near
    Sx = 2.0*near / (range * aspectRatio + range * aspectRatio)
    Sy = near / range
    Sz = -(far + near) / (far - near)
    Pz = -(2.0*far*near) / (far - near)
    return @SMatrix GLfloat[ Sx  0.0  0.0  0.0;
                            0.0   Sy  0.0  0.0;
                            0.0  0.0   Sz   Pz;
                            0.0  0.0 -1.0  0.0]
end
get_projective_matrix(window, camera) = get_projective_matrix(camera, get_current_aspectratio(window))

function get_projective_matrix(window::GLFW.Window, camera::OrthographicCamera)
    width, height = GLFW.GetFramebufferSize(window)
    right = width/2
    left = -width/2
    top = height/2
    bottom = -height/2
    near = camera.near
    far = camera.far
    Sx = 2 / (right - left)
    Sy = 2 / (top - bottom)
    Sz = -2 / (far - near)
    Px = -(right + left) / (right-left)
    Py = -(top + bottom) / (top - bottom)
    Pz = -(far + near) / (far - near)
    return @SMatrix GLfloat[ Sx  0.0  0.0   Px;
                            0.0   Sy  0.0   Py;
                            0.0  0.0   Sz   Pz;
                            0.0  0.0  0.0  1.0]
end

"""
    updatecamera!(window, camera)
Incrementally update camera's position/orientation.
### Keymap
- `A`: slide left
- `D`: slide right
- `W`: move forward
- `S`: move backward
- `Q`: move upward
- `E`: move downward
- `ARROW LEFT`: yaw left
- `ARROW RIGHT`: yaw right
- `ARROW UP`: pitch up
- `ARROW DOWN`: pitch down
- `Z`: roll left
- `C`: roll right
"""
function updatecamera! end

let
    previousTime = time_ns()/1e9
    global function updatecamera!(window::GLFW.Window, camera::AbstractCamera)
        currentTime = time_ns()/1e9
        elapsedTime = currentTime - previousTime
        previousTime = currentTime
        moveFlag = false
        singleFrameMove = zeros(GLfloat, 3)
        GLFW.GetKey(window, GLFW.KEY_A) && (singleFrameMove[1] -= camera.movingSpeed * elapsedTime; moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_D) && (singleFrameMove[1] += camera.movingSpeed * elapsedTime; moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_Q) && (singleFrameMove[2] -= camera.movingSpeed * elapsedTime; moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_E) && (singleFrameMove[2] += camera.movingSpeed * elapsedTime; moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_S) && (singleFrameMove[3] -= camera.movingSpeed * elapsedTime; moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_W) && (singleFrameMove[3] += camera.movingSpeed * elapsedTime; moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_LEFT)  && (rotate!(camera, camera.up, camera.headingSpeed * elapsedTime); moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_RIGHT) && (rotate!(camera, camera.up, -camera.headingSpeed * elapsedTime); moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_UP)    && (rotate!(camera, camera.right, camera.headingSpeed * elapsedTime); moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_DOWN)  && (rotate!(camera, camera.right, -camera.headingSpeed * elapsedTime); moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_Z)     && (rotate!(camera, camera.forward, camera.headingSpeed * elapsedTime); moveFlag=true)
        GLFW.GetKey(window, GLFW.KEY_C)     && (rotate!(camera, camera.forward, -camera.headingSpeed * elapsedTime); moveFlag=true)
        moveFlag && return update_view_matrix!(camera, singleFrameMove)
        return get_view_matrix(camera)
    end
end

function setposition!(camera::AbstractCamera, x::AbstractVector)
     camera.position = x
     camera.viewMatrix = inv(vcat([camera.rotationMatrix camera.position], SMatrix{1,4,GLfloat,4}(0,0,0,1)))
     return camera
end

function setrotation!(camera::AbstractCamera, x::Quaternion; fwd=[0,0,-1], rgt=[1,0,0], up=[0,1,0])
    camera.quaternion = x
    camera.rotationMatrix = rotationmatrix(camera.quaternion)
    camera.forward = camera.rotationMatrix * fwd
    camera.right = camera.rotationMatrix * rgt
    camera.up = camera.rotationMatrix * up
    camera.viewMatrix = inv(vcat([camera.rotationMatrix camera.position], SMatrix{1,4,GLfloat,4}(0,0,0,1)))
    return camera
end

function resetcamera!(camera::AbstractCamera)
    setposition!(camera, [0,0,0])
    setrotation!(camera, qrotation([1,0,0], 0))
    return camera
end

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_path: OpenGL 4 Tutorials/19_fog/camera.jl

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
# Tests of intstats.jl

using MLBase
using Base.Test

# repeach

@test repeach(2:4, 2) == [2, 2, 3, 3, 4, 4]
@test repeach(2:4, 1:3) == [2, 3, 3, 4, 4, 4]

# repeachcol

a = rand(4, 3)

@test repeachcol(a, 2) == a[:, [1,1,2,2,3,3]]
@test repeachcol(a, 1:3) == a[:, [1,2,2,3,3,3]]

# repeachrow

a = rand(3, 4)
@test repeachrow(a, 2) == a[[1,1,2,2,3,3], :]
@test repeachrow(a, 1:3) == a[[1,2,2,3,3,3], :]

# unique_inverse
a = [:a, :a, :b, :c, :b, :a]
ui = MLBase.unique_inverse(a)
@test isa(ui, Tuple{Vector{Symbol}, Vector{Vector{Int}}})
b = Array{Symbol}(mapreduce(length, +, ui[2]))
for (obj, idx) in zip(ui...) b[idx] = obj end
@test a == b


text:
struct BSONObjectId
    data::NTuple{12, UInt8}
end

Base.isless(x::BSONObjectId, y::BSONObjectId) = x.data < y.data

function BSONObjectId(x::AbstractVector{UInt8})
    length(x) == 12 || throw(ArgumentError("ObjectId bytes must be 12 long"))
    BSONObjectId(NTuple{12, UInt8}(x))
end

function BSONObjectId(s::AbstractString)
    length(s) == 24 || throw(ArgumentError("ObjectId hex string must be 24 characters"))
    BSONObjectId(hex2bytes(s))
end

function Base.show(io::IO, x::BSONObjectId)
    print(io, bytes2hex(collect(x.data)))
    nothing
end

time_part_(x::BSONObjectId) = (Int(x.data[1]) << 24) | (Int(x.data[2]) << 16) | (Int(x.data[3]) << 8) | Int(x.data[4])

Dates.DateTime(x::BSONObjectId) = DateTime(Dates.UTM(Dates.UNIXEPOCH + time_part_(x) * 1000))

Base.time(x::BSONObjectId) = Float64(time_part_(x))

mutable struct BSONObjectIdGenerator
    cache_pad1::NTuple{4, UInt128}
    seq::Threads.Atomic{UInt32}
    rnd::NTuple{5, UInt8}
    cache_pad2::NTuple{4, UInt128}

    function BSONObjectIdGenerator()
        seq = Threads.Atomic{UInt32}(rand(UInt32))
        rnd = rand(UInt64)
        new(
            (UInt128(0), UInt128(0), UInt128(0), UInt128(0)),
            seq,
            (
                (rnd >> 32) % UInt8,
                (rnd >> 24) % UInt8,
                (rnd >> 16) % UInt8,
                (rnd >> 8) % UInt8,
                rnd % UInt8,
            ),
            (UInt128(0), UInt128(0), UInt128(0), UInt128(0)),
        )
    end
end

function Base.show(io::IO, x::BSONObjectIdGenerator)
    print(io, "BSONObjectIdGenerator($(x.seq), $(x.rnd))")
end

if Sys.islinux()
    struct LinuxTimespec
        seconds::Clong
        nanoseconds::Cuint
    end
    @inline function seconds_from_epoch_()
        ts = Ref{LinuxTimespec}()
        ccall(:clock_gettime, Cint, (Cint, Ref{LinuxTimespec}), 0, ts)
        x = ts[]
        x.seconds % UInt32
    end
else
    @inline function seconds_from_epoch_()
        tv = Libc.TimeVal()
        tv.sec % UInt32
    end
end

@inline id_from_parts_(t, r, s) = BSONObjectId((
    (t >> 24) % UInt8,
    (t >> 16) % UInt8,
    (t >> 8) % UInt8,
    t % UInt8,
    r[1],
    r[2],
    r[3],
    r[4],
    r[5],
    (s >> 16) % UInt8,
    (s >> 8) % UInt8,
    s  % UInt8,
))

@inline Base.getindex(x::BSONObjectIdGenerator) = id_from_parts_(
    seconds_from_epoch_(),
    x.rnd,
    Threads.atomic_add!(x.seq, UInt32(1))
)

struct BSONObjectIdIterator
    t::UInt32
    r::NTuple{5, UInt8}
    first::UInt32
    past_last::UInt32
end

Base.eltype(::BSONObjectIdIterator) = BSONObjectId

@inline Base.length(x::BSONObjectIdIterator) = (x.past_last - x.first) % Int

@inline function Base.iterate(x::BSONObjectIdIterator, cur::UInt32 = x.first)
    if cur == x.past_last
        nothing
    else
        id_from_parts_(x.t, x.r, cur), cur + UInt32(1)
    end
end

@inline function Base.getindex(x::BSONObjectIdGenerator, range)
    n = length(range) % UInt32
    first = Threads.atomic_add!(x.seq, n)
    BSONObjectIdIterator(seconds_from_epoch_(), x.rnd, first, first + n)
end

const default_object_id_generator = BSONObjectIdGenerator()

@inline BSONObjectId() = default_object_id_generator[]

@inline bson_object_id_range(n::Integer) = default_object_id_generator[1:n]


text:
import LinearAlgebra: conj!

"""
    Initialise SamplePath on time grid t by copying x into each value of the field yy
"""
#initSamplePath(t,x) = Bridge.samplepath(t,  x)
initSamplePath(t,x) = SamplePath(t, [copy(x) for s in t])

function Bridge.sample!(W::SamplePath{Vector{T}}, P::Wiener{Vector{T}}, y1 = W.yy[1]) where {T}
    y = copy(y1)
    copyto!(W.yy[1], y)

    for i = 2:length(W.tt)
        rootdt = sqrt(W.tt[i]-W.tt[i-1])
        for k in eachindex(y)
            y[k] =  y[k] + rootdt*randn(T)
        end
        copyto!(W.yy[i], y)
    end
    #println(W.yy[1])
    W
end

function Bridge.solve(solver::EulerMaruyama!, u, W::SamplePath, P::Bridge.ProcessOrCoefficients)
    N = length(W)

    tt = W.tt
    yy = [copy(u)]

    i = 1
    dw = W.yy[i+1] - W.yy[i]
    t¯ = tt[i]
    dt = tt[i+1] - t¯

    tmp1 = Bridge._b!((i,t¯), u, 0*u, P)
    tmp2 = Bridge.σ(t¯, u, dw, P)

    y = u + tmp1*dt + tmp2

    #dump(y)
    #error("here")

    for i in 2:N-1
        t¯ = tt[i]
        dt = tt[i+1] - t¯
        push!(yy, y)
        if dw isa Number
            dw = W.yy[i+1] - W.yy[i]
        else
            for k in eachindex(dw)
                dw[k] = W.yy[i+1][k] - W.yy[i][k]
            end
        end

        Bridge._b!((i,t¯), y, tmp1, P)
        Bridge.σ!(t¯, y, dw, tmp2, P)

        for k in eachindex(y)
            y[k] = y[k] + tmp1[k]*dt + tmp2[k]
        end
    end
    copyto!(yy[end], Bridge.endpoint(y, P))
    SamplePath(tt, yy)
end

function Bridge.solve!(solver::EulerMaruyama!, Y::SamplePath, u, W::SamplePath, P::Bridge.ProcessOrCoefficients)
    N = length(W)

    tt = W.tt
    yy = Y.yy
    copyto!(yy[1], u)

    i = 1
    dw = W.yy[i+1] - W.yy[i]
    t¯ = tt[i]
    dt = tt[i+1] - t¯

    tmp1 = Bridge._b!((i,t¯), u, 0*u, P)
    tmp2 = Bridge.σ(t¯, u, dw, P)

    y = u + tmp1*dt + tmp2

    for i in 2:N-1
        t¯ = tt[i]
        dt = tt[i+1] - t¯
        copyto!(yy[i], y)
        if dw isa Number
            dw = W.yy[i+1] - W.yy[i]
        else
            for k in eachindex(dw)
                dw[k] = W.yy[i+1][k] - W.yy[i][k]
            end
        end

        Bridge._b!((i,t¯), y, tmp1, P)
        Bridge.σ!(t¯, y, dw, tmp2, P)

        for k in eachindex(y)
            y[k] = y[k] + tmp1[k]*dt + tmp2[k]
        end
    end
    copyto!(yy[end], Bridge.endpoint(y, P))
    SamplePath(tt, yy)
end

struct StratonovichHeun! <: Bridge.SDESolver
end

function Bridge.solve!(solver::StratonovichHeun!, Y::SamplePath, u, W::SamplePath, P::Bridge.ProcessOrCoefficients)
    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    tt = Y.tt
    tt[:] = W.tt
    yy = Y.yy
    y = copy(u)
    ȳ = copy(u)

    tmp1 = copy(y)
    tmp2 = copy(y)
    tmp3 = copy(y)
    tmp4 = copy(y)

    dw = copy(W.yy[1])
    for i in 1:N-1
        t¯ = tt[i]
        dt = tt[i+1] - t¯
        copyto!(yy[i], y)
        if dw isa Number
            dw = W.yy[i+1] - W.yy[i]
        else
            for k in eachindex(dw)
                dw[k] = W.yy[i+1][k] - W.yy[i][k]
            end
        end

        Bridge._b!((i,t¯), y, tmp1, P)
        Bridge.σ!(t¯, y, dw, tmp2, P)

        for k in eachindex(y)
            ȳ[k] = y[k] + tmp1[k]*dt + tmp2[k] # Euler prediction
        end

        Bridge._b!((i + 1,t¯ + dt), ȳ, tmp3, P) # coefficients at ȳ
        #Bridge.σ!(t¯ + dt, ȳ, dw2, tmp4, P)  # original implementation
        Bridge.σ!(t¯ + dt, ȳ, dw, tmp4, P)

        for k in eachindex(y)
            y[k] = y[k] + 0.5*((tmp1[k] + tmp3[k])*dt + tmp2[k] + tmp4[k])
        end
    end
    copyto!(yy[end], Bridge.endpoint(y, P))
    Y
end

function LinearAlgebra.naivesub!(At::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector = b)
    A = At.parent
    n = size(A, 2)
    if !(n == length(b) == length(x))
        throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal"))
    end
    @inbounds for j in n:-1:1
        iszero(A.data[j,j]) && throw(SingularException(j))
        xj = x[j] = A.data[j,j] \ b[j]
        for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
            b[i] -= A.data[j,i] * xj
        end
    end
    x
end

function LinearAlgebra.naivesub!(At::Adjoint{<:Any,<:LowerTriangular}, B::AbstractMatrix, X::AbstractMatrix = B)
    A = At.parent
    n = size(A, 2)
    if !(n == size(B,1) == size(B,2) == size(X,1) == size(X,2))
        throw(DimensionMismatch())
    end
    @inbounds for k in 1:n
        for j in n:-1:1
            iszero(A.data[j,j]) && throw(SingularException(j))
            xjk = X[j,k] = A.data[j,j] \ B[j,k]
            for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                B[i,k] -= A.data[j,i] * xjk
            end
        end
    end
    X
end

function LinearAlgebra.naivesub!(A::LowerTriangular, B::AbstractMatrix, X::AbstractMatrix = B)
    n = size(A,2)
    if !(n == size(B,1) == size(X,1))
        throw(DimensionMismatch())
    end
    if !(size(B,2) == size(X,2))
        throw(DimensionMismatch())
    end

    @inbounds for k in 1:size(B,2)
        for j in 1:n
            iszero(A.data[j,j]) && throw(SingularException(j))
            xjk = X[j,k] = A.data[j,j] \ B[j,k]
            for i in j+1:n
                B[i,k] -= A.data[i,j] * xjk
            end
        end
    end
    X
end

function LinearAlgebra.naivesub!(A::UpperTriangular, B::AbstractMatrix, X::AbstractMatrix = B)
    n = size(A, 2)
    if !(n == size(B,1) == size(X,1))
        throw(DimensionMismatch())
    end
    if !(size(B,2) == size(X,2))
        throw(DimensionMismatch())
    end

    @inbounds for k in 1:size(B, 2)
        for j in n:-1:1
            iszero(A.data[j,j]) && throw(SingularException(j))
            xjk = X[j,k] = A.data[j,j] \ B[j,k]
            for i in j-1:-1:1 # counterintuitively 1:j-1 performs slightly better
                B[i,k] -= A.data[i,j] * xjk
            end
        end
    end
    X
end

function lyapunovpsdbackward_step!(t, dt, Paux,Hend⁺,H⁺)
    B = Matrix(Bridge.B(t - dt/2, Paux))
    ϕ = (I + 1/2*dt*B)\(I - 1/2*dt*B)
    #Ht .= ϕ *(Hend⁺ + 1/2*dt*Bridge.a(t - dt, Paux))* ϕ' + 1/2*dt*Bridge.a(t, Paux)
    H⁺ .= ϕ *(Hend⁺ + 1/2*dt*Bridge.a(t - dt, Paux))* conj!(copy(ϕ)) + 1/2*dt*Bridge.a(t, Paux)
    H⁺
end

"""
Version where B̃ and ã do not depend on try
"""
function lyapunovpsdbackward_step!(t, dt, Paux,Hend⁺,H⁺,B̃, ã)
    ϕ = (I + 1/2*dt*B̃)\(I - 1/2*dt*B̃)
    #Ht .= ϕ *(Hend⁺ + 1/2*dt*Bridge.a(t - dt, Paux))* ϕ' + 1/2*dt*Bridge.a(t, Paux)
    H⁺ .= ϕ *(Hend⁺ + 1/2*dt*ã)* conj!(copy(ϕ)) + 1/2*dt*ã
    H⁺
end

"""
Compute transpose of square matrix of Unc matrices

A = reshape([Unc(1:4), Unc(5:8), Unc(9:12), Unc(13:16)],2,2)
B = copy(A)
A
conj!(B)
"""
function conj!(A::Array{<:Unc,2}) ## FIX, no copy of A wanted
    m, n = size(A)
    B = reshape(copy(A), n, m)
    for i in 1:m
        for j in 1:n
                B[j, i] = A[i, j]'
        end
    end
    B
end

function conj2(A::Array{T,2}) where {T<:Unc}
    At =  Matrix{T}(undef,size(A,2),size(A,1))
    for i in 1:size(A,2)
        for j in 1:size(A,1)
            At[i,j] = A[j,i]'
        end
    end
    At
end

"""
    Forward simulate landmarks process specified by P on grid t.
    Returns driving motion W and landmarks process X
    t: time grid
    x0: starting point
    P: landmarks specification
"""
function landmarksforward(t, x0::State{Pnt}, P) where Pnt
    W = initSamplePath(t,  zeros(Pnt, dimwiener(P)))
    sample!(W, Wiener{Vector{Pnt}}())
    # forward simulate landmarks
    X = initSamplePath(t,x0)
    #println("Solve for forward process:")
    solve!(EulerMaruyama!(), X, x0, W, P)
    W, X
end

# """
#     Forward simulate landmarks process specified by P using Wiener process W.
#     Writes into X
#     x0: starting point
#     P: landmarks specification
#     W: driving Wiener process
#     X: samplepath written into
#
#     Need that W.tt == X.tt
# """
# function landmarksforward!(x0::State{Pnt}, P, W, X) where Pnt
#     N = length(W)
#     N != length(X) && error("X and W differ in length.")
#
#     solve!(EulerMaruyama!(), X, x0, W, P)
# end

tc(t,T) = t.* (2 .-t/T)
extractcomp(v,i) = map(x->x[i], v)

"""
    Adapting Radford Neal's R implementation of Hamiltonian Monte Carlo with
    stepsize ϵ and L steps
"""
function HMC(U, ∇U, ϵ, L, current_q)
    q = current_q
    p = randn(length(q)) # independent standard normal variates
    current_p = p
    # Make a half step for momentum at the beginning
    p = p - ϵ * ∇U(q) / 2
    # Alternate full steps for position and momentum
    for i in 1:L
        # Make a full step for the position
        q = q + ϵ * p
        # Make a full step for the momentum, except at end of trajectory
        if !(i==L)
            p = p - ϵ * ∇U(q)
        end
    end
    # Make a half step for momentum at the end.
    p = p - ϵ * ∇U(q) / 2
    # Negate momentum at end of trajectory to make the proposal symmetric
    p = -p
    # Evaluate potential and kinetic energies at start and end of trajectory
    current_U = U(current_q)
    current_K = sum(current_pˆ2) / 2
    proposed_U = U(q)
    proposed_K = sum(pˆ2) / 2
    # Accept or reject the state at end of trajectory, returning either
    # the position at the end of the trajectory or the initial position
    if rand() < exp(current_U-proposed_U+current_K-proposed_K)
        return q # accept
    else
        return current_q # reject
    end
end

function deepcopyto!(dest::AbstractArray{T1,N}, src::AbstractArray{T2,N}) where {T1,T2,N}
    checkbounds(dest, axes(src)...)
    src′ = Base.unalias(dest, src)
    for I in eachindex(IndexStyle(src′,dest), src′)
        @inbounds dest[I] = deepcopy(src′[I])
    end
    dest
end

struct FlatPrior end
import Distributions.logpdf
logpdf(::FlatPrior, _x) = 0.0


text:
################################################################################
#
#  fmpz_mod_poly.jl : Flint fmpz_mod_poly (polynomials over Z/nZ, large modulus)
#
################################################################################

export FmpzModPolyRing, fmpz_mod_poly, factor

################################################################################
#
#  Type and parent object methods
#
################################################################################

parent(a::fmpz_mod_poly) = a.parent

base_ring(R::FmpzModPolyRing) = R.base_ring

base_ring(a::fmpz_mod_poly) = base_ring(parent(a))

elem_type(::Type{fmpz_mod_poly}) = fmpz_mod_poly

elem_type(::Type{FmpzModPolyRing}) = fmpz_mod_poly

parent_type(::Type{fmpz_mod_poly}) = FmpzModPolyRing

function check_parent(x::fmpz_mod_poly, y::fmpz_mod_poly)
  parent(x) != parent(y) && error("Parents must coincide")
  nothing
end

################################################################################
#
#  Basic manipulation
#
################################################################################

length(x::fmpz_mod_poly) = ccall((:fmpz_mod_poly_length, :libflint), Int,
                               (Ref{fmpz_mod_poly}, ), x)

degree(x::fmpz_mod_poly) = ccall((:fmpz_mod_poly_degree, :libflint), Int,
                               (Ref{fmpz_mod_poly}, ), x)

function coeff(x::fmpz_mod_poly, n::Int)
  n < 0 && throw(DomainError())
  z = fmpz()
  ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz_mod_poly}, Int), z, x, n)
  return base_ring(x)(z)
end

zero(R::FmpzModPolyRing) = R(0)

one(R::FmpzModPolyRing) = R(1)

gen(R::FmpzModPolyRing) = R([fmpz(0), fmpz(1)])

isgen(a::fmpz_mod_poly) = (degree(a) == 1 &&
                              iszero(coeff(a,0)) && isone(coeff(a,1)))

iszero(a::fmpz_mod_poly) = Bool(ccall((:fmpz_mod_poly_is_zero, :libflint), Int32,
                              (Ref{fmpz_mod_poly}, ), a))

var(R::FmpzModPolyRing) = R.S

modulus(a::fmpz_mod_poly) = a.parent.n

modulus(R::FmpzModPolyRing) = R.n

function deepcopy_internal(a::fmpz_mod_poly, dict::ObjectIdDict)
  z = fmpz_mod_poly(modulus(a), a)
  z.parent = a.parent
  return z
end

################################################################################
#
#  AbstractString I/O
#
################################################################################

function show(io::IO, x::fmpz_mod_poly)
   len = length(x)
   S = var(parent(x))

   if len == 0
      print(io, base_ring(x)(0))
   else
      for i = 1:len - 1
         c = coeff(x, len - i)
         bracket = needs_parentheses(c)
         if !iszero(c)
            if i != 1 && !isnegative(c)
               print(io, "+")
            end
            if !isone(c) && (c != -1 || show_minus_one(typeof(c)))
               if bracket
                  print(io, "(")
               end
               show(io, c)
               if bracket
                  print(io, ")")
               end
               print(io, "*")
            end
            if c == -1 && !show_minus_one(typeof(c))
               print(io, "-")
            end
            print(io, string(S))
            if len - i != 1
               print(io, "^")
               print(io, len - i)
            end
         end
      end
      c = coeff(x, 0)
      bracket = needs_parentheses(c)
      if !iszero(c)
         if len != 1 && !isnegative(c)
            print(io, "+")
         end
         if bracket
            print(io, "(")
         end
         show(io, c)
         if bracket
            print(io, ")")
         end
      end
   end
end

function show(io::IO, R::FmpzModPolyRing)
  print(io, "Univariate Polynomial Ring in ")
  print(io, string(var(R)))
  print(io, " over ")
  print(io, base_ring(R))
end

show_minus_one(::Type{fmpz_mod_poly}) = true

################################################################################
#
#  Canonicalization
#
################################################################################

canonical_unit(a::fmpz_mod_poly) = canonical_unit(lead(a))

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::fmpz_mod_poly)
  z = parent(x)()
  ccall((:fmpz_mod_poly_neg, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}), z, x)
  return z
end

################################################################################
#
#   Binary operations
#
################################################################################

function +(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:fmpz_mod_poly_add, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},  Ref{fmpz_mod_poly}),
               z, x, y)
  return z
end

function -(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:fmpz_mod_poly_sub, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
              z, x, y)
  return z
end

function *(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:fmpz_mod_poly_mul, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},  Ref{fmpz_mod_poly}),
              z, x, y)
  return z
end

###############################################################################
#
#  Ad hoc binary operations
#
###############################################################################

function *(x::fmpz_mod_poly, y::fmpz)
  z = parent(x)()
  ccall((:fmpz_mod_poly_scalar_mul_fmpz, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}), z, x, y)
  return z
end

*(x::fmpz, y::fmpz_mod_poly) = y*x

*(x::fmpz_mod_poly, y::Integer) = x*fmpz(y)

*(x::Integer, y::fmpz_mod_poly) = y*x

function *(x::fmpz_mod_poly, y::Generic.Res{fmpz})
  (base_ring(x) != parent(y)) && error("Must have same parent")
  return x*y.data
end

*(x::Generic.Res{fmpz}, y::fmpz_mod_poly) = y*x

function +(x::fmpz_mod_poly, y::Int)
  z = parent(x)()
  ccall((:fmpz_mod_poly_add_si, :libflint), Void,
    (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int), z, x, y)
  return z
end

+(x::Int, y::fmpz_mod_poly) = +(y, x)

function +(x::fmpz_mod_poly, y::fmpz)
  z = parent(x)()
  ccall((:fmpz_mod_poly_add_fmpz, :libflint), Void,
    (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}), z, x, y)
  return z
end

+(x::fmpz, y::fmpz_mod_poly) = y + x

+(x::fmpz_mod_poly, y::Integer) = x + fmpz(y)

+(x::Integer, y::fmpz_mod_poly) = fmpz(y) + x

function +(x::fmpz_mod_poly, y::Generic.Res{fmpz})
  (base_ring(x) != parent(y)) && error("Elements must have same parent")
  return x + y.data
end

+(x::Generic.Res{fmpz}, y::fmpz_mod_poly) = y + x

function -(x::fmpz_mod_poly, y::Int)
  z = parent(x)()
  ccall((:fmpz_mod_poly_sub_si, :libflint), Void,
    (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int), z, x, y)
  return z
end

function -(x::Int, y::fmpz_mod_poly)
  z = parent(y)()
  ccall((:fmpz_mod_poly_si_sub, :libflint), Void,
    (Ref{fmpz_mod_poly}, Int, Ref{fmpz_mod_poly}), z, x, y)
  return z
end

function -(x::fmpz_mod_poly, y::fmpz)
  z = parent(x)()
  ccall((:fmpz_mod_poly_sub_fmpz, :libflint), Void,
    (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}), z, x, y)
  return z
end

function -(x::fmpz, y::fmpz_mod_poly)
  z = parent(y)()
  ccall((:fmpz_mod_poly_fmpz_sub, :libflint), Void,
    (Ref{fmpz_mod_poly}, Ref{fmpz}, Ref{fmpz_mod_poly}), z, x, y)
  return z
end

-(x::fmpz_mod_poly, y::Integer) = x - fmpz(y)

-(x::Integer, y::fmpz_mod_poly) = fmpz(x) - y

function -(x::fmpz_mod_poly, y::Generic.Res{fmpz})
  (base_ring(x) != parent(y)) && error("Elements must have same parent")
  return x - y.data
end

function -(x::Generic.Res{fmpz}, y::fmpz_mod_poly)
   (parent(x) != base_ring(y)) && error("Elements must have same parent")
   return x.data - y
end

################################################################################
#
#  Powering
#
################################################################################

function ^(x::fmpz_mod_poly, y::Int)
  y < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:fmpz_mod_poly_pow, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int), z, x, y)
  return z
end

################################################################################
#
#  Comparison
#
################################################################################

function ==(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x, y)
  return Bool(ccall((:fmpz_mod_poly_equal, :libflint), Int32,
             (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}), x, y))
end

################################################################################
#
#  Ad hoc comparisons
#
################################################################################

function ==(x::fmpz_mod_poly, y::Generic.Res{fmpz})
  base_ring(x) != parent(y) && error("Incompatible base rings in comparison")
  if length(x) > 1
     return false
  elseif length(x) == 1
     u = fmpz()
     ccall((:fmpz_mod_poly_get_coeff_fmpz, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz_mod_poly}, Int), u, x, 0)
     return u == y
  else
    return iszero(y)
  end
end

==(x::Generic.Res{fmpz}, y::fmpz_mod_poly) = y == x

################################################################################
#
#  Truncation
#
################################################################################

function truncate(a::fmpz_mod_poly, n::Int)
  n < 0 && throw(DomainError())

  z = deepcopy(a)

  if length(z) <= n
    return z
  end

  ccall((:fmpz_mod_poly_truncate, :libflint), Void,
          (Ref{fmpz_mod_poly}, Int), z, n)
  return z
end

function mullow(x::fmpz_mod_poly, y::fmpz_mod_poly, n::Int)
  check_parent(x, y)
  n < 0 && throw(DomainError())

  z = parent(x)()
  ccall((:fmpz_mod_poly_mullow, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int),
                z, x, y, n)
  return z
end

###############################################################################
#
#   Reversal
#
###############################################################################

function reverse(x::fmpz_mod_poly, len::Int)
  len < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:fmpz_mod_poly_reverse, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int), z, x, len)
  return z
end

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpz_mod_poly, len::Int)
  len < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:fmpz_mod_poly_shift_left, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int), z, x, len)
  return z
end

function shift_right(x::fmpz_mod_poly, len::Int)
  len < 0 && throw(DomainError())
  z = parent(x)()
  ccall((:fmpz_mod_poly_shift_right, :libflint), Void,
            (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int), z, x, len)
  return z
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x, y)
  iszero(y) && throw(DivideError())
  d = fmpz()
  q = parent(x)()
  r = parent(x)()
  ccall((:fmpz_mod_poly_divrem_f, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
               d, q, r, x, y)
  d != 1 && error("Impossible inverse in divexact")
  return q
end

div(x::fmpz_mod_poly, y::fmpz_mod_poly) = divexact(x,y)

################################################################################
#
#  Ad hoc exact division
#
################################################################################

function divexact(x::fmpz_mod_poly, y::Generic.Res{fmpz})
  base_ring(x) != parent(y) && error("Elements must have same parent")
  iszero(y) && throw(DivideError())
  q = parent(x)()
  ccall((:fmpz_mod_poly_scalar_div_fmpz, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}),
               q, x, y.data)
  return q
end

function divexact(x::fmpz_mod_poly, y::fmpz)
  y == 0 && throw(DivideError())
  q = parent(x)()
  ccall((:fmpz_mod_poly_scalar_div_fmpz, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}),
               q, x, y)
  return q
end

function divexact(x::fmpz_mod_poly, y::Int)
  y == 0 && throw(DivideError())
  q = parent(x)()
  ccall((:fmpz_mod_poly_scalar_div_fmpz, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}),
               q, x, fmpz(y))
  return q
end

################################################################################
#
#  Division with remainder
#
################################################################################

function divrem(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  iszero(y) && throw(DivideError())
  q = parent(x)()
  r = parent(x)()
  d = fmpz()
  ccall((:fmpz_mod_poly_divrem_f, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
          d, q, r, x, y)
  d > 1 && error("Impossible inverse in divrem")
  return q, r
end

################################################################################
#
#  Remainder
#
################################################################################

function rem(x::fmpz_mod_poly, y::fmpz_mod_poly)
  q, r = divrem(x, y)
  return r
end

mod(x::fmpz_mod_poly, y::fmpz_mod_poly) = rem(x, y)

################################################################################
#
#  Removal and valuation
#
################################################################################

function divides(z::fmpz_mod_poly, x::fmpz_mod_poly)
   q, r = divrem(z, x)
   return iszero(r), q
end

################################################################################
#
#  GCD
#
################################################################################

function gcd(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x, y)
  z = parent(x)()
  f = fmpz()
  ccall((:fmpz_mod_poly_gcd_f, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
              f, z, x, y)
  f > 1 && error("Impossible inverse: $(f) divides modulus")
  return z
end

function gcdx(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x, y)
  g = parent(x)()
  s = parent(x)()
  t = parent(x)()
  f = fmpz()
  ccall((:fmpz_mod_poly_xgcd_f, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},
         Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}), f, g, s, t, x, y)
  f > 1 && error("Impossible inverse: $(f) divides modulus")
  return g, s, t
end

function gcdinv(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  length(y) <= 1 && error("Length of second argument must be >= 2")
  g = parent(x)()
  s = parent(x)()
  f = fmpz()
  ccall((:fmpz_mod_poly_gcdinv_f, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
          f, g, s, x, y)
  f > 1 && error("Impossible inverse: $(f) divides modulus")
  return g, s
end

################################################################################
#
#  Modular arithmetic
#
################################################################################

function invmod(x::fmpz_mod_poly, y::fmpz_mod_poly)
  length(y) == 0 && error("Second argument must not be 0")
  check_parent(x, y)
  if length(y) == 1
    return parent(x)(inv(eval(x, coeff(y, 0))))
  end
  z = parent(x)()
  r = ccall((:fmpz_mod_poly_invmod, :libflint), Int32,
            (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
               z, x, y)
  r == 0 ? error("Impossible inverse in invmod") : return z
end

function mulmod(x::fmpz_mod_poly, y::fmpz_mod_poly, z::fmpz_mod_poly)
  check_parent(x, y)
  check_parent(y, z)
  w = parent(x)()
  ccall((:fmpz_mod_poly_mulmod, :libflint), Void,
        (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
        w, x, y, z)
  return w
end

function powmod(x::fmpz_mod_poly, e::Int, y::fmpz_mod_poly)
  check_parent(x, y)
  z = parent(x)()

  if e < 0
    g, x = gcdinv(x, y)
    if g != 1
      error("Element not invertible")
    end
    e = -e
  end

  ccall((:fmpz_mod_poly_powmod_ui_binexp, :libflint), Void,
        (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Int, Ref{fmpz_mod_poly}), z, x, e, y)

  return z
end

doc"""
    powmod(x::fmpz_mod_poly, e::fmpz, y::fmpz_mod_poly)
> Return $x^e \pmod{y}$.
"""
function powmod(x::fmpz_mod_poly, e::fmpz, y::fmpz_mod_poly)
  z = parent(x)()

  if e < 0
    g, x = gcdinv(x, y)
    if g != 1
      error("Element not invertible")
    end
    e = -e
  end

  ccall((:fmpz_mod_poly_powmod_fmpz_binexp, :libflint), Void,
        (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz}, Ref{fmpz_mod_poly}),
            z, x, e, y)
  return z
end

################################################################################
#
#  Resultant
#
################################################################################

function resultant(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  z = parent(x)()
  !isprobabprime(modulus(x)) && error("Modulus not prime in resultant")
  r = fmpz()
  ccall((:fmpz_mod_poly_resultant, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}), r, x, y)
  return base_ring(x)(r)
end

################################################################################
#
#  Evaluation
#
################################################################################

function evaluate(x::fmpz_mod_poly, y::Generic.Res{fmpz})
  base_ring(x) != parent(y) && error("Elements must have same parent")
  z = fmpz()
  ccall((:fmpz_mod_poly_evaluate_fmpz, :libflint), Void,
              (Ref{fmpz}, Ref{fmpz_mod_poly}, Ref{fmpz}), z, x, (y.data))
  return parent(y)(z)
end

################################################################################
#
#  Derivative
#
################################################################################

function derivative(x::fmpz_mod_poly)
  z = parent(x)()
  ccall((:fmpz_mod_poly_derivative, :libflint), Void,
        (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}), z, x)
  return z
end

###############################################################################
#
#   Integral
#
###############################################################################

function integral(x::fmpz_mod_poly)
   len = length(x)
   v = Array{Generic.Res{fmpz}}(len + 1)
   v[1] = zero(base_ring(x))
   for i = 1:len
      v[i + 1] = divexact(coeff(x, i - 1), base_ring(x)(i))
   end
   return parent(x)(v)
end

################################################################################
#
#  Composition
#
################################################################################

function compose(x::fmpz_mod_poly, y::fmpz_mod_poly)
  check_parent(x,y)
  z = parent(x)()
  ccall((:fmpz_mod_poly_compose, :libflint), Void,
          (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}),
               z, x, y)
  return z
end

################################################################################
#
#  Lifting
#
################################################################################

doc"""
    function lift(R::FmpzPolyRing, y::fmpz_mod_poly)
> Lift from a polynomial over $\mathbb{Z}/n\mathbb{Z}$ to a polynomial over
> $\mathbb{Z}$ with minimal reduced nonnegative coefficients. The ring `R`
> specifies the ring to lift into.
"""
function lift(R::FmpzPolyRing, y::fmpz_mod_poly)
   z = fmpz_poly()
   ccall((:fmpz_mod_poly_get_fmpz_poly, :libflint), Void,
          (Ref{fmpz_poly}, Ref{fmpz_mod_poly}), z, y)
   z.parent = R
  return z
end

################################################################################
#
#  Irreducibility
#
################################################################################

doc"""
    isirreducible(x::fmpz_mod_poly)
> Return `true` if $x$ is irreducible, otherwise return `false`.
"""
function isirreducible(x::fmpz_mod_poly)
  !isprobabprime(modulus(x)) && error("Modulus not prime in isirreducible")
  return Bool(ccall((:fmpz_mod_poly_is_irreducible, :libflint), Int32,
          (Ref{fmpz_mod_poly}, ), x))
end

################################################################################
#
#  Squarefree testing
#
################################################################################

doc"""
    issquarefree(x::fmpz_mod_poly)
> Return `true` if $x$ is squarefree, otherwise return `false`.
"""
function issquarefree(x::fmpz_mod_poly)
   !isprobabprime(modulus(x)) && error("Modulus not prime in issquarefree")
   return Bool(ccall((:fmpz_mod_poly_is_squarefree, :libflint), Int32,
      (Ref{fmpz_mod_poly}, ), x))
end

################################################################################
#
#  Factorization
#
################################################################################

doc"""
    factor(x::fmpz_mod_poly)
> Return the factorisation of $x$.
"""
function factor(x::fmpz_mod_poly)
  !isprobabprime(modulus(x)) && error("Modulus not prime in factor")
  fac = _factor(x)
  return Fac(parent(x)(lead(x)), fac)
end

function _factor(x::fmpz_mod_poly)
  fac = fmpz_mod_poly_factor(parent(x).n)
  ccall((:fmpz_mod_poly_factor, :libflint), UInt,
          (Ref{fmpz_mod_poly_factor}, Ref{fmpz_mod_poly}), fac, x)
  res = Dict{fmpz_mod_poly, Int}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:fmpz_mod_poly_factor_get_fmpz_mod_poly, :libflint), Void,
         (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly_factor}, Int), f, fac, i - 1)
    e = unsafe_load(fac.exp, i)
    res[f] = e
  end
  return res
end

doc"""
    factor_squarefree(x::fmpz_mod_poly)
> Return the squarefree factorisation of $x$.
"""
function factor_squarefree(x::fmpz_mod_poly)
  !isprobabprime(modulus(x)) && error("Modulus not prime in factor_squarefree")
  fac = _factor_squarefree(x)
  return Fac(parent(x)(lead(x)), fac)
end

function _factor_squarefree(x::fmpz_mod_poly)
  fac = fmpz_mod_poly_factor(parent(x).n)
  ccall((:fmpz_mod_poly_factor_squarefree, :libflint), UInt,
          (Ref{fmpz_mod_poly_factor}, Ref{fmpz_mod_poly}), fac, x)
  res = Dict{fmpz_mod_poly, Int}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:fmpz_mod_poly_factor_get_fmpz_mod_poly, :libflint), Void,
         (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly_factor}, Int), f, fac, i - 1)
    e = unsafe_load(fac.exp, i)
    res[f] = e
  end
  return res
end

doc"""
    factor_distinct_deg(x::fmpz_mod_poly)
> Return the distinct degree factorisation of a squarefree polynomial $x$.
"""
function factor_distinct_deg(x::fmpz_mod_poly)
  !issquarefree(x) && error("Polynomial must be squarefree")
  !isprobabprime(modulus(x)) && error("Modulus not prime in factor_distinct_deg")
  degs = Array{Int}(degree(x))
  degss = [ pointer(degs) ]
  fac = fmpz_mod_poly_factor(parent(x).n)
  ccall((:fmpz_mod_poly_factor_distinct_deg, :libflint), UInt,
          (Ref{fmpz_mod_poly_factor}, Ref{fmpz_mod_poly}, Ptr{Ref{Int}}),
          fac, x, degss)
  res = Dict{Int, fmpz_mod_poly}()
  for i in 1:fac.num
    f = parent(x)()
    ccall((:fmpz_mod_poly_factor_get_fmpz_mod_poly, :libflint), Void,
         (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly_factor}, Int), f, fac, i - 1)
    res[degs[i]] = f
  end
  return res
end

################################################################################
#
#  Unsafe functions
#
################################################################################

function zero!(x::fmpz_mod_poly)
  ccall((:fmpz_mod_poly_zero, :libflint), Void,
                   (Ref{fmpz_mod_poly}, ), x)
  return x
end

function fit!(x::fmpz_mod_poly, n::Int)
  ccall((:fmpz_mod_poly_fit_length, :libflint), Void,
                   (Ref{fmpz_mod_poly}, Int), x, n)
  return nothing
end

function setcoeff!(x::fmpz_mod_poly, n::Int, y::UInt)
  ccall((:fmpz_mod_poly_set_coeff_ui, :libflint), Void,
                   (Ref{fmpz_mod_poly}, Int, UInt), x, n, y)
  return x
end

function setcoeff!(x::fmpz_mod_poly, n::Int, y::Int)
  ccall((:fmpz_mod_poly_set_coeff_ui, :libflint), Void,
                   (Ref{fmpz_mod_poly}, Int, UInt), x, n, mod(y, x.n))
  return x
end

function setcoeff!(x::fmpz_mod_poly, n::Int, y::fmpz)
  ccall((:fmpz_mod_poly_set_coeff_fmpz, :libflint), Void,
                   (Ref{fmpz_mod_poly}, Int, Ref{fmpz}), x, n, y)
  return x
end

setcoeff!(x::fmpz_mod_poly, n::Int, y::Integer) = setcoeff!(x, n, fmpz(y))

setcoeff!(x::fmpz_mod_poly, n::Int, y::Generic.Res{fmpz}) = setcoeff!(x, n, y.data)

function add!(z::fmpz_mod_poly, x::fmpz_mod_poly, y::fmpz_mod_poly)
  ccall((:fmpz_mod_poly_add, :libflint), Void,
     (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},  Ref{fmpz_mod_poly}), z, x, y)
  return z
end

function addeq!(z::fmpz_mod_poly, y::fmpz_mod_poly)
  ccall((:fmpz_mod_poly_add, :libflint), Void,
     (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},  Ref{fmpz_mod_poly}), z, z, y)
  return z
end

function sub!(z::fmpz_mod_poly, x::fmpz_mod_poly, y::fmpz_mod_poly)
  ccall((:fmpz_mod_poly_sub, :libflint), Void,
     (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},  Ref{fmpz_mod_poly}), z, x, y)
  return z
end

function mul!(z::fmpz_mod_poly, x::fmpz_mod_poly, y::fmpz_mod_poly)
  ccall((:fmpz_mod_poly_mul, :libflint), Void,
     (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly},  Ref{fmpz_mod_poly}), z, x, y)
  return z
end

################################################################################
#
#  Promotion rules
#
################################################################################

promote_rule(::Type{fmpz_mod_poly}, ::Type{V}) where {V <: Integer} = fmpz_mod_poly

promote_rule(::Type{fmpz_mod_poly}, ::Type{fmpz}) = fmpz_mod_poly

promote_rule(::Type{fmpz_mod_poly}, ::Type{Generic.Res{fmpz}}) = fmpz_mod_poly

###############################################################################
#
#   Polynomial substitution
#
###############################################################################

function (f::fmpz_mod_poly)(a::Generic.Res{fmpz})
   if parent(a) != base_ring(f)
      return subst(f, a)
   end
   return evaluate(f, a)
end

################################################################################
#
#  Parent object call overloads
#
################################################################################

function (R::FmpzModPolyRing)()
  z = fmpz_mod_poly(R.n)
  z.parent = R
  return z
end

function (R::FmpzModPolyRing)(x::fmpz)
  z = fmpz_mod_poly(R.n, x)
  z.parent = R
  return z
end

function (R::FmpzModPolyRing)(x::Integer)
  z = fmpz_mod_poly(R.n, fmpz(x))
  z.parent = R
  return z
end

function (R::FmpzModPolyRing)(x::Generic.Res{fmpz})
  base_ring(R) != parent(x) && error("Wrong parents")
  z = fmpz_mod_poly(R.n, x.data)
  z.parent = R
  return z
end

function (R::FmpzModPolyRing)(arr::Array{fmpz, 1})
  z = fmpz_mod_poly(R.n, arr)
  z.parent = R
  return z
end

function (R::FmpzModPolyRing)(arr::Array{Generic.Res{fmpz}, 1})
  if length(arr) > 0
     (base_ring(R) != parent(arr[1])) && error("Wrong parents")
  end
  z = fmpz_mod_poly(R.n, arr)
  z.parent = R
  return z
end

(R::FmpzModPolyRing)(arr::Array{T, 1}) where {T <: Integer} = R(map(base_ring(R), arr))

function (R::FmpzModPolyRing)(x::fmpz_poly)
  z = fmpz_mod_poly(R.n, x)
  z.parent = R
  return z
end

function (R::FmpzModPolyRing)(f::fmpz_mod_poly)
   parent(f) != R && error("Unable to coerce polynomial")
   return f
end

################################################################################
#
#  Polynomial ring constructor
#
################################################################################

function PolynomialRing(R::Generic.ResRing{fmpz}, s::AbstractString; cached=true)
   parent_obj = FmpzModPolyRing(R, Symbol(s), cached)

   return parent_obj, parent_obj([R(0), R(1)])
end