text:
module CURAND

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, libraryPropertyType, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK
using ..CUDA: libcurand, @retry_reclaim

using CEnum

# core library
include("libcurand_common.jl")
include("error.jl")
include("libcurand.jl")

# low-level wrappers
include("wrappers.jl")

# high-level integrations
include("random.jl")

# thread cache for task-local library handles
const CURAND_THREAD_RNGs = Vector{Union{Nothing,RNG}}()
const GPUARRAY_THREAD_RNGs = Vector{Union{Nothing,GPUArrays.RNG}}()

function default_rng()
    tid = Threads.threadid()
    if @inbounds CURAND_THREAD_RNGs[tid] === nothing
        ctx = context()
        CURAND_THREAD_RNGs[tid] = get!(task_local_storage(), (:CURAND, ctx)) do
            rng = RNG()
            Random.seed!(rng)
            rng
        end
    end
    something(@inbounds CURAND_THREAD_RNGs[tid])
end

function GPUArrays.default_rng(::Type{<:CuArray})
    tid = Threads.threadid()
    if @inbounds GPUARRAY_THREAD_RNGs[tid] === nothing
        ctx = context()
        GPUARRAY_THREAD_RNGs[tid] = get!(task_local_storage(), (:GPUArraysRNG, ctx)) do
            dev = device()
            N = attribute(dev, DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
            state = CuArray{NTuple{4, UInt32}}(undef, N)
            rng = GPUArrays.RNG(state)
            Random.seed!(rng)
            rng
        end
    end
    something(@inbounds GPUARRAY_THREAD_RNGs[tid])
end

function __init__()
    resize!(CURAND_THREAD_RNGs, Threads.nthreads())
    fill!(CURAND_THREAD_RNGs, nothing)

    resize!(GPUARRAY_THREAD_RNGs, Threads.nthreads())
    fill!(GPUARRAY_THREAD_RNGs, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        CURAND_THREAD_RNGs[tid] = nothing
        GPUARRAY_THREAD_RNGs[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        CURAND_THREAD_RNGs[tid] = nothing
        GPUARRAY_THREAD_RNGs[tid] = nothing
    end
end

@deprecate seed!() CUDA.seed!()
@deprecate seed!(seed) CUDA.seed!(seed)

end


text:
module BerryPhase

export inputfiles, interall, eigenvec, berryphase
export Model

abstract Model
include("../util.jl/bond.jl")

function header_interall(io, m::Model)
  println(io,
"""
=== header
NInterAll $(length(m.coeff) * num_bonds(m))
=== reserved
=== reserved
=== end of header """)
end

include("dimer.jl")
include("alt-chain.jl")
include("alt-ladder.jl")

include("mkparams.jl")
include("eigenvec.jl")

function berryphase(model::Model, M::Integer=10)
  generate_defs(model)
  interall(model, 0.0)
  run(pipeline(`../HPhi -e namelist.def`, stdout="std.out", stderr="std.err", append=true))

  bp = 0.0
  vec0 = eigenvec()
  N = length(vec0)

  ref = [complex(float(i)) for i in 1:N]

  #  ref = zeros(Complex128, N)
  #  ref[end] = one(Complex128)

  # ref = rand(Complex128, N)

  normalize!(ref)
  gauge_fix!(vec0, ref)
  vec_old = deepcopy(vec0)
  dt = 2pi/M
  for i in 1:(M-1)
    t = i*dt
    interall(model, t)
    run(pipeline(`../HPhi -e namelist.def`, stdout="std.out", stderr="std.err", append=true))
    vec_new = eigenvec()
    gauge_fix!(vec_new, ref)
    ip = zero(Complex128)
    @inbounds for j in 1:N
      ip += vec_new[j]' * vec_old[j]
      vec_old[j] = vec_new[j]
    end
    bp += angle(ip)
  end
  ip = dot(vec0, vec_old)
  bp += angle(ip)
  return mod(bp, 2pi)
end

end ## of module BerryPhase

ext: jl

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
using .GeneticVariation

export bcf, vcf

"""
    openvcf(::String)
Open VCF file (`.vcf(.gz)`, or `.bcf(.gz)`) and return an `IO` stream in reading mode `"r"`.
"""
function openvcf(infile::String)
    if endswith(infile, ".gz")
        throw(ArgumentError("Please load in GZip.jl with \`using GZip\`"))
    else
        return open(infile, "r")
    end
end

function bcf(infile::String; rename_loci::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false)
    bases = (A = Int8(1), T = Int8(2), C = Int8(3), G = Int8(4), miss = Int8(0))
    stream = BCF.Reader(openvcf(infile))
    nmarkers = countlines(openvcf(infile)) - length(BCF.header(stream)) - 1
    sample_ID = header(stream).sampleID
    nsamples = length(sample_ID)
    loci_names = fill("marker", nmarkers)
    geno_df = DataFrame(:name => sample_ID, :population =>  "missing")
    if silent == false
        @info "\n$(abspath(infile))\n$nsamples samples detected\n$nmarkers markers detected\npopulation info must be added <---"
    end
        for record in stream
        ref_alt = Dict(-1 => "miss", 0 => BCF.ref(record), [i => j for (i,j) in enumerate(BCF.alt(record))]...)
        raw_geno = BCF.genotype(record, 1:nsamples, "GT")
        conv_geno = map(raw_geno) do rg
            tmp = replace.(rg, "." => "-1")
            ig = collect(parse.(Int8, split(tmp, r"\/|\|")))
            [bases[Symbol(ref_alt[i])] for i in ig] |> sort |> Tuple
        end
        insertcols!(geno_df, Symbol(BCF.chrom(record) * "_" * string(BCF.pos(record))) => conv_geno)
    end
    close(stream)
    if rename_snp
        rnm = append!([:name, :population], [Symbol.("snp_" * i) for i in string.(1:nmarkers)])
        rename!(geno_df, rnm)
    end
    stacked_geno_df = DataFrames.stack(geno_df, DataFrames.Not(1:2))
    rename!(stacked_geno_df, [:name, :population, :locus, :genotype])
    # set columns as PooledArrays
    select!(
        stacked_geno_df,
        :name => PooledArray => :name,
        :population => PooledArray => :population,
        :locus => (i -> PooledArray(i |> Vector{String})) => :locus,
        :genotype
    )
    # replace missing genotypes as missing
    stacked_geno_df.genotype = map(stacked_geno_df.genotype) do geno
        if all(0 .== geno)
            return missing
        else
            return geno
        end
    end
    sort!(stacked_geno_df, [:name, :locus])
    #meta_df = generate_meta(stacked_geno_df)
    # ploidy finding
    meta_df = DataFrames.combine(DataFrames.groupby(stacked_geno_df, :name),
        :genotype => (i -> Int8(length(first(skipmissing(i))))) => :ploidy
    )
    insertcols!(meta_df, 2, :population => "missing")
    insertcols!(meta_df, 4, :longitude => Vector{Union{Missing, Float32}}(undef, nsamples), :latitude => Vector{Union{Missing, Float32}}(undef, nsamples))

    pd_out = PopData(meta_df, stacked_geno_df)
    !allow_monomorphic && drop_monomorphic!(pd_out)
    return pd_out
end

### VCF parsing ###

function vcf(infile::String; rename_snp::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false)
    bases = (A = Int8(1), T = Int8(2), C = Int8(3), G = Int8(4), miss = Int8(0))
    stream = VCF.Reader(openvcf(infile))
    nmarkers = countlines(openvcf(infile)) - length(VCF.header(stream)) - 1
    sample_ID = header(stream).sampleID
    nsamples = length(sample_ID)
    loci_names = fill("marker", nmarkers)
    geno_df = DataFrame(:name => sample_ID, :population =>  "missing")
    if silent == false
        @info "\n$(abspath(infile))\n$nsamples samples detected\n$nmarkers markers detected\npopulation info must be added <---"
    end
    for record in stream
        ref_alt = Dict(-1 => "miss", 0 => VCF.ref(record), [i => j for (i,j) in enumerate(VCF.alt(record))]...)
        raw_geno = VCF.genotype(record, 1:nsamples, "GT")
        conv_geno = map(raw_geno) do rg
            tmp = replace.(rg, "." => "-1")
            ig = collect(parse.(Int8, split(tmp, r"\/|\|")))
            [bases[Symbol(ref_alt[i])] for i in ig] |> sort |> Tuple
        end
        insertcols!(geno_df, Symbol(VCF.chrom(record) * "_" * string(VCF.pos(record))) => conv_geno)
    end
    close(stream)
    if rename_snp
        rnm = append!([:name, :population], [Symbol.("snp_" * i) for i in string.(1:nmarkers)])
        rename!(geno_df, rnm)
    end
    stacked_geno_df = DataFrames.stack(geno_df, DataFrames.Not(1:2))
    rename!(stacked_geno_df, [:name, :population, :locus, :genotype])
    # set columns as PooledArrays
    select!(
        stacked_geno_df,
        :name => PooledArray => :name,
        :population => PooledArray => :population,
        :locus => (i -> PooledArray(i |> Vector{String})) => :locus,
        :genotype
    )
    # replace missing genotypes as missing
    stacked_geno_df.genotype = map(stacked_geno_df.genotype) do geno
        if all(0 .== geno)
            return missing
        else
            return geno
        end
    end
    sort!(stacked_geno_df, [:name, :locus])
    # ploidy finding
    #meta_df = generate_meta(stacked_geno_df)
    meta_df = DataFrames.combine(DataFrames.groupby(stacked_geno_df, :name),
        :genotype => (i -> Int8(length(first(skipmissing(i))))) => :ploidy
    )
    insertcols!(meta_df, 2, :population => "missing")
    insertcols!(meta_df, 4, :longitude => Vector{Union{Missing, Float32}}(undef, nsamples), :latitude => Vector{Union{Missing, Float32}}(undef, nsamples))

    pd_out = PopData(meta_df, stacked_geno_df)
    !allow_monomorphic && drop_monomorphic!(pd_out)
    return pd_out
end


text:
include("../utils.jl")

using RecursiveArrayTools: ArrayPartition

@testset "Product manifold" begin
    @test_throws MethodError ProductManifold()
    M1 = Sphere(2)
    M2 = Euclidean(2)
    @test (@inferred ProductManifold(M1, M2)) isa ProductManifold
    Mse = ProductManifold(M1, M2)
    @test Mse == M1 × M2
    @test Mse == ProductManifold(M1) × M2
    @test Mse == ProductManifold(M1) × ProductManifold(M2)
    @test Mse == M1 × ProductManifold(M2)
    @test Mse[1] == M1
    @test Mse[2] == M2
    @test injectivity_radius(Mse) ≈ π
    @test injectivity_radius(
        Mse,
        ProductRetraction(ExponentialRetraction(), ExponentialRetraction()),
    ) ≈ π
    @test injectivity_radius(Mse, ExponentialRetraction()) ≈ π
    @test injectivity_radius(
        Mse,
        ProductRepr([0.0, 1.0, 0.0], [0.0, 0.0]),
        ProductRetraction(ExponentialRetraction(), ExponentialRetraction()),
    ) ≈ π
    @test injectivity_radius(
        Mse,
        ProductRepr([0.0, 1.0, 0.0], [0.0, 0.0]),
        ExponentialRetraction(),
    ) ≈ π
    @test is_default_metric(Mse, ProductMetric())

    @test Manifolds.number_of_components(Mse) == 2
    # test that arrays are not points
    @test_throws DomainError is_point(Mse, [1, 2], true)
    @test check_point(Mse, [1, 2]) isa DomainError
    @test_throws DomainError is_vector(Mse, 1, [1, 2], true; check_base_point=false)
    @test check_vector(Mse, 1, [1, 2]; check_base_point=false) isa DomainError
    #default fallbacks for check_size, Product not working with Arrays
    @test Manifolds.check_size(Mse, zeros(2)) isa DomainError
    @test Manifolds.check_size(Mse, zeros(2), zeros(3)) isa DomainError
    types = [Vector{Float64}]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_STATIC_SIZED && push!(types, MVector{5,Float64})

    retraction_methods = [
        Manifolds.ProductRetraction(
            ManifoldsBase.ExponentialRetraction(),
            ManifoldsBase.ExponentialRetraction(),
        ),
    ]
    inverse_retraction_methods = [
        Manifolds.InverseProductRetraction(
            ManifoldsBase.LogarithmicInverseRetraction(),
            ManifoldsBase.LogarithmicInverseRetraction(),
        ),
    ]

    @testset "get_component, set_component!, getindex and setindex!" begin
        p1 = ProductRepr([0.0, 1.0, 0.0], [0.0, 0.0])
        @test get_component(Mse, p1, 1) == p1.parts[1]
        @test get_component(Mse, p1, Val(1)) == p1.parts[1]
        @test p1[Mse, 1] == p1.parts[1]
        @test p1[Mse, Val(1)] == p1.parts[1]
        @test p1[Mse, 1] isa Vector
        @test p1[Mse, Val(1)] isa Vector
        p2 = [10.0, 12.0]
        set_component!(Mse, p1, p2, 2)
        @test get_component(Mse, p1, 2) == p2
        p1[Mse, 2] = 2 * p2
        @test p1[Mse, 2] == 2 * p2
        p3 = [11.0, 15.0]
        set_component!(Mse, p1, p3, Val(2))
        @test get_component(Mse, p1, Val(2)) == p3
        p1[Mse, Val(2)] = 2 * p3
        @test p1[Mse, Val(2)] == 2 * p3

        p1ap = ArrayPartition([0.0, 1.0, 0.0], [0.0, 0.0])
        @test get_component(Mse, p1ap, 1) == p1ap.x[1]
        @test get_component(Mse, p1ap, Val(1)) == p1ap.x[1]
        @test p1ap[Mse, 1] == p1ap.x[1]
        @test p1ap[Mse, Val(1)] == p1ap.x[1]
        @test p1ap[Mse, 1] isa Vector
        @test p1ap[Mse, Val(1)] isa Vector
        set_component!(Mse, p1ap, p2, 2)
        @test get_component(Mse, p1ap, 2) == p2
        p1ap[Mse, 2] = 2 * p2
        @test p1ap[Mse, 2] == 2 * p2
        p3 = [11.0, 15.0]
        set_component!(Mse, p1ap, p3, Val(2))
        @test get_component(Mse, p1ap, Val(2)) == p3
        p1ap[Mse, Val(2)] = 2 * p3
        @test p1ap[Mse, Val(2)] == 2 * p3

        p1c = copy(p1)
        p1c.parts[1][1] = -123.0
        @test p1c.parts[1][1] == -123.0
        @test p1.parts[1][1] == 0.0
        copyto!(p1c, p1)
        @test p1c.parts[1][1] == 0.0
    end

    @testset "copyto!" begin
        p = ProductRepr([0.0, 1.0, 0.0], [0.0, 0.0])
        X = ProductRepr([1.0, 0.0, 0.0], [1.0, 0.0])
        q = allocate(p)
        copyto!(Mse, q, p)
        @test p.parts == q.parts
        Y = allocate(X)
        copyto!(Mse, Y, p, X)
        @test Y.parts == X.parts
    end

    @testset "Broadcasting" begin
        p1 = ProductRepr([0.0, 1.0, 0.0], [0.0, 1.0])
        p2 = ProductRepr([3.0, 4.0, 5.0], [2.0, 5.0])
        br_result = p1 .+ 2.0 .* p2
        @test br_result isa ProductRepr
        @test br_result.parts[1] ≈ [6.0, 9.0, 10.0]
        @test br_result.parts[2] ≈ [4.0, 11.0]

        br_result .= 2.0 .* p1 .+ p2
        @test br_result.parts[1] ≈ [3.0, 6.0, 5.0]
        @test br_result.parts[2] ≈ [2.0, 7.0]

        br_result .= p1
        @test br_result.parts[1] ≈ [0.0, 1.0, 0.0]
        @test br_result.parts[2] ≈ [0.0, 1.0]

        @test axes(p1) == (Base.OneTo(2),)

        # errors
        p3 = ProductRepr([3.0, 4.0, 5.0], [2.0, 5.0], [3.0, 2.0])
        @test_throws DimensionMismatch p1 .+ p3
        @test_throws DimensionMismatch p1 .= p3
    end

    @testset "CompositeManifoldError" begin
        Mpr = Sphere(2) × Sphere(2)
        p1 = [1.0, 0.0, 0.0]
        p2 = [0.0, 1.0, 0.0]
        X1 = [0.0, 1.0, 0.2]
        X2 = [1.0, 0.0, 0.2]
        p = ProductRepr(p1, p2)
        X = ProductRepr(X1, X2)
        pf = ProductRepr(p1, X1)
        Xf = ProductRepr(X1, p2)
        @test is_point(Mpr, p, true)
        @test_throws CompositeManifoldError is_point(Mpr, X, true)
        @test_throws ComponentManifoldError is_vector(Mpr, pf, X, true)
        @test_throws ComponentManifoldError is_vector(Mpr, p, Xf, true)
    end

    @testset "arithmetic" begin
        Mee = ProductManifold(Euclidean(3), Euclidean(2))
        p1 = ProductRepr([0.0, 1.0, 0.0], [0.0, 1.0])
        p2 = ProductRepr([1.0, 2.0, 0.0], [2.0, 3.0])

        @test isapprox(Mee, p1 + p2, ProductRepr([1.0, 3.0, 0.0], [2.0, 4.0]))
        @test isapprox(Mee, p1 - p2, ProductRepr([-1.0, -1.0, 0.0], [-2.0, -2.0]))
        @test isapprox(Mee, -p1, ProductRepr([0.0, -1.0, 0.0], [0.0, -1.0]))
        @test isapprox(Mee, p1 * 2, ProductRepr([0.0, 2.0, 0.0], [0.0, 2.0]))
        @test isapprox(Mee, 2 * p1, ProductRepr([0.0, 2.0, 0.0], [0.0, 2.0]))
        @test isapprox(Mee, p1 / 2, ProductRepr([0.0, 0.5, 0.0], [0.0, 0.5]))
    end

    @testset "Show methods" begin
        Mse2 = ProductManifold(M1, M1, M2, M2)
        @test sprint(show, Mse2) == "ProductManifold($(M1), $(M1), $(M2), $(M2))"
        withenv("LINES" => 10, "COLUMNS" => 100) do
            @test sprint(show, "text/plain", ProductManifold(M1)) ==
                  "ProductManifold with 1 submanifold:\n $(M1)"
            @test sprint(show, "text/plain", Mse2) ==
                  "ProductManifold with 4 submanifolds:\n $(M1)\n $(M1)\n $(M2)\n $(M2)"
            return nothing
        end
        withenv("LINES" => 7, "COLUMNS" => 100) do
            @test sprint(show, "text/plain", Mse2) ==
                  "ProductManifold with 4 submanifolds:\n $(M1)\n ⋮\n $(M2)"
            return nothing
        end

        @test sprint(show, "text/plain", ProductManifold(Mse, Mse)) == """
        ProductManifold with 2 submanifolds:
         ProductManifold(Sphere(2, ℝ), Euclidean(2; field = ℝ))
         ProductManifold(Sphere(2, ℝ), Euclidean(2; field = ℝ))"""

        p = Manifolds.ProductRepr(Float64[1, 0, 0])
        @test sprint(show, "text/plain", p) == """
        ProductRepr with 1 submanifold component:
         Component 1 =
          3-element $(sprint(show, Vector{Float64})):
           1.0
           0.0
           0.0"""

        p = Manifolds.ProductRepr(
            Float64[1, 0, 0],
            Float64[0, 1, 0],
            Float64[1, 2],
            Float64[3, 4],
        )
        @test sprint(show, "text/plain", p) == """
        ProductRepr with 4 submanifold components:
         Component 1 =
          3-element $(sprint(show, Vector{Float64})):
           1.0
           0.0
           0.0
         Component 2 =
          3-element $(sprint(show, Vector{Float64})):
           0.0
           1.0
           0.0
         Component 3 =
          2-element $(sprint(show, Vector{Float64})):
           1.0
           2.0
         Component 4 =
          2-element $(sprint(show, Vector{Float64})):
           3.0
           4.0"""

        p = Manifolds.ProductRepr(
            Float64[1, 0, 0],
            Float64[0, 1, 0],
            Float64[1, 2],
            Float64[3, 4],
            Float64[5, 6],
        )
        @test sprint(show, "text/plain", p) == """
        ProductRepr with 5 submanifold components:
         Component 1 =
          3-element $(sprint(show, Vector{Float64})):
           1.0
           0.0
           0.0
         Component 2 =
          3-element $(sprint(show, Vector{Float64})):
           0.0
           1.0
           0.0
         ⋮
         Component 4 =
          2-element $(sprint(show, Vector{Float64})):
           3.0
           4.0
         Component 5 =
          2-element $(sprint(show, Vector{Float64})):
           5.0
           6.0"""
    end

    M3 = Manifolds.Rotations(2)
    Mser = ProductManifold(M1, M2, M3)

    @test submanifold(Mser, 2) == M2
    @test (@inferred submanifold(Mser, Val((1, 3)))) == M1 × M3
    @test submanifold(Mser, 2:3) == M2 × M3
    @test submanifold(Mser, [1, 3]) == M1 × M3

    pts_sphere = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    pts_r2 = [[0.0, 0.0], [1.0, 0.0], [0.0, 0.1]]
    angles = (0.0, π / 2, 2π / 3)
    pts_rot = [[cos(ϕ) sin(ϕ); -sin(ϕ) cos(ϕ)] for ϕ in angles]
    pts = [ProductRepr(p[1], p[2], p[3]) for p in zip(pts_sphere, pts_r2, pts_rot)]
    test_manifold(
        Mser,
        pts,
        test_injectivity_radius=false,
        is_tangent_atol_multiplier=1,
        exp_log_atol_multiplier=1,
        test_inplace=true,
        test_rand_point=true,
        test_rand_tvector=true,
    )

    @testset "product vector transport" begin
        p = ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0])
        q = ProductRepr([0.0, 1.0, 0.0], [2.0, 0.0])
        X = log(Mse, p, q)
        m = ProductVectorTransport(ParallelTransport(), ParallelTransport())
        Y = vector_transport_to(Mse, p, X, q, m)
        Z = -log(Mse, q, p)
        @test isapprox(Mse, q, Y, Z)
    end

    @testset "Implicit product vector transport" begin
        p = ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0])
        q = ProductRepr([0.0, 1.0, 0.0], [2.0, 0.0])
        X = log(Mse, p, q)
        for m in [ParallelTransport(), SchildsLadderTransport(), PoleLadderTransport()]
            Y = vector_transport_to(Mse, p, X, q, m)
            Z1 = vector_transport_to(
                Mse.manifolds[1],
                submanifold_component.([p, X, q], Ref(1))...,
                m,
            )
            Z2 = vector_transport_to(
                Mse.manifolds[2],
                submanifold_component.([p, X, q], Ref(2))...,
                m,
            )
            Z = ProductRepr(Z1, Z2)
            @test isapprox(Mse, q, Y, Z)
        end
        for m in [ParallelTransport(), SchildsLadderTransport(), PoleLadderTransport()]
            Y = vector_transport_direction(Mse, p, X, X, m)
            Z1 = vector_transport_direction(
                Mse.manifolds[1],
                submanifold_component.([p, X, X], Ref(1))...,
                m,
            )
            Z2 = vector_transport_direction(
                Mse.manifolds[2],
                submanifold_component.([p, X, X], Ref(2))...,
                m,
            )
            Z = ProductRepr(Z1, Z2)
            @test isapprox(Mse, q, Y, Z)
        end
    end
    @testset "Parallel transport" begin
        p = ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0])
        q = ProductRepr([0.0, 1.0, 0.0], [2.0, 0.0])
        X = log(Mse, p, q)
        # to
        Y = parallel_transport_to(Mse, p, X, q)
        Z1 = parallel_transport_to(
            Mse.manifolds[1],
            submanifold_component.([p, X, q], Ref(1))...,
        )
        Z2 = parallel_transport_to(
            Mse.manifolds[2],
            submanifold_component.([p, X, q], Ref(2))...,
        )
        Z = ProductRepr(Z1, Z2)
        @test isapprox(Mse, q, Y, Z)
        Ym = allocate(Y)
        parallel_transport_to!(Mse, Ym, p, X, q)
        @test isapprox(Mse, q, Y, Z)

        # direction
        Y = parallel_transport_direction(Mse, p, X, X)
        Z1 = parallel_transport_direction(
            Mse.manifolds[1],
            submanifold_component.([p, X, X], Ref(1))...,
        )
        Z2 = parallel_transport_direction(
            Mse.manifolds[2],
            submanifold_component.([p, X, X], Ref(2))...,
        )
        Z = ProductRepr(Z1, Z2)
        @test isapprox(Mse, q, Y, Z)
        Ym = allocate(Y)
        parallel_transport_direction!(Mse, Ym, p, X, X)
        @test isapprox(Mse, q, Ym, Z)
    end

    @testset "ProductRepr" begin
        @test (@inferred convert(
            ProductRepr{Tuple{T,Float64,T} where T},
            ProductRepr(9, 10, 11),
        )) == ProductRepr(9, 10.0, 11)
        @test (@inferred convert(ProductRepr, ProductRepr(9, 10, 11))) ===
              ProductRepr(9, 10, 11)

        p = ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0])
        @test p == ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0])
        @test submanifold_component(Mse, p, 1) === p.parts[1]
        @test submanifold_component(Mse, p, Val(1)) === p.parts[1]
        @test submanifold_component(p, 1) === p.parts[1]
        @test submanifold_component(p, Val(1)) === p.parts[1]
        @test submanifold_components(Mse, p) === p.parts
        @test submanifold_components(p) === p.parts
        @test allocate(p, Int, 10) isa Vector{Int}
        @test length(allocate(p, Int, 10)) == 10
    end

    @testset "ArrayPartition" begin
        p = ArrayPartition([1.0, 0.0, 0.0], [0.0, 0.0])
        @test submanifold_component(Mse, p, 1) === p.x[1]
        @test submanifold_component(Mse, p, Val(1)) === p.x[1]
        @test submanifold_component(p, 1) === p.x[1]
        @test submanifold_component(p, Val(1)) === p.x[1]
        @test submanifold_components(Mse, p) === p.x
        @test submanifold_components(p) === p.x
    end

    for TP in [ProductRepr, ArrayPartition]
        @testset "TP=$TP" begin
            Ts = SizedVector{3,Float64}
            Tr2 = SizedVector{2,Float64}
            pts_sphere = [
                convert(Ts, [1.0, 0.0, 0.0]),
                convert(Ts, [0.0, 1.0, 0.0]),
                convert(Ts, [0.0, 0.0, 1.0]),
            ]
            pts_r2 = [
                convert(Tr2, [0.0, 0.0]),
                convert(Tr2, [1.0, 0.0]),
                convert(Tr2, [0.0, 0.1]),
            ]

            pts = [TP(p[1], p[2]) for p in zip(pts_sphere, pts_r2)]
            basis_types = (
                DefaultOrthonormalBasis(),
                ProjectedOrthonormalBasis(:svd),
                get_basis(Mse, pts[1], DefaultOrthonormalBasis()),
                DiagonalizingOrthonormalBasis(
                    TP(SizedVector{3}([0.0, 1.0, 0.0]), SizedVector{2}([1.0, 0.0])),
                ),
            )
            distr_M1 = Manifolds.uniform_distribution(M1, pts_sphere[1])
            distr_M2 = Manifolds.projected_distribution(
                M2,
                Distributions.MvNormal(zero(pts_r2[1]), 1.0),
            )
            distr_tv_M1 = Manifolds.normal_tvector_distribution(M1, pts_sphere[1], 1.0)
            distr_tv_M2 = Manifolds.normal_tvector_distribution(M2, pts_r2[1], 1.0)
            @test injectivity_radius(Mse, pts[1]) ≈ π
            @test injectivity_radius(Mse) ≈ π
            @test injectivity_radius(Mse, pts[1], ExponentialRetraction()) ≈ π
            @test injectivity_radius(Mse, ExponentialRetraction()) ≈ π

            @test ManifoldsBase.allocate_coordinates(
                Mse,
                pts[1],
                Float64,
                number_of_coordinates(Mse, DefaultOrthogonalBasis()),
            ) isa Vector{Float64}

            test_manifold(
                Mse,
                pts;
                point_distributions=[
                    Manifolds.ProductPointDistribution(distr_M1, distr_M2),
                ],
                tvector_distributions=[
                    Manifolds.ProductFVectorDistribution(distr_tv_M1, distr_tv_M2),
                ],
                test_injectivity_radius=true,
                test_musical_isomorphisms=true,
                musical_isomorphism_bases=[DefaultOrthonormalBasis()],
                test_tangent_vector_broadcasting=true,
                test_project_tangent=true,
                test_project_point=true,
                test_mutating_rand=false,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                test_riesz_representer=true,
                test_default_vector_transport=true,
                test_rand_point=true,
                test_rand_tvector=true,
                vector_transport_methods=[
                    ProductVectorTransport(ParallelTransport(), ParallelTransport()),
                    ProductVectorTransport(
                        SchildsLadderTransport(),
                        SchildsLadderTransport(),
                    ),
                    ProductVectorTransport(PoleLadderTransport(), PoleLadderTransport()),
                ],
                basis_types_vecs=(basis_types[1], basis_types[3], basis_types[4]),
                basis_types_to_from=basis_types,
                is_tangent_atol_multiplier=1,
                exp_log_atol_multiplier=1,
            )
            @test number_eltype(pts[1]) === Float64

            @test (@inferred ManifoldsBase._get_vector_cache_broadcast(pts[1])) ===
                  Val(false)
        end
    end

    @testset "vee/hat" begin
        M1 = Rotations(3)
        M2 = Euclidean(3)
        M = M1 × M2

        e = Matrix{Float64}(I, 3, 3)
        p = ProductRepr(exp(M1, e, hat(M1, e, [1.0, 2.0, 3.0])), [1.0, 2.0, 3.0])
        X = [0.1, 0.2, 0.3, -1.0, 2.0, -3.0]

        Xc = hat(M, p, X)
        X2 = vee(M, p, Xc)
        @test isapprox(X, X2)
    end

    @testset "get_coordinates" begin
        # make sure `get_coordinates` does not return an `ArrayPartition`
        p1 = ProductRepr([0.0, 1.0, 0.0], [0.0, 0.0])
        X1 = ProductRepr([1.0, 0.0, -1.0], [1.0, 0.0])
        Tp1Mse = TangentSpaceAtPoint(Mse, p1)
        c = get_coordinates(Tp1Mse, p1, X1, DefaultOrthonormalBasis())
        @test c isa Vector

        p1ap = ArrayPartition([0.0, 1.0, 0.0], [0.0, 0.0])
        X1ap = ArrayPartition([1.0, 0.0, -1.0], [1.0, 0.0])
        Tp1apMse = TangentSpaceAtPoint(Mse, p1ap)
        cap = get_coordinates(Tp1apMse, p1ap, X1ap, DefaultOrthonormalBasis())
        @test cap isa Vector
    end

    @testset "Basis printing" begin
        p = ProductRepr([1.0, 0.0, 0.0], [1.0, 0.0])
        B = DefaultOrthonormalBasis()
        Bc = get_basis(Mse, p, B)
        Bc_components_s = sprint.(show, "text/plain", Bc.data.parts)
        @test sprint(show, "text/plain", Bc) == """
        $(typeof(B)) for a product manifold
        Basis for component 1:
        $(Bc_components_s[1])
        Basis for component 2:
        $(Bc_components_s[2])
        """
    end

    @testset "Basis-related errors" begin
        a = ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0])
        B = CachedBasis(DefaultOrthonormalBasis(), ProductBasisData(([],)))
        @test_throws AssertionError get_vector!(
            Mse,
            a,
            ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0]),
            [1.0, 2.0, 3.0, 4.0, 5.0], # this is one element too long, hence assertionerror
            B,
        )
        @test_throws MethodError get_vector!(
            Mse,
            a,
            ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0]),
            [1.0, 2.0, 3.0, 4.0],
            B, # empty elements yield a submanifold MethodError
        )
    end

    @testset "allocation promotion" begin
        M2c = Euclidean(2; field=ℂ)
        Msec = ProductManifold(M1, M2c)
        @test Manifolds.allocation_promotion_function(Msec, get_vector, ()) === complex
        @test Manifolds.allocation_promotion_function(Mse, get_vector, ()) === identity
    end

    @testset "empty allocation" begin
        p = allocate_result(Mse, uniform_distribution)
        @test isa(p, ProductRepr)
        @test size(p[Mse, 1]) == (3,)
        @test size(p[Mse, 2]) == (2,)
    end

    @testset "Uniform distribution" begin
        Mss = ProductManifold(Sphere(2), Sphere(2))
        p = rand(uniform_distribution(Mss))
        @test is_point(Mss, p)
        @test is_point(Mss, rand(uniform_distribution(Mss, p)))
    end

    @testset "Atlas & Induced Basis" begin
        M = ProductManifold(Euclidean(2), Euclidean(2))
        p = ProductRepr(zeros(2), ones(2))
        X = ProductRepr(ones(2), 2 .* ones(2))
        A = RetractionAtlas()
        a = get_parameters(M, A, p, p)
        p2 = get_point(M, A, p, a)
        @test all(p2.parts .== p.parts)
    end

    @testset "metric conversion" begin
        M = SymmetricPositiveDefinite(3)
        N = ProductManifold(M, M)
        e = EuclideanMetric()
        p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]
        q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]
        P = ProductRepr(p, q)
        X = ProductRepr(log(M, p, q), log(M, q, p))
        Y = change_metric(N, e, P, X)
        Yc = ProductRepr(
            change_metric(M, e, p, log(M, p, q)),
            change_metric(M, e, q, log(M, q, p)),
        )
        @test norm(N, P, Y - Yc) ≈ 0
        Z = change_representer(N, e, P, X)
        Zc = ProductRepr(
            change_representer(M, e, p, log(M, p, q)),
            change_representer(M, e, q, log(M, q, p)),
        )
        @test norm(N, P, Z - Zc) ≈ 0
    end
end


text:
using Base.Meta

"""
Removes boxing of native types:
    (:call, :(Base.box), Type, val)
    (:call, :(Base.Math.box), Type, val)
    (:call, :(Base.FastMath.box), Type, val)
    => val
"""
function unbox(tree)
    depthfirst(tree) do e
        if isexpr(e, :call)
            if e.args[1] in (GlobalRef(Base,:box),
                             GlobalRef(Base.Math,:box),
                             GlobalRef(Base.FastMath,:box))
                return e.args[3]
            end
        end
        return e
    end
end


text:
using Alea, FactCheck

facts("Distributions") do

  context("Uniform") do
    U=Uniform(0, 4)
    @fact mean(U) --> 2 "mean"
    @fact var(U) --> roughly(4//3) "variance"
  end

  @fact Statistics._dummy_func_for_code_coverage(1) --> 1

end


text:
using LinearAlgebra, Elemental, Test

m, n = 100, 80
A = Elemental.Matrix(Float64)
Elemental.gaussian!(A, m, n)
U_Elemental, s_Elemental, V_Elemental = svd(A)
U_LAPACK, s_LAPACK, V_LAPACK = svd(Matrix(A))
@test abs.(Matrix(U_Elemental)'*U_LAPACK) ≈ Matrix(I, n, n)
@test Matrix(s_Elemental) ≈ s_LAPACK
@test abs.(Matrix(V_Elemental)'*V_LAPACK) ≈ Matrix(I, n, n)


text:
int_rules_4_1_12 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.12*(e*x)^m*(a+b*sin(c+d*x^n))^p =#
    @apply_utils Antiderivative(sin(~(d') * (~(e') + ~(f') * ~x) ^ 2), ~x) => (sqrt(pi / 2) / (~f * Rt(~d, 2))) * FresnelS(sqrt(2 / pi) * Rt(~d, 2) * (~e + ~f * ~x)) <-- FreeQ([~d, ~e, ~f], ~x)
    @apply_utils Antiderivative(cos(~(d') * (~(e') + ~(f') * ~x) ^ 2), ~x) => (sqrt(pi / 2) / (~f * Rt(~d, 2))) * FresnelC(sqrt(2 / pi) * Rt(~d, 2) * (~e + ~f * ~x)) <-- FreeQ([~d, ~e, ~f], ~x)
    @apply_utils Antiderivative(sin(~c + ~(d') * (~(e') + ~(f') * ~x) ^ 2), ~x) => sin(~c) * Antiderivative(cos(~d * (~e + ~f * ~x) ^ 2), ~x) + cos(~c) * Antiderivative(sin(~d * (~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x)
    @apply_utils Antiderivative(cos(~c + ~(d') * (~(e') + ~(f') * ~x) ^ 2), ~x) => cos(~c) * Antiderivative(cos(~d * (~e + ~f * ~x) ^ 2), ~x) - sin(~c) * Antiderivative(sin(~d * (~e + ~f * ~x) ^ 2), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n), ~x) => (im / 2) * Antiderivative(ℯ ^ (-(~c) * im - ~d * im * (~e + ~f * ~x) ^ ~n), ~x) - (im / 2) * Antiderivative(ℯ ^ (~c * im + ~d * im * (~e + ~f * ~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && IGtQ(~n, 2)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n), ~x) => (1 / 2) * Antiderivative(ℯ ^ (-(~c) * im - ~d * im * (~e + ~f * ~x) ^ ~n), ~x) + (1 / 2) * Antiderivative(ℯ ^ (~c * im + ~d * im * (~e + ~f * ~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~f], ~x) && IGtQ(~n, 2)
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * sin(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 1))
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * cos(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 1))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (-1 / ~f) * Subst(Antiderivative((~a + ~b * sin(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ 2, ~x), ~x, 1 / (~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && EqQ(~n, -2)))
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (-1 / ~f) * Subst(Antiderivative((~a + ~b * cos(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ 2, ~x), ~x, 1 / (~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && EqQ(~n, -2)))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~n * ~f)) * Subst(Antiderivative((~x) ^ (1 / ~n - 1) * (~a + ~b * sin(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~e + ~f * ~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && IntegerQ(1 / ~n))
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~n * ~f)) * Subst(Antiderivative((~x) ^ (1 / ~n - 1) * (~a + ~b * cos(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~e + ~f * ~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && IntegerQ(1 / ~n))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], (k / ~f) * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * sin(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~e + ~f * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && FractionQ(~n))
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], (k / ~f) * Subst(Antiderivative((~x) ^ (k - 1) * (~a + ~b * cos(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~e + ~f * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f], ~x) && (IGtQ(~p, 0) && FractionQ(~n))
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n), ~x) => (im / 2) * Antiderivative(ℯ ^ (-(~c) * im - ~d * im * (~e + ~f * ~x) ^ ~n), ~x) - (im / 2) * Antiderivative(ℯ ^ (~c * im + ~d * im * (~e + ~f * ~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~n], ~x)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n), ~x) => (1 / 2) * Antiderivative(ℯ ^ (-(~c) * im - ~d * im * (~e + ~f * ~x) ^ ~n), ~x) + (1 / 2) * Antiderivative(ℯ ^ (~c * im + ~d * im * (~e + ~f * ~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~f, ~n], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * sin(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && IGtQ(~p, 1)
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~a + ~b * cos(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n], ~x) && IGtQ(~p, 1)
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~p, ~x) => Unintegrable((~a + ~b * sin(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~p, ~x) => Unintegrable((~a + ~b * cos(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * sin(~c + ~d * ExpandToSum(~u, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && Not(LinearMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * cos(~c + ~d * ExpandToSum(~u, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n, ~p], ~x) && (LinearQ(~u, ~x) && Not(LinearMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * sin(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * sin(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(a') + ~(b') * cos(~u)) ^ ~(p'), ~x) => Antiderivative((~a + ~b * cos(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative(sin(~(d') * (~x) ^ ~n) / ~x, ~x) => SinIntegral(~d * (~x) ^ ~n) / ~n <-- FreeQ([~d, ~n], ~x)
    @apply_utils Antiderivative(cos(~(d') * (~x) ^ ~n) / ~x, ~x) => CosIntegral(~d * (~x) ^ ~n) / ~n <-- FreeQ([~d, ~n], ~x)
    @apply_utils Antiderivative(sin(~c + ~(d') * (~x) ^ ~n) / ~x, ~x) => sin(~c) * Antiderivative(cos(~d * (~x) ^ ~n) / ~x, ~x) + cos(~c) * Antiderivative(sin(~d * (~x) ^ ~n) / ~x, ~x) <-- FreeQ([~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(cos(~c + ~(d') * (~x) ^ ~n) / ~x, ~x) => cos(~c) * Antiderivative(cos(~d * (~x) ^ ~n) / ~x, ~x) - sin(~c) * Antiderivative(sin(~d * (~x) ^ ~n) / ~x, ~x) <-- FreeQ([~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * sin(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IntegerQ(Simplify((~m + 1) / ~n)) && (EqQ(~p, 1) || (EqQ(~m, ~n - 1) || IntegerQ(~p) && GtQ(Simplify((~m + 1) / ~n), 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a + ~b * cos(~c + ~d * ~x)) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (IntegerQ(Simplify((~m + 1) / ~n)) && (EqQ(~p, 1) || (EqQ(~m, ~n - 1) || IntegerQ(~p) && GtQ(Simplify((~m + 1) / ~n), 0))))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x) && IntegerQ(Simplify((~m + 1) / ~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n), ~x) => (2 / ~n) * Subst(Antiderivative(sin(~a + ~b * (~x) ^ 2), ~x), ~x, (~x) ^ (~n / 2)) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && EqQ(~m, ~n / 2 - 1)
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n), ~x) => (2 / ~n) * Subst(Antiderivative(cos(~a + ~b * (~x) ^ 2), ~x), ~x, (~x) ^ (~n / 2)) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && EqQ(~m, ~n / 2 - 1)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sin(~(c') + ~(d') * (~x) ^ ~n), ~x) => (-((~e) ^ (~n - 1)) * (~e * ~x) ^ ((~m - ~n) + 1) * cos(~c + ~d * (~x) ^ ~n)) / (~d * ~n) + (((~e) ^ ~n * ((~m - ~n) + 1)) / (~d * ~n)) * Antiderivative((~e * ~x) ^ (~m - ~n) * cos(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(~n, ~m + 1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cos(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e) ^ (~n - 1) * (~e * ~x) ^ ((~m - ~n) + 1) * sin(~c + ~d * (~x) ^ ~n)) / (~d * ~n) - (((~e) ^ ~n * ((~m - ~n) + 1)) / (~d * ~n)) * Antiderivative((~e * ~x) ^ (~m - ~n) * sin(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(~n, ~m + 1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * sin(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e * ~x) ^ (~m + 1) * sin(~c + ~d * (~x) ^ ~n)) / (~e * (~m + 1)) - ((~d * ~n) / ((~e) ^ ~n * (~m + 1))) * Antiderivative((~e * ~x) ^ (~m + ~n) * cos(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * cos(~(c') + ~(d') * (~x) ^ ~n), ~x) => ((~e * ~x) ^ (~m + 1) * cos(~c + ~d * (~x) ^ ~n)) / (~e * (~m + 1)) + ((~d * ~n) / ((~e) ^ ~n * (~m + 1))) * Antiderivative((~e * ~x) ^ (~m + ~n) * sin(~c + ~d * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e], ~x) && (IGtQ(~n, 0) && LtQ(~m, -1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sin(~(c') + ~(d') * (~x) ^ ~n), ~x) => (im / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) * im - ~d * im * (~x) ^ ~n), ~x) - (im / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c * im + ~d * im * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cos(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) * im - ~d * im * (~x) ^ ~n), ~x) + (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c * im + ~d * im * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + (~(b') * (~x) ^ ~n) / 2) ^ 2, ~x) => (1 / 2) * Antiderivative((~x) ^ ~m, ~x) - (1 / 2) * Antiderivative((~x) ^ ~m * cos(2 * ~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + (~(b') * (~x) ^ ~n) / 2) ^ 2, ~x) => (1 / 2) * Antiderivative((~x) ^ ~m, ~x) + (1 / 2) * Antiderivative((~x) ^ ~m * cos(2 * ~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((~x) ^ (~m + 1) * sin(~a + ~b * (~x) ^ ~n) ^ ~p) / (~m + 1) - ((~b * ~n * ~p) / (~m + 1)) * Antiderivative(sin(~a + ~b * (~x) ^ ~n) ^ (~p - 1) * cos(~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b], ~x) && (IGtQ(~p, 1) && (EqQ(~m + ~n, 0) && (NeQ(~n, 1) && IntegerQ(~n))))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((~x) ^ (~m + 1) * cos(~a + ~b * (~x) ^ ~n) ^ ~p) / (~m + 1) + ((~b * ~n * ~p) / (~m + 1)) * Antiderivative(cos(~a + ~b * (~x) ^ ~n) ^ (~p - 1) * sin(~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b], ~x) && (IGtQ(~p, 1) && (EqQ(~m + ~n, 0) && (NeQ(~n, 1) && IntegerQ(~n))))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((~n * sin(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) - ((~x) ^ ~n * cos(~a + ~b * (~x) ^ ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p)) + ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * sin(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1, 0) && GtQ(~p, 1))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (~n * cos(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) + ((~x) ^ ~n * sin(~a + ~b * (~x) ^ ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p) + ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * cos(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1, 0) && GtQ(~p, 1))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((((~m - ~n) + 1) * (~x) ^ ((~m - 2 * ~n) + 1) * sin(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) - ((~x) ^ ((~m - ~n) + 1) * cos(~a + ~b * (~x) ^ ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p)) + ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * sin(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x)) - ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2)) * Antiderivative((~x) ^ (~m - 2 * ~n) * sin(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b], ~x) && (GtQ(~p, 1) && (IGtQ(~n, 0) && IGtQ(~m, 2 * ~n - 1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((((~m - ~n) + 1) * (~x) ^ ((~m - 2 * ~n) + 1) * cos(~a + ~b * (~x) ^ ~n) ^ ~p) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) + ((~x) ^ ((~m - ~n) + 1) * sin(~a + ~b * (~x) ^ ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / (~b * ~n * ~p) + ((~p - 1) / ~p) * Antiderivative((~x) ^ ~m * cos(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x)) - ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2)) * Antiderivative((~x) ^ (~m - 2 * ~n) * cos(~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b], ~x) && (GtQ(~p, 1) && (IGtQ(~n, 0) && IGtQ(~m, 2 * ~n - 1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((((~x) ^ (~m + 1) * sin(~a + ~b * (~x) ^ ~n) ^ ~p) / (~m + 1) - (~b * ~n * ~p * (~x) ^ (~m + ~n + 1) * cos(~a + ~b * (~x) ^ ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) - (((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * sin(~a + ~b * (~x) ^ ~n) ^ ~p, ~x)) + (((~b) ^ 2 * (~n) ^ 2 * ~p * (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b], ~x) && (GtQ(~p, 1) && (IGtQ(~n, 0) && (ILtQ(~m, -2 * ~n + 1) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((((~x) ^ (~m + 1) * cos(~a + ~b * (~x) ^ ~n) ^ ~p) / (~m + 1) + (~b * ~n * ~p * (~x) ^ (~m + ~n + 1) * sin(~a + ~b * (~x) ^ ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) - (((~b) ^ 2 * (~n) ^ 2 * (~p) ^ 2) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * cos(~a + ~b * (~x) ^ ~n) ^ ~p, ~x)) + (((~b) ^ 2 * (~n) ^ 2 * ~p * (~p - 1)) / ((~m + 1) * (~m + ~n + 1))) * Antiderivative((~x) ^ (~m + 2 * ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p - 2), ~x) <-- FreeQ([~a, ~b], ~x) && (GtQ(~p, 1) && (IGtQ(~n, 0) && (ILtQ(~m, -2 * ~n + 1) && NeQ(~m + ~n + 1, 0))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (k / ~e) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * sin(~c + (~d * (~x) ^ (k * ~n)) / (~e) ^ ~n)) ^ ~p, ~x), ~x, (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~p) && (IGtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (k / ~e) * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * cos(~c + (~d * (~x) ^ (k * ~n)) / (~e) ^ ~n)) ^ ~p, ~x), ~x, (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IntegerQ(~p) && (IGtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~p, 1) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~x) ^ ~n * cos(~a + ~b * (~x) ^ ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (~n * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) + ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1, 0) && (LtQ(~p, -1) && NeQ(~p, -2)))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((-((~x) ^ ~n) * sin(~a + ~b * (~x) ^ ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (~n * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) + ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (EqQ((~m - 2 * ~n) + 1, 0) && (LtQ(~p, -1) && NeQ(~p, -2)))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~x) ^ ((~m - ~n) + 1) * cos(~a + ~b * (~x) ^ ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) * (~x) ^ ((~m - 2 * ~n) + 1) * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) + ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) + ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) * Antiderivative((~x) ^ (~m - 2 * ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b], ~x) && (LtQ(~p, -1) && (NeQ(~p, -2) && (IGtQ(~n, 0) && IGtQ(~m, 2 * ~n - 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => ((-((~x) ^ ((~m - ~n) + 1)) * sin(~a + ~b * (~x) ^ ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) * (~x) ^ ((~m - 2 * ~n) + 1) * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 2)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) + ((~p + 2) / (~p + 1)) * Antiderivative((~x) ^ ~m * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) + ((((~m - ~n) + 1) * ((~m - 2 * ~n) + 1)) / ((~b) ^ 2 * (~n) ^ 2 * (~p + 1) * (~p + 2))) * Antiderivative((~x) ^ (~m - 2 * ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 2), ~x) <-- FreeQ([~a, ~b], ~x) && (LtQ(~p, -1) && (NeQ(~p, -2) && (IGtQ(~n, 0) && IGtQ(~m, 2 * ~n - 1))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~a + ~b * sin(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && (IntegerQ(~m) && EqQ(~n, -2))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -(Subst(Antiderivative((~a + ~b * cos(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && (IntegerQ(~m) && EqQ(~n, -2))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (-k / ~e) * Subst(Antiderivative((~a + ~b * sin(~c + ~d / ((~e) ^ ~n * (~x) ^ (k * ~n)))) ^ ~p / (~x) ^ (k * (~m + 1) + 1), ~x), ~x, 1 / (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => With([k = Denominator(~m)], (-k / ~e) * Subst(Antiderivative((~a + ~b * cos(~c + ~d / ((~e) ^ ~n * (~x) ^ (k * ~n)))) ^ ~p / (~x) ^ (k * (~m + 1) + 1), ~x), ~x, 1 / (~e * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -((~e * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative((~a + ~b * sin(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && Not(RationalQ(~m))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~m * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => -((~e * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative((~a + ~b * cos(~c + ~d * (~x) ^ -(~n))) ^ ~p / (~x) ^ (~m + 2), ~x), ~x, 1 / ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IGtQ(~p, 0) && (ILtQ(~n, 0) && Not(RationalQ(~m))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * sin(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = Denominator(~n)], k * Subst(Antiderivative((~x) ^ (k * (~m + 1) - 1) * (~a + ~b * cos(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, ~x), ~x, (~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && (IntegerQ(~p) && FractionQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a + ~b * sin(~c + ~d * (~x) ^ Simplify(~n / (~m + 1)))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a + ~b * cos(~c + ~d * (~x) ^ Simplify(~n / (~m + 1)))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~e * ~x) ^ ~m * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && (IntegerQ(~p) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && Not(IntegerQ(~n)))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sin(~(c') + ~(d') * (~x) ^ ~n), ~x) => (im / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) * im - ~d * im * (~x) ^ ~n), ~x) - (im / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c * im + ~d * im * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cos(~(c') + ~(d') * (~x) ^ ~n), ~x) => (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (-(~c) * im - ~d * im * (~x) ^ ~n), ~x) + (1 / 2) * Antiderivative((~e * ~x) ^ ~m * ℯ ^ (~c * im + ~d * im * (~x) ^ ~n), ~x) <-- FreeQ([~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~p, ~x) => Antiderivative(ExpandTrigReduce((~e * ~x) ^ ~m, (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && IGtQ(~p, 0)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~e * ~x) ^ ~m * (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~e * ~x) ^ ~m * (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * sin(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~u)) ^ ~(p'), ~x) => Antiderivative((~e * ~x) ^ ~m * (~a + ~b * cos(ExpandToSum(~u, ~x))) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~m, ~p], ~x) && (BinomialQ(~u, ~x) && Not(BinomialMatchQ(~u, ~x)))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~n * ~f)) * Subst(Antiderivative(ExpandIntegrand((~a + ~b * sin(~c + ~d * ~x)) ^ ~p, (~x) ^ (1 / ~n - 1) * ((~g - (~e * ~h) / ~f) + (~h * (~x) ^ (1 / ~n)) / ~f) ^ ~m, ~x), ~x), ~x, (~e + ~f * ~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IGtQ(~p, 0) && IntegerQ(1 / ~n))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (1 / (~n * ~f)) * Subst(Antiderivative(ExpandIntegrand((~a + ~b * cos(~c + ~d * ~x)) ^ ~p, (~x) ^ (1 / ~n - 1) * ((~g - (~e * ~h) / ~f) + (~h * (~x) ^ (1 / ~n)) / ~f) ^ ~m, ~x), ~x), ~x, (~e + ~f * ~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IGtQ(~p, 0) && IntegerQ(1 / ~n))
    #= Antiderivative(((~g')+(~h')*(~x))^(~m')*((~a')+(~b')*sin((~c')+(~d')*((~e')+(~f')*(~x))^(~n)))^(~p'),(~x)*Symbol) :=  1/(n*f^(m+1))*Subst(Antiderivative(ExpandIntegrand((a+b*sin(c+d*x))^p,x^(1/n-1)*( f*g-e*h+h*x^(1/n))^m,x),x),x,(e+f*x)^n) <-- FreeQ([a,b,c,d,e,f,g,h],x)  && IGtQ(p,0) && IntegerQ(m) && IntegerQ(1/n) =#
    #= Antiderivative(((~g')+(~h')*(~x))^(~m')*((~a')+(~b')*cos((~c')+(~d')*((~e')+(~f')*(~x))^(~n)))^(~p'),(~x)*Symbol) :=  1/(n*f^(m+1))*Subst(Antiderivative(ExpandIntegrand((a+b*cos(c+d*x))^p,x^(1/n-1)*( f*g-e*h+h*x^(1/n))^m,x),x),x,(e+f*x)^n) <-- FreeQ([a,b,c,d,e,f,g,h],x)  && IGtQ(p,0) && IntegerQ(m) && IntegerQ(1/n) =#
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = If(FractionQ(~n), Denominator(~n), 1)], (k / (~f) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand((~a + ~b * sin(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, (~x) ^ (k - 1) * ((~f * ~g - ~e * ~h) + ~h * (~x) ^ k) ^ ~m, ~x), ~x), ~x, (~e + ~f * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~p, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => Module([k = If(FractionQ(~n), Denominator(~n), 1)], (k / (~f) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand((~a + ~b * cos(~c + ~d * (~x) ^ (k * ~n))) ^ ~p, (~x) ^ (k - 1) * ((~f * ~g - ~e * ~h) + ~h * (~x) ^ k) ^ ~m, ~x), ~x), ~x, (~e + ~f * ~x) ^ (1 / k))) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h], ~x) && (IGtQ(~p, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~f) * Subst(Antiderivative(((~h * ~x) / ~f) ^ ~m * (~a + ~b * sin(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~e + ~f * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IGtQ(~p, 0) && EqQ(~f * ~g - ~e * ~h, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => (1 / ~f) * Subst(Antiderivative(((~h * ~x) / ~f) ^ ~m * (~a + ~b * cos(~c + ~d * (~x) ^ ~n)) ^ ~p, ~x), ~x, ~e + ~f * ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m], ~x) && (IGtQ(~p, 0) && EqQ(~f * ~g - ~e * ~h, 0))
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~g + ~h * ~x) ^ ~m * (~a + ~b * sin(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~(g') + ~(h') * ~x) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~(e') + ~(f') * ~x) ^ ~n)) ^ ~(p'), ~x) => Unintegrable((~g + ~h * ~x) ^ ~m * (~a + ~b * cos(~c + ~d * (~e + ~f * ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p], ~x)
    @apply_utils Antiderivative((~v) ^ ~(m') * (~(a') + ~(b') * sin(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~v, ~x) ^ ~m * (~a + ~b * sin(~c + ~d * ExpandToSum(~u, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (LinearQ(~u, ~x) && (LinearQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && LinearMatchQ(~v, ~x))))
    @apply_utils Antiderivative((~v) ^ ~(m') * (~(a') + ~(b') * cos(~(c') + ~(d') * (~u) ^ ~n)) ^ ~(p'), ~x) => Antiderivative(ExpandToSum(~v, ~x) ^ ~m * (~a + ~b * cos(~c + ~d * ExpandToSum(~u, ~x) ^ ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n, ~p], ~x) && (LinearQ(~u, ~x) && (LinearQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && LinearMatchQ(~v, ~x))))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * cos(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => sin(~a + ~b * (~x) ^ ~n) ^ (~p + 1) / (~b * ~n * (~p + 1)) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (EqQ(~m, ~n - 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * sin(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => -(cos(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (EqQ(~m, ~n - 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * sin(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * cos(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => ((~x) ^ ((~m - ~n) + 1) * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (((~m - ~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (~m - ~n) * sin(~a + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (LtQ(0, ~n, ~m + 1) && NeQ(~p, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * cos(~(a') + ~(b') * (~x) ^ ~(n')) ^ ~(p') * sin(~(a') + ~(b') * (~x) ^ ~(n')), ~x) => (-((~x) ^ ((~m - ~n) + 1)) * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) + (((~m - ~n) + 1) / (~b * ~n * (~p + 1))) * Antiderivative((~x) ^ (~m - ~n) * cos(~a + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~p], ~x) && (LtQ(0, ~n, ~m + 1) && NeQ(~p, -1))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.12 (e x)^m (a+b sin(c+d x^n))^p.jl


text:
count_nonzeros(IJ, mn) = [count(isequal(ij), IJ) for ij in 1:mn]

"""
    nonzero_minors(N, overlap) -> (ks,ls)

Find all (distinct) minor determinants of order `N` of the orbital
`overlap` matrix that do not vanish, i.e. all non-vanishing minors are
guaranteed to be present, but not all of the returned minors are
guaranteed to be non-zero. Vanishing minors returned arise when the
overlap matrix is rank deficient, which is unlikely to happen when
computing energy expressions, but must still be guarded against. This
is most easily checked by actually calculating the [`cofactor`](@ref),
which is most likely desired anyway.
"""
function nonzero_minors(N::Integer, overlap::AbstractSparseMatrix)
    ks = Vector{Vector{Int}}()
    ls = Vector{Vector{Int}}()

    (I,J,V) = findnz(overlap)
    m,n = size(overlap)
    m == n || throw(ArgumentError("Overlap matrix must be square"))

    # First deal with trivial special cases.
    N > m && return ks, ls
    if N == m
        push!(ks, collect(1:m))
        push!(ls, collect(1:m))
        return ks, ls
    end

    # Number of non-zero matrix elements for each row/column.
    nzrowcount = count_nonzeros(I, m)
    nzcolcount = count_nonzeros(J, n)
    # Vanishing rows/columns.
    zrows = findall(iszero, nzrowcount)
    zcols = findall(iszero, nzcolcount)
    # Minimum order of Γ⁽ᵖ⁾
    prmin = length(zrows)
    pcmin = length(zcols)
    pmin = max(prmin, pcmin)

    # Then deal with next set of trivial special cases.
    N < pmin && return ks, ls
    if N == prmin && N == pcmin
        push!(ks, zrows)
        push!(ls, zcols)
        return ks, ls
    end

    # The general case.
    Nr = N - prmin
    Nc = N - pcmin
    # Find all rows supporting a specific column j.
    I′ = [I[findall(isequal(j), J)] for j in 1:n]
    # Find all columns supporting a specific row i.
    J′ = [J[findall(isequal(i), I)] for i in 1:m]

    for k in combinations(unique(I),Nr)
        # Find all columns supported by any of the rows in the
        # combination k.
        cols = unique(vcat([J′[i] for i in k]...))
        # For each combination k of rows, we can at most strike
        # out Nc columns.
        colbudget = Nc
        colsmuststrike = Vector{Int}()

        for j in cols
            # Test if column j is solely supported by the row
            # combination k.
            colsupport = setdiff(I′[j], k)
            if isempty(colsupport)
                colbudget -= 1
                colbudget < 0 && break
                push!(colsmuststrike, j)
            end
        end
        colbudget < 0 && continue

        Ncm = length(colsmuststrike)

        # The columns we must strike out fill the column budget.
        if Ncm == Nc
            push!(ks, sort(vcat(k,zrows)))
            push!(ls, sort(vcat(colsmuststrike,zcols)))
            continue
        end

        # Find all combinations of the remaining columns.
        colchoices = combinations(setdiff(unique(J), colsmuststrike), Nc-Ncm)
        for l in colchoices
            # Find all rows that would be affected if column
            # combination l is stricken out.
            colsupport = unique(vcat([I′[j] for j in l]...))
            supported = true
            for i in colsupport
                # If the row i supporting one of the columns in the
                # column combination l is already a candidate to be
                # stricken out, it does not matter if its support
                # vanishes.
                i ∈ k && continue
                # Otherwise, if its support vanishes, then l is an
                # unviable column combination.
                rowsupport = setdiff(J′[i], l)
                if isempty(rowsupport)
                    supported = false
                    break
                end
            end
            if supported
                push!(ks, sort(vcat(k,zrows)))
                push!(ls, sort(vcat(l,colsmuststrike,zcols)))
            end
        end
    end

    ks, ls
end


text:
infile = "./src/Day13/input.txt"
f = readlines(infile)
arrival = parse(Int, f[1])
buses = parse.(Int,filter(x -> x != "x", split(f[2], ",")))

function find_next_bus(arrival, buses)
    minutes_since_departure = mod.(arrival, buses)
    if any(minutes_since_departure .== 0)
        idx = findfirst(==(0), minutes_since_departure)
        return buses[idx], 0
    end
    time_to_next = buses .- minutes_since_departure
    idx = findmin(time_to_next)
    return buses[idx[2]], idx[1]
end
id, wait =  find_next_bus(arrival, buses)
id*wait

starttimes = [(parse(Int,step), t - 1) for (t, step) in enumerate(split(f[2],",")) if step != "x"]

function find_match(start_t, multiple, nextbus)
    step, offset = nextbus
    check = start_t
    while mod(check + offset, step) != 0
        check +=  multiple
    end
    multiple = lcm(multiple, step)
    return check, multiple
end

testlist1 = "17,x,13,19"
testlist2 = "67,7,59,61"

find_match(0,17,(13,2))
find_match(102, 221, (19,3))

function find_schedule(starttimes)
    multiple, check = starttimes[1]
    for nextbus in starttimes[2:end]
        check, multiple = find_match(check, multiple, nextbus)
    end
    return check
end

a = find_schedule(starttimes)


text:
using FlowFarm
using Test

@testset "FlowFarm.jl" begin
    # Write your own tests here.
end


text:
"""
Routes fileformat

Example:
#;ROUTES;3;446;20181206;20190102;8;EAR_P
L;ABIRO1AGMTT;AP;999999999999;000000000000;ABIRO;SP;1
L;ABIRO1AGMTT;AP;999999999999;000000000000;*2TNG;DBP;2
L;ABIRO1DGMTT;DP;999999999999;000000000000;*TAN3;DBP;1
L;ABIRO1DGMTT;DP;999999999999;000000000000;ABIRO;SP;2
L;ABIRO2AGMTT;AP;999999999999;000000000000;ABIRO;SP;1
L;ABIRO2AGMTT;AP;999999999999;000000000000;*2TNG;DBP;2
L;ABIRO2DGMTT;DP;999999999999;000000000000;*1TNG;DBP;1
L;ABIRO2DGMTT;DP;999999999999;000000000000;ABIRO;SP;2
"""
module Routes

export read

struct RoutePoint
    wp::String
    location_type::String
end

struct Route
    type::String
    route::Vector{RoutePoint}
end

function read(file)
    routes = Dict{String, Route}()
    name = ""
    type = ""
    route = Vector{RoutePoint}()
    for route_line in eachline(file)
        split_route = split(route_line, ';')
        if split_route[1] == "L"
            if split_route[2] != name
                # Create route
                if split_route[2] != ""
                    routes[name] = Route(type, route)
                end
                name = split_route[2]
                type = split_route[3]
                route = Vector{RoutePoint}()
            end
            route = vcat(route, RoutePoint(split_route[6], split_route[7]))
        else
            # Initial line
        end
    end
    # Create route
    routes[name] = Route(type, route)
    return routes
end

end # module


text:
using Test

@testset "MCMCChains" begin
    # run plotting tests
    println("Plotting")
    @time include("plot_test.jl")

    # run function tests
    println("Diagnostics")
    @time include("diagnostic_tests.jl")

    # run tests for missing values
    println("Missing values")
    @time include("missing_tests.jl")

    # run tests for describing sections
    println("Sections")
    @time include("sections_tests.jl")

    # run tests for accessing parameters
    println("Accessing parameters")
    @time include("get_tests.jl")

    # run tests for serialization
    println("Serialization")
    @time include("serialization_tests.jl")

    # run tests for sampling api
    println("Sampling")
    @time include("sampling_tests.jl")

    # run tests for array constructor
    println("Array")
    @time include("arrayconstructor_tests.jl")

    # run tests for dataframe constructor
    println("Dataframe")
    @time include("dfconstructor_tests.jl")

    # run tests for dataframe summary
    println("Summary")
    @time include("summarize_tests.jl")

    # run tests for posterior stats
    println("Model statistics")
    @time include("modelstats_test.jl")

    # run tests for concatenation
    println("Concatenation")
    @time include("concatenation_tests.jl")
end


text:
# using Gadget2 standard units:
#     velocity is canonical momentum, in km/s
#     distance is in kpc/h
#     mass is in 1e10 Ms/h
const G = 43007.1
const H0 = 0.1
const C = 299792.458

# using Dark Sky cosmology:
const ΩΛ = 0.7048
const Ω0 = 1 - ΩΛ
const h = 0.6881
const σ8 = 0.8344

const RHO_CRIT_TODAY = 3*H0^2/8/π/G

# Hubble as function of scale factor a, or z
H(z) = H0*sqrt(Ω0.*(z+1.0).^3+ΩΛ)

# see here: http://mnras.oxfordjournals.org/content/322/2/419.full
OM(z) = Ω0 / (Ω0 + ΩΛ./(z+1.0).^3)
OL(z) = ΩΛ / (Ω0.*(z+1.0).^3 + ΩΛ)

function unnormalizedDa(z)
    o0 = OMa(z)
    ol = OLa(z)
    up = 5o0
    down = o0^(4/7) - ol + (1+o0/2)*(1+ol/70)
    a*up/2down
end
D(z) = unnormalizedDa(z) / unnormalizedDa(1.0)
D2(z) = -3/7*D(z)^2

F(z) = OMa(z)^(5/9)
F2(z) = 2*OMa(z)^(6/11)

# mean/critical density
ρc(z) = ρ0(z)./OM(z)

A(z) = 1.0./(z+1.0)

ρ0(z) = Ω0.*RHO_CRIT_TODAY.*(1+z).^3 # mean rho
function RXb(z, m, X)
    vX = m./ρ0(z)./X
    (vX.*3./4./π).^(1/3)
end
R200b(z, m) = RXb(z, m, 200.0)
function RXc(z, m, X)
    vX = m./ρc(z)./X
    (vX.*3./4./π).^(1/3)
end
R200c(z, m) = RXc(z, m, 200.0)

RVIR(z, m) = RXc(z, m, Δcvir(z))

# line of sight comoving distance
# (see: https://arxiv.org/pdf/astro-ph/9905116v4.pdf)
E(z) = sqrt(Ω0.*(1+z).^3+ΩΛ)
const DH = C/H0
Dc(z) = DH.*quadgk(z->1.0./E(z), 0, z)[1]

# angular distance, in arcseconds
rad_to_as(Θ) = Θ*206265
as_to_rad(Θ) = Θ/206265
Da_rad(z, dθ) = Dc(z)./(1+z).*dθ
Da_as(z, dθ) = Da_rad(z, as_to_rad(dθ))

#

function Δcvir(z)
    x = OM(z)-1.0
    18.*π.*π+82.*x-39.*x.*x
end

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
using Documenter, PGEN

makedocs(
    format = Documenter.HTML(),
    sitename = "PGEN.jl",
    authors = "Seyoon Ko",
    clean = true,
    debug = true,
    pages = [
        "PGEN.jl Tutorial" => "index.md",
        "PGEN format description" => "PGEN_description.md"
    ]
)

deploydocs(
    repo   = "github.com/OpenMendel/PGEN.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    devbranch = "main"
)


text:
module DynamicRemoteCall

using Distributed

export dynamicremotecall!

mutable struct RepeatCounter
    i::Int
    lock::Bool
    RepeatCounter() = new(0, true)
end

try_(c::RepeatCounter) = c.lock

lock_(c::RepeatCounter) = c.lock = false

unlock_(c::RepeatCounter) = c.lock = true

function readandincrement(c::RepeatCounter)
    while true
        if try_(c)
            lock_(c)
            i = c.i += 1
            unlock_(c)
            return i
        end
    end
end

function remotecall_subprocess!(f, results::Vector, counter::RepeatCounter, pid::Int, repeattimes, args...; kwargs...)
    while true
        i = readandincrement(counter)
        if i <= repeattimes
            results[i] = remotecall_fetch(f, pid, args...; kwargs...)
        else
            break
        end
    end
end

function dynamicremotecall!(f, results::Vector, workers::AbstractArray{Int}, args...; kwargs...)
    counter = RepeatCounter()
    repeattimes = length(results)
    @sync for pid in workers
        @async remotecall_subprocess!(f, results, counter, pid, repeattimes, args...; kwargs...)
    end
end

dynamicremotecall!(f, results::Vector, args...; kwargs...) = dynamicremotecall!(f, results, workers(), args...; kwargs...)

end # module

max_stars_repo_name: wangl-cc/DynamicRemoteCall-jl

max_issues_repo_name: wangl-cc/DynamicRemoteCall-jl

max_forks_repo_name: wangl-cc/DynamicRemoteCall-jl


text:
module ChowLiuTrees

include("information.jl")
include("chow_liu_trees.jl")

end # module


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    Arrow.Primitive

An `ArrowVector` where each element is a "fixed size" scalar of some kind, like an integer, float, decimal, or time type.
"""
struct Primitive{T, A} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    validity::ValidityBitmap
    data::A
    ℓ::Int64
    metadata::Union{Nothing, Dict{String, String}}
end

Primitive(::Type{T}, b::Vector{UInt8}, v::ValidityBitmap, data::A, l, meta) where {T, A} =
    Primitive{T, A}(b, v, data, l, meta)

Base.size(p::Primitive) = (p.ℓ,)

function Base.copy(p::Primitive{T, A}) where {T, A}
    if nullcount(p) == 0 && T === eltype(A)
        return copy(p.data)
    else
        return convert(Array, p)
    end
end

@propagate_inbounds function Base.getindex(p::Primitive{T}, i::Integer) where {T}
    @boundscheck checkbounds(p, i)
    if T >: Missing
        return @inbounds (p.validity[i] ? ArrowTypes.fromarrow(T, p.data[i]) : missing)
    else
        return @inbounds ArrowTypes.fromarrow(T, p.data[i])
    end
end

@propagate_inbounds function Base.setindex!(p::Primitive{T}, v, i::Integer) where {T}
    @boundscheck checkbounds(p, i)
    if T >: Missing
        if v === missing
            @inbounds p.validity[i] = false
        else
            @inbounds p.data[i] = convert(Base.nonmissingtype(T), v)
        end
    else
        @inbounds p.data[i] = convert(Base.nonmissingtype(T), v)
    end
    return v
end

arrowvector(::PrimitiveKind, x::Primitive, i, nl, fi, de, ded, meta; kw...) = x

function arrowvector(::PrimitiveKind, x, i, nl, fi, de, ded, meta; kw...)
    validity = ValidityBitmap(x)
    return Primitive(eltype(x), UInt8[], validity, x, length(x), meta)
end

function compress(Z::Meta.CompressionType, comp, p::P) where {P <: Primitive}
    len = length(p)
    nc = nullcount(p)
    validity = compress(Z, comp, p.validity)
    data = compress(Z, comp, p.data)
    return Compressed{Z, P}(p, [validity, data], len, nc, Compressed[])
end

function makenodesbuffers!(col::Primitive{T}, fieldnodes, fieldbuffers, bufferoffset, alignment) where {T}
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    # adjust buffer offset, make primitive array buffer
    bufferoffset += blen
    blen = len * sizeof(Base.nonmissingtype(T))
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    return bufferoffset + padding(blen, alignment)
end

function writebuffer(io, col::Primitive{T}, alignment) where {T}
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    writebitmap(io, col, alignment)
    n = writearray(io, Base.nonmissingtype(T), col.data)
    @debug 1 "writing array: col = $(typeof(col.data)), n = $n, padded = $(padding(n, alignment))"
    writezeros(io, paddinglength(n, alignment))
    return
end


text:
module FastArithmetic

using Nemo

export monomialToDual

"""
    monomialToDual{T <: RingElem}(a::Array{T,1},P::PolyElem{T})

Convert monomial (canonical) coordinates "a" to dual coordinates
with respect to the trace generated by P.

# Remark
* a and P are over a field k which must be a perfect field ;
* P must be monic and squarefree ;
* k[x]/(P) is not necessarily a field since P is not necessarily irreducible.
"""
function monomialToDual{T <: RingElem}(a::Array{T,1},P::PolyElem{T})

  # We first check if the elements of b and coefficients
  # of P belong to the same field
  k::Nemo.Ring=parent(a[1])
  @assert k==base_ring(P)

  # We set some constants, P is of degree m and belongs to the ring R = k[t]
  m::Int64=degree(P)
  R::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(R)

  # We set Q to 1/rev(P,m+1) mod t^m
  Q::PolyElem{T}=reverse(P,m+1)
  Q=gcdinv(Q,t^m)[2]

  # Finally we compute b = rev(P'a mod P, m)Q mod t^m
  A::PolyElem{T}=R(a)
  b::PolyElem{T}=(reverse((derivative(P)*A)%P,m)*Q)%(t^m)
  return T[coeff(b,i) for i in 0:(m-1)]
end

export monomialToDual_pre
function monomialToDual_pre{T <: RingElem}(a::Array{T,1},P::PolyElem{T},TP::PolyElem{T})

  # We first check if the elements of b and coefficients
  # of P belong to the same field
  k::Nemo.Ring=parent(a[1])
  @assert k==base_ring(P)

  # We set some constants, P is of degree m and belongs to the ring R = k[t]
  m::Int64=degree(P)
  R::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(R)

  # Finally we compute b = rev(P'a mod P, m)TP mod t^m
  A::PolyElem{T}=R(a)
  b::PolyElem{T}=(reverse((derivative(P)*A)%P,m)*TP)%(t^m)
  return T[coeff(b,i) for i in 0:(m-1)]
end

export dualToMonomial

"""
    dualToMonomial{T}(b::Array{T,1},P::PolyElem{T})

Convert dual coordinates "b" to monomial coordinates with respect to the trace
generated by P.

# Remark
* a and P are over a field k which must be a perfect field ;
* P must be monic and squarefree ;
* k[x]/(P) is not necessarily a field since P is not necessarily irreducible.
"""
function dualToMonomial{T}(b::Array{T,1},P::PolyElem{T})

  # We first check if the elements of b and coefficients
  # of P belong to the same field
  k::Nemo.Ring=parent(b[1])
  @assert k==base_ring(P)

  # We set some constants, P is of degree m and belongs to the ring R = k[t]
  m::Int64=degree(P)
  R::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(R)

  # We compute S = 1/P' mod P, c = rev(P, m+1)b, d = cS mod P
  S::PolyElem{T}=gcdinv(derivative(P),P)[2]
  c::PolyElem{T}=(reverse(P,m+1)*R(b))%(t^m)
  c=reverse(c,m)
  d::PolyElem{T}=(c*S)%P
  return T[coeff(d,i) for i in 0:(m-1)]
end

export dualToMonomial_pre
function dualToMonomial_pre{T}(b::Array{T,1},P::PolyElem{T}, S::PolyElem{T})

  # We first check if the elements of b and coefficients
  # of P belong to the same field
  k::Nemo.Ring=parent(b[1])
  @assert k==base_ring(P)

  # We set some constants, P is of degree m and belongs to the ring R = k[t]
  m::Int64=degree(P)
  R::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(R)

  # We compute S = 1/P' mod P, c = rev(P, m+1)b, d = cS mod P
  c::PolyElem{T}=(reverse(P,m+1)*R(b))%(t^m)
  c=reverse(c,m)
  d::PolyElem{T}=(c*S)%P
  return T[coeff(d,i) for i in 0:(m-1)]
end

export mulT

"""
    mulT{T}(c::Array{T,1},P::PolyElem{T},n::Int64)

The tranposition of the naive algorithm of multiplication by P.

# Arguments
* c::Array{T,1} must have length m (degree of P) + n.

# Remark
* closely linked with the middle product [1] ;
* I don't really see why n is an argument, since we could obtain it by
computing n = length(c) - degree(P).

# References
* [1] : G. Hanrot, M. Quercia, and P. Zimmerman. The middle product algorithm I.
Appl. Algebra Engrg. Comm. Comput., 14(6):415-438, 2004.
"""
function mulT{T}(c::Array{T,1},P::PolyElem{T},n::Int64)

  # We first check if the elements of b and coefficients
  # of P belong to the same field
  m::Int64=degree(P)
  k::Nemo.Ring=base_ring(P)
  @assert k==parent(c[1])

  # We compute R = k[t]
  R::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(R)

  # And we transpose all the basic operations
  p::Array{T,1}=T[k(coeff(P,j)) for j in 0:m]
  b::Array{T,1}=Array(T,n+1) # array filled with "#undef"
  for i in 1:(n+1)
    b[i]=k(0)
  end
  for i in range(m+n,-1,m+n+1)
    for j in range(min(m,i),-1,min(m,i)-max(0,i-n)+1)
      b[i-j+1]=b[i-j+1]+p[j+1]*c[i+1]
    end
  end
  return R(T[b[i] for i in 1:(n+1)])
end

export mulTmid

"""
    mulTmid{T}(c::Array{T,1},P::PolyElem{T},n::Int64)

The tranposition of the algorithm of multiplication by P. Using middle product.

# Arguments
* c::Array{T,1} must have length m (degree of P) + n + 1 since
  it represents a polynomial in R_{n+m}[X].

# Remark
* the middle product formula is in [1]
* mulTmid(⋅,P,n) : R_{n+m}[X] ⟶ R_n[X]

# References
* [1] : G. Hanrot, M. Quercia, and P. Zimmerman. The middle product algorithm I.
Appl. Algebra Engrg. Comm. Comput., 14(6):415-438, 2004.
"""
function mulTmid{T}(c::Array{T,1},P::PolyElem{T},n::Int64)

  # We first check if the elements of b and coefficients
  # of P belong to the same ring
  m::Int64=degree(P)
  k::Nemo.Ring=base_ring(P)
  @assert k==parent(c[1])

  # And we compute the middle product of P and c
  R::Nemo.Ring=parent(P)
  C::PolyElem{T}=R(c)
  Q::PolyElem{T}=reverse(P,m+1)
  return shift_right(mullow(Q,C,n+m+1),m)
end

export remT

"""
    remT{T}(r::Array{T,1},P::PolyElem{T},n::Int64)

Transposition of the remainder by P algorithm.

An other linear extension algorithm. Take the r first values of a linear
recurring sequence of minimal polynomial P and compute the n first values.
"""
function remT{T}(r::Array{T,1},P::PolyElem{T},n::Int64)

  # We first check if the elements of b and coefficients
  # of P belong to the same ring
  m::Int64=degree(P)
  k::Nemo.Ring=base_ring(P)
  @assert k==parent(r[1])

  # We compute the ring K = k[t]
  K::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(K)

  # We compute α = 1/rev(P, m+1) mod t^(n-m)
  R::PolyElem{T}=K(T[r[i] for i in 1:m]) # useless creation of a list... memory !
  α::PolyElem{T}=reverse(P,m+1)
  α=gcdinv(α,t^(n-m))[2]

  # We add zeros to match the awaited size
  b::Array{T,1}=copy(r) # copy + push! : cost ?
  while length(b)<n
    push!(b,k(0))
  end

  # We compute the final result
  ans = R-shift_left((mullow(α,mulTmid(b,P,n-m),n-m)),m)
  return T[coeff(ans,j) for j in 0:n-1]
end

export remT_pre
function remT_pre{T}(r::Array{T,1},P::PolyElem{T},n::Int64,TP::PolyElem{T})

  # We first check if the elements of b and coefficients
  # of P belong to the same ring
  m::Int64=degree(P)
  k::Nemo.Ring=base_ring(P)
  @assert k==parent(r[1])

  # We compute the ring K = k[t]
  K::Nemo.Ring=parent(P)
  t::PolyElem{T}=gen(K)

  # We compute α = 1/rev(P, m+1) mod t^(n-m)
  R::PolyElem{T}=K(T[r[i] for i in 1:m]) # useless creation of a list... memory !
  α::PolyElem{T}=TP % t^(n-m)

  # We add zeros to match the awaited size
  b::Array{T,1}=copy(r) # copy + push! : cost ?
  while length(b)<n
    push!(b,k(0))
  end

  # We compute the final result
  ans = R-shift_left((mullow(α,mulTmid(b,P,n-m),n-m)),m)
  return T[coeff(ans,j) for j in 0:n-1]
end

export remTnaif

"""
    remTnaif{T}(r::Array{T,1},P::PolyElem{T},n::Int64)

Linear extension algorithm.

Take the m first elements of a Linear recurring sequence with minimal
polynomial P and compute the n first.
"""
function remTnaif{T}(r::Array{T,1},P::PolyElem{T},n::Int64)
  m::Int64=degree(P)
  p::Array{T,1}=T[coeff(P,j) for j in 0:m]
  b::Array{T,1}=copy(r)
  while length(b)<n
      s=sum([-1*p[j]*b[end-m+j] for j in 1:m])
      push!(b,p[end]^(-1)*s)
  end
  return b
end

export mulModT

"""
    mulModT{T}(P::Array{T,1},Q::PolyElem{T},R::PolyElem{T},n::Int64)

Transposition of the modular multiplication.

# Arguments
* P::Array{T,1} represents a polynomial of degree r - 1, so is of size r (where r=deg(R)) ;
* Q::PolyElem{T} is the fixed polynomial used in mul(⋅,Q) : P ⟼ PQ ;
* R::PolyElem{T} is the polynomial used in the remainder rem(⋅,R) : P ⟼ P mod R ;
* n::Int64 is the degree of the output (at most).

# Remark
* mulModT(⋅,Q,R,n) : R_{r-1}[X] ⟶ R_n[X].
"""
function mulModT{T}(P::Array{T,1},Q::PolyElem{T},R::PolyElem{T},n::Int64)
  q,r=degree(Q),degree(R)
  a = remT(P,R,n+q+1)
  return mulTmid(a,Q,n)
end

export mulModT_pre
function mulModT_pre{T}(P::Array{T,1},Q::PolyElem{T},R::PolyElem{T},n::Int64,TR::PolyElem{T})
  q,r=degree(Q),degree(R)
  a = remT_pre(P,R,n+q+1,TR)
  return mulTmid(a,Q,n)
end

export embed

"""
    embed{T}(b::Array{T,1},P::PolyElem{T},c::Array{T,1},Q::PolyElem{T},r::Int64=0)

Compute the embeding of Π={bc | b ∈ k[x]/(P) , c ∈ k[y]/(Q)} ⊂ k[x,y]/(P,Q) in k[z]/(R).
"""
function embed{T}(b::Array{T,1},P::PolyElem{T},c::Array{T,1},Q::PolyElem{T},r::Int64=0)
  if r==0
    r=length(b)*length(c)
  end
  t::Array{T,1}=remT(b,P,r)
  u::Array{T,1}=remT(c,Q,r)
  return T[t[j]*u[j] for j in 1:r]
end

export berlekampMassey

"""
    berlekampMassey{T <: RingElem}(a::Array{T,1},n::Int64,S=0)

Compute the minimal polynomial of a linear recurring sequence.
"""
function berlekampMassey{T <: RingElem}(a::Array{T,1},n::Int64,S=0)
  if S==0
    k::Nemo.Ring=parent(a[1])
    S::Nemo.Ring,x::PolyElem{T}=PolynomialRing(k,"x")
  else
    x=gen(S)
  end
  m::Int64=2*n-1
  R0::PolyElem{T}=x^(2*n)
  R1::PolyElem{T}=S(reverse(a))
  V0::PolyElem{T}=S(0)
  V1::PolyElem{T}=S(1)
  while n<=degree(R1)
    Q::PolyElem{T},R::PolyElem{T}=divrem(R0,R1)
    V::PolyElem{T}=V0-Q*V1
    V0=V1
    V1=V
    R0=R1
    R1=R
  end
  return V1*lead(V1)^(-1)
end

export computeR

"""
    computeR{T}(P::PolyElem{T},Q::PolyElem{T})

Compute the composed product R of P and Q.
"""
function computeR{T}(P::PolyElem{T},Q::PolyElem{T})
  m::Int64=degree(P)
  n::Int64=degree(Q)
  k::Nemo.Ring=base_ring(P)

  up::Array{T,1}=monomialToDual([k(1)],P)
  uq::Array{T,1}=monomialToDual([k(1)],Q)

  t::Array{T,1}=embed(up,P,uq,Q,2*m*n)

  return berlekampMassey(t,m*n,parent(P))
end

export project

"""
    project{T}(a::Array{T,1},P::PolyElem{T},Q::PolyElem{T})

Compute the section of the embedding k[x]/(P) ⟶ k[z]/(R), where R = P ⊙ Q.
"""
function project{T}(a::Array{T,1},P::PolyElem{T},Q::PolyElem{T})
  n::Int64=degree(Q)
  m::Int64=degree(P)
  c::Array{T,1}=Array(T,n)
  k::Nemo.Ring=base_ring(Q)
  c[1]=k(1)
  for j in 2:n
    c[j]=k(0)
  end
  u::Array{T,1}=remT(c,Q,m*n) # it seems to be the only thing expensive
  K::Nemo.Ring=parent(P)
  d::PolyElem{T}=K([a[j]*u[j] for j in 1:(m*n)])%P
  return T[coeff(d,j) for j in 0:(m-1)]
end

export phi1

"""
    phi1{T}(b::Array{T,2},P::PolyElem{T},Q::PolyElem{T})

Compute the isomorphism Φ : k[x,y]/(P,Q) ⟶ k[z]/(R).

# Arguments
* b::Array{T,2} is a 2d array reprenting (b_ji) where 0 <= i < m and 0 <= j < n, it is
the transpose of the b_ij of the paper because we prefer to extract columns than lines for
types reasons.
"""
function phi1{T}(b::Array{T,2},P::PolyElem{T},Q::PolyElem{T})
  k::Nemo.Ring=parent(b[1,1])
  @assert k==base_ring(Q)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  u::Array{T,1}=remT(monomialToDual(T[k(1)],P),P,m*(n+1)-1)
  a::Array{T,1}=Array(T,m*n)
  for j in 1:m*n
    a[j]=k(0)
  end
  for i in 1:m
    t::Array{T,1}=remT(b[:,i],Q,m*n)
    for j in 1:m*n
      a[j]=a[j]+t[j]*u[i+j-1] # /!\ indices
    end
  end
  return a
end

export phi1_pre
function phi1_pre{T}(b::Array{T,2},P::PolyElem{T},Q::PolyElem{T}, up::Array{T,1})
  k::Nemo.Ring=parent(b[1,1])
  @assert k==base_ring(Q)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  u::Array{T,1}=remT(up,P,m*(n+1)-1)
  a::Array{T,1}=Array(T,m*n)
  for j in 1:m*n
    a[j]=k(0)
  end
  for i in 1:m
    t::Array{T,1}=remT(b[:,i],Q,m*n)
    for j in 1:m*n
      a[j]=a[j]+t[j]*u[i+j-1] # /!\ indices
    end
  end
  return a
end

export inversePhi1

"""
    inversePhi1{T}(a::Array{T,1},P::PolyElem{T},Q::PolyElem{T})

Compute Φ^(-1) : k[z]/(R) ⟶ k[x,y]/(P,Q).

# Output
* b::Array{T,2} is a 2d array reprenting (b_ji) where 0 <= i < m and 0 <= j < n, it is
the transpose of the b_ij of the paper because we prefer to extract columns than lines for
types reasons.
"""
function inversePhi1{T}(a::Array{T,1},P::PolyElem{T},Q::PolyElem{T})
  k::Nemo.Ring=parent(a[1,1])
  @assert k==base_ring(Q)
  K::Nemo.Ring=parent(Q)
  y::PolyElem{T}=gen(K)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  b::Array{T,2}=Array(T,(n,m))
  u::Array{T,1}=remT(monomialToDual(T[k(1)],P),P,m*(n+1)-1)
  for i in m:-1:1
    d::PolyElem{T}=K([a[j]*u[i+j-1] for j in 1:(m*n)])%Q
    for j in 1:n
    b[j,i]=coeff(d,j-1)
    end
  end
  return b
end

export inversePhi1_pre
function inversePhi1_pre{T}(a::Array{T,1},P::PolyElem{T},Q::PolyElem{T}, up::Array{T,1})
  k::Nemo.Ring=parent(a[1,1])
  @assert k==base_ring(Q)
  K::Nemo.Ring=parent(Q)
  y::PolyElem{T}=gen(K)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  b::Array{T,2}=Array(T,(n,m))
  u::Array{T,1}=remT(up,P,m*(n+1)-1)
  for i in m:-1:1
    d::PolyElem{T}=K([a[j]*u[i+j-1] for j in 1:(m*n)])%Q
    for j in 1:n
    b[j,i]=coeff(d,j-1)
    end
  end
  return b
end

export phi2

"""
    phi2{T}(b::Array{T,2},P::PolyElem{T},Q::PolyElem{T})

Compute the isomorphism Φ : k[x,y]/(P,Q) ⟶ k[z]/(R).

# Argument
* This time b::Array{T,2} is the same as the one in the text.
"""
function phi2{T,Y}(b::Array{T,2},P::Y,Q::Y,R::Y)
  k::Nemo.Ring=parent(b[1,1])
  @assert k==base_ring(P)
  K::Nemo.Ring=parent(P)
  z=gen(K)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  N::Int64=n+m-1
  p::Int64=ceil(sqrt(N))
  q::Int64=ceil(N/p)
  y::Array{T,1}=monomialToDual(T[k(0),k(1)],Q)
  up::Array{T,1}=monomialToDual(T[k(1)],P)
  S::Y=K(dualToMonomial(embed(up,P,y,Q),R))
  U::Y=gcdinv(S,R)[2]

  Sprime::Array{Y,1}=Array(Y,q+1)
  Sprime[1]=K(1)

  for i in 2:(q+1)
    Sprime[i]=mulmod(Sprime[i-1],S,R)
  end

  MT::Nemo.Ring=MatrixSpace(K,q,n)
  mt::MatElem=MT()

  for i in 1:q
    c::Array{T,1}=T[coeff(Sprime[i],h) for h in 0:(m*n-1)]
    for j in 1:n
      mt[i,j]=K(c[(j-1)*m+1:j*m])
    end
  end

  MC::Nemo.Ring=MatrixSpace(K,p,q)
  mc::MatElem=MC()

  for i in 1:p
    for j in 1:q
      mc[i,j]=K(T[h+i*q+j-m-2 < 1 ? k(0) : h+i*q+j-m-2 > n ? k(0) : b[h,h+i*q+j-m-2] for h in 1:m])
    end
  end

  Mv::MatElem=mc*mt

  V::Array{Y,1}=Array(Y,p)
  for i in 1:p
    V[i]=sum([shift_left(Mv[i,j],(j-1)*m) for j in 1:n])%R
  end

  a::Y=K()

  for i in p:-1:1
    a=(Sprime[q+1]*a+V[i])%R
  end

  a=mulmod(a,U^(m-1),R)

  return T[coeff(a,i) for i in 0:(m*n-1)]
end

export phi2_pre
function phi2_pre{T,Y}(b::Array{T,2},P::Y,Q::Y,R::Y, up::Array{T,1})
  k::Nemo.Ring=parent(b[1,1])
  @assert k==base_ring(P)
  K::Nemo.Ring=parent(P)
  z=gen(K)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  N::Int64=n+m-1
  p::Int64=ceil(sqrt(N))
  q::Int64=ceil(N/p)
  y::Array{T,1}=monomialToDual(T[k(0),k(1)],Q)
  S::Y=K(dualToMonomial(embed(up,P,y,Q),R))
  U::Y=gcdinv(S,R)[2]

  Sprime::Array{Y,1}=Array(Y,q+1)
  Sprime[1]=K(1)

  for i in 2:(q+1)
    Sprime[i]=mulmod(Sprime[i-1],S,R)
  end

  MT::Nemo.Ring=MatrixSpace(K,q,n)
  mt::MatElem=MT()

  for i in 1:q
    c::Array{T,1}=T[coeff(Sprime[i],h) for h in 0:(m*n-1)]
    for j in 1:n
      mt[i,j]=K(c[(j-1)*m+1:j*m])
    end
  end

  MC::Nemo.Ring=MatrixSpace(K,p,q)
  mc::MatElem=MC()

  for i in 1:p
    for j in 1:q
      mc[i,j]=K(T[h+i*q+j-m-2 < 1 ? k(0) : h+i*q+j-m-2 > n ? k(0) : b[h,h+i*q+j-m-2] for h in 1:m])
    end
  end

  Mv::MatElem=mc*mt

  V::Array{Y,1}=Array(Y,p)
  for i in 1:p
    V[i]=sum([shift_left(Mv[i,j],(j-1)*m) for j in 1:n])%R
  end

  a::Y=K()

  for i in p:-1:1
    a=(Sprime[q+1]*a+V[i])%R
  end

  a=mulmod(a,U^(m-1),R)

  return T[coeff(a,i) for i in 0:(m*n-1)]
end

export inversePhi2

"""
    inversePhi2{T}(a::Array{T,1},P::Y,Q::Y)

Compute Φ^(-1) : k[z]/(R) ⟶ k[x,y]/(P,Q).

# Output
* This time b::Array{T,2} is the same as the one in the text.
"""
function inversePhi2{T,Y}(a::Array{T,1},P::Y,Q::Y,R::Y)
  k::Nemo.Ring=parent(a[1])
  @assert k==base_ring(P)
  K::Nemo.Ring=parent(P)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  N::Int64=n+m-1
  p::Int64=ceil(sqrt(N))
  q::Int64=ceil(N/p)
  y::Array{T,1}=monomialToDual(T[k(0),k(1)],Q)
  up::Array{T,1}=monomialToDual(T[k(1)],P)
  S::Y=K(dualToMonomial(embed(up,P,y,Q),R))
  U::Y=gcdinv(S,R)[2]

  Sprime::Array{Y,1}=Array(Y,q+1)
  Sprime[1]=K(1)

  for i in 2:(q+1)
    Sprime[i]=mulmod(Sprime[i-1],S,R)
  end

  MT::Nemo.Ring=MatrixSpace(K,n,q)
  mt::MatElem=MT()

  for i in 1:q
    c::Array{T,1}=T[coeff(Sprime[i],h) for h in 0:(m*n-1)]
    for j in 1:n
      mt[j,i]=reverse(K(c[(j-1)*m+1:j*m]),m) # reverse in order to do transposed product
    end
  end

  u::Y=powmod(U,m-1,R)
  W::Y=mulTmid(remT(a,R,2*m*n-1),u,m*n-1)
  a=T[coeff(W,j) for j in 0:(m*n-1)]

  V::Array{Array{T,1},1}=Array(Array{T,1},p)
  A::Y=K()

  for i in 1:p
    V[i]=remT(a,R,m*n+m-1)
    A=mulTmid(remT(a,R,2*m*n-1),Sprime[q+1],m*n-1) # That's the thing taking time !!!
    a=T[coeff(A,j) for j in 0:(m*n-1)]
  end

  MV::Nemo.Ring=MatrixSpace(K,p,n)
  mv::MatElem=MV()

  for i in 1:p, j in 1:n
    mv[i,j]=K(V[i][(j-1)*m+1:(j-1)*m+2*m-1])
  end

  mc::MatElem=mv*mt
  cc::Array{Y,1}=Array(Y,p*q)

  for i in 1:p, j in 1:q
    cc[(i-1)*q+j]=shift_right(truncate(mc[i,j],2*m-1),m-1)
  end

  return T[coeff(cc[j-i+m],i-1) for i in 1:m, j in 1:n]
end

export inversePhi2_pre
function inversePhi2_pre{T,Y}(a::Array{T,1},P::Y,Q::Y,R::Y,up::Array{T,1})
  k::Nemo.Ring=parent(a[1])
  @assert k==base_ring(P)
  K::Nemo.Ring=parent(P)
  m::Int64=degree(P)
  n::Int64=degree(Q)
  N::Int64=n+m-1
  p::Int64=ceil(sqrt(N))
  q::Int64=ceil(N/p)
  y::Array{T,1}=monomialToDual(T[k(0),k(1)],Q)
  S::Y=K(dualToMonomial(embed(up,P,y,Q),R))
  U::Y=gcdinv(S,R)[2]

  Sprime::Array{Y,1}=Array(Y,q+1)
  Sprime[1]=K(1)

  for i in 2:(q+1)
    Sprime[i]=mulmod(Sprime[i-1],S,R)
  end

  MT::Nemo.Ring=MatrixSpace(K,n,q)
  mt::MatElem=MT()

  for i in 1:q
    c::Array{T,1}=T[coeff(Sprime[i],h) for h in 0:(m*n-1)]
    for j in 1:n
      mt[j,i]=reverse(K(c[(j-1)*m+1:j*m]),m) # reverse in order to do transposed product
    end
  end

  u::Y=powmod(U,m-1,R)
  W::Y=mulTmid(remT(a,R,2*m*n-1),u,m*n-1)
  a=T[coeff(W,j) for j in 0:(m*n-1)]

  V::Array{Array{T,1},1}=Array(Array{T,1},p)
  A::Y=K()

  for i in 1:p
    V[i]=remT(a,R,m*n+m-1)
    A=mulTmid(remT(a,R,2*m*n-1),Sprime[q+1],m*n-1) # That's the thing taking time !!!
    a=T[coeff(A,j) for j in 0:(m*n-1)]
  end

  MV::Nemo.Ring=MatrixSpace(K,p,n)
  mv::MatElem=MV()

  for i in 1:p, j in 1:n
    mv[i,j]=K(V[i][(j-1)*m+1:(j-1)*m+2*m-1])
  end

  mc::MatElem=mv*mt
  cc::Array{Y,1}=Array(Y,p*q)

  for i in 1:p, j in 1:q
    cc[(i-1)*q+j]=shift_right(truncate(mc[i,j],2*m-1),m-1)
  end

  return T[coeff(cc[j-i+m],i-1) for i in 1:m, j in 1:n]
end

println("FastArithmetic comes with even less warranty\n")

end

max_stars_repo_name: artichaut/fastarithmetic

max_issues_repo_name: artichaut/fastarithmetic

max_forks_repo_name: artichaut/fastarithmetic