text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbcoons(u1::NURBS1D{I,F},u2::NURBS1D{I,F},
             v1::NURBS1D{I,F},v2::NURBS1D{I,F}
            )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}

Construction of a Coons patch.

Construction of a bilinearly blended Coons surface patch from four NURBS
curves that define the boundary.

The orientation of the four NURBS boundary curves.

       ^ V direction
       |
       |      u2
       ------->--------
       |              |
       |              |
    v1 ^   Surface    ^ v2
       |              |
       |              |
       ------->-----------> U direction
              u1

# Arguments:
- `u1`: NURBS curve defining the bottom U direction boundary of the constructed
  NURBS surface.
- `u2`: NURBS curve defining the top U direction boundary of the constructed
  NURBS surface.
- `v1`: NURBS curve defining the bottom V direction boundary of the constructed
  NURBS surface.
- `v1`: NURBS curve defining the top V direction boundary of the constructed
  NURBS surface.

# Output:
- `srf`: Coons NURBS surface patch.

# Examples:
```julia
julia> srf = nrbcoons(u1, u2, v1, v2)
```

Define four NURBS curves and construct a Coons surface patch.
```julia
julia> pnts = [0.0  3.0  4.5  6.5 8.0 10.0;
               0.0  0.0  0.0  0.0 0.0  0.0;
               2.0  2.0  7.0  4.0 7.0  9.0];
julia> crv1 = nrbmak(pnts, vec([0 0 0 1/3 0.5 2/3 1 1 1]));
julia> pnts= [ 0.0  3.0  5.0  8.0 10.0;
              10.0 10.0 10.0 10.0 10.0;
               3.0  5.0  8.0  6.0 10.0];
julia> crv2 = nrbmak(pnts, vec([0 0 0 1/3 2/3 1 1 1]));
julia> pnts= [0.0 0.0 0.0  0.0;
              0.0 3.0 8.0 10.0;
              2.0 0.0 5.0  3.0];
julia> crv3 = nrbmak(pnts, vec([0 0 0 0.5 1 1 1]));
julia> pnts= [10.0 10.0 10.0 10.0 10.0;
               0.0   3.0  5.0  8.0 10.0;
               9.0   7.0  7.0 10.0 10.0];
julia> crv4 = nrbmak(pnts, vec([0 0 0 0.25 0.75 1 1 1]));
julia> srf = nrbcoons(crv1, crv2, crv3, crv4);
julia> nrbplot(srf,[20;20]);
```

"""
function nrbcoons(u1::NURBS1D{I,F},u2::NURBS1D{I,F},
                  v1::NURBS1D{I,F},v2::NURBS1D{I,F}
                 )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}

r1 = nrbruled(u1, u2);
r2 = nrbtransp(nrbruled(v1, v2));
t  = nrb4surf(u1.coefs[:,1], u1.coefs[:,end], u2.coefs[:,1], u2.coefs[:,end]);

# raise all surfaces to a common degree
du = max(r1.order[1], r2.order[1], t.order[1]);
dv = max(r1.order[2], r2.order[2], t.order[2]);
r1 = nrbdegelev(r1, [du - r1.order[1]; dv - r1.order[2]]);
r2 = nrbdegelev(r2, [du - r2.order[1]; dv - r2.order[2]]);
t  = nrbdegelev(t,  [du -  t.order[1]; dv -  t.order[2]]);

# merge the knot vectors, to obtain a common knot vector

# U knots
k1 = r1.knots[1];
k2 = r2.knots[1];
k3 =  t.knots[1];
k = unique([k1;k2;k3]);
n = length(k);
# TODO this is bad, increasing the size without allocating
kua = Vector{F}();
kub = Vector{F}();
kuc = Vector{F}();
for i in 1:n
  i1 = sum(x-> x == k[i],k1);
  i2 = sum(x-> x == k[i],k2);
  i3 = sum(x-> x == k[i],k3);
  m = max(i1, i2, i3);
  if m-i1>0; append!(kua, fill(k[i],m-i1)); end #if
  if m-i2>0; append!(kub, fill(k[i],m-i2)); end #if
  if m-i3>0; append!(kuc, fill(k[i],m-i3)); end #if
end # for i

# U knots
k1 = r1.knots[2];
k2 = r2.knots[2];
k3 =  t.knots[2];
k = unique([k1;k2;k3]);
n = length(k);
# TODO this is bad, increasing the size without allocating
kva = Vector{F}();
kvb = Vector{F}();
kvc = Vector{F}();
for i in 1:n
  i1 = sum(x-> x == k[i],k1);
  i2 = sum(x-> x == k[i],k2);
  i3 = sum(x-> x == k[i],k3);
  m = max(i1, i2, i3);
  if m-i1>0; append!(kva, fill(k[i],m-i1)); end #if
  if m-i2>0; append!(kvb, fill(k[i],m-i2)); end #if
  if m-i3>0; append!(kvc, fill(k[i],m-i3)); end #if
end # for i

if !(isempty(kua)) || !(isempty(kva)); r1 = nrbkntins(r1, [kua,kva]); end # if
if !(isempty(kub)) || !(isempty(kvb)); r2 = nrbkntins(r2, [kub,kvb]); end # if
if !(isempty(kuc)) || !(isempty(kvc)); t  = nrbkntins(t , [kuc,kvc]); end # if

coefs = r1.coefs + r2.coefs - t.coefs;

return nrbmak(coefs, r1.knots);
end # nrbcoons


text:
#=
   My faster rewrites of `min`, `max`, `minmax`, `maxmin`.

   These are processing isomorphs for all normal uses.
   Julia's NaN, NaN32 are handled properly.
   Negatively signed NaNs, which cannot appear unless
   intentionally generated, are not IEEE754 compliant.

   For use with `Float32`, `Float64`.
=#

const FastFloat = Union{Float32, Float64}

function fastermin(a::T, b::T) where {T<:FastFloat}
    a_b = a - b
    (signbit(a_b) || isnan(a)) ? a : b
end

function fastermax(a::T, b::T) where {T<:FastFloat}
    b_a = b - a
    (signbit(b_a) || isnan(a)) ? a : b
end

fasterminmax(a::T, b::T) where {T<:FastFloat} =
    (min(a, b), max(a, b))

fastermaxmin(a::T, b::T) where {T<:FastFloat} =
    (max(a, b), min(a, b))

#=
   their utilzation
=#

min(x::Robust32, y::Robust32) = Rob32(fastermin(value64(x), value64(y)))
max(x::Robust32, y::Robust32) = Rob32(fastermax(value64(x), value64(y)))
minmax(x::Robust32, y::Robust32) = Rob32(fasterminmax(value64(x), value64(y)))
maxmin(x::Robust32, y::Robust32) = Rob32(fastmaxmin(value64(x), value64(y)))

min(x::Robust32, y::Float64) = Rob32(fastermin(value64(x), y))
max(x::Robust32, y::Float64) = Rob32(fastermax(value64(x), y))
minmax(x::Robust32, y::Float64) = Rob32(fasterminmax(value64(x), y))
maxmin(x::Robust32, y::Float64) = Rob32(fasterminmax(value64(x), y))

min(x::Float64, y::Robust32) = Rob32(fastermin(x, value64(y)))
max(x::Float64, y::Robust32) = Rob32(fastermax(x, value64(y)))
minmax(x::Float64, y::Robust32) = Rob32(fasterminmax(x, value64(y)))
maxmin(x::Float64, y::Robust32) = Rob32(fasterminmax(x, value64(y)))

min(x::Robust32, y::Real) = Rob32(fastermin(value64(x), Float64(y)))
max(x::Robust32, y::Real) = Rob32(fastermax(value64(x), Float64(y)))
minmax(x::Robust32, y::Real) = Rob32(fasterminmax(value64(x), Float64(y)))
maxmin(x::Robust32, y::Real) = Rob32(fasterminmax(value64(x), Float64(y)))

min(x::Real, y::Robust32) = Rob32(fastermin(Float64(x), value64(y)))
max(x::Real, y::Robust32) = Rob32(fastermax(Float64(x), value64(y)))
minmax(x::Real, y::Robust32) = Rob32(fasterminmax(Float64(x), value64(y)))
maxmin(x::Real, y::Robust32) = Rob32(fasterminmax(Float64(x), value64(y)))

ext: jl


text:
#!/usr/bin/julia

# A simple implementation of the continued fraction factorization method.
# (variation of the Brillhart-Morrison algorithm).

# See also:
#   https://en.wikipedia.org/wiki/Pell%27s_equation
#   https://en.wikipedia.org/wiki/Continued_fraction_factorization
#   https://trizenx.blogspot.com/2018/10/continued-fraction-factorization-method.html

# "Gaussian elimination" algorithm from:
#    https://github.com/martani/Quadratic-Sieve

using Primes

const BIG_ONE  = big"1"
const BIG_ZERO = big"0"

function issquare(n)
    isqrt(n)^2 == n
end

function issquarefree(n)
    for (p,e) in factor(n)
        (e >= 2) && return false
    end
    return true
end

function next_multiplier(n, k)

    k += 2

    while (!issquarefree(k) || gcd(k,n) != 1)
        k += 1
    end

    return k
end

function gaussian_elimination(A, n)

    m = length(A)
    I = [BIG_ONE << k for k in 0:m]

    nrow = 0

    for col in (1 : min(m, n))

        npivot = 0

        for row in (nrow+1 : m)
            if (((A[row] >> (col-1)) & 1) == 1)
                npivot = row
                nrow += 1
                break
            end
        end

        npivot == 0 && continue

        if (npivot != nrow)
            A[[npivot,nrow]] = A[[nrow, npivot]]
            I[[npivot,nrow]] = I[[nrow, npivot]]
        end

        for row in (nrow+1 : m)
            if (((A[row] >> (col-1)) & 1) == 1)
                A[row] = xor(A[row], A[nrow])
                I[row] = xor(I[row], I[nrow])
            end
        end
    end

    return I
end

function check_factor(n, g, factors)

    while (rem(n,g) == 0)

        n = div(n,g)
        push!(factors, g)

        if (isprime(n))
            push!(factors, n)
            return 1
        end
    end

    return n
end

function is_smooth_over_prod(n, k)

    g = gcd(n, k)

    while (g > 1)
        n = div(n, g)
        while (rem(n, g) == 0)
            n = div(n, g)
        end
        n == 1 && return true
        g = gcd(n, g)
    end

    n == 1
end

function jacobi(a, n)
    a %= n
    result = 1
    while a != 0
        while iseven(a)
            a >>= 1
            ((n % 8) in [3, 5]) && (result *= -1)
        end
        a, n = n, a
        (a % 4 == n % 4 == 3) && (result *= -1)
        a %= n
    end
    return n == 1 ? result : 0
end

function cffm(n, multiplier = 1)

    n <= 1     && return []
    isprime(n) && return [n]

    if (iseven(n))

        v = 0
        while (iseven(n))
            v += 1
            n >>= 1
        end

        arr1 = [big"2" for i in 1:v]
        arr2 = cffm(n)

        return append!(arr1, arr2)
    end

    if (issquare(n))
        f = cffm(isqrt(n))
        append!(f, f)
        return sort(f)
    end

    N = n*multiplier

    x = isqrt(N)
    y = x
    z = 1
    w = x+x
    r = w

    nf = round(Int64, exp(sqrt(log(n) * log(log(n))))^(sqrt(2) / 4) / 1.25)
    factor_base = [2]

    begin
        p = 3
        while (length(factor_base) < nf)
            if (jacobi(N, p) >= 0)
                push!(factor_base, p)
            end
            p = nextprime(p+1)
        end
    end

    # ~ B = round(Int64, exp(sqrt(log(n) * log(log(n))) / 2))
    # ~ factor_base = []

    # ~ for p in (primes(B))
        # ~ if (jacobi(N,p) >= 0)
            # ~ push!(factor_base, p)
        # ~ end
    # ~ end

    factor_prod  = prod([big(k) for k in factor_base])
    factor_index = Dict{Int64, Int64}()

    for k in (1:length(factor_base))
        factor_index[factor_base[k]] = k-1
    end

    function exponent_signature(factors)
        sig = BIG_ZERO

        for (p,e) in factors
            if (isodd(e))
                sig |= (BIG_ONE << factor_index[p])
            end
        end

        return sig
    end

    L = length(factor_base)+1

    Q = []
    A = []

    (f1, f2) = (BIG_ONE, x)

    while (length(A) < L)

        y = (r*z - y)
        z = div((N - y*y), z)
        r = div((x + y), z)

        (f1, f2) = (f2, rem((r*f2 + f1), n))

        if (issquare(z))
            g = gcd(f1 - isqrt(z), n)

            if (g > 1 && g < n)
                arr1 = cffm(g)
                arr2 = cffm(div(n,g))
                return sort(append!(arr1, arr2))
            end
        end

        if (z > 1 && is_smooth_over_prod(z, factor_prod))
            push!(A, exponent_signature(factor(z)))
            push!(Q, [f1, z])
        end

        if (z == 1)
            println("z == 1 -> trying again with a multiplier...")
            return cffm(n, next_multiplier(n, multiplier))
        end
    end

    while (length(A) < L)
        push!(A, 0)
    end

    I = gaussian_elimination(A, length(A))

    LR = 0
    for k in (length(A):-1:1)
        if (A[k] != 0)
            LR = k+1
            break
        end
    end

    remainder = n
    factors   = []

    function cfrac_find_factors(solution)

        X = BIG_ONE
        Y = BIG_ONE

        for i in 0:length(Q)-1
            if (((solution >> i) & 1) == 1)

                X *= Q[i+1][1]
                Y *= Q[i+1][2]

                g = gcd(X - isqrt(Y), remainder)

                if (g > 1 && g < remainder)

                    remainder = check_factor(remainder, g, factors)

                    if (remainder == 1)
                        return true
                    end
                end
            end
        end

        return false
    end

    for k in LR:length(I)
        cfrac_find_factors(I[k]) && break
    end

    final_factors = []

    for f in (factors)
        if (isprime(f))
            push!(final_factors, f)
        else
            append!(final_factors, cffm(f))
        end
    end

    if (remainder != 1)
        if (remainder != n)
            append!(final_factors, cffm(remainder))
        else
            push!(final_factors, remainder)
        end
    end

    # Failed to factorize n (try again with a multiplier)
    if (remainder == n)
        println("Trying again with a multiplier...")
        return cffm(n, next_multiplier(n, multiplier))
    end

    # Return all prime factors of n
    return sort(final_factors)
end

if (length(ARGS) >= 1)
    for n in (ARGS)
        println(n, " = ", cffm(parse(BigInt, n)))
    end
else
    @time println("2^128 + 1 = ", cffm((big"2")^128 + 1))
end

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
"""
    ReadonlyTiffDiskArray(mappedtype, rawtype, ifds, dims) -> ReadonlyTiffDiskArray

A lazy representation of a OMETIFF file. This custom type is needed since TIFF
files are laid out noncontiguously and nonregularly. It uses an internal index
to determine the mapping from indices to the locations of data slices on disk.
These slices are generally XY slices and are usually loaded in all at once so it
is quickly loaded into an internal cache to speed up the process. Externally,
this type should behave very similarly to an in-memory array, albeit with a
higher cost of accessing an element.

$(FIELDS)
"""
mutable struct ReadonlyTiffDiskArray{T <: Gray, R, N1, N2} <: AbstractArray{T, N2}
    """
    A map of dimensions (sans XY) to the corresponding [`IFD`](@ref)
    """
    ifds::OrderedDict{NTuple{N1, Int}, IFD}

    """
    The full set of dimensions of the TIFF file, including XY
    """
    dims::NTuple{N2, Int}

    """
    An internal cache to fill when reading from disk
    """
    cache::Array{R, 2}

    """
    The dimension indices corresponding to the slice currently in the cache
    """
    cache_index::NTuple{N1, Int}

    function ReadonlyTiffDiskArray(::Type{T}, ::Type{R}, ifds::OrderedDict{NTuple{N1, Int}, IFD}, dims::NTuple{N2, Int}) where {T, R, N1, N2}
        if N2 - 2 != N1
            error("$N2 dimensions given, but the IFDs are indexed on $N1 dimensions instead of "*
                  "expected $(N2-2).")
        end
        new{T, R, N1, N2}(ifds, dims, Array{R}(undef, dims[1], dims[2]), (-1, -1, -1, -1))
    end
end

Base.size(A::ReadonlyTiffDiskArray) = A.dims

function Base.getindex(A::ReadonlyTiffDiskArray{Gray{T}, R, N1, N2}, i1::Int, i2::Int, i::Vararg{Int, N1}) where {T, R, N1, N2}
    # check the loaded cache is already the correct slice
    if A.cache_index == i
        return Gray(reinterpret(T, A.cache[i2, i1]))
    end

    ifd = A.ifds[i]

    # if the file isn't open, lets open a handle and update it
    if !isopen(ifd.file.io)
        path = ifd.file.filepath
        ifd.file.io = getstream(open(path))
    end

    n_strips = length(ifd.strip_offsets)
    strip_len = floor(Int, (size(A.cache, 1) * size(A.cache, 2)) / n_strips)

    # if the data is striped then we need to change the buffer shape so that we
    # can read into it. This should be replaced with a view of cache in Julia
    # >1.4, see https://github.com/JuliaLang/julia/pull/33046
    if n_strips > 1 && size(tmp) != (strip_len, )
        tmp = Array{R}(undef, strip_len)
    else
        tmp = A.cache
    end

    _read_ifd_data!(ifd, A.cache, tmp)

    A.cache_index = i

    return Gray(reinterpret(T, A.cache[i2, i1]))
end

function Base.setindex!(A::ReadonlyTiffDiskArray{Gray{T}, R, N1, N2}, X, I...) where {T, R, N1, N2}
    error("This array is on disk and is read only. Convert to a mutable in-memory version by running "*
          "`copy(arr)`. \n\n𝗡𝗼𝘁𝗲: For large files this can be quite expensive. A future PR will add "*
          "support for reading and writing to/from disk.")
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_m8_opt_rho13_5(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_m8_opt_rho13_5!(A_copy)
end

@inline function exp_mono_m8_opt_rho13_5!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 Ba9 Bb9 B9 T2k11
    # Computing B2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.003584018414055765
    coeff2=0.16376260367565243
    coeff3=0.04917059370539749
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.30528232877254763
    coeff2=0.23746768414771152
    coeff3=0.0023273769915510135
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots3,memslots4)
    # Deallocating Ba3 in slot 3
    # Deallocating Bb3 in slot 4
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0005515398938947843
    coeff2=0.02533967378579724
    coeff3=0.017067262331527198
    coeff4=0.9638185168573307
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.008872840493334803
    coeff2=0.25571769129339744
    coeff3=0.0037256194747560907
    coeff4=0.002052240135641998
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots3,memslots4)
    # Deallocating Ba4 in slot 3
    # Deallocating Bb4 in slot 4
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0012630062267903097
    coeff2=-0.0010962696861478022
    coeff3=0.00022839563887310042
    coeff4=0.013673005502563394
    coeff5=0.8216219421030205
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.8293727566445452
    coeff2=0.1675033881979259
    coeff3=0.0030839993650848517
    coeff4=0.0031897875289976277
    coeff5=3.234759791150232e-5
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots3,memslots4)
    # Deallocating Ba5 in slot 3
    # Deallocating Bb5 in slot 4
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=4.437908455542293e-5
    coeff2=0.0003258145596094315
    coeff3=0.00022081927999800417
    coeff4=-0.13039438638797307
    coeff5=0.8475672106088131
    coeff6=0.8145974867536132
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.10275230374603168
    coeff2=0.25648676493721656
    coeff3=0.005628536962635935
    coeff4=0.008329230127484285
    coeff5=0.00021583757875408033
    coeff6=-1.0968604718569746e-7
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots3,memslots4)
    # Deallocating Ba6 in slot 3
    # Deallocating Bb6 in slot 4
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-2.8201037195155804e-6
    coeff2=4.8728643050646066e-6
    coeff3=2.2723646529748017e-6
    coeff4=-0.0061294802470251075
    coeff5=0.0688513777099278
    coeff6=-0.12038470670165373
    coeff7=0.9948034699390058
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.004650329563627054
    coeff2=0.2598813154789629
    coeff3=0.006290014251994485
    coeff4=0.010576806709956765
    coeff5=0.0003949349245182783
    coeff6=8.856309878298438e-6
    coeff7=7.553245375703866e-9
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots3,memslots4)
    # Deallocating Ba7 in slot 3
    # Deallocating Bb7 in slot 4
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=1.0002604649982476
    coeff2=0.2422828857323093
    coeff3=0.022001533431922154
    coeff4=0.15933029209747282
    coeff5=0.03184541272812379
    coeff6=0.007635550247315368
    coeff7=0.0008615033049494069
    coeff8=3.134369854398072e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=1.0002604649982476
    coeff2=0.2422828857323093
    coeff3=0.022001533431922154
    coeff4=0.15933029209747282
    coeff5=0.03184541272812379
    coeff6=0.007635550247315368
    coeff7=0.0008615033049494069
    coeff8=3.134369854398072e-5
    memslots4 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots4,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots3,memslots4)
    # Deallocating Ba8 in slot 3
    # Deallocating Bb8 in slot 4
    # Computing Ba9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.1346603332119568
    coeff2=-5.9198518169123926e-5
    coeff3=-8.660244928553117e-6
    coeff4=-8.731780915688725e-5
    coeff5=-2.1711606510903095e-5
    coeff6=-5.9723490061415135e-6
    coeff7=-7.70252692547925e-7
    coeff8=-3.1855167508881976e-8
    coeff9=1.0081137407871135
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing Bb9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=-0.06540541282818779
    coeff2=6.846710342127941e-5
    coeff3=2.738382227960167e-6
    coeff4=-6.3730674335546956e-6
    coeff5=-5.739005263969626e-6
    coeff6=-2.157568588204766e-6
    coeff7=-3.3636703984819904e-7
    coeff8=-1.5592924041261544e-8
    coeff9=0.9946441317159572
    # Smart lincomb recycle B7
    memslots9 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots9,true,I*coeff1,true,true)
    # Computing B9 with operation: mult
    mul!(memslots4,memslots3,memslots9)
    # Deallocating Ba9 in slot 3
    # Deallocating Bb9 in slot 9
    # Computing T2k11 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B8+x*B9
    coeff1=0.008776543046053789
    coeff2=-1.3078408314728238e-5
    coeff3=-9.364829337290504e-7
    coeff4=-4.881187067247219e-6
    coeff5=-6.561388347938557e-7
    coeff6=-1.0279110807751132e-7
    coeff7=-5.458270135394687e-9
    coeff8=-0.06726697857924235
    coeff9=0.9964850858948753
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots10 .+ coeff9.*memslots4
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 2
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B8 in slot 10
    # Deallocating B9 in slot 4
    return memslots1 # Returning T2k11
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using ApproximateVI
using Test

include("testentropy.jl")

@testset "ApproximateVI.jl" begin
    testentropy()
end


text:
include("maximin_core.jl")
include("cc_maximin_core.jl")
include("robust_maximin_core.jl")
include("minimax_core.jl")


text:
@discord_object struct AuditLog
    webhooks::Vector{Webhook}
    users::Vector{User}
    audit_log_entries::Vector{AuditLogEntry}
    integrations::Vector{Integration}
end


text:
# # Database Interface
# This file describes the methods that should or can be implemented
# for subtypes of `AbstractDB`.

# First of all, we make it so that any `AbstractResult` is broadcasted wholly:
Broadcast.broadcastable( db :: AbstractDB ) = Ref( db );

# ## Mandatory Methods

# A database should be constructed using the `init_db` method:
"Constructor for empty database of type `T`."
function init_db( :: T, ::Type{ <: AbstractResult },
    :: Type{<:NothingOrSaveable}) :: T where T<:Type{<:AbstractDB}
    nothing
end

# Internally the variables might be scaled (or *transformed*).
# The `is_transformed` method gives an indication to whether or
# not the site vectors of the stored results are transformed or not.
# `set_transformed!` is used to set the flag.
"Bool indicating if the database data been transformed."
is_transformed( :: AbstractDB ) :: Bool = false
"Set the flag indicating whether the database data has been transformed or not."
set_transformed!( :: AbstractDB, :: Bool ) :: Nothing = nothing

# The results are indexed with integers and `get_ids` should
# return a vector or iterator of all result ids:
"List of all `id :: Int` belonging to the stored results."
get_ids( db :: AbstractDB  ) = Int[]

# An `id` can then be used to retrieve a result from the database:
"Get result with `id` from database `db`."
function get_result( db :: AbstractDB{R,I}, id :: Int ) :: R where {R <: AbstractResult,I}
    return R()
end

# When a new result is added to the database, `next_id` is called to get
# its id:
"Return an id for the next result to be added to `db`."
function next_id( db :: AbstractDB ) :: Int
	return -1
end

"Add result `res` to database `db`."
_add_result!(db :: AbstractDB, res :: AbstractResult) = nothing

# There is only one simple method to put an `AbstractIterSaveable` into a database:
"""
    stamp!(db, ids)

Put the saveable `ids` into the database `db`.
"""
function stamp!( db :: AbstractDB, ids :: NothingOrSaveable) :: Nothing
	return nothing
end

# ## Derived methods

# ### Getters
# There are getters for the types …
"Return type of results stored in database."
get_res_type( :: AbstractDB{R,I} ) where {R,I} = R
"Return type of `AbstractIterSaveable`s stored in database."
get_saveable_type( :: AbstractDB{R,I}) where {R,I} = I

# … and for the values:
"Return the evaluation site vector for result with `id` in database `db`."
get_value( db :: AbstractDB, id :: Int ) = get_value( get_result( db, id ) )
"Return the evaluation value vector for result with `id` in database `db`."
get_site( db :: AbstractDB, id :: Int ) = get_site( get_result( db, id) )

"Return a vector of all evaluation site vectors stored in database."
function get_sites( db :: AbstractDB)
	return [ get_site( get_result( db, id ) ) for id = get_ids( db ) ]
end

"Return a vector of all evaluation value vectors stored in database."
function get_values( db :: AbstractDB  )
	return [ get_value( get_result( db, id ) ) for id = get_ids( db ) ]
end

# For getting the number of results in a database we owerwrite
# `Base.length`.
"Number of entries in database."
Base.length( db :: AbstractDB ) :: Int = length(get_sites(db))

# There is a default implmentation to get the ids of results in
# a database that don't have valid value vectors.
# I recommend to owerwrite this:
"Return `true` if the result with `id` in `db` has a valid evaluation vector."
function get_evaluated_flag( db, id ) :: Bool
    res = get_result(db, id)
    has_valid_value( res ) && return true
    return false
end

"Return vector of ids of database `db` that are not evaluated yet."
function _missing_ids( db :: AbstractDB )
    missing_ids = Int[]
    for id = get_ids( db )
        if !get_evaluated_flag(db, id)
            push!(missing_ids, id)
        end
    end
    return _missing_ids
end

# ### Setters

# `new_result!` is meant to construct a new result and return its data base index.
# It is implemented so that if `id` is provided as a positive integer, the new result
# has that id. Else, `next_id` should be called internally.
"Add a new result to the database, return its id of type Int."
function new_result!( db :: AbstractDB{R,I}, x :: Vec, y :: Vec, id :: Int = - 1 ) where{R,I}
    new_id = id < 0 ? next_id(db) : id
	new_result = init_res( R, x, y, new_id )
	push!(db.res, new_result)
	if !has_valid_value(new_result)
		set_evaluated_flag!( db, new_id, false )
	end
	db.num_entries += 1
	return new_id
end

# These setters are used in (un)transforming the database.
# They are based on the setters for `AbstractResult`:
"Set site of result with `id` in database `db` to `x`."
function set_site!(db, id, x) :: Nothing
    set_site!( get_result(db,id), x )
end

"Set value of result with `id` in database `db` to `x`."
function set_value!(db, id, y) :: Nothing
    set_value!( get_result(db,id), y )
end

# If you overwrite `get_evaluated_flag` you might want to
# overwrite `set_evaluated_flag!` too:
"Set the evaluation status for result with `id` to `state`."
set_evaluated_flag!( db :: AbstractDB, id :: Int, state = true) = nothing

# ### Miscellaneous

# To find a result by its values we have:
"""
    find_result(db, x, y)
Return id of a result in `db` that has site `x` and value `y` or return -1
if there is no such result.
"""
function find_result( db :: AbstractDB, x :: Vec, y :: Vec  ) :: Int
    for id ∈ get_ids( db )
        if get_site( db, id ) == x && get_value( db, id ) == y
            return id
        end
    end
    return -1
end

# The above function is utilized in `ensure_contains_values!`:
"Return id of result in `db` with site `x` and values `y`. Create if necessary."
function ensure_contains_values!( db :: AbstractDB, x :: Vec, y :: Vec ) :: Int
    x_pos = find_result(db, x,y);
    if x_pos < 0
        x_pos = new_result!(db, x, y);
    end
    return x_pos
end

# The `eval_missing` method is important for the new two stage model
# construction process and called after "preparing" the models for
# updates but before calling the `update_model` methods:

"Evaluate all unevaluated results in `db` using objectives of `mop`."
function eval_missing!( db :: AbstractDB, mop :: AbstractMOP ) :: Nothing

    missing_ids = _missing_ids(db)

    @logmsg loglevel2 "Performing $(length(missing_ids)) objective evaluations into the database."
    ## evaluate everything in one go to exploit parallelism
    eval_sites = [ get_site( db, id ) for id in missing_ids ]
    eval_values = eval_all_objectives.(mop, eval_sites)

    @assert length(eval_sites) == length(eval_values) == length(missing_ids) "Number of evaluation results does not match."
    for (i,id) in enumerate(missing_ids)
        set_value!( db, id, eval_values[i] )
    end

    for id in missing_ids
        set_evaluated_flag!(db, id)
    end
    return nothing
end

#src NOTE the modifying methods require the `get_site` and `get_value`
#src to return mutable arrays for result. Can we ensure this?

# Internally the variables might be scaled.
# For conversion we offer the `scale!` and `unscale!` defaults:
"Scale the site of result with `id` in database `db` using bounds of `mop`."
function scale!( db :: AbstractDB, id :: Int, mop :: AbstractMOP ) :: Nothing
    set_site!( db, id, scale( get_site( db, id ), mop ) )
    return nothing
end

"Unscale the site of result with `id` in database `db` using bounds of `mop`."
function unscale!( db :: AbstractDB, id :: Int, mop :: AbstractMOP ) :: Nothing
    set_site!( db, id, unscale( get_site( db, id), mop ) )
    return nothing
end

# Also, the value vectors might be resorted internally.
# This is taken care of by the following methods:

"Apply internal objective sorting to result with `id` in `db`."
function apply_internal_sorting!( db :: AbstractDB, id :: Int, mop :: AbstractMOP ) :: Nothing
	set_value!( db, id, apply_internal_sorting( get_value(db, id), mop ) )
    return nothing
end

"Reverse internal sorting of objectives for the result with `id` in `db`."
function reverse_internal_sorting!( db :: AbstractDB, id :: Int, mop :: AbstractMOP ) :: Nothing
	set_value!( db, id, reverse_internal_sorting( get_value(db, id), mop ) )
    nothing
end

# Both, variable scaling and objective sorting, is combined in the
# `(un)transform!` methods:

"Apply scaling and objectives sorting to each result in database `db`."
function transform!( db :: AbstractDB, mop :: AbstractMOP ) :: Nothing
    if !is_transformed(db)
        for id in get_ids(db)
            scale!( db, id, mop)
            apply_internal_sorting!( db, id, mop )
        end
        set_transformed!(db, true)
    end
    nothing
end

"Undo scaling and objectives sorting to each result in database `db`."
function untransform!( db :: AbstractDB, mop :: AbstractMOP ) :: Nothing
    if is_transformed(db)
        for id in get_ids(db)
            unscale!( db, id, mop)
            reverse_internal_sorting!( db, id, mop )
        end
        set_transformed!(db, false)
    end
    nothing
end

# We have a generic copy function that returns a new database containing
# results from the old database:
"""
    copy_db( old_db, result_type, saveable_type )

Return a new database of same 'base' type but possibly with different result and
saveable type.
"""
function copy_db( old_db :: DBT; res_type = Nothing, saveable_type :: Type = Nothing ) where DBT <: AbstractDB
    try
        base_type = @eval $(DBT.name.name)
        _res_type = res_type <: Nothing ? get_res_type( old_db ) : res_type
        _saveable_type = saveable_type <: Nothing ? get_saveable_type( old_db ) : saveable_type
        new_db = init_db(base_type, _res_type, _saveable_type )
        for id = get_ids(old_db)
            res = get_result( old_db, id )
            new_result!(new_db, get_site(res), get_value(res) )
        end
        @logmsg loglevel2 "Copied database with new saveable type."
        return new_db
    catch e
        @error "Failed to copy database with new saveable type." exception=(e, catch_backtrace())
        return old_db
    end
end

# Finally, this little helper returns ids of database results that
# conform to some variable bounds and is used in the model construction:
"Return indices of results in `db` that lie in a box with corners `lb` and `ub`."
function results_in_box_indices(db, lb, ub, exclude_indices = Int[] )
	return [ id for id = get_ids( db ) if
		id ∉ exclude_indices && all(lb .<= get_site(db,id) .<= ub ) ]
end


text:
export cache

## CachedOperator

mutable struct CachedOperator{T<:Number,DM<:AbstractMatrix,M<:Operator,DS,RS,BI} <: Operator{T}
    op::M
    data::DM
    datasize::Tuple{Int,Int}
    domainspace::DS
    rangespace::RS
    bandwidths::BI
    padding::Bool   # records whether the operator should be padded for QR operations
                    #  note that BandedMatrix/AlmostBandedMatrix don't need this
                    #  as the padding is done in the bandwidths
end

CachedOperator(op::Operator,data::AbstractMatrix,sz::Tuple{Int,Int},ds,rs,bi,pd=false) =
    CachedOperator{eltype(data),typeof(data),typeof(op),
                    typeof(ds),typeof(rs),typeof(bi)}(op,data,sz,ds,rs,bi,pd)

CachedOperator(op::Operator,data::AbstractMatrix,sz::Tuple{Int,Int},pd=false) =
    CachedOperator(op,data,sz,domainspace(op),rangespace(op),bandwidths(op),pd)
CachedOperator(op::Operator,data::AbstractMatrix,padding=false) = CachedOperator(op,data,size(data),padding)

function default_CachedOperator(op::Operator;padding::Bool=false)
    if isbanded(op)
        CachedOperator(BandedMatrix,op;padding=padding)
    elseif isbandedblockbanded(op) && !padding
        CachedOperator(BandedBlockBandedMatrix,op)
    elseif isblockbanded(op)
        CachedOperator(BlockBandedMatrix,op;padding=padding)
    elseif israggedbelow(op)
        CachedOperator(RaggedMatrix,op;padding=padding)
    else
        CachedOperator(Matrix,op;padding=padding)
    end
end

CachedOperator(op::Operator;padding::Bool=false) = default_CachedOperator(op;padding=padding)

"""
    cache(op::Operator)

Caches the entries of an operator, to speed up multiplying a Fun by the operator.
"""
cache(O::Operator;kwds...) = CachedOperator(O;kwds...)
cache(::Type{MT},O::Operator;kwds...) where {MT<:AbstractMatrix} = CachedOperator(MT,O;kwds...)

convert(::Type{Operator{T}},S::CachedOperator{T}) where {T} = S
convert(::Type{Operator{T}},S::CachedOperator) where {T} =
    CachedOperator(convert(Operator{T}, S.op),convert(AbstractMatrix{T}, S.data),
                    S.datasize,S.domainspace,S.rangespace,S.bandwidths)

domainspace(C::CachedOperator) = C.domainspace
rangespace(C::CachedOperator) = C.rangespace
bandwidths(C::CachedOperator{T,BM,M}) where {T<:Number,BM<:BandedMatrix,M<:Operator} = C.bandwidths
blockbandwidths(C::CachedOperator{T,BM,M}) where {T<:Number,BM<:BandedMatrix,M<:Operator} = C.bandwidths

# TODO: cache this information as well
@wrapperstructure CachedOperator

@propagate_inbounds function Base.getindex(B::CachedOperator,k::Integer,j::Integer)
    resizedata!(B,k,j)
    if k ≤ size(B.data,1) && j ≤ size(B.data,2)
        B.data[k,j]
    else
        zero(eltype(B))
    end
end

function Base.getindex(B::CachedOperator,k::AbstractRange,j::AbstractRange)
    if !isempty(k) && !isempty(j)
        resizedata!(B,maximum(k),maximum(j))
        B.data[k,j]
    else
        Matrix{eltype(B)}(length(k),length(j))
    end
end

function Base.getindex(B::CachedOperator,k::Integer)
    if size(B,1)==1
        B[1,k]
    elseif size(B,2)==1
        B[k,1]
    else
        error("Not implemented")
    end
end

@propagate_inbounds function Base.setindex!(B::CachedOperator,v,k,j)
    resizedata!(B,k,j)
    B.data[k,j] = v
end

@propagate_inbounds function Base.setindex!(B::CachedOperator,v,k::Integer)
    if size(B,1)==1
        B[1,k] = v
    elseif size(B,2)==1
        B[k,1] = v
    else
        error("Not implemented")
    end
 end

resizedata!(B::CachedOperator,::Colon,m::Integer) = resizedata!(B,size(B,1),m)
resizedata!(B::CachedOperator,n::Integer,::Colon) = resizedata!(B,n,size(B,2))

function mul_coefficients(B::CachedOperator,v::AbstractVector{T}) where T<:Number
    resizedata!(B,:,length(v))

    B.data*pad(v,size(B.data,2))
end


text:
"""
    BiComplex{T <: Real} = Complexification{Cmplex{T}}

A bi-complex number is a complexification of the complex numbers.
"""
const BiComplex{T <: Real} = Complexification{Cmplex{T}}

function BiComplex(z::Cmplex{T}) where T <: Real
    BiComplex{T}(z, zero(Cmplex{T}))
end

function BiComplex(a::T, b::T, c::T, d::T) where T <: Real
    BiComplex{T}(Cmplex{T}(a, b), Cmplex{T}(c, d))
end

function BiComplex(a::T, b::T, c::T) where T <: Real
    BiComplex{T}(Cmplex{T}(a, b), Cmplex{T}(c, zero(T)))
end

function BiComplex(a::T, b::T) where T <: Real
    BiComplex{T}(Cmplex{T}(a, b), zero(Cmplex{T}))
end

function BiComplex(a::T) where T <: Real
    BiComplex{T}(Cmplex{T}(a, zero(T)), zero(Cmplex{T}))
end

function BiComplex(a::Real, b::Real)
    BiComplex(promote(a, b)...)
end

function BiComplex(a::Real, b::Real, c::Real)
    BiComplex(promote(a, b, c)...)
end

function BiComplex(a::Real, b::Real, c::Real, d::Real)
    BiComplex(promote(a, b, c, d)...)
end

"""
    unreal(z::BiComplex)

The unreal part of a bi-complex number is a 3-vector.
"""
function unreal(z::BiComplex)
    vcat(unreal(z.l), array(z.r))
end

function show(io::IO, z::BiComplex)
    print(io, "[1: ")
    print(io, z.l.l)
    print(io, ", H: ")
    print(io, z.l.r)
    print(io, ", G: ")
    print(io, z.r.l)
    print(io, ", GH: ")
    print(io, z.r.r)
    print(io, "]")
end

function random(::Type{BiComplex{T}}) where T <: Real
    BiComplex{T}(random(Cmplex{T}), random(Cmplex{T}))
end

function barH(z::BiComplex{T}) where T <: Real
    BiComplex{T}(barH(z.l), barH(z.r))
end

function barG(z::BiComplex{T}) where T <: Real
    BiComplex{T}(z.l, -(z.r))
end

function quadrance(z::BiComplex)
    (z * barG(z)).l
end

function norm(z::BiComplex)
    abs2(quadrance(z))
end

function iszerodivisor(z::BiComplex)
    iszero(quadrance(z))
end

function inv(z::BiComplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    q = quadrance(z)
    barG(z) * barH(q) / abs2(q)
end

function lzng(w::BiComplex, x::BiComplex, y::BiComplex, z::BiComplex)
    w * barH(x) * barG(y) * barH(barG(z))
end

function lozenge(z0::BiComplex, z1::BiComplex, z2::BiComplex, z3::BiComplex)
    # Even permutations
    total = lzng(z0, z1, z2, z3)
    total += lzng(z2, z3, z0, z1)
    total += lzng(z1, z2, z0, z3)
    total += lzng(z0, z3, z1, z2)
    total += lzng(z2, z0, z1, z3)
    total += lzng(z1, z3, z2, z0)
    total += lzng(z0, z2, z3, z1)
    total += lzng(z3, z1, z0, z2)
    total += lzng(z2, z1, z3, z0)
    total += lzng(z3, z0, z2, z1)
    total += lzng(z1, z0, z3, z2)
    total += lzng(z3, z2, z1, z0)

    # Odd permutations
    total += lzng(z1, z2, z3, z0)
    total += lzng(z3, z0, z1, z2)
    total += lzng(z2, z0, z3, z1)
    total += lzng(z3, z1, z2, z0)
    total += lzng(z0, z1, z3, z2)
    total += lzng(z3, z2, z0, z1)
    total += lzng(z1, z0, z2, z3)
    total += lzng(z2, z3, z1, z0)
    total += lzng(z0, z2, z1, z3)
    total += lzng(z1, z3, z0, z2)
    total += lzng(z2, z1, z0, z3)
    total += lzng(z0, z3, z2, z1)

    total / 24
end

function anti_lozenge(z0::BiComplex, z1::BiComplex, z2::BiComplex, z3::BiComplex)
    # Even permutations
    total = lzng(z0, z1, z2, z3)
    total += lzng(z2, z3, z0, z1)
    total += lzng(z1, z2, z0, z3)
    total += lzng(z0, z3, z1, z2)
    total += lzng(z2, z0, z1, z3)
    total += lzng(z1, z3, z2, z0)
    total += lzng(z0, z2, z3, z1)
    total += lzng(z3, z1, z0, z2)
    total += lzng(z2, z1, z3, z0)
    total += lzng(z3, z0, z2, z1)
    total += lzng(z1, z0, z3, z2)
    total += lzng(z3, z2, z1, z0)

    # Odd permutations
    total -= lzng(z1, z2, z3, z0)
    total -= lzng(z3, z0, z1, z2)
    total -= lzng(z2, z0, z3, z1)
    total -= lzng(z3, z1, z2, z0)
    total -= lzng(z0, z1, z3, z2)
    total -= lzng(z3, z2, z0, z1)
    total -= lzng(z1, z0, z2, z3)
    total -= lzng(z2, z3, z1, z0)
    total -= lzng(z0, z2, z1, z3)
    total -= lzng(z1, z3, z0, z2)
    total -= lzng(z2, z1, z0, z3)
    total -= lzng(z0, z3, z2, z1)

    total
end

function exp(z::BiComplex)
    a, b, c, d = z.l.l, z.l.r, z.r.l, z.r.r
    e_a = exp(a)
    c_b, s_b = cos(b), sin(b)
    c_c, s_c = cos(c), sin(c)
    c_d, s_d = cosh(d), sinh(d)

    e_a * BiComplex(
        (c_b * c_c * c_d) + (s_b * s_c * s_d),
        (s_b * c_c * c_d) - (c_b * s_c * s_d),
        (c_b * s_c * c_d) - (s_b * c_c * s_d),
        (s_b * s_c * c_d) + (c_b * c_c * s_d)
    )
end


text:
"""
    Repo{T}

one-dimensional array, containing classifiable elements of type T.

Each element of Repo can be assigned to a class, identified by an integer.
All elements are intially in class one. Importantly, the time needed for
accessing elements, reassigning an element to a different class and
counting the number of elements is independent of the number of elements in the
class. This can be used for example to write highly-performant event-driven
simulations.

# Example
```jldoctest
julia> using Repos

julia> r=Repo(["Apple","Carrot","Cherry","Pepper"],2)
Repository of 4 objects in 2 classes

julia> setclass!(r,1,2)

julia> setclass!(r,3,2)

julia> print_repo(r)
Repository of 2 classes
  Class 1
    1 - 4: Pepper
    2 - 2: Carrot
  Class 2
    1 - 3: Cherry
    2 - 1: Apple

julia> fruits=class(r,2)
Class of 2 objects

julia> print_repo(fruits)
Class of 2 objects
    1 - 3: Cherry
    2 - 1: Apple
```
"""
mutable struct Repo{T} <: AbstractVector{T}
    starts::Array{Int,1}
    counts::Array{Int,1}
    ids::Array{Int,1}
    pos::Array{Int,1}
    cargo::Array{T,1}
    N::Int
    C::Int
end

"""
    RepoClass{T}

A class of elements in a Repo structure.

The RepoClass is an interface for the Repo structure. This means that instances
of RepoClass automatically update if the corresponding Repo is changed.

# Example
```jldoctest
julia> using Repos

julia> r=Repo(["Peter","Paul","Mary"],2)
Repository of 3 objects in 2 classes

julia> infected=class(r,2)
Class of 0 objects

julia> setclass!(r,1,2)

julia> infected
Class of 1 objects

julia> print_repo(infected)
Class of 1 objects
    1 - 1: Peter
```
"""
struct RepoClass{T} <: AbstractVector{T}
    repo::Repo{T}
    c::Int
end


text:
using Signals
using Test


text:
"""
    FullGenerator{T} <: Generator{T}

An Matrix or SparseMatrixCSC representation of the generator of a generator of a DiscretisedFluidQueue.

Higher memory requirements than LazyGenerator (as blocks are duplicated) but much faster matrix arithmetic.
"""
struct FullGenerator{T} <: Generator{T}
    B::Union{Array{Float64,2}, SparseArrays.SparseMatrixCSC{Float64,Int}}
    dq::DiscretisedFluidQueue{T}
end
FullGenerator(B,dq) = FullGenerator{typeof(dq.mesh)}(B,dq)

size(B::FullGenerator) = size(B.B)
# getindex(B::FullGenerator,i) = B.B[i]
getindex(B::FullGenerator,i,j) = B.B[i,j]
getindex(B::FullGenerator,t::Tuple,s::Tuple) = _getindex_from_tuple(B,t,s)
setindex!(B::FullGenerator,x::Float64,i,j) = (B.B[i,j]=x)
sum(B::FullGenerator; kwargs...) = sum(B.B; kwargs...)

fast_mul(A::AbstractArray{<:Real,2}, B::FullGenerator) = A*B.B
fast_mul(B::FullGenerator, A::AbstractArray{<:Real,2}) = fast_mul(A,B)
fast_mul(A::FullGenerator, B::FullGenerator) = A.B*B.B
fast_mul(A::FullGenerator, x::Real) = A.B*x
fast_mul(x::Real,A::FullGenerator) = fast_mul(A, x)

function show(io::IO, mime::MIME"text/plain", B::FullGenerator)
    if VERSION >= v"1.6"
        show(io, mime, B.B)
    else
        show(io, mime, Matrix(B.B))
    end
end

"""
    build_full_generator(dq::DiscretisedFluidQueue; v::Bool = false)

Returns a SparseMatrixCSC generator of a DiscretisedFluidQueue.
"""
function build_full_generator(dq::DiscretisedFluidQueue; v::Bool=false)
    lazy = build_lazy_generator(dq; v=v)
    return build_full_generator(lazy)
end

"""
    build_full_generator(lzB::LazyGenerator)

Returns a SparseMatrixCSC representation of a LazyGenerator.
"""
function build_full_generator(lzB::LazyGenerator)
    B = fast_mul(SparseArrays.SparseMatrixCSC{Float64,Int}(LinearAlgebra.I(size(lzB,1))),lzB)
    return FullGenerator(B,lzB.dq)
end


text:
@testset "Utilities" begin

    T = rand_kruskal3(2, (10, 20, 30), true)
    @test size(T) == (10, 20, 30)

    @testset "khatrirao()" begin
        A = randn(5,4)
        @test_throws DimensionMismatch TensorDecompositions.khatrirao!(randn(5,5), A)
        @test_throws DimensionMismatch TensorDecompositions.khatrirao!(randn(4,5), A)
        @test TensorDecompositions.khatrirao!(randn(5,4), A) == A
        @test TensorDecompositions.khatrirao(A) == A

        @test_throws DimensionMismatch TensorDecompositions.khatrirao!(randn(30,4), rand(2,4), rand(5, 4), rand(3, 5))
        @test_throws DimensionMismatch TensorDecompositions.khatrirao!(randn(30,5), rand(2,4), rand(5, 4), rand(3, 4))
        @test_throws DimensionMismatch TensorDecompositions.khatrirao!(randn(29,4), rand(2,4), rand(5, 4), rand(3, 4))

        @test @inferred(TensorDecompositions.khatrirao!(randn(30,4), randn(2,4), randn(5, 4), randn(3, 4))) isa Matrix{Float64}

        @test_throws DimensionMismatch TensorDecompositions.khatrirao(rand(2,4), rand(5, 4), rand(3, 5))
        X = @inferred TensorDecompositions.khatrirao(rand(2,4), rand(5, 4), rand(3, 4))
        @test X isa Matrix{Float64}
        @test size(X) == (30, 4)
    end

    @testset "_row_unfold()" begin
        res = @inferred TensorDecompositions._row_unfold(T, 1)
        @test size(res) == (10, 600)

        res = @inferred TensorDecompositions._row_unfold(T, 2)
        @test size(res) == (20, 300)

        res = @inferred TensorDecompositions._row_unfold(T, 3)
        @test size(res) == (30, 200)
    end

    @testset "_col_unfold()" begin
        res = @inferred TensorDecompositions._col_unfold(T, 1)
        @test size(res) == (600, 10)

        res = @inferred TensorDecompositions._col_unfold(T, 2)
        @test size(res) == (300, 20)

        res = @inferred TensorDecompositions._col_unfold(T, 3)
        @test size(res) == (200, 30)
    end

    @testset "tensorcontractmatrices()" begin
        factors = TensorDecompositions._random_factors(size(T), (5, 2, 6))
        res = @inferred tensorcontractmatrices(T, factors)
        @test size(res) == (5, 2, 6)
    end

end


text:
function eig_min!(linear_solver::abstract_linear_system_solver, K::SparseMatrixCSC{Float64,Int64}, x_0::Array{Float64,1}, n::Int64, m::Int64, tol::Float64)
	# computes smallest magnitude eigenvalue and eigenvector viz inverse iteration
	try
		@assert(n + m == size(K,1))
		@assert(size(K,1) == size(K,2))
		@assert(n == length(x_0))

		eigen_value = 0.0;
		x = x_0;
		sol = zeros(n + m);
		MAX_IT = 20;
		err = Inf;
		for i = 1:MAX_IT

			rhs = [x; zeros(m)];
			linear_solver.ls_solve!(rhs, sol);
			x_old = x;
			x = sol[1:n];
			eigen_value = 1.0/norm(x,2)
			x = x*eigen_value;

			err = norm(x - x_old, 1)

			if err < tol
				println(eigen_value)
				println("converged!")
				break
			end
		end

		return eigen_value, x, err
	catch e
		println("ERROR in eig_min")
		throw(e)
	end
end

function is_B_matrix(B::SparseMatrixCSC{Float64,Int64})
	# checks that matrix is made of 1x1 or 2x2 blocks

end

function eigsigns_of_B(B::SparseMatrixCSC{Float64,Int64}, toler::Float64)
	#

	@assert(size(B,1) == size(B,2))
	n = size(B,1);
	trans_eigenvals = zeros(n);
	i = 1;
	while i <= n
		if i == n || B[i, i+1] == 0
			trans_eigenvals[i] = B[i, i];
			i = i + 1;
		else
			trans_eigenvals[i:(i+1)] = eigvals(full(B[i:i+1,i:i+1]));
			i = i + 2;
		end
	end

	pos_eigs = sum(trans_eigenvals .> toler);
	neg_eigs = sum(trans_eigenvals .< toler);
	zero_eigs = length(trans_eigenvals) - pos_eigs - neg_eigs;

	return pos_eigs, neg_eigs, zero_eigs
end

if false
	H = spdiagm([1.1,-1.0]);
	A = sparse([[1.0, 1.0] ])
	K = [ [H A]; [A' spzeros(1,1)] ];

	linear_solver = linear_solver_JULIA();
	linear_solver.initialize(K);
	linear_solver.ls_factor();

	x_0 = randn(2);

	eig_min!(linear_solver, K, x_0, 2, 1, 1e-6)
end

max_forks_repo_name: ohinder/Everest


text:
# Binding query parameters.

mutable struct BindNode <: AbstractSQLNode
    over::Union{SQLNode, Nothing}
    list::Vector{SQLNode}
    label_map::OrderedDict{Symbol, Int}

    function BindNode(; over = nothing, list, label_map = nothing)
        if label_map !== nothing
            return new(over, list, label_map)
        end
        n = new(over, list, OrderedDict{Symbol, Int}())
        for (i, l) in enumerate(n.list)
            name = label(l)
            if name in keys(n.label_map)
                err = DuplicateLabelError(name, path = [l, n])
                throw(err)
            end
            n.label_map[name] = i
        end
        n
    end
end

BindNode(list...; over = nothing) =
    BindNode(over = over, list = SQLNode[list...])

"""
    Bind(; over = nothing; list)
    Bind(list...; over = nothing)

The `Bind` node binds the query parameters in an inner query to make it
a correlated subquery.

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id]);

julia> visit_occurrence = SQLTable(:visit_occurrence, columns = [:visit_occurrence_id, :person_id]);

julia> q = From(person) |>
           Where(Fun.exists(From(visit_occurrence) |>
                            Where(Get.person_id .== Var.person_id) |>
                            Bind(Get.person_id)));

julia> print(render(q))
SELECT "person_1"."person_id"
FROM "person" AS "person_1"
WHERE (EXISTS (
  SELECT NULL
  FROM "visit_occurrence" AS "visit_occurrence_1"
  WHERE ("visit_occurrence_1"."person_id" = "person_1"."person_id")
))
```
"""
Bind(args...; kws...) =
    BindNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Bind), pats::Vector{Any}) =
    dissect(scr, BindNode, pats)

function PrettyPrinting.quoteof(n::BindNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Bind))
    if isempty(n.list)
        push!(ex.args, Expr(:kw, :list, Expr(:vect)))
    else
        append!(ex.args, quoteof(n.list, ctx))
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::BindNode) =
    label(n.over)

rebase(n::BindNode, n′) =
    BindNode(over = rebase(n.over, n′), list = n.list)


text:
solve(i::BudgetedMSetInstance, ::LinearProgramming; kwargs...) = budgeted_msets_lp(i; kwargs...)

function _budgeted_msets_lp_sub(i::BudgetedMSetInstance, solver)
  model = Model(solver)
  @variable(model, x[1:length(values(i))], Bin)
  @objective(model, Max, dot(x, values(i)))
  @constraint(model, sum(x) <= m(i))
  @constraint(model, c, dot(x, weights(i)) >= 0)

  set_silent(model)

  return model, x, c
end

function budgeted_msets_lp(i::BudgetedMSetInstance; solver=nothing, β::Int=budget(i))
  # Solve for all budgets at once, even though it is not really more efficient than looping outside this function.
  model, x, c = _budgeted_msets_lp_sub(i, solver)
  set_normalized_rhs(c, β)
  optimize!(model)

  V = Array{Float64, 3}(undef, m(i), dimension(i) + 1, budget(i) + 1)
  S = Dict{Tuple{Int, Int, Int}, Vector{Int}}()

  if termination_status(model) == MOI.OPTIMAL
    sol = findall(JuMP.value.(x) .>= 0.5)
    V[m(i), 0 + 1, β + 1] = objective_value(model)
    S[m(i), 0, β] = sol
  else
    V[m(i), 0 + 1, β + 1] = -Inf
    S[m(i), 0, β] = Int[-1]
  end

  return BudgetedMSetSolution(i, S[m(i), 0, β], V, S)
end

function budgeted_msets_lp_select(i::BudgetedMSetInstance, budgets; solver=nothing)
  model, x, c = _budgeted_msets_lp_sub(i, solver)

  V = Array{Float64, 3}(undef, m(i), dimension(i) + 1, budget(i) + 1)
  S = Dict{Tuple{Int, Int, Int}, Vector{Int}}()
  for budget in budgets
    set_normalized_rhs(c, budget)
    optimize!(model)

    if termination_status(model) == MOI.OPTIMAL
      sol = findall(JuMP.value.(x) .>= 0.5)
      V[m(i), 0 + 1, budget + 1] = objective_value(model)
      S[m(i), 0, budget] = sol
    else
      V[m(i), 0 + 1, budget + 1] = -Inf
      S[m(i), 0, budget] = Int[-1]
    end
  end

  return BudgetedMSetSolution(i, S[m(i), 0, maximum(budgets)], V, S)
end

function budgeted_msets_lp_all(i::BudgetedMSetInstance; solver=nothing, max_budget::Int=budget(i))
  # Solve for all budgets at once, even though it is not really more efficient than looping outside this function.
  return budgeted_msets_lp_select(i, 0:max_budget, solver=solver)
end


text:
# Combined 1D/2D plots across multiple frames.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, PyPlot, Glob, Printf, LaTeXStrings
using Vlasiator: set_args, prep2d, set_colorbar, set_plot

@assert matplotlib.__version__ >= "3.4" "Matplotlib version too old: subfigure is not supported!"

const files = glob("bulk*.vlsv", ".")
const nfile = length(files)

meta = load(files[1])

Re = Vlasiator.Re # Earth radii
const x1, x2 = 8.0, 29.0
point1 = [x1, 0, 0] .* Re
point2 = [x2, 0, 0] .* Re

cellids, distances, coords = getcellinline(meta, point1, point2)

const loc = range(x1, x2, length=length(cellids))

# Upper/lower limits for each variable
const ρmin, ρmax = 0.0, 10.0     # [amu/cc]
const vmin, vmax = -640.0, 0.0   # [km/s]
const pmin, pmax = 0.0, 1.82     # [nPa]
const bmin, bmax = -25.0, 60.0   # [nT]
const vamin, vamax = 0.0, 250.0  # [km/s]
const vsmin, vsmax = 0.0, 400.0  # [km/s]

const fontsize = 14

function init_figure()
   fig = plt.figure(constrained_layout=true, figsize=(12, 12))
   subfigs = fig.subfigures(1, 2, wspace=0.05)

   axsL = subfigs[1].subplots(4, 1, sharex=true)
   axsR = subfigs[2].subplots(2, 1, sharex=true, sharey=true)

   # Set line plots' axes
   axsL[end].set_xlim(x1, x2)

   axsL[1].set_ylim(ρmin, ρmax)
   axsL[2].set_ylim(vmin, vmax)
   axsL[3].set_ylim(pmin, pmax)
   axsL[4].set_ylim(bmin, bmax)
   for ax in axsL
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.grid(true)
   end

   axsL[end].set_xlabel(L"x [$R_E$]"; fontsize)
   axsL[1].set_ylabel("Density [amu/cc]"; fontsize)
   axsL[2].set_ylabel("Velocity [km/s]"; fontsize)
   axsL[3].set_ylabel("Pressure [nPa]"; fontsize)
   axsL[4].set_ylabel("Magnetic field [nT]"; fontsize)

   # Set contour plots' axes and colorbars
   cmap = matplotlib.cm.turbo
   colorscale = Linear
   axisunit = RE

   pArgs1 = set_args(meta, "VA", axisunit; normal=:none)
   cnorm1, cticks1 = set_colorbar(colorscale, vamin, vamax)
   cnorm2, cticks2 = set_colorbar(colorscale, vsmin, vsmax)

   for ax in axsR
      ax.set_aspect("equal")

      # Set border line widths
      for loc in ("left", "bottom", "right", "top")
         edge = get(ax.spines, loc, nothing)
         edge.set_linewidth(2.0)
      end
      ax.xaxis.set_tick_params(width=2.0, length=3)
      ax.yaxis.set_tick_params(width=2.0, length=3)
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())

      ax.set_ylabel(L"Y [$R_E$]"; fontsize)
   end

   axsR[2].set_xlabel(L"X [$R_E$]"; fontsize)

   axsR[1].set_title("Alfven speed", fontsize="x-large")
   axsR[2].set_title("Sound speed", fontsize="x-large")

   fig.suptitle("Density Pulse Run", fontsize="xx-large")

   return fig, axsL, axsR, cmap, cnorm1, cticks1, cnorm2, cticks2
end

"Update frame."
function process(i, file)
   println("i = $i/$nfile, file = $file")
   local meta = load(file)

   p_extract = readvariable(meta, "vg_pressure", cellids) .* 1e9 |> vec # [nPa]
   rho_extract = readvariable(meta, "proton/vg_rho", cellids) |> vec
   v_extract = readvariable(meta, "proton/vg_v", cellids)
   vmag2_extract = sum(v_extract.^2, dims=1) |> vec
   pdyn_extract = rho_extract .* Vlasiator.mᵢ .* vmag2_extract .* 1e9 # [nPa]

   bz = readvariable(meta, "vg_b_vol", cellids)[3,:] .* 1e9 #[nT]

   imagnetopause_ = findfirst(<(0.0), bz)

   axsL[1].plot(loc, rho_extract ./ 1e6, label="Proton density", color="#1f77b4")
   vl1 = axsL[1].vlines(loc[imagnetopause_], ρmin, ρmax;
      colors="r", linestyle="dashed", alpha=0.5)

   axsL[2].plot(loc, v_extract[1,:] ./ 1e3, label="Vx", color="#1f77b4")
   vl2 = axsL[2].vlines(loc[imagnetopause_], vmin, vmax;
      colors="r", linestyle="dashed", alpha=0.5)

   axsL[3].plot(loc, pdyn_extract, label="Dynamic", color="#1f77b4")
   axsL[3].plot(loc, p_extract, label="Thermal", color="#ff7f0e")
   vl3 = axsL[3].vlines(loc[imagnetopause_], pmin, pmax;
      colors="r", linestyle="dashed", alpha=0.5)

   axsL[4].plot(loc, bz, label="Bz", color="#1f77b4")
   hl4 = axsL[4].hlines(0.0, loc[1], loc[end]; colors="k", linestyle="dashed", alpha=0.2)
   vl4 = axsL[4].vlines(loc[imagnetopause_], bmin, bmax;
       colors="r", linestyle="dashed", alpha=0.5)

   str_title = @sprintf "Sun-Earth line, t= %4.1fs" meta.time
   subfigs[1].suptitle(str_title, fontsize="x-large")

   axsL[2].legend(;loc="upper right", fontsize)
   axsL[3].legend(;loc="lower right", fontsize)
   axsL[4].legend(;loc="upper right", fontsize)

   x, y = Vlasiator.get_axis(pArgs1)

   data = prep2d(meta, "VA")'
   c1 = axsR[1].pcolormesh(x, y, data ./ 1e3, norm=cnorm1, cmap=cmap, shading="nearest")

   data = prep2d(meta, "VS")'
   c2 = axsR[2].pcolormesh(x, y, data ./ 1e3, norm=cnorm2, cmap=cmap, shading="nearest")

   if i == 1
      cb1 = colorbar(c1; ax=axsR[1], ticks=cticks1, fraction=0.046, pad=0.04)
      cb1.ax.set_ylabel("[km/s]"; fontsize)
      cb1.outline.set_linewidth(1.0)

      cb2 = colorbar(c2; ax=axsR[2], ticks=cticks2, fraction=0.046, pad=0.04)
      cb2.ax.set_ylabel("[km/s]"; fontsize)
      cb2.outline.set_linewidth(1.0)
   end

   savefig("out/"*lpad(i, 4, '0')*".png", bbox_inches="tight")

   for ax in axsL
      for line in ax.get_lines()
         line.remove()
      end
   end
   vl1.remove()
   vl2.remove()
   vl3.remove()
   vl4.remove()
   hl4.remove()
end

# Loop over snapshots
for (i, file) in enumerate(files)
   process(i, file)
end

close(fig)
println("Finished!")


text:
function HaltonSequenceGenerator(arg0::jint)
    return HaltonSequenceGenerator((jint,), arg0)
end

function HaltonSequenceGenerator(arg0::jint, arg1::Vector{jint}, arg2::Vector{jint})
    return HaltonSequenceGenerator((jint, Vector{jint}, Vector{jint}), arg0, arg1, arg2)
end

function get_next_index(obj::HaltonSequenceGenerator)
    return jcall(obj, "getNextIndex", jint, ())
end

function next_vector(obj::HaltonSequenceGenerator)
    return jcall(obj, "nextVector", Vector{jdouble}, ())
end

function skip_to(obj::HaltonSequenceGenerator, arg0::jint)
    return jcall(obj, "skipTo", Vector{jdouble}, (jint,), arg0)
end


text:
module CoveragePackage

abstract type AbstractExtractor end

struct ExtractKeyAsField{S,V} <: AbstractExtractor
	key::S
	item::V
end

Base.hash(e::ExtractKeyAsField, h::UInt) = hash((e.key, e.item), h)
Base.:(==)(e1::ExtractKeyAsField, e2::ExtractKeyAsField) = e1.key == e2.key && e1.item == e2.item

end

max_stars_repo_name: racinmat/CoveragePackage

max_issues_repo_name: racinmat/CoveragePackage

max_forks_repo_name: racinmat/CoveragePackage


text:
################################################################################
#
#  FakeFmpqMat.jl : Model fmpq_mat's as fmpz_mat's with denominator
#
################################################################################

export num, den

numerator(x::FakeFmpqMat) = deepcopy(x.num)

denominator(x::FakeFmpqMat) = deepcopy(x.den)

rows(x::FakeFmpqMat) = rows(x.num)

cols(x::FakeFmpqMat) = cols(x.num)

function simplify_content!(x::FakeFmpqMat)
  c = content(x.num)
  c = gcd(c, x.den)
  if c != 1
    divexact!(x.num, x.num, c)
    divexact!(x.den, x.den, c)
  end
  y = canonical_unit(x.den)
  if !isone(y)
    mul!(x.den, x.den, y)
    mul!(x.num, x.num, y)
  end
end

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

function hash(a::FakeFmpqMat, b::UInt)
  h = xor(Base.hash(a.num, b), Base.hash(a.den, b))
  h = xor(h, Base.hash(b))
  h = (h << 1) | (h >> (sizeof(Int)*8 - 1))
  return h
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FakeFmpqMat)
  print(io, "FakeFmpqMat with numerator\n", x.num, "\nand denominator ", x.den)
end

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

function -(x::FakeFmpqMat)
  return FakeFmpqMat(-x.num, x.den, true)
end

function inv(x::FakeFmpqMat)
  i, d_i = pseudo_inv(x.num)
  #TODO gcd d_i and x.den 1st!!!
  i *= x.den
  z = FakeFmpqMat(i,d_i)
  simplify_content!(z)
  return z
end

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

function +(x::FakeFmpqMat, y::FakeFmpqMat)
  t = y.den*x.num + x.den*y.num
  d = x.den*y.den
  z = FakeFmpqMat(t,d)
  return z
end

function *(x::FakeFmpqMat, y::FakeFmpqMat)
  t = x.num*y.num
  d = x.den * y.den
  z = FakeFmpqMat(t, d)
  return z
end

function mul!(z::FakeFmpqMat, x::FakeFmpqMat, y::FakeFmpqMat)
  mul!(z.num, x.num, y.num)
  mul!(z.den, x.den, y.den)
  simplify_content!(z)
  return z
end

################################################################################
#
#  Adhoc binary operations
#
################################################################################

function *(x::FakeFmpqMat, y::fmpz_mat)
  t = x.num*y
  z = FakeFmpqMat(t, denominator(x))
  simplify_content!(z)
  return z
end

function *(x::fmpz_mat, y::FakeFmpqMat)
  t = x*y.num
  z = FakeFmpqMat(t, denominator(y))
  simplify_content!(z)
  return z
end

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

==(x::FakeFmpqMat, y::FakeFmpqMat) = (x.num == y.num) && (x.den == y.den)

################################################################################
#
#  Conversion
#
################################################################################

to_array(x::FakeFmpqMat) = (x.num, x.den)

function to_fmpz_mat(x::FakeFmpqMat)
  x.den != 1 && error("Denominator has to be 1")
  return numerator(x)
end

################################################################################
#
#  Hermite normal form for numerator
#
################################################################################

function hnf!(x::FakeFmpqMat)
  x.num = hnf(x.num)
end

function hnf(x::FakeFmpqMat, shape = :lowerleft)
  h = _hnf(x.num, shape)
  return FakeFmpqMat(h, denominator(x))
end

################################################################################
#
#  Sub
#
################################################################################

function sub(x::FakeFmpqMat, r::UnitRange{Int}, c::UnitRange{Int})
  z = FakeFmpqMat(sub(x.num, r, c), x.den)
  return z
end

function Base.deepcopy_internal(x::FakeFmpqMat, dict::ObjectIdDict)
  z = FakeFmpqMat()
  z.num = Base.deepcopy_internal(x.num, dict)
  z.den = Base.deepcopy_internal(x.den, dict)
  z.rows = rows(x)
  z.cols = cols(x)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
  return z
end

################################################################################
#
#  Determinant
#
################################################################################

function det(x::FakeFmpqMat)
  rows(x) != cols(x) && error("Matrix must be square")

  return det(x.num)//(x.den)^rows(x)
end