text:
doc"""
`CoulombPotential(Z;R=(0.,0.,0.))`

represents the Coulomb potential $v(r) = -Z/|r-R|$
"""
type CoulombPotential <: ScalarFunc3d{Float64}
    Z::Float64
    R::Point3{Float64}
end

CoulombPotential(Z;R=(0.,0.,0.)) = CoulombPotential(Z,R)

function Base.call(f::CoulombPotential, x::Float64, y::Float64, z::Float64)
    x0,y0,z0 = f.R
    r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2)
    -f.Z/r
end

Base.call(f::CoulombPotential, p::Point3{Float64}) = f(p[1],p[2],p[3])
@vectorize2 Base.call CoulombPotential Point3{Float64} Float64

function setposition(f::CoulombPotential, R::Point3{Float64})
    f.R = R
end

max_stars_repo_name: johncsnyder/Element


text:
@enum IndexType::UInt8 NotAnIndex=0 LoopIndex=1 ComputedIndex=2 SymbolicIndex=3

Base.:|(u::Unsigned, it::IndexType) = u | UInt8(it)
Base.:(==)(u::Unsigned, it::IndexType) = (u % UInt8) == UInt8(it)

function _append_fields!(t::Expr, body::Expr, sym::Symbol, ::Type{T}) where {T}
  gf = GlobalRef(Core,:getfield)
  for f ∈ 1:fieldcount(T)
    TF = fieldtype(T, f)
    Base.issingletontype(TF) && continue
    gfcall = Expr(:call, gf, sym, f)
    if fieldcount(TF) ≡ 0
      push!(t.args, gfcall)
    elseif TF <: DataType
      push!(t.args, Expr(:call, Expr(:curly, lv(:StaticType), gfcall)))
    else
      newsym = gensym(sym)
      push!(body.args, Expr(:(=), newsym, gfcall))
      _append_fields!(t, body, newsym, TF)
    end
  end
  return nothing
end
@generated function flatten_to_tuple(r::T) where {T}
  body = Expr(:block, Expr(:meta,:inline))
  t = Expr(:tuple)
  if Base.issingletontype(T)
    nothing
  elseif fieldcount(T) ≡ 0
    push!(t.args, :r)
  elseif T <: DataType
    push!(t.args, Expr(:call, Expr(:curly, lv(:StaticType), :r)))
  else
    _append_fields!(t, body, :r, T)
  end
  push!(body.args, t)
  body
end
function rebuild_fields(offset::Int, ::Type{T}) where {T}
  gf = GlobalRef(Core,:getfield)
  call = (T <: Tuple) ? Expr(:tuple) : Expr(:new, T)
  for f ∈ 1:fieldcount(T)
    TF = fieldtype(T, f)
    if Base.issingletontype(TF)
      push!(call.args, TF.instance)
    elseif fieldcount(TF) ≡ 0
      push!(call.args, Expr(:call, gf, :t, (offset += 1), false))
    elseif TF <: DataType
      push!(call.args, Expr(:call, lv(:gettype), Expr(:call, gf, :t, (offset += 1), false)))
    else
      arg, offset = rebuild_fields(offset, TF)
      push!(call.args, arg)
    end
  end
  return call, offset
end
@generated function reassemble_tuple(::Type{T}, t::Tuple) where {T}
  if Base.issingletontype(T)
    return T.instance
  elseif fieldcount(T) ≡ 0
    call = Expr(:call, GlobalRef(Core,:getfield), :t, 1, false)
  elseif T <: DataType
    call = Expr(:call, lv(:gettype), Expr(:call, GlobalRef(Core,:getfield), :t, 1, false))
  else
    call, _ = rebuild_fields(0, T)
  end
  Expr(:block, Expr(:meta,:inline), call)
end

"""
    ArrayRefStruct

A condensed representation of an [`ArrayReference`](@ref).
It supports array-references with up to 8 indexes, where the data for each consecutive index is packed into corresponding 8-bit fields
of `index_types` (storing the enum `IndexType`), `indices` (the `id` for each index symbol), and `offsets` (currently unused).
"""
struct ArrayRefStruct{array,ptr}
  index_types::UInt128
  indices::UInt128
  offsets::UInt128
  strides::UInt128
end
array_and_ptr(@nospecialize(ar::ArrayRefStruct{a,p})) where {a,p} = (a::Symbol,p::Symbol)
# array(@nospecialize(ar::ArrayRefStruct{a,p})) where {a,p} = a::Symbol
# ptr(@nospecialize(ar::ArrayRefStruct{a,p})) where {a,p}   = p::Symbol

function findindoradd!(v::Vector{T}, s::T) where {T}
  ind = findfirst(==(s), v)
  ind === nothing || return ind
  push!(v, s)
  length(v)
end
function ArrayRefStruct(ls::LoopSet, mref::ArrayReferenceMeta, arraysymbolinds::Vector{Symbol}, ids::Vector{Int})
  index_types = zero(UInt128)
  indices = zero(UInt128)
  offsets = zero(UInt128)
  strides = zero(UInt128)
  @unpack loopedindex, ref = mref
  indv = ref.indices
  offv = ref.offsets
  strv = ref.strides
  # we can discard that the array was considered discontiguous, as it should be recovered from type information
  start = 1 + (first(indv) === DISCONTIGUOUS)
  for (n,ind) ∈ enumerate(@view(indv[start:end]))
    index_types <<= 8
    indices <<= 8
    offsets <<= 8
    offsets |= (offv[n] % UInt8)
    strides <<= 8
    strides |= (strv[n] % UInt8)
    if loopedindex[n]
      index_types |= LoopIndex
      if strv[n] ≠ 0
        indices |= getloopid(ls, ind)
      end
    else
      parent = get(ls.opdict, ind, nothing)
      @assert !(parent === nothing) "Index $ind not found in array."
      # if parent === nothing
      #     index_types |= SymbolicIndex
      #     indices |= findindoradd!(arraysymbolinds, ind)
      # else
      index_types |= ComputedIndex
      indices |= ids[identifier(parent)]
      # end
    end
  end
  ArrayRefStruct{mref.ref.array,mref.ptr}( index_types, indices, offsets, strides )
end

"""
    OperationStruct

A condensed representation of an [`Operation`](@ref).
"""
struct OperationStruct <: AbstractLoopOperation
    # instruction::Instruction
  loopdeps::UInt128
  reduceddeps::UInt128
  childdeps::UInt128
  parents₀::UInt128
  parents₁::UInt128
  node_type::OperationType
  symid::UInt16
  array::UInt8
end
optype(os) = os.node_type

function findmatchingarray(ls::LoopSet, mref::ArrayReferenceMeta)
  id = 0x01
  for r ∈ ls.refs_aliasing_syms
    r == mref && return id
    id += 0x01
  end
  0x00
end
filled_8byte_chunks(u::T) where {T<:Unsigned} = sizeof(T) - (leading_zeros(u) >>> 3)

function shifted_loopset(ls::LoopSet, loopsyms::Vector{Symbol})
  ld = zero(UInt128) # leading_zeros(ld) >> 2 yields the number of loopdeps
  for d ∈ loopsyms
    ld <<= 4
    ld |= getloopid(ls, d)::Int
  end
  ld
end
loopdeps_uint(ls::LoopSet, op::Operation) = shifted_loopset(ls, loopdependencies(op))
reduceddeps_uint(ls::LoopSet, op::Operation) = shifted_loopset(ls, reduceddependencies(op))
childdeps_uint(ls::LoopSet, op::Operation) = shifted_loopset(ls, reducedchildren(op))
function parents_uint(oppv::AbstractVector{Operation})
  p = zero(UInt128)
  for parent ∈ oppv
    p <<= 16
    p |= identifier(parent)
  end
  p
end
function parents_uint(op::Operation)
  opv = parents(op)
  N = length(opv)
  @assert N ≤ 16
  p0 = parents_uint(view(opv, 1:min(8,N)))
  p1 = N > 8 ? parents_uint(view(opv, 9:N)) : zero(p0)
  p0, p1
end
function recursively_set_parents_true!(x::Vector{Bool}, op::Operation)
  x[identifier(op)] && return nothing # don't redescend
  x[identifier(op)] = true
  for opp ∈ parents(op)
    recursively_set_parents_true!(x, opp)
  end
  return nothing
end
function getroots(ls::LoopSet)::Vector{Bool}
  rooted = Vector{Bool}(undef, length(operations(ls)))
  getroots!(rooted, ls)
end
function getroots!(rooted::Vector{Bool}, ls::LoopSet)
  fill!(rooted, false)
  ops = operations(ls)
  for or ∈ ls.outer_reductions
    recursively_set_parents_true!(rooted, ops[or])
  end
  for op ∈ ops
    isstore(op) && recursively_set_parents_true!(rooted, op)
  end
  remove_outer_reducts!(rooted, ls)
  return rooted
end
function OperationStruct!(varnames::Vector{Symbol}, ids::Vector{Int}, ls::LoopSet, op::Operation)
  instr = instruction(op)
  ld = loopdeps_uint(ls, op)
  rd = reduceddeps_uint(ls, op)
  cd = childdeps_uint(ls, op)
  p0, p1 = parents_uint(op)
  array = accesses_memory(op) ? findmatchingarray(ls, op.ref) : 0x00
  ids[identifier(op)] = id = findindoradd!(varnames, name(op))
  OperationStruct(
    ld, rd, cd, p0, p1, op.node_type, id, array
  )
end
## turn a LoopSet into a type object which can be used to reconstruct the LoopSet.

@inline zerorangestart(r::Base.OneTo) = CloseOpen(maybestaticlast(r))
@inline zerorangestart(r::CloseOpen) = CloseOpen(length(r))
@inline zerorangestart(r::CloseOpen{Zero}) = r
@inline zerorangestart(r::AbstractUnitRange) = Zero():One():(maybestaticlast(r)-maybestaticfirst(r))
@inline zerorangestart(r::AbstractRange) = Zero():static_step(r):(maybestaticlast(r)-maybestaticfirst(r))
@inline zerorangestart(r::CartesianIndices) = CartesianIndices(map(zerorangestart, r.indices))
@inline zerorangestart(r::ArrayInterface.OptionallyStaticUnitRange{StaticInt{1}}) = CloseOpen(maybestaticlast(r))

function loop_boundary!(q::Expr, ls::LoopSet, loop::Loop, shouldindbyind::Bool)
  if isstaticloop(loop) || loop.rangesym === Symbol("")
    call = Expr(:call, :(:))
    f = gethint(first(loop))
    s = gethint(step(loop))
    l = gethint(last(loop))
    if !shouldindbyind
      l -= f
      f = 0
    end
    pushexpr!(call, staticexpr(f))
    isone(s) || pushexpr!(call, staticexpr(s))
    pushexpr!(call, staticexpr(l))
    push!(q.args, call)
  elseif shouldindbyind
    push!(q.args, loop.rangesym)
  else
    push!(q.args, Expr(:call, lv(:zerorangestart), loop.rangesym))
  end
end

function loop_boundaries(ls::LoopSet, shouldindbyind::Vector{Bool})
    lbd = Expr(:tuple)
    for (ibi, loop) ∈ zip(shouldindbyind, ls.loops)
        loop_boundary!(lbd, ls, loop, ibi)
    end
    lbd
end

tuple_expr(v) = tuple_expr(identity, v)
function tuple_expr(f, v)
    t = Expr(:tuple)
    for vᵢ ∈ v
        push!(t.args, f(vᵢ))
    end
    t
end

function argmeta_and_consts_description(ls::LoopSet, arraysymbolinds)
  Expr(
    :tuple,
    length(arraysymbolinds),
    tuple_expr(ls.outer_reductions),
    tuple_expr(first, ls.preamble_symsym),
    tuple_expr(ls.preamble_symint),
    tuple_expr(ls.preamble_symfloat),
    tuple_expr(ls.preamble_zeros),
    tuple_expr(ls.preamble_funcofeltypes)
  )
end
@inline vdata(v::Vec) = getfield(v, :data)
@inline vdata(v::VecUnroll) = getfield(v, :data)
@inline vdata(x) = x

function loopset_return_value(ls::LoopSet, ::Val{extract}) where {extract}
    @assert !iszero(length(ls.outer_reductions))
    if isone(length(ls.outer_reductions))
        op = getop(ls, ls.outer_reductions[1])
        if extract
            # if (isu₁unrolled(op) | isu₂unrolled(op))
                Expr(:call, :vdata, Symbol(mangledvar(op), "##onevec##"))
            # else
                # Expr(:call, :data, mangledvar(op))
            # end
        else
            Symbol(mangledvar(op), "##onevec##")
        end
    else#if length(ls.outer_reductions) > 1
        ret = Expr(:tuple)
        ops = operations(ls)
        for or ∈ ls.outer_reductions
            op = ops[or]
            if extract
                push!(ret.args, Expr(:call, :vdata, Symbol(mangledvar(op), "##onevec##")))
            else
                push!(ret.args, Symbol(mangledvar(ops[or]), "##onevec##"))
            end
        end
        ret
    end
end
const DROPPEDCONSTANT = Instruction(Symbol("##DROPPED#CONSTANT##"),Symbol("##DROPPED#CONSTANT##"))
function skip_constant(instr::Instruction)
  (((instr == LOOPCONSTANT) || (instr.mod === :numericconstant)) || (instr == DROPPEDCONSTANT)) || instr.mod === GLOBALCONSTANT
end

function add_reassigned_syms!(q::Expr, ls::LoopSet)
    for op ∈ operations(ls)
        if isconstant(op)
            instr = instruction(op)
            skip_constant(instr) || push!(q.args, instr.instr)
        end
    end
end
function add_external_functions!(q::Expr, ls::LoopSet)
    for op ∈ operations(ls)
        if iscompute(op)
            instr = instruction(op)
            if instr.mod !== :LoopVectorization
                push!(q.args, instr.instr)
            end
        end
    end
end

function check_if_empty(ls::LoopSet, q::Expr)
    lb = loop_boundaries(ls, fill(false, length(ls.loops)))
    Expr(:if, Expr(:call, :!, Expr(:call, :any, :isempty, lb)), q)
end

val(x) = Expr(:call, Expr(:curly, :Val, x))

@inline gespf1(x, i) = gesp(x, i)
@inline gespf1(x::StridedPointer{T,1}, i::Tuple{I}) where {T,I<:Integer} = gesp(x, i)
@inline gespf1(x::StridedBitPointer{T,1}, i::Tuple{I}) where {T,I<:Integer} = gesp(x, i)
@inline gespf1(x::StridedPointer{T,1}, i::Tuple{Zero}) where {T} = x
@inline gespf1(x::StridedBitPointer{T,1}, i::Tuple{Zero}) where {T} = x
@generated function gespf1(x::AbstractStridedPointer{T,N,C,B,R}, i::Tuple{I}) where {T,N,I<:Integer,C,B,R}
  ri = argmin(R)
  quote
    $(Expr(:meta,:inline))
    p, li = VectorizationBase.tdot(x, (vsub_nsw(getfield(i,1,false),one($I)),), strides(x))
    ptr = gep(p, li)
    si = ArrayInterface.StrideIndex{1,$(R[ri],),$(C===1 ? 1 : 0)}(
      (getfield(strides(x), $ri, false),), (Zero(),)
    )
    stridedpointer(ptr, si, StaticInt{$(B===1 ? 1 : 0)}())
  end
end
@generated function gespf1(x::AbstractStridedPointer{T,N,C,B,R}, ::Tuple{VectorizationBase.NullStep}) where {T,N,C,B,R}
  ri = argmin(R)
  quote
    $(Expr(:meta,:inline))
    si = ArrayInterface.StrideIndex{1,$(R[ri],),$(C===1 ? 1 : 0)}(
      (getfield(strides(x), $ri, false),),
      (getfield(offsets(x), $ri, false),)
    )
    stridedpointer(pointer(x), si, StaticInt{$(B==1 ? 1 : 0)}())
  end
end
function findfirstcontaining(ref, ind)
  for (i,indr) ∈ enumerate(getindicesonly(ref))
    ind === indr && return i
  end
  0
end
function should_zerorangestart(ls::LoopSet, allarrayrefs::Vector{ArrayReferenceMeta}, name_to_array_map::Vector{Vector{Int}}, isrooted::Vector{Bool})
  loops = ls.loops
  shouldindbyind = fill(false, length(loops))
  for (i,loop) ∈ enumerate(loops)
    ind = loop.itersymbol
    if isloopvalue(ls, ind, isrooted)
      # we don't zero the range if it is used as a loopvalue
      shouldindbyind[i] = true
      continue
    end
    # otherwise, we need
    for namev ∈ name_to_array_map
      baseref = allarrayrefs[first(namev)]
      # firstcontainsind relies on stripping of duplicate inds in parsing
      firstcontainsind = findfirstcontaining(baseref, ind)
      basestride = firstcontainsind == 0 ? 0 : getstrides(baseref)[firstcontainsind]
      allsame = true
      # The idea here is that if any ref to the same array doesn't have `ind`,
      # we can't offset that dimension because different inds will clash.
      # Because offseting the array means counter-offseting the range, we need
      # to be consistent, and check that all arrays are valid first.
      for j ∈ @view(namev[2:end])
        ref = allarrayrefs[j]
        if (firstcontainsind ≠ findfirstcontaining(ref, ind)) || ((firstcontainsind ≠ 0) && (basestride ≠ getstrides(ref)[firstcontainsind]))
          allsame = false
          break
        end
      end
      if !allsame
        shouldindbyind[i] = true
        break
      end
    end
  end
  return shouldindbyind
end
function check_shouldindbyind(ls::LoopSet, ind::Symbol, shouldindbyind::Vector{Bool})
  for (i,loop) ∈ enumerate(ls.loops)
    loop.itersymbol === ind && return shouldindbyind[i]
  end
  true
end

@inline densewrapper(sp, A) = sp
# @inline dummy_ptrarray(sp::AbstractStridedPointer{T,N}, A::AbstractArray{T,N}) where {T,N} = PtrArray(sp, VectorizationBase.zerotuple(Val{N}()), VectorizationBase.val_dense_dims(A))
@inline densewrapper(sp::AbstractStridedPointer{T,N}, A::AbstractArray{T,N}) where {T,N} = _densewrapper(sp, VectorizationBase.val_dense_dims(A))
@inline _densewrapper(sp, ::Nothing) = sp
@inline _densewrapper(sp::AbstractStridedPointer, ::Val{D}) where {D} = VectorizationBase.DensePointerWrapper{D}(sp)

# write a "check_loops_safe_to_zerorangestart
# that will be used to
# 1) decide whether to zerorangestart
# 2) decide whether to gesp that loopstart inside `add_grouped_strided_pointer`
function add_grouped_strided_pointer!(extra_args::Expr, ls::LoopSet)
  allarrayrefs, name_to_array_map, unique_to_name_and_op_map = uniquearrayrefs_csesummary(ls)
  gsp = Expr(:call, lv(:grouped_strided_pointer))
  tgarrays = Expr(:tuple)
  # refs_to_gesp = ArrayReferenceMeta[]
  gespsummaries = Tuple{Int,Vector{Symbol}}[]
  i = 0
  preserve_assignment = Expr(:tuple); preserve = Symbol[];
  @unpack equalarraydims, refs_aliasing_syms = ls
    # duplicate_map = collect(1:length(refs_aliasing_syms))
  duplicate_map = Vector{Int}(undef, length(refs_aliasing_syms))

  # for (i,j) ∈ enumerate(array_refs_with_same_name) # iterate over unique names
  #   ar = allarrayrefs[j]
  #   gespinds = cse_constant_offsets!(ls, allarrayrefs, j, array_refs_with_same_name, arrayref_to_name_op_collection)
  # end
  for (j,ref) ∈ enumerate(refs_aliasing_syms)
    vpref = vptr(ref)
    duplicate = false
    for k ∈ 1:j-1 # quadratic, but should be short enough so that this is faster than O(1) algs
      if vptr(refs_aliasing_syms[k]) === vpref
        duplicate = true
        break
      end
    end
    duplicate && continue
    duplicate_map[j] = (i += 1)
    found = false
    for k ∈ eachindex(allarrayrefs)
      if sameref(allarrayrefs[k], ref)
        gespindsummary = cse_constant_offsets!(ls, allarrayrefs, k, name_to_array_map, unique_to_name_and_op_map)
        push!(gespsummaries, (k, gespindsummary))
        found = true
        break
      end
    end
    @assert found
    push!(preserve, presbufsym(ref.ref.array))
  end
  roots = getroots(ls)
  shouldindbyind = should_zerorangestart(ls, allarrayrefs, name_to_array_map, roots)
  for (k,gespindsummary) ∈ gespsummaries
    ref = allarrayrefs[k]
    gespinds = calcgespinds(ls, ref, gespindsummary, shouldindbyind, name_to_array_map[first(first(unique_to_name_and_op_map[k]))], unique_to_name_and_op_map)
    push!(tgarrays.args, Expr(:call, lv(:densewrapper), Expr(:call, lv(:gespf1), vptr(ref), gespinds), name(ref)))
  end
  push!(gsp.args, tgarrays)
  matcheddims = Expr(:tuple)
  for (vptrs,dims) ∈ equalarraydims
    t = Expr(:tuple)
    for (vp,d) ∈ zip(vptrs,dims)
      _id = findfirst(Base.Fix2(===,vp) ∘ vptr, refs_aliasing_syms)
      _id === nothing && continue
      push!(t.args, Expr(:tuple, duplicate_map[_id], d))
    end
    length(t.args) > 1 && push!(matcheddims.args, t)
  end
  push!(gsp.args, val(matcheddims))
  gsps = gensym!(ls, "#grouped#strided#pointer#")
  push!(extra_args.args, gsps)
  pushpreamble!(ls, Expr(:(=), gsps, Expr(:call, GlobalRef(Core,:getfield), gsp, 1)))
  preserve, shouldindbyind, roots
end

# first_cache() = ifelse(gt(num_cache_levels(), StaticInt{2}()), StaticInt{2}(), StaticInt{1}())
# function _first_cache_size(::StaticInt{FCS}) where {FCS}
#     L1inclusive = StaticInt{FCS}() - VectorizationBase.cache_size(One())
#     ifelse(eq(first_cache(), StaticInt(2)) & VectorizationBase.cache_inclusive(StaticInt(2)), L1inclusive, StaticInt{FCS}())
# end
# _first_cache_size(::Nothing) = StaticInt(262144)
# first_cache_size() = _first_cache_size(cache_size(first_cache()))

@generated function _turbo_config_val(
    ::Val{CNFARG}, ::StaticInt{W}, ::StaticInt{RS}, ::StaticInt{AR}, ::StaticInt{NT},
    ::StaticInt{CLS}, ::StaticInt{L1}, ::StaticInt{L2}, ::StaticInt{L3}
) where {CNFARG,W,RS,AR,CLS,L1,L2,L3,NT}
  inline,u₁,u₂,v,BROADCAST,thread = CNFARG
  nt = min(thread % UInt, NT % UInt)
  t = Expr(:tuple, inline, u₁, u₂, v, BROADCAST, W, RS, AR, CLS, L1, L2, L3, nt)
  length(CNFARG) == 7 && push!(t.args, CNFARG[7])
  Expr(:call, Expr(:curly, :Val, t))
end
@inline function avx_config_val(
    ::Val{CNFARG}, ::StaticInt{W}
) where {CNFARG,W}
    _turbo_config_val(
        Val{CNFARG}(), StaticInt{W}(), register_size(), available_registers(), lv_max_num_threads(),
        cache_linesize(), cache_size(StaticInt(1)), cache_size(StaticInt(2)), cache_size(StaticInt(3))
    )
end
function find_samename_constparent(op::Operation, opname::Symbol)
    for opp ∈ parents(op)
        (((isconstant(opp) && instruction(opp) == LOOPCONSTANT) && (name(opp) === opname))) && return opp
        opptemp = find_samename_constparent(opp, opname)
        opptemp === opp || return opptemp
    end
    op
end
function remove_outer_reducts!(roots::Vector{Bool}, ls::LoopSet)
    ops = operations(ls)
    for or ∈ ls.outer_reductions
        op = ops[or]
        optemp = find_samename_constparent(op, name(op))
        if optemp ≢ op
            roots[identifier(optemp)] = false
        end
    end
end

function split_ifelse!(
  ls::LoopSet, preserve::Vector{Symbol}, shouldindbyind::Vector{Bool}, roots::Vector{Bool}, extra_args::Expr, k::Int,
  inlineu₁u₂::Tuple{Bool,Int8,Int8,Int8}, thread::UInt, debug::Bool
)
  roots[k] = false
  op = operations(ls)[k]
  op.instruction = DROPPEDCONSTANT
  op.node_type = constant
  # we want to eliminate
  parents_op = parents(op)
  condop = first(parents_op)
  # create one loop where `opp` is true, and a second where it is `false`
  prepre = ls.prepreamble; append!(prepre.args, ls.preamble.args)
  ls.prepreamble = Expr(:block); ls.preamble = Expr(:block);
  ls_true = deepcopy(ls)
  lsfalse = ls
  true_ops = operations(ls_true)
  falseops = operations(lsfalse)
  true_op = parents(true_ops[k])[2]
  falseop = parents_op[3]
  true_op.dependencies = loopdependencies(op)
  falseop.dependencies = loopdependencies(op)
  true_op.reduced_children = reducedchildren(op)
  falseop.reduced_children = reducedchildren(op)
  condop_count = 0
  for i ∈ eachindex(falseops)
    fop = falseops[i]
    parents_false = parents(fop)
    for (j,opp) ∈ enumerate(parents_false)
      if opp === op # then ops[i]'s jth parent is the ifelse
        # These reduction to scalar instructions are added for non-outer reductions initialized with non-constant ops
        # So we check if now
        # if (j == 2) && (Base.sym_in(instruction(fop).instr, (:reduced_add, :reduced_prod, :reduced_max, :reduced_min, :reduced_all, :reduced_any)))
          # if isconstantop(true_op)
          #   (true_ops[i]).instruction = Instruction(:identity)
          # end
          # if isconstantop(falseop)
          #   fop.instruction = Instruction(:identity)
          # end
        # end
        parents(true_ops[i])[j] = true_op
        parents_false[j] = falseop
      end
      condop_count += roots[i] & (condop === opp)
    end
  end
  roots[identifier(condop)] &= condop_count > 0
  q = :(if $(name(condop))
      $(generate_call_split(ls_true, preserve, shouldindbyind, roots, copy(extra_args), inlineu₁u₂, thread, debug))
    else
      $(generate_call_split(lsfalse, preserve, shouldindbyind, roots, extra_args, inlineu₁u₂, thread, debug))
    end)
  push!(prepre.args, q)
  prepre
end

function generate_call(ls::LoopSet, inlineu₁u₂::Tuple{Bool,Int8,Int8,Int8}, thread::UInt, debug::Bool)
  extra_args = Expr(:tuple)
  fill_children!(ls)
  preserve, shouldindbyind, roots = add_grouped_strided_pointer!(extra_args, ls)
  generate_call_split(ls, preserve, shouldindbyind, roots, extra_args, inlineu₁u₂, thread, debug)
end
function generate_call_split(
  ls::LoopSet, preserve::Vector{Symbol}, shouldindbyind::Vector{Bool}, roots::Vector{Bool}, extra_args::Expr,
  inlineu₁u₂::Tuple{Bool,Int8,Int8,Int8}, thread::UInt, debug::Bool
)
  for (k,op) ∈ enumerate(operations(ls))
    parents_op = parents(op)
    if (iscompute(op) && (instruction(op).instr === :ifelse)) && (length(parents_op) == 3) && isconstantop(first(parents_op))
      return split_ifelse!(ls, preserve, shouldindbyind, roots, extra_args, k, inlineu₁u₂, thread, debug)
    end
  end
  return generate_call_types(ls, preserve, shouldindbyind, roots, extra_args, inlineu₁u₂, thread, debug)
end

# Try to condense in type stable manner
function generate_call_types(
  ls::LoopSet, preserve::Vector{Symbol}, shouldindbyind::Vector{Bool}, roots::Vector{Bool}, extra_args::Expr,
  (inline,u₁,u₂,v)::Tuple{Bool,Int8,Int8,Int8}, thread::UInt, debug::Bool
)
  # good place to check for split
  operation_descriptions = Expr(:tuple)
  varnames = Symbol[]; ids = Vector{Int}(undef, length(operations(ls)))
  ops = operations(ls)
  for op ∈ ops
    instr::Instruction = instruction(op)
    if (isconstant(op) && (instr == LOOPCONSTANT)) && (!roots[identifier(op)])
      instr = op.instruction = DROPPEDCONSTANT
    end
    push!(operation_descriptions.args, QuoteNode(instr.mod))
    push!(operation_descriptions.args, QuoteNode(instr.instr))
    push!(operation_descriptions.args, OperationStruct!(varnames, ids, ls, op))
  end
  arraysymbolinds = Symbol[]
  arrayref_descriptions = Expr(:tuple)
  duplicate_ref = fill(false, length(ls.refs_aliasing_syms))
  for (j,ref) ∈ enumerate(ls.refs_aliasing_syms)
    vpref = vptr(ref)
    # duplicate_ref[j] ≠ 0 && continue
    duplicate_ref[j] && continue
    push!(arrayref_descriptions.args, ArrayRefStruct(ls, ref, arraysymbolinds, ids))
  end
  argmeta = argmeta_and_consts_description(ls, arraysymbolinds)
  loop_bounds = loop_boundaries(ls, shouldindbyind)
  loop_syms = tuple_expr(QuoteNode, ls.loopsymbols)
  func = debug ? lv(:_turbo_loopset_debug) : lv(:_turbo_!)
  lbarg = debug ? Expr(:call, :typeof, loop_bounds) : loop_bounds
  configarg = (inline,u₁,u₂,v,ls.isbroadcast,thread)
  unroll_param_tup = Expr(:call, lv(:avx_config_val), :(Val{$configarg}()), VECTORWIDTHSYMBOL)
  q = Expr(:call, func, unroll_param_tup, val(operation_descriptions), val(arrayref_descriptions), val(argmeta), val(loop_syms))

  add_reassigned_syms!(extra_args, ls) # counterpart to `add_ops!` constants
  for (opid,sym) ∈ ls.preamble_symsym # counterpart to process_metadata! symsym extraction
    if instruction(ops[opid]) ≠ DROPPEDCONSTANT
      push!(extra_args.args, sym)
    end
  end
  append!(extra_args.args, arraysymbolinds) # add_array_symbols!
  add_external_functions!(extra_args, ls) # extract_external_functions!
  add_outerreduct_types!(extra_args, ls) # extract_outerreduct_types!
  if debug
    vecwidthdefq = Expr(:block)
    push!(q.args, Expr(:tuple, lbarg, extra_args))
  else
    vargsym = gensym(:vargsym)
    vecwidthdefq = Expr(:block, Expr(:(=), vargsym, Expr(:tuple, lbarg, extra_args)))
    push!(q.args, Expr(:call, GlobalRef(Base,:Val), Expr(:call, GlobalRef(Base,:typeof), vargsym)), Expr(:(...), Expr(:call, lv(:flatten_to_tuple), vargsym)))
  end
  define_eltype_vec_width!(vecwidthdefq, ls, nothing, true)
  push!(vecwidthdefq.args, q)
  if debug
    pushpreamble!(ls,vecwidthdefq)
    Expr(:block, ls.prepreamble, ls.preamble)
  else
    setup_call_final(ls, setup_outerreduct_preserve(ls, vecwidthdefq, preserve))
  end
end
# @inline reductinittype(::T) where {T} = StaticType{T}()
typeof_expr(op::Operation) = Expr(:call, GlobalRef(Base,:typeof), name(op))
eltype_expr(op::Operation) = Expr(:call, GlobalRef(Base,:eltype), name(op))
function add_outerreduct_types!(extra_args::Expr, ls::LoopSet) # extract_outerreduct_types!
  for or ∈ ls.outer_reductions
    op = operations(ls)[or]
    if instruction(op).instr ≢ :ifelse
      push!(extra_args.args, eltype_expr(op))
    else
      push!(extra_args.args, name(op))
    end
  end
end
"""
    check_args(::Vararg{AbstractArray})

LoopVectorization will optimize an `@turbo` loop if `check_args` on each on the indexed abstract arrays returns true.
It returns true for `AbstractArray{T}`s when `check_type(T) == true` and the array or its parent is a `StridedArray` or `AbstractRange`.

To provide support for a custom array type, ensure that `check_args` returns true, either through overloading it or subtyping `DenseArray`.
Additionally, define `pointer` and `stride` methods.
"""
@inline function check_args(A::AbstractArray{T}) where {T}
    check_type(T) && check_device(ArrayInterface.device(A))
end
@inline check_args(A::BitVector) = true
@inline check_args(A::BitArray) = iszero(size(A,1) & 7)
@inline check_args(::VectorizationBase.AbstractStridedPointer) = true
@inline function check_args(x)
    # @info "`LoopVectorization.check_args(::$(typeof(x))) == false`, therefore compiling a probably slow `@inbounds @fastmath` fallback loop." maxlog=1
    false
end
@inline check_args(A, B, C::Vararg{Any,K}) where {K} = check_args(A) && check_args(B, C...)
@inline check_args(::AbstractRange{T}) where {T} = check_type(T)
@inline check_args(::UpTri) = false
@inline check_args(::LoTri) = false
@inline check_args(::Diagonal) = false
@inline check_args(::Type{T}) where {T} = check_type(T)
"""
    check_type(::Type{T}) where {T}

Returns true if the element type is supported.
"""
@inline check_type(::Type{T}) where {T <: NativeTypes} = true
@inline check_type(::Type{T}) where {T} = false
@inline check_type(::Type{T}) where {T <: AbstractSIMD} = true
@inline check_device(::ArrayInterface.CPUPointer) = true
@inline check_device(::ArrayInterface.CPUTuple) = true
@inline check_device(x) = false

function check_args_call(ls::LoopSet)
    q = Expr(:call, lv(:check_args))
    append!(q.args, ls.includedactualarrays)
    for r ∈ ls.outer_reductions
        push!(q.args, Expr(:call, :typeof, name(ls.operations[r])))
    end
    q
end

make_fast(q) = Expr(:macrocall, Symbol("@fastmath"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), q)
make_crashy(q) = Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__,Symbol(@__FILE__)), q)

@inline vecmemaybe(x::NativeTypes) = x
@inline vecmemaybe(x::VectorizationBase._Vec) = Vec(x)
@inline vecmemaybe(x::Tuple) = VectorizationBase.VecUnroll(x)
@inline vecmemaybe(x::Mask) = x

function gc_preserve(call::Expr, preserve::Vector{Symbol})
    q = Expr(:gc_preserve, call)
    append!(q.args, preserve)
    q
end

# @generated function ifelse_reduce(f::F, x::Vec{W,T}) where {F,T,W}
#   Wt = W >>> 1
#   # uw = Symbol(:x_, Wt)
#   # lw = Symbol(:x_, Wt)
#   xw = Symbol(:x_, Wt)
#   q = Expr(:block, Expr(:meta,:inline), :($xw = f($(VectorizationBase.uppervector)(x), $(VectorizationBase.lowervector)(x))))
#   while Wt > 2
#     Wt >>>= 1
#     ex = :(f($(VectorizationBase.uppervector)($xw), $(VectorizationBase.lowervector)($xw)))
#     xw = Symbol(:x_, Wt)
#     push!(q.args, :($xw = $ex))
#   end

# end

# function setup_call_inline(ls::LoopSet, inline::Bool, u₁::Int8, u₂::Int8, thread::Int)
#   call, preserve = generate_call_split(ls, (inline,u₁,u₂), thread % UInt, false)
#   setup_call_ret!(ls, call, preserve)
# end
setup_outerreduct_preserve_mangler(op::Operation) = Symbol(mangledvar(op), "##onevec##")
function setup_outerreduct_preserve(ls::LoopSet, call::Expr, preserve::Vector{Symbol})
  iszero(length(ls.outer_reductions)) && return gc_preserve(call, preserve)
  retv = loopset_return_value(ls, Val(false))
  q = Expr(:block, gc_preserve(Expr(:(=), retv, call), preserve))
  for or ∈ ls.outer_reductions
    op = ls.operations[or]
    var = name(op)
    # push!(call.args, Symbol("##TYPEOF##", var))
    instr = instruction(op)
    out = setup_outerreduct_preserve_mangler(op)
    reducq = if instr.instr ≢ :ifelse
      Expr(:call, reduction_scalar_combine(op), Expr(:call, lv(:vecmemaybe), out), var)
    else
      opinstr = ifelse_reduction(:IfElseReduced, op) do opv
        opvname = name(opv)
        oporig = gensym(opvname)
        pushfirst!(q.args, Expr(:(=), oporig, opvname))
        Expr(:call, lv(:vecmemaybe), setup_outerreduct_preserve_mangler(opv)), (oporig,)
      end
      Expr(:call, opinstr, Expr(:call, lv(:vecmemaybe), out), var)
    end
    push!(q.args, Expr(:(=), var, reducq))
  end
  q
end
function setup_call_final(ls::LoopSet, q::Expr)
  pushpreamble!(ls, q)
  push!(ls.preamble.args, nothing)
  return ls.preamble
end
function setup_call_debug(ls::LoopSet)
  generate_call(ls, (false,zero(Int8),zero(Int8),zero(Int8)), zero(UInt), true)
end
function setup_call(
  ls::LoopSet, q::Expr, source::LineNumberNode, inline::Bool, check_empty::Bool, u₁::Int8, u₂::Int8, v::Int8, thread::Int, warncheckarg::Int
)
  # We outline/inline at the macro level by creating/not creating an anonymous function.
  # The old API instead was based on inlining or not inline the generated function, but
  # the generated function must be inlined into the initial loop preamble for performance reasons.
  # Creating an anonymous function and calling it also achieves the outlining, while still
  # inlining the generated function into the loop preamble.
  lnns = extract_all_lnns(q)
  pushfirst!(lnns, source)
  call = generate_call(ls, (inline, u₁, u₂, v), thread%UInt, false)
  call = check_empty ? check_if_empty(ls, call) : call
  argfailure = make_crashy(make_fast(q))
  if warncheckarg ≠ 0
    warnstring = "$(first(lnns)):\n`LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead.\nUse `warn_check_args=false`, e.g. `@turbo warn_check_args=false ...`, to disable this warning."
    warning = :(@warn $warnstring)
    warncheckarg > 0 && push!(warning.args, :(maxlog=$warncheckarg))
    argfailure = Expr(:block,  warning, argfailure)
  end
  pushprepreamble!(ls, Expr(:if, check_args_call(ls), call, argfailure))
  prepend_lnns!(ls.prepreamble, lnns)
  return ls.prepreamble
end

ext: jl


text:
"""
This  package deals with cyclotomic numbers,  the complex numbers which are
linear  combinations  of  roots  of  unity  with  rational coefficients. It
depends on the packages `ModuleElt` and `Primes`. It can also handle linear
combinations  of roots of  unity with coefficients  in an arbitrary ring of
real  numbers like  `Floats`. The  cyclotomic numbers  can be  converted to
`Complex` numbers.

The  cyclotomic numbers form a field, the  cyclotomic field. It is also the
maximal  extension of the rationals which  has an abelian Galois group. Its
ring of integers are the sums of roots of unity with integer coefficients.

Cyclotomic  numbers are very  important for finite  groups, since character
values of finite groups are cyclotomic integers.

This package is a port of the GAP implementation of cyclotomics, which uses
a  normal form  given by  writing them  in the  Zumbroich basis.  This form
allows to find the smallest Cyclotomic field which contains a given number,
and  decide in particular  if a cyclotomic  is zero. Let ζₙ=exp(2iπ/n). The
Zumbroich  basis is  a particular  subset of  size φ(n) of 1,ζₙ,ζₙ²,…,ζₙⁿ⁻¹
which forms a basis of ℚ (ζₙ).  The reference is

T. Breuer, Integral bases for subfields of cyclotomic fields AAECC 8 (1997)

I  started  this  file  by  porting  Christian  Stump's Sage code, which is
simpler  to understand than GAP's C  code.

As GAP does, I lower automatically numbers after each computation (that is,
reduce  them to the smallest cyclotomic field where they belong). Currently
the  code is somewhat  slower (depending on  the operation it  has the same
speed or is slower up to 50%) than the C code in GAP but there are probably
opportunities to optimize that I missed.

What GAP does which I do not do is convert automatically a Cyclotomic which
is  rational to a Rational,  a Rational which is  integral to an Integer, a
BigInt  which is  small to  an Int,  etc… This  is a tremendously important
optimization  but because of type stability in Julia it needs a new type of
number to be added to Julia, which I am not competent enough to try.

This package is similar (and mostly compatible) with Marek Kaluba's package
`Cyclotomics`,  whose existence I discovered after writing this package. We
discussed merging them but concluded it would be a lot of work for benefits
which are not clear currently. Some differences are:

  - I define two types in this  package: `Root1` represents a root of unity,
    and  `Cyc` a cyclotomic number. The advantage of having a separate type
    for  roots of  unity is  that computations  are very  fast for them, of
    which  I take advantage  in the package  `CycPol` for polynomials whose
    zeros are roots of unity.

  - In Kaluba's package  numbers are not  systematically lowered but only on
    demand  (like  for  printing).  this  speeds  up some computations by a
    factor  approaching 2, but it also makes  some computations I had to do
    infeasible,  like the following one (which if not lowering involves too
    large fields); the answer is `-36ζ₃²`:

```julia-rep1
julia> prod(x->1-x,[E(3),E(3,2),E(6),E(6,5),E(8),E(8),E(8,5),E(8,7),E(9,2),
E(9,5),E(9,8),E(12,7),E(12,11),E(16),E(16,3),E(16,5),E(16,9),E(16,11),E(16,13),
E(18,5),E(18,5),E(18,11),E(18,11),E(18,17),E(18,17),E(21,2),E(21,5),E(21,8),
E(21,11),E(21,17),E(21,20),
E(27,2),E(27,5),E(27,8),E(27,11),E(27,14),E(27,17),E(27,20),E(27,23),E(27,26),
E(32,7),E(32,15),E(32,23),E(32,31),E(39),E(39,4),E(39,7),E(39,10),E(39,16),
E(39,19),E(39,22),E(39,25),E(39,28),E(39,31),E(39,34),E(39,37),E(42),E(42,13),
E(42,19),E(42,25),E(42,31),E(42,37),E(48,11),E(48,19),E(48,35),E(48,43),E(60,7),
E(60,19),E(60,31),E(60,43),E(78,5),E(78,11),E(78,17),E(78,23),E(78,29),E(78,35),
E(78,41),E(78,47),E(78,53),E(78,59),E(78,71),E(78,77),E(80,7),E(80,23),E(80,31),
E(80,39),E(80,47),E(80,63),E(80,71),E(80,79),E(88,3),E(88,19),E(88,27),E(88,35),
E(88,43),E(88,51),E(88,59),E(88,67),E(88,75),E(88,83),E(90),E(90,7),E(90,13),
E(90,19),E(90,31),E(90,37),E(90,43),E(90,49),E(90,61),E(90,67),E(90,73),
E(90,79),E(96,5),E(96,13),E(96,29),E(96,37),E(96,53),E(96,61),E(96,77),E(96,85),
E(104),E(104,9),E(104,17),E(104,25),E(104,33),E(104,41),E(104,49),E(104,57),
E(104,73),E(104,81),E(104,89),E(104,97),E(144),E(144,17),E(144,25),E(144,41),
E(144,49),E(144,65),E(144,73),E(144,89),E(144,97),E(144,113),E(144,121),
E(144,137),E(152,5),E(152,13),E(152,21),E(152,29),E(152,37),E(152,45),
E(152,53),E(152,61),E(152,69),E(152,77),E(152,85),E(152,93),E(152,101),
E(152,109),E(152,117),E(152,125),E(152,141),E(152,149),E(204,11),E(204,23),
E(204,35),E(204,47),E(204,59),E(204,71),E(204,83),E(204,95),E(204,107),
E(204,131),E(204,143),E(204,155),E(204,167),E(204,179),E(204,191),E(204,203)])
```
If you `develop` my package it is easy to use the strategy of lowering only
on  demand or to use alternate implementations like dense vectors or sparse
vectors (like `Cyclotomics`) --- I prepared boolean flags to choose various
implementations in the code. I choosed the implementation with `ModuleElts`
and systematic lowering as giving the best results.

The main way to build a Cyclotomic number is to use the function `E(n,k=1)`
which   constructs  the  `Root1`  equal  to   `ζₙᵏ`,  and  to  make  linear
combinations of such numbers.

# Examples
```julia-repl
julia> E(3,2) # a root of unity
Root1: ζ₃²

julia> E(3)+E(4) # nice display at the repl
Cyc{Int64}: ζ₁₂⁴-ζ₁₂⁷-ζ₁₂¹¹
```
```julia-rep1
julia> print(E(3)+E(4)) # otherwise give output which can be parsed back
E(12,4)-E(12,7)-E(12,11)
```
```julia-repl
julia> E(12,11)-E(12,7) # square roots of integers are recognized on output
Cyc{Int64}: √3

# but you can prevent that recognition
julia> repr(E(12,11)-E(12,7),context=(:limit=>true,:quadratic=>false))
"-ζ₁₂⁷+ζ₁₂¹¹"

julia> a=E(3)+E(3,2)
Cyc{Int64}: -1

julia> conductor(a) # a has been lowered to ℚ (ζ₁)=ℚ
1

julia> typeof(Int(a))
Int64
```
```julia-rep1
julia> Int(E(4))
ERROR: InexactError: convert(Int64, E(4))
```
```julia-repl
julia> inv(1+E(4)) # like for Integers inverse involves floats
Cyc{Float64}: 0.5-0.5ζ₄

julia> 1//(1+E(4))  # but not written this way
Cyc{Rational{Int64}}: (1-ζ₄)/2

julia> Cyc(1//2+im) # we can convert Gaussian rationals to cyclotomics
Cyc{Rational{Int64}}: (1+2ζ₄)/2

julia> conj(1+E(4)) # complex conjugate
Cyc{Int64}: 1-ζ₄

julia> real(E(3))  # real part
Cyc{Rational{Int64}}: -1/2

julia> Rational{Int}(real(E(3)))
-1//2

julia> imag(E(3))  # imaginary part
Cyc{Rational{Int64}}: √3/2

julia> c=Cyc(E(9))   # the normal form in the Zumbroich basis has two terms
Cyc{Int64}: -ζ₉⁴-ζ₉⁷

julia> Root1(c) #  but you can convert back to Root1 if possible
Root1: ζ₉

julia> Root1(1+E(4)) # the constructor Root1 returns nothing for a non-root
```

The  group of  roots of  unity is  isomorphic to  ℚ /ℤ  , thus  `Root1` are
represented internally by a rational number in `[0,1[`.

```julia-repl
julia> Root1(;r=1//4) # this constructor ensures the fraction is in [0,1[
Root1: ζ₄

julia> c=E(4)*E(3) # fast computation if staying inside roots of unity
Root1: ζ₁₂⁷
```
`Root1` have the same operations as `Cyc`, but are first converted to `Cyc`
for  any  operation  other  than  `one,  isone,  *, ^, inv, conj, /, //`. A
`Root1`  can be  raised to  a `Rational`  power, which  extracts a  root if
needed.  A `Root1` can be  taken apart using `order`  and `exponent` --- if
`a`  and `b` are prime to each other,  `a` is the order of `E(a,b)` and `b`
(taken  mod `a`) is the exponent. Note  that the order is not the conductor
since `E(6)==-E(3)` has order 6 and conductor 3.

```julia-repl
julia> c=Complex{Float64}(E(3))  # convert to Complex{float} is sometimes useful
-0.4999999999999999 + 0.8660254037844387im
```
In  presence  of  a  `Cyc`  a  number  `<:Real`  or  `<:Complex{<:Real}` is
converted to a `Cyc`.

```julia-repl
julia> 0.0+E(3)
Cyc{Float64}: 1.0ζ₃

julia> E(3)+1//2
Cyc{Rational{Int64}}: √-3/2

julia> E(3)+im
Cyc{Int64}: ζ₁₂⁴-ζ₁₂⁷-ζ₁₂¹¹
```
The  function  `complex`  converts  a  `Cyc{T}`  to  a  `Complex{T}` if the
conductor is 1 or 4, to a `Complex{float(T)}` otherwise.

```julia-repl
julia> complex(E(4))
0 + 1im

julia> complex(E(3))
-0.4999999999999999 + 0.8660254037844387im
```
`Cyc`s have methods `copy, hash, ==, cmp, isless` (total order) so they can
be  keys in hashes, elements of sets,  and can be sorted. Cyclotomics which
are  integers  or  rationals  compare  correctly  to `Real`s (this does not
extend to irrational real `Cyc`s):

```julia-repl
julia> -1<Cyc(0)<1
true
```

`Cyc`s have the operations `+, -, *, /, //, inv, ^, conj, abs2, abs, image,
real, reim, isinteger, isreal, one, isone, zero, iszero, complex, adjoint`.
Cyclotomics   with  rational  or  integer  coefficients  have  the  methods
`numerator`  and  `denominator`:  a  `Cyc`  `x`  is a Cyclotomic integer if
`denominator(x)==1`   and  then  `numerator(x)`   gives  the  corresponding
`Cyc{<:Integer}`.

You  can pick apart a cyclotomic in various ways. The fastest is to use the
iterator  `pairs` which, for a cyclotomic  `a` of conductor `e` iterates on
the  pairs `(i,c)` such that  `a` has a non-zero  coefficient `c` on `ζₑⁱ`.
You  can also get  the coefficient of  `ζₑⁱ` as `a[i]`  but it is slower to
iterate  on coefficients  this way.  Finally you  can get (efficiently) the
vector of all coefficients by `coefficients`.

```julia-repl
julia> a=E(3)+E(4)
Cyc{Int64}: ζ₁₂⁴-ζ₁₂⁷-ζ₁₂¹¹

julia> collect(pairs(a))
3-element Vector{Pair{Int64, Int64}}:
  4 => 1
  7 => -1
 11 => -1

julia> a[6],a[7]
(0, -1)

julia> coefficients(a)
12-element Vector{Int64}:
  0
  0
  0
  0
  1
  0
  0
 -1
  0
  0
  0
 -1

julia> valtype(a) # the type of the coefficients of a
Int64
```

For more information see the docstring for the methods Quadratic,
galois, root, conjugates.

Finally, a benchmark:

```benchmark
julia> function testmat(p)
         ss=[[i,j] for i in 0:p-1 for j in i+1:p-1]
         [(E(p,i'*reverse(j))-E(p,i'*j))//p for i in ss,j in ss]
       end
testmat (generic function with 1 method)

julia> @btime CyclotomicNumbers.testmat(12)^2;  # on Julia 1.7.2
  287.781 ms (3774785 allocations: 279.66 MiB)
```
The equivalent in GAP:

```
testmat:=function(p)local ss;ss:=Combinations([0..p-1],2);
  return List(ss,i->List(ss,j->(E(p)^(i*Reversed(j))-E(p)^(i*j))/p));
end;
```
testmat(12)^2 takes 0.31s in GAP3, 0.22s in GAP4
"""
module CyclotomicNumbers
export coefficients, root, E, Cyc, conductor, galois, Root1, Quadratic,
       order, conjugates, modZ

#---- formatting utilities duplicated here to avoid dependency ----------
const ok="([^-+*/]|√-|{-)*"
const par="(\\([^()]*\\))"
const nobf=Regex("^[-+]?$ok$par*$ok(/+)?[0-9]*\$")
const nob=Regex("^[-+]?$ok$par*$ok\$")

function bracket_if_needed(c::String;allow_frac=false)
  e=allow_frac ? nobf : nob
  if match(e,c)!==nothing c
  else "("*c*")"
  end
end

function format_coefficient(c::String;allow_frac=false)
  if c=="1" ""
  elseif c=="-1" "-"
  else bracket_if_needed(c;allow_frac)
  end
end

function stringexp(io::IO,n::Integer)
  if isone(n) ""
  elseif get(io,:TeX,false)
    "^"*(n in 0:9 ? string(n) : "{"*string(n)*"}")
  elseif get(io,:limit,false)
    if n<0 res=['⁻']; n=-n else res=Char[] end
    for i in reverse(digits(n))
      push!(res,['⁰','¹','²','³','⁴','⁵','⁶','⁷','⁸','⁹'][i+1])
    end
    String(res)
  else "^"*string(n)
  end
end

function stringind(io::IO,n::Integer)
  if get(io,:TeX,false)
    n in 0:9 ? "_"*string(n) : "_{"*string(n)*"}"
  elseif get(io,:limit,false)
    if n<0 res=['₋']; n=-n else res=Char[] end
    for i in reverse(digits(n)) push!(res,Char(0x2080+i)) end
    String(res)
  else "_"*string(n)
  end
end

#---- number theory utilities duplicated here to avoid dependency ----------
"`CyclotomicNumbers.prime_residues(n)` the numbers less than `n` and prime to `n` "
function prime_residues(n)
  if n==1 return [0] end
  pp=trues(n-1)
  for i in 2:div(n,2)
    if !pp[i] || n%i!=0 continue end
    pp[i:i:n-1].=false
  end
  (1:n-1)[pp]
end

using Primes: Primes

if isdefined(Primes,:eachfactor)
  const EF=Primes.eachfactor
else
  const fact_dict=Dict(2=>Primes.factor(2))
  EF(n)=get!(()->Primes.factor(n),fact_dict,n)
end

#------------------------ type Root1 ----------------------------------
struct Root1 <: Number # E(c,n)
  r::Rational{Int}
  global Root1_(x)=new(x)
end

"`modZ(x::Rational{<:Integer})` returns `x mod ℤ ` as a rational in `[0,1[`"
modZ(x::Rational{<:Integer})=Base.unsafe_rational(mod(numerator(x),
                                    denominator(x)),denominator(x))

" `E(n,p=1)` returns the `Root1` equal to `ζₙᵖ`"
E(c,n=1)=Root1_(modZ(n//c))

Base.exponent(a::Root1)=numerator(a.r)
order(a::Root1)=denominator(a.r)
conductor(a::Root1)=order(a)%4==2 ? div(order(a),2) : order(a)

Root1(;r::Rational{<:Integer}=0//1)=Root1_(modZ(r))

function Root1(c::Real)
  if c==1 Root1_(0//1)
  elseif c==-1 Root1_(1//2)
  elseif c==0 Cyc(0)
  else nothing
  end
end

Base.broadcastable(r::Root1)=Ref(r)

function Base.show(io::IO, ::MIME"text/plain", r::Root1)
  if !haskey(io,:typeinfo) print(io,"Root1: ") end
  show(io,r)
end

function Base.show(io::IO, r::Root1)
  repl=get(io,:limit,false)
  TeX=get(io,:TeX,false)
  d=exponent(r)
  c=order(r)
  if repl || TeX
    if c==1 print(io,"1")
    elseif c==2 print(io,"-1")
    else print(io,TeX ? "\\zeta" : "ζ",stringind(io,c),stringexp(io,d))
    end
  else
    print(io,(d==1 || d==0) ? "E($c)" : "E($c,$d)")
  end
end

function Base.cmp(a::Root1,b::Root1)
  r=cmp(conductor(a),conductor(b))
  if !iszero(r) return r end
  cmp(Cyc(a),Cyc(b))
end

Base.isless(a::Root1,b::Root1)=cmp(a,b)==-1
Base.isless(c::Root1,d::Real)=Cyc(c)<Cyc(d)
Base.isless(d::Real,c::Root1)=Cyc(d)<Cyc(c)
Base.:(==)(a::Root1,b::Root1)=a.r==b.r
Base.one(a::Root1)=Root1_(0//1)
Base.zero(a::Root1)=zero(Cyc{Int})
Base.isone(a::Root1)=iszero(a.r)
Base.:*(a::Root1,b::Root1)=Root1_(modZ(a.r+b.r))

Base.:^(a::Root1,n::Integer)=Root1_(modZ(n*a.r))
Base.:^(a::Root1,r::Rational{<:Integer})=root(a,denominator(r))^numerator(r)
Base.inv(a::Root1)=Root1_(modZ(-a.r))
Base.conj(a::Root1)=inv(a)
Base.:/(a::Root1,b::Root1)=a*inv(b)
Base.://(a::Root1,b::Root1)=a/b

#------------------------ type Cyc ----------------------------------
const impl=:MM # I tried 4 different implementations. For testmat(12)^2
    # MM=ModuleElt is fastest
    # MM=HModuleElt is 50% slower than ModuleElt
    # :svec is 25% slower than ModuleElt
    # :vec is 45% slower than ModuleElt
const lazy=false # whether to lower all the time or on demand

if impl==:vec
struct Cyc{T <: Real}<: Number   # a cyclotomic number
  d::Vector{T} # the i-th element is the coefficient on ζⁱ⁻¹
  global function Cyc_(d::Vector{T}) where T<:Real
    new{T}(d)
  end
end
function Cyc(c::Integer,v::AbstractVector)
  if c!=length(v) error("c=$c but length(v)=$(length(v))\n") end
  Cyc_(v)
end
Base.pairs(c::Cyc)=Iterators.map(x->x[1]-1=>x[2],
                              Iterators.filter(x->x[2]!=0,enumerate(c.d)))
elseif impl==:svec
using SparseArrays
mutable struct Cyc{T <: Real}<:Number    # a cyclotomic number
  d::SparseVector{T,Int} # d[i]==coeff on ζⁱ⁻¹ (where i∈ zumbroich_basis(n))
  global function Cyc_(d::SparseVector{T}) where T<:Real
#   for debugging you may uncomment the following line
#   if !issorted(d.nzind) || any(iszero,d.nzval) error(d) end
    new{T}(d)
  end
end
Base.pairs(c::Cyc)=Iterators.map(x->x[1]-1=>x[2],zip(c.d.nzind,c.d.nzval))
function Cyc(c::Integer,v::SparseVector)
  if c!=length(v) error("c=$c but length(v)=$(length(v))\n") end
  Cyc_(v)
end

elseif impl==:MM
using ModuleElts
const MM=ModuleElt # you can try with HModuleElt
mutable struct Cyc{T <: Real}<: Number   # a cyclotomic number
  n::Int
  d::MM{Int,T} # list of pairs: i=>coeff on ζⁱ (where i∈ zumbroich_basis(n))
end
Base.pairs(c::Cyc)=c.d
end

if impl==:MM
  conductor(c::Cyc)=c.n
  Base.getindex(c::Cyc,i::Integer)=c.d[mod(i,conductor(c))]
else
  conductor(c::Cyc)=length(c.d)
  Base.getindex(c::Cyc,i::Integer)=c.d[mod(i,conductor(c))+1]
end

Base.valtype(c::Cyc{T}) where T =T # how to recover T from c

"""
   `conductor(c::Cyc)`
   `conductor(a::AbstractArray)`

returns the smallest positive integer  n such that `c∈ ℚ (ζₙ)` (resp. all
elements of `a` are in `ℚ (ζₙ)`).

```julia-repl
julia> conductor(E(9))
9

julia> conductor([E(3),1//2,E(4)])
12
```
"""
conductor(a::AbstractArray)=lcm(conductor.(a))
conductor(i::Integer)=1 # for convenience
conductor(i::Rational{<:Integer})=1 # for convenience

"""
  CyclotomicNumbers.zumbroich_basis(n::Int)

  returns  the Zumbroich basis of  ℚ (ζₙ) as the  vector of i in 0:n-1 such
  that `ζₙⁱ` is in the basis.
"""
function zumbroich_basis(n::Int)
#  Note that this function is not used in the rest of the code. We use Elist.
  if n==1 return [0] end
  function J(k::Int, p::Int) # see [Breuer] Rem. 1 p. 283
    if k==0 if p==2 return 0:0 else return 1:p-1 end
    elseif p==2 return 0:1
    else return div(1-p,2):div(p-1,2)
    end
  end
  nfact=EF(n)
  res=
  let J=J
    [[div(n*i,p^k) for i in J(k-1,p)] for (p,np) in nfact for k in 1:np]
  end
  v=sum.(vec(collect(Iterators.product(res...))))
  sort(v.%n)
end

"""
`coefficients(c::Cyc)`

for  a cyclotomic `c` of conductor `n`,  returns a vector `v` of length `n`
such that `c==∑ᵢ vᵢ₋₁ ζⁱ`.

```julia-repl
julia> coefficients(Cyc(E(9)))
9-element Vector{Int64}:
  0
  0
  0
  0
 -1
  0
  0
 -1
  0
```
"""
function coefficients(c::Cyc)
if impl==:svec return Array(c.d)
elseif impl==:vec return copy(c.d)
else
  res=zeros(valtype(c),conductor(c))
  @inbounds for (i,v) in pairs(c) res[i+1]=v end
  res
end
end

"""
`denominator(c::Cyc{Rational})`

returns the smallest `d` such that `d*c` has integral coefficients (thus is
an algebraic integer).
"""
Base.denominator(c::Cyc)=lcm(denominator.(values(c.d)))

"""
`numerator(c::Cyc{Rational})`

returns `denominator(c)*c` as a cyclotomic over the integers.
"""
Base.numerator(c::Cyc{<:Union{T,Rational{T}}}) where T<:Integer=Cyc{T}(c*denominator(c))

if impl==:vec
Cyc(i::Real)=Cyc_([i])
elseif impl==:svec
Cyc(i::Real)=Cyc_(iszero(i) ? spzeros(typeof(i),1) : SparseVector(1,[1],[i]))
else
Cyc(i::Real)=Cyc(1,MM(i==0 ? Pair{Int,typeof(i)}[] : [0=>i];check=false))
end
Cyc{T}(i::Real) where T<:Real=Cyc(T(i))

Cyc{T}(a::Root1) where T<:Real=Cyc{T}(Cyc(a))

Base.zero(c::Cyc{T}) where T=Cyc{T}(0)

if impl==:svec obviouslyzero(c::Cyc)=nnz(c.d)==0 # faster than iszero(c.d)
else           obviouslyzero(c::Cyc)=iszero(c.d)
end
if lazy Base.iszero(c::Cyc)=obviouslyzero(lower!(c))
else    Base.iszero(c::Cyc)=obviouslyzero(c)
end
Base.zero(::Type{Cyc{T}}) where T=Cyc{T}(0)
Base.one(c::Cyc{T}) where T =Cyc{T}(1)
Base.isone(c::Cyc)=isone(conductor(c)) && isone(num(c))

function Cyc(c::Complex)
  if iszero(imag(c)) return Cyc(real(c)) end
if impl==:vec
  Cyc(4,[real(c),imag(c),zero(real(c)),zero(real(c))])
elseif impl==:svec
  Cyc_(dropzeros!(SparseVector(4,[1,2],[real(c),imag(c)])))
else
  if iszero(real(c)) Cyc(4,MM(1=>imag(c);check=false))
  else Cyc(4,MM(0=>real(c),1=>imag(c);check=false))
  end
end
end

Cyc{T}(c::Complex) where T=T(real(c))+E(4)*T(imag(c))

function Cyc{T}(c::Cyc{T1}) where {T,T1}
  if T==T1 return c end
  if impl==:MM Cyc(conductor(c),convert(MM{Int,T},c.d))
  else Cyc(conductor(c),T.(c.d))
  end
end

# num(c): value of c as a real when conductor(c)==1
if impl==:MM
  num(c::Cyc{T}) where T =iszero(c) ? zero(T) : last(first(c.d))
else
  num(c::Cyc)=c.d[1]
end

function (::Type{T})(c::Cyc)where T<:Union{Integer,Rational}
  @static if lazy lower!(c) end
  if conductor(c)==1 return T(num(c)) end
  throw(InexactError(:convert,T,c))
end

function Base.convert(::Type{T},c::Cyc;check=true)where T<:AbstractFloat
  if check && !isreal(c) throw(InexactError(:convert,T,c)) end
  real(convert(Complex{T},c))
end

function (::Type{T})(c::Cyc)where T<:AbstractFloat
  convert(T,c;check=true)
end

(::Type{T})(a::Root1) where T<:Number = T(Cyc(a))

function Complex{T}(c::Cyc)where T<:AbstractFloat
  sum(((k,v),)->T(v)*cispi(2*T(k)/conductor(c)),pairs(c);init=Complex{T}(0.0))
end

function Complex{T}(c::Cyc)where T<:Union{Integer,Rational}
  @static if lazy lower!(c) end
  if conductor(c)==1 return Complex{T}(num(c)) end
  if conductor(c)==4
    res=Complex{T}(0)
    for (k,v) in pairs(c)
      res+=k==0 ? v : im*v
    end
    return res
  end
  throw(InexactError(:convert,Complex{T},c))
end

Complex{T}(a::Root1) where T=Complex{T}(Cyc(a))
function Base.complex(c::Cyc{T}) where T
  @static if lazy lower!(c) end
  (conductor(c)==1||conductor(c)==4) ? Complex{T}(c) : Complex{float(T)}(c)
end

Base.complex(a::Root1)=complex(Cyc(a))

function Base.isinteger(c::Cyc)
  @static if lazy lower!(c) end
  conductor(c)==1 && isinteger(num(c))
end

Base.isinteger(a::Root1)=a.r==0//1 || a.r==1//2

Base.isreal(c::Cyc)=conductor(c)==1 || c==conj(c)
Base.isreal(a::Root1)=isinteger(a)

function Base.real(c::Cyc)
  @static if lazy lower!(c) end
  if conductor(c)==1 return num(c) end
  (c+conj(c))*1//2
end

Base.real(a::Root1)=real(Cyc(a))

function Base.imag(c::Cyc)
  @static if lazy lower!(c) end
  if conductor(c)==1 return 0 end
  E(4)*(conj(c)-c)*1//2
end

Base.imag(a::Root1)=imag(Cyc(a))
Base.reim(c::Cyc)=(real(c),imag(c))
Base.abs2(a::Root1)=abs2(Cyc(a))
Base.abs(a::Root1)=one(a)

# memoize Elist
const Elist_dict=Dict{Tuple{Int,Int},Pair{Bool,Vector{Int}}}((1,0)=>(true=>[0]))
"""
  CyclotomicNumbers.Elist(n,i)

  expresses  ζₙⁱ  in  zumbroich_basis(n):  it  is  a  sum  of some ζₙʲ with
  coefficients all 1 or all -1. The result is a Pair sgn=>inds where sgn is
  true  if coefficients are all 1 and false otherwise, and inds is the list
  of i in 0:n-1 such that ζₙⁱ occurs with a non-zero coefficient.

  Should only be called with i∈ 0:n-1
"""
function Elist(n::Int,i::Int)
  get!(Elist_dict,(n,i)) do
    mp=Int[]
    j=i
    for (p,np) in EF(n)
      f=p^np
      m=div(n,f)
      cnt=mod(j*invmod(m,f),f)
      j-=cnt*m
      if p==2
        if 1==div(cnt,div(f,p)) push!(mp,p) end
      else
        tmp=zeros(Int,np)
        for k in 1:np
          f=div(f,p)
          tmp[k],cnt=divrem(cnt,f)
        end
        for k in np-1:-1:1
          if tmp[k+1]>div(p-1,2)
            tmp[k]+=1
            if k==1 && tmp[k]==p tmp[k]=0 end
          end
        end
        if tmp[1]==0 push!(mp,p) end
      end
    end
    if isempty(mp) return true=>[i] end
    v=vec(sum.(Iterators.product((div(n,p)*(1:p-1) for p in mp)...)))
    v=sort((i.+v).%n)
    iseven(length(mp))=>v
  end
end

# p iterator of pairs i=>c meaning c*E(n,i)
# This constructor is not in API since the result may or may not need lowering
function Cyc(n::Integer,::Type{T},p) where T
  res=if     impl==:vec  zeros(T,n)
      elseif impl==:svec spzeros(T,n)
      elseif impl==:MM   Pair{Int,T}[]
      end
  for (i,c) in p
    (s,v)=Elist(n,mod(i,n))
    if !s c=-c end
    if impl==:MM for k in v push!(res,k=>c) end
    else @inbounds view(res,v.+1).+=c
    end
  end
if impl==:MM Cyc(n,MM(res))
elseif impl==:svec Cyc(n,dropzeros!(res))
else Cyc(n,res)
end
end

function Cyc(a::Root1) # the result is guaranteed lowered
  n=order(a)
  e=exponent(a)
  if n%4==2 return -E(div(n,2),div(e,2)+div(n+2,4)) end
  Cyc(n,Int,(e=>1,))
end

function Base.promote_rule(a::Type{Cyc{T1}},b::Type{T2})where {T1,T2<:Real}
  Cyc{promote_type(T1,T2)}
end

function Base.promote_rule(a::Type{Cyc{T1}},b::Type{Complex{T2}})where {T1,T2<:Real}
  Cyc{promote_type(T1,T2)}
end

function Base.promote_rule(a::Type{Cyc{T1}},b::Type{Cyc{T2}})where {T1,T2}
  Cyc{promote_type(T1,T2)}
end

# total order is necessary to put Cycs in a sorted list
# for conductor(c)==1  a<b is as expected
function Base.cmp(a::Cyc,b::Cyc)
  @static if lazy lower!(a);lower!(b) end
  t=cmp(conductor(a),conductor(b))
  if !iszero(t) return t end
  conductor(a)==1 ? cmp(num(a),num(b)) : cmp(a.d,b.d)
end

if lazy
function Base.:(==)(a::Cyc,b::Cyc)
  if conductor(a)==conductor(b) return a.d==b.d end
  lower!(a);lower!(b)
  conductor(a)==conductor(b) && a.d==b.d
end
else
Base.:(==)(a::Cyc,b::Cyc)=conductor(a)==conductor(b) && a.d==b.d
end
Base.isless(a::Cyc,b::Cyc)=cmp(a,b)==-1
Base.isless(c::Cyc,d::Real)=c<Cyc(d)
Base.isless(d::Real,c::Cyc)=Cyc(d)<c
Base.isless(a::Root1,b::Cyc)=Cyc(a)<b
Base.isless(a::Cyc,b::Root1)=a<Cyc(b)

# hash is necessary to put Cycs as keys of a Dict or make a Set
function Base.hash(a::Cyc, h::UInt)
  @static if lazy lower!(a) end
  hash(a.d, hash(conductor(a), h))
end

function Base.show(io::IO, ::MIME"text/html", a::Cyc)
  print(io, "\$")
  show(IOContext(io,:TeX=>true),a)
  print(io, "\$")
end

function Base.show(io::IO, ::MIME"text/plain", r::Cyc)
  if !haskey(io,:typeinfo) print(io,typeof(r),": ") end
  show(io,r)
end

function normal_show(io::IO,p::Cyc)
  repl=get(io,:limit,false)
  TeX=get(io,:TeX,false)
  if valtype(p)<:Rational{<:Integer}
    den=denominator(p)
    p=Cyc{typeof(den)}(p*den)
  else
    den=1
  end
  res=""
  for (deg,v) in pairs(p)
    if deg==0 t=string(v)
    else
      t=format_coefficient(string(v))
      if repl || TeX
        r=(TeX ? "\\zeta" : "ζ")*stringind(io,conductor(p))*stringexp(io,deg)
      else
        r=(deg==1 ? "E($(conductor(p)))" : "E($(conductor(p)),$deg)")
      end
      t*=r
    end
    if t[1]!='-' t="+"*t end
    res*=t
  end
  if res[1]=='+' res=res[2:end] end
  if !isone(den)
    res=bracket_if_needed(res)
    res*=repl||TeX ? "/$den" : "//$den"
  end
  res
end

function Base.show(io::IO, p::Cyc)
  quadratic=get(io,:quadratic,true)
  repl=get(io,:limit,false)
  TeX=get(io,:TeX,false)
  @static if lazy lower!(p) end
  if conductor(p)==1
    n=num(p)
    if n isa Integer || n isa Rational{<:Integer}
      if denominator(n)==1
        if repl||TeX||haskey(io,:typeinfo) print(io,numerator(n))
        else print(io,"Cyc{",typeof(n),"}(",numerator(n),")")
        end
        return
      end
    else
      if repl||TeX||haskey(io,:typeinfo) print(io,n)
      else print(io,"Cyc{",typeof(n),"}(",n,")")
      end
      return
    end
  end
  rqq=[normal_show(io,p)]
  if quadratic && (valtype(p)<:Integer || valtype(p)<:Rational{<:Integer})
    q=Quadratic(p)
    if !isnothing(q) push!(rqq,repr(q;context=io)) end
    for test in [1-E(4),1+E(4),E(3),E(3,2),1-E(3),1-E(3,2),1+E(3),1+E(3,2),root(-3)]
      if !iszero(conductor(p)%conductor(test)) continue end
      q=Quadratic(p*1//test)
      if isnothing(q) continue end
      rq=repr(q;context=io)
      rq=format_coefficient(rq;allow_frac=true)
      t=format_coefficient(normal_show(io,test))
      if !isempty(rq) && rq[1]=='-' rq="-"*t*rq[2:end] else rq=t*rq end
      push!(rqq,rq)
    end
  end
  print(io,rqq[argmin(length.(rqq))])
end

# write a,b in common field Q(ζ_n)
function promote_conductor(a::Cyc{T},b::Cyc{T})where T
  if conductor(a)==conductor(b) return (a, b) end
  n=lcm(conductor(a),conductor(b))
  if n!=conductor(a)
    m=div(n,conductor(a))
    let m=m
      a=Cyc(n,valtype(a),i*m=>v for (i,v) in pairs(a))
    end
  end
  if n!=conductor(b)
    m=div(n,conductor(b))
    let m=m # the eltype of next iterator is Any!
      b=Cyc(n,valtype(b),i*m=>v for (i,v) in pairs(b))
    end
  end
  (a,b)
end

function Base.:+(x::Cyc,y::Cyc)
  a,b=promote(x,y)
  if obviouslyzero(a) return b end
  if obviouslyzero(b) return a end
  a,b=promote_conductor(a,b)
if impl==:svec res=Cyc_(dropzeros!(a.d+b.d))
else           res=Cyc(conductor(a),a.d+b.d)
end
  @static if !lazy lower!(res) end
  res
end

Base.:-(a::Cyc)=Cyc(conductor(a),-a.d)

function Base.:-(x::Cyc,y::Cyc)
  a,b=promote(x,y)
  if obviouslyzero(a) return -b end
  if obviouslyzero(b) return a end
  a,b=promote_conductor(a,b)
if impl==:svec res=Cyc_(dropzeros!(a.d-b.d))
else           res=Cyc(conductor(a),a.d-b.d)
end
  @static if !lazy lower!(res) end
  res
end

function Base.div(c::Cyc,a::Real)
if impl==:MM
  n=merge(div,c.d,a)
  Cyc(iszero(n) ? 1 : conductor(c),n)
else
  res=div.(c.d,a)
  iszero(res) ? Cyc{eltype(res)}(0) : Cyc(conductor(c),res)
end
end

function Base.:*(c::Cyc,a::Real)
  # the check for one saves a lot of time in some applications
  if isone(a) return Cyc{promote_type(valtype(c),typeof(a))}(c) end
  if iszero(a) return zero(Cyc{promote_type(valtype(c),typeof(a))}) end
  Cyc(conductor(c),c.d*a)
end
Base.:*(a::Real,c::Cyc)=c*a

function Base.:*(a::Cyc,b::Cyc;reduce=!lazy)
  a,b=promote(a,b)
  if obviouslyzero(a) return a end
  if obviouslyzero(b) return b end
  if conductor(a)==1 return b*num(a)
  elseif conductor(b)==1 return a*num(b)
  end
  n=lcm(conductor(a),conductor(b))
  na=div(n,conductor(a))
  nb=div(n,conductor(b))
  c=Cyc(n,valtype(a),na*i+nb*j=>ai*bj for (i,ai) in pairs(a),(j,bj) in pairs(b))
  reduce ? lower!(c) : c
end

Base.://(c::Cyc,a::Real)=Cyc(conductor(c),impl==:MM ? c.d//a : c.d.//a)
Base.://(a::Cyc,c::Cyc)=a*inv(Cyc{promote_type(valtype(c),Rational{Int})}(c))
Base.://(a::Real,c::Cyc)=a*inv(Cyc{promote_type(valtype(c),Rational{Int})}(c))
Base.://(c::Root1,a::Real)=Cyc(c)//a
Base.://(c::Root1,a::Cyc)=Cyc(c)//a
Base.://(a::Real,c::Root1)=a//Cyc(c)
Base.://(a::Cyc,c::Root1)=a//Cyc(c)
Base.:/(c::Cyc,a::Real)=c*inv(a)
Base.:/(a::Cyc,c::Cyc)=a*inv(c)
Base.div(c::Root1,a)=div(Cyc(c),a)

# change c to have data n,v
function Cyc!(c,n,v)
  if impl==:svec
    c.d=v
  elseif impl==:vec
    resize!(c.d,length(v))
    c.d.=v
  else
    c.n=n
    c.d=v
  end
  c
end

#---- defined here to avoid julia1.8 dependency ----------
if VERSION.major==1 && VERSION.minor<8
function allequal(a)
  if isempty(a) return true end
  o=first(a)
  all(==(o),a)
end
end

function lower!(c::Cyc) # write c in smallest Q(ζ_n) where it sits
  n=conductor(c)
 # println("lowering $(conductor(c)):$(c.d)")
  if n==1 return c end
if impl==:MM
  if obviouslyzero(c) return Cyc!(c,1,zero(c.d)) end
elseif impl==:vec
  if obviouslyzero(c) return Cyc!(c,1,zeros(valtype(c),1)) end
elseif impl==:svec
  if obviouslyzero(c) return Cyc!(c,1,spzeros(valtype(c),1)) end
end
  for (p,np) in EF(n)
    m=div(n,p)
if impl==:vec
    kk=filter(i->c.d[i]!=0,eachindex(c.d))
    val=c.d[kk]
    if np>1
      if all(k->(k-1)%p==0,kk)
        v=zeros(valtype(c),m)
        view(v,map(x->1+div(x-1,p),kk)).=val
        return lower!(Cyc!(c,m,v))
      end
    elseif iszero(length(kk)%(p-1))
      kk=kk.-1
      cnt=zeros(Int,m)
      for k in kk cnt[1+(k%m)]+=1 end
      if !all(x->iszero(x) || x==p-1,cnt) continue end
      u=(0:length(cnt)-1)[cnt.!=0]
      kk=@. div(u+m*mod(-u,p)*invmod(m,p),p)%m
      if !issorted(kk) sort!(kk) end
      v=zeros(valtype(c),m)
      if p==2
#       view(v,kk.+1).=[c.d[1+(k*p)%n] for k in kk]
        @views @. v[kk+1].=c.d[1+(kk*p)%n]
        return lower!(Cyc!(c,m,v))
      elseif all(k->allequal(map(i->c.d[1+(m*i+k*p)%n],1:p-1)),kk)
#       view(v,kk.+1).=[-c.d[1+(m+k*p)%n] for k in kk]
        @views @. v[kk+1].=-c.d[1+(m+kk*p)%n]
        return lower!(Cyc!(c,m,v))
      end
    end
elseif impl==:svec
    kk=c.d.nzind
    val=c.d.nzval
    if np>1
      if all(k->(k-1)%p==0,kk)
        return lower!(Cyc!(c,m,SparseVector(m,map(x->1+div(x-1,p),kk),val)))
      end
    elseif iszero(length(kk)%(p-1))
      cnt=zeros(Int,m)
      for k in kk cnt[1+(k-1)%m]+=1 end
      if !all(x->iszero(x) || x==p-1,cnt) continue end
      u=(0:length(cnt)-1)[cnt.!=0]
      kk=@. div(u+m*mod(-u,p)*invmod(m,p),p)%m
      if !issorted(kk) sort!(kk) end
      let c=c, kk=kk, p=p, m=m
      if p==2
        return lower!(Cyc!(c,m,SparseVector(m,kk.+1,[c.d[1+(k*p)%n] for k in kk])))
      elseif all(k->allequal(c.d[1+(m*i+k*p)%n] for i in 1:p-1),kk)
        return lower!(Cyc!(c,m,SparseVector(m,kk.+1,[-c.d[1+(m+k*p)%n] for k in kk])))
      end
      end
    end
elseif impl==:MM
    if np>1
      if all(k->first(k)%p==0,c.d)
        return lower!(Cyc!(c,m,MM(div(k,p)=>v for (k,v) in c.d;check=false)))
      end
    elseif iszero(length(c.d)%(p-1))
      cnt=zeros(Int,m)
      for (k,v) in c.d cnt[1+(k%m)]+=1 end
      if !all(x->iszero(x) || x==p-1,cnt) continue end
  #   u=findall(!iszero,cnt).-1 # next line is quite faster
  #   u=(0:length(cnt)-1)[cnt.!=0]
      u=[i-1 for (i,x) in pairs(cnt) if !iszero(x)]
      kk=@. div(u+m*mod(-u,p)*invmod(m,p),p)%m
      if p==2
        return lower!(Cyc!(c,m,MM(k=>c.d[(k*p)%n] for k in kk)))
      elseif all(k->allequal(c.d[(m*i+k*p)%n] for i in 1:p-1),kk)
        return lower!(Cyc!(c,m,MM(k=>-c.d[(m+k*p)%n] for k in kk)))
      end
    end
end
  end
  c
end

galois(c::Rational,n::Int)=c

galois(c::Integer,n::Int)=c

"""
  galois(c::Cyc,n::Int) applies to c the galois automorphism
  of Q(ζ_conductor(c)) raising all roots of unity to the n-th power.
  n should be prime to conductor(c).
# Examples
```julia-repl
julia> galois(1+E(4),-1) # galois(c,-1) is the same as conj(c)
Cyc{Int64}: 1-ζ₄

julia> galois(root(5),2)==-root(5)
true
```
"""
function galois(c::Cyc,n::Int)
  if gcd(n,conductor(c))!=1
    throw(DomainError(n,"should be prime to conductor $(conductor(c))"))
  end
  if obviouslyzero(c) return c end
  Cyc(conductor(c),valtype(c),e*n=>p for (e,p) in pairs(c))
end

function galois(c::Root1,n::Int)
  d=order(c)
  if gcd(n,d)!=1 error("$n should be prime to order($c)=$d") end
  Root1(;r=n*c.r)
end

Base.conj(c::Cyc)=galois(c,-1)

"""
`conjugates(c)`

returns the list of distinct galois conjugates of `c` (over the Rationals),
starting with c

```julia-repl
julia> conjugates(1+root(5))
2-element Vector{Cyc{Int64}}:
 1+√5
 1-√5
```
"""
function conjugates(c::Union{Cyc,Root1}) # Root1 or Cyc
  res=[c]
  for i in prime_residues(c isa Cyc ? conductor(c) : order(c))[2:end]
    c1=galois(c,i)
    if !(c1 in res) push!(res,c1) end
  end
  res
end

conjugates(c::Union{Rational{<:Integer},Integer})=[c]

function Base.inv(c::Cyc)
  if conductor(c)==1 return Cyc(inv(num(c))) end
  l=conjugates(c)
  r=l[2]
  for t in l[3:end] r=*(r,t;reduce=false) end
  r/num(*(c,r;reduce=true))
end

Base.:^(a::Cyc, n::Integer)=n>=0 ? Base.power_by_squaring(a,n) :
                                   Base.power_by_squaring(inv(a),-n)

Base.abs2(c::Cyc)=c*conj(c)
Base.abs(c::Cyc)=abs(complex(c))
Base.adjoint(c::Cyc)=conj(c)

"""
`Root1(c)`

`c`  should be a cyclotomic number (a  `Cyc`), or a `Real`. `Root1` returns
`E(n,e)` if `c==E(n,e)`, and `nothing` if `c` is not a root of unity.

```julia-repl
julia> r=Root1(-E(9,2)-E(9,5))
Root1: ζ₉⁸

julia> order(r)
9

julia> exponent(r)
8

julia> Cyc(r) # the Zumbroich normal form is a sum of 2 roots of unity
Cyc{Int64}: -ζ₉²-ζ₉⁵

julia> Root1(-E(9,4)-E(9,5)) # nothing
```
"""
function Root1(c::Cyc)
  @static if lazy lower!(c) end
  if !(all(x->last(x)==1,pairs(c)) || all(x->last(x)==-1,pairs(c)))
    return nothing
  end
  n=conductor(c)
  for i in prime_residues(n)
    if c==E(n,i) return Root1_(i//n) end
    if -c==E(n,i)
      if n%2==0 return E(n,div(n,2)+i)
      else return E(2n,n+2*i)
      end
    end
  end
  # return nothing
end

Base.:(==)(a::Root1,b::Number)=Cyc(a)==b # too expensive in lazy case

function Base.:*(a::Cyc,b::Root1)
  n=lcm(conductor(a),order(b))
  na=div(n,conductor(a))
  nb=div(n,order(b))
  res=Cyc(n,valtype(a),na*i+nb*exponent(b)=>va for (i,va) in pairs(a))
  @static if !lazy lower!(res) end
  res
end
Base.:*(b::Root1,a::Cyc)=a*b

Base.:+(a::Root1,b::Root1)=Cyc(a)+Cyc(b)
Base.:-(a::Root1,b::Root1)=Cyc(a)-Cyc(b)
Base.:-(r::Root1)=-Cyc(r)
Base.promote_rule(a::Type{Root1},b::Type{Cyc{T}}) where T =b
Base.promote_rule(a::Type{Root1},b::Type{<:Real})=Cyc{b}
Base.promote_rule(a::Type{Root1},b::Type{Complex{T}}) where T =Cyc{promote_type(T,Int)}

struct Quadratic
  a
  b
  root
  den
end

"""
  `Quadratic(c::Cyc)`

determines  if  `c`  lives  in  a  quadratic  extension  of  `ℚ `. The call
`q=Quadratic(c)`  returns a  struct `Quadratic`  with fields  `q.a`, `q.b`,
`q.root`,  `q.den` representing `c` as `(q.a + q.b root(q.root))//q.den` if
such a representation is possible or returns `q===nothing` otherwise.

# Examples
```julia-repl
julia> Quadratic(E(3,2)-2E(3))
(1-3√-3)/2

julia> Quadratic(1+E(5))

```
"""
function Quadratic(c::Cyc)
  den=denominator(c)
  c=numerator(c)
  if conductor(c)==1 return Quadratic(num(c),0,1,den) end
  f=Primes.factor(conductor(c))
  v2=get(f,2,0)
  if v2>3 || (v2==2 && any(p->p[1]!=2 && p[2]!=1,f)) ||
     (v2<2 && any(x->x!=1,values(f)))
    return nothing
  end
  f=keys(f)
  if v2==0
    sqr=conductor(c)
    if sqr%4==3 sqr=-sqr end
    gal=conjugates(c)
    if length(gal)!=2 return nothing end
    a=numerator(convert(valtype(c),sum(gal)))  # trace of 'c' over the rationals
    if iseven(length(f)) b=2*c[1]-a
    else b=2*c[1]+a
    end
    if iseven(a) && iseven(b) a>>=1; b>>=1; d=1
    else d=2
    end
  elseif v2==2
    sqr=conductor(c)>>2
    if sqr==1 a=c[0];b=-c[1]
    else
      a=c[4]
      if iseven(length(f)) a=-a end
      b=-c[sqr+4]
    end
    if sqr%4==1 sqr=-sqr; b=-b end
    d=1
  else		# v2 = 1 or 3
    sqr=conductor(c)>>2
    if sqr==2
      a=c[0];b=c[1]
      if b==c[3] sqr=-2 end
    else
      a=c[8]
      if iseven(length(f)) a=-a end
      b=c[(sqr>>1)+8]
      if b!=-c[3*(sqr>>1)-8] sqr=-sqr
      elseif (sqr>>1)%4==3 b=-b
      end
    end
    d=1
  end
  if d*c!=a+b*root(sqr) return nothing end
  return Quadratic(a,b,sqr,den*d)
end

function Base.show(io::IO,q::Quadratic)
  repl=get(io,:limit,false)
  TeX=get(io,:TeX,false)
  rq=string(q.a)
  if q.b!=0
    if iszero(q.a) rq=""
    elseif q.b>0 rq*="+" end
    rq*=q.b==1 ? "" : q.b==-1 ? "-" : string(q.b)
    r=string(q.root)
    rq*=TeX ? "\\sqrt{$r}" : repl ? "√$r" : "root($r)"
    if !iszero(q.a) && q.den!=1 rq="("*rq*")" end
  end
  print(io,rq)
  if q.den!=1 && rq!="0" print(io,(repl||TeX) ? "/" : "//",q.den) end
end

const inforoot=Ref(false)
xprint(x...;p...)=print(IOContext(stdout,:limit=>true,p...),x...)
function proot(x,n,r)
  xprint("root(",x;quadratic=false)
  if n!=2 xprint(",",n) end
  xprint(")=",r,"\n";quadratic=false)
end
const Irootdict=Dict{Tuple{Int,Int},Union{Cyc{Int},Int}}()

"""
`root(x,n=2)`

computes  the `n`-th root of `x` when we know  how to do it. We know how to
compute  `n`-th  roots  for  roots  of  unity, square roots of integers and
`n`-th  roots  of  perfect  `n`-th  powers  of  integers or square roots of
integers.

```julia-repl
julia> root(-1)
Cyc{Int64}: ζ₄

julia> root(E(4))
Root1: ζ₈

julia> root(27,6)
Cyc{Int64}: √3
```
"""
function root(x::Integer,n=2)
  if isone(n) || (!lazy && isone(x)) return x end
  get!(Irootdict,(n,x)) do
    if x==1 || (x==-1 && isodd(n)) return x end
    if x<0 && n==2 return E(4)*root(-x) end
    l=Primes.factor(x)
    if any(y->(2y)%n!=0,values(l))
      if x==-1 return root(E(2),n) end
      error("root($x,$n) not implemented")
    end
    a=prod(p^div(pow,n) for (p,pow) in l)
    b=[p for (p,pow) in l if pow%n!=0]
    for p in b
      if p%4==1 a*=sum(k->E(p,k^2),1:p)
      elseif p==2 a*=E(8)-E(8,3)
      elseif p%4==3 a*=E(4,3)*sum(k->E(p,k^2),1:p)
      end
    end
    if inforoot[] proot(x,n,a) end
    a
  end
end

root(x::Rational,n=2)=root(numerator(x),n)//root(denominator(x),n)

# find the "canonical" best of the n possible roots
function root(r::Root1,n=2)
  if iszero(n) return one(r) end
  d=order(r)
  j=1
  n1=n
  while true
    k=gcd(n1,d)
    n1=div(n1,k)
    j*=k
    if k==1 break end
  end
  E(j*d,exponent(r)*gcdx(n1,d)[2])
end

const Crootdict=Dict{Tuple{Int,Cyc},Cyc}()
function root(x::Cyc,n=2)
  if isone(n) || isone(x) return x end
  get!(Crootdict,(n,x)) do
  d=denominator(x)
  if d!=1 return root(numerator(x),n)//root(d,n) end
  d=gcd(coefficients(x))
  if d!=1 return root(div(x,d),n)*root(d,n) end
    r=Root1(x)
    if isnothing(r)
      if conductor(x)>1 error("cannot compute root($x,$n)") end
      return root(num(x),n)
    end
    res=Cyc(root(r,n))
    if inforoot[] proot(x,n,res) end
    res
  end
end

Base.gcd(v::Vector{<:Cyc})=one(v[1])
Base.gcd(a::Cyc,b::Cyc)=one(a)
Base.gcd(a::Cyc,b::Number)=one(a)
Base.gcd(b::Number,a::Cyc)=one(a)

# testmat(12)^2
# 347.534 ms (4367402 allocations: 366.17 MiB) in 1.5.3
# 565.431 ms (5861810 allocations: 775.28 MiB) in 1.5.3, HModuleElts
function testmat(p)
  ss=[[i,j] for i in 0:p-1 for j in i+1:p-1]
  [(E(p,i'*reverse(j))-E(p,i'*j))//p for i in ss,j in ss]
end
end


text:
struct ImmutableX
    x
end

mutable struct MutableX
    x
end

counter = 0

struct CounterX
    x
    id
    function CounterX(x)
        global counter
        counter += 1
        new(x, counter)
    end
end

a = ImmutableX(1)
b = ImmutableX(1)
a_mutable = MutableX(1)
b_mutable = MutableX(1)
a_array = ImmutableX([1])
b_array = ImmutableX([1])
a_counter = CounterX(1)
b_counter = CounterX(1)

println(a === b) # True
println(a_mutable === b_mutable) #False
println(a_array === b_array) # False
println(a_counter === b_counter) # False
println(a_counter == b_counter) # False

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
# Practical Machine learning
# Artificial Neural Network
# Chapter 11
function sigmoidGradient(z)
  g = SIGMOIDGRADIENT(z)
  return g
end


text:
# Script to demonstrate convergence analysis of Ecut values for DFTK.jl

using Atomistic
using AtomsBase
using DFTK
using NBodySimulator
using Unitful
using UnitfulAtomic

setup_threading(n_blas = 4)

N = 8
element = :Ar
box_size = 1.5u"nm" # this number was chosen arbitrarily
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = DynamicSystem(initial_bodies, box_size)
eq_simulator = NBSimulator(
    Δt = Δt,
    steps = 20000,
    t₀ = 0.0u"s",
    thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
)
potential = LJPotential(
    ϵ = 1.657e-21u"J",
    σ = 0.34u"nm",
    R = 0.765u"nm"
)
eq_result = @time simulate(initial_system, eq_simulator, potential)

display(@time plot_rdf(eq_result, potential.σ, 0.5))

dftk_potential = DFTKPotential(
    psp = ElementPsp(:Ar, psp = load_psp(list_psp(:Ar, functional = "lda")[1].identifier)),
    lattice = box_size * [[1.0 0 0]; [0 1.0 0]; [0 0 1.0]],
    Ecut = 5u"hartree",
    kgrid = [1, 1, 1],
    damping = 0.7,
    mixing = LdosMixing()
)

display(@time analyze_convergence(get_system(eq_result), dftk_potential, (5:2.5:25) .* u"hartree"))

;

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
#TODO: download https://cdn.plot.ly/plotly-1.54.2.min.js to deps/ if it doesn't exist
file_path = ""
if get(ENV, "PLOTS_HOST_DEPENDENCY_LOCAL", "false") == "true"
    global file_path
    local_fn = joinpath(dirname(@__FILE__), "plotly-1.54.2.min.js")
    if !isfile(local_fn)
        @info("Cannot find deps/plotly-1.54.2.min.js... downloading latest version.")
        download("https://cdn.plot.ly/plotly-1.54.2.min.js", local_fn)
        isfile(local_fn) && (file_path = local_fn)
    else
        file_path = local_fn
    end
end

open("deps.jl", "w") do io
    println(io, "const plotly_local_file_path = $(repr(file_path))")
end


text:
module Commands

using ..PanGraph: panic

export Arg, Command
export flags, usage, arg

Maybe{T} = Union{Missing,T}

mutable struct Arg
    type::Type
    meta::AbstractString
    flag::NamedTuple{(:short,:long), Tuple{Maybe{AbstractString},Maybe{AbstractString}}}
    help::AbstractString

    value
end

flags(arg::Arg) =
    if !ismissing(arg.flag.long) && !ismissing(arg.flag.short)
        "$(arg.flag.long), $(arg.flag.short)"
    elseif !ismissing(arg.flag.long) && ismissing(arg.flag.short)
        "$(arg.flag.long)"
    elseif ismissing(arg.flag.long) && !ismissing(arg.flag.short)
        "$(arg.flag.short)"
    else
        ""
    end

mutable struct Command
    cmd::AbstractString
    use::AbstractString
    meta::AbstractString
    foot::AbstractString

    sub::Array{Command,1}
    arg::Array{Arg,1}

    run::Maybe{Function}
end

Command(cmd,help,meta,foot,arg::Array{Arg})     = Command(cmd,help,meta,foot,Command[],arg,missing)
Command(cmd,help,meta,foot,sub::Array{Command}) = Command(cmd,help,meta,foot,sub,Arg[],missing)

Command(cmd,help,meta,foot,arg::Array{Arg},fn::Function)     = Command(cmd,help,meta,foot,Command[],arg,fn)
Command(cmd,help,meta,foot,sub::Array{Command},fn::Function) = Command(cmd,help,meta,foot,sub,Arg[],fn)

function usage(io::IO, cmd::Command)
    print(io,"""
    usage: $(cmd.use)

        $(cmd.meta)

    """)

    if length(cmd.sub) > 0
        println(io,"The available commands are:\n")

        maxlen = maximum(length(sub.cmd) for sub in cmd.sub)
        for sub in cmd.sub
            println(io, "\t$(sub.cmd)$(" "^(maxlen-length(sub.cmd)))\t$(sub.meta)")
        end

		println(io,"")
    end

    if length(cmd.arg) > 0
        println(io,"The available options are:\n")

        maxlen = maximum(length(flags(arg)) for arg in cmd.arg)
        for (i,arg) in enumerate(cmd.arg)
            flag = flags(arg)
            println(io, "\t$(arg.help)")
            println(io, "\t$(flag)$(" "^(maxlen-length(flag)))\t$(arg.meta)")
            println(io, "\t$(" "^maxlen)\t$(arg.type) default: $(arg.value)")
            if i < length(cmd.arg)
                println("")
            end
        end

		println(io,"")
    end

    if length(cmd.foot) > 0
        println(io,"The expected arguments are:\n")
        println(io,"\t", join(split(cmd.foot,"\n"), "\n\t"))
    end

	return 2
end

usage(cmd::Command) = usage(stderr, cmd)

function Base.parse(cmd::Command, args)
    for sub in cmd.sub
        if sub.cmd == args[1]
            cmd.run = sub.run
            return args[2:end]
        end
    end

    if length(cmd.sub) > 0
        if length(cmd.sub) ≥ 1 && args[1] == "-h"
            usage(cmd)
            exit(2)
        end
        panic("""
        pangraph $(args[1]): unknown command
        Run 'pangraph help' for usage.
        """)
    end

	length(cmd.arg) > 0 || return length(args) > 0 ? args : nothing

    let arg = ""
        itr = Iterators.Stateful(args)
        if isempty(itr)
            usage(cmd)
            return nothing
        end

        arguments = String[]
        @label ARGLOOP #-------------------------
        !isempty(itr) || @goto ENDLOOP

        word = peek(itr)
        if !startswith(word, "-")
            push!(arguments,popfirst!(itr))
            @goto ARGLOOP
        end

        popfirst!(itr) # advance iterator
        for arg ∈ cmd.arg
            if arg.flag.long == word || arg.flag.short == word
                arg.value =
                if arg.type == Bool
                    true
                else
                    if isempty(itr)
                        panic("""flag '$(word)' requires non-missing input data of type '$(arg.type)'
                        Run 'pangraph help' for usage.
                        """)
                    end

                    val = popfirst!(itr)
                    try
                        arg.type == String ? val : Base.parse(arg.type, val)
                    catch e
                        if isa(e, ArgumentError)
                            panic("""flag '$(word)' requires input data of type '$(arg.type)'. recieved '$(val)'
                            Run 'pangraph help' for usage.
                            """)
                        else
                            throw(e)
                        end
                    end
                end

                @goto ARGLOOP
            end
        end

        if word == "-h"
            usage(cmd)
            exit(2)
        end

        panic("""provided flag '$(word)' not expected
        Run 'pangraph help' for usage.
        """)

        @label ENDLOOP #-------------------------

        return arguments
    end
end

function Base.run(cmd::Command, args)
    if ismissing(cmd.run)
        error("command $(cmd.cmd) does not have a valid execution function")
    end
    cmd.run(args)
end

function arg(cmd::Command, short::AbstractString)
    i = 1
    while i ≤ length(cmd.arg) && cmd.arg[i].flag.short != short
        i += 1
    end
    cmd.arg[i].flag.short != short && error("flag $(short) not found")

    return cmd.arg[i].value
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using Base.Test

julia = Base.julia_cmd().exec[1]
build_app_jl = joinpath(@__DIR__, "..", "build_app.jl")
examples_hello = joinpath(@__DIR__, "..", "examples", "hello.jl")

builddir = mktempdir()
@assert isdir(builddir)

@testset "Basic file resource args" begin
# Test the build_app.jl script CLI args.
res1 = @__FILE__ # haha copy this file itself as a "resource"!
res2 = joinpath(@__DIR__, "runtests.jl") # lol sure this is a resource, why not.
NEWARGS = Base.shell_split("""--verbose
        -R $res1 --resource $res2 -L $res1 --lib $res2
        $examples_hello "HelloWorld" $builddir""")
eval(:(ARGS = $NEWARGS))
@test 0 == include("$build_app_jl")
@test isdir("$builddir/HelloWorld.app")
@test isfile("$builddir/HelloWorld.app/Contents/MacOS/hello")

# Make sure the specified resources and libs were copied:
@test isfile("$builddir/HelloWorld.app/Contents/Resources/$(basename(res1))")
@test isfile("$builddir/HelloWorld.app/Contents/Resources/$(basename(res2))")
@test isfile("$builddir/HelloWorld.app/Contents/Libraries/$(basename(res1))")
end

@testset "Exits without juliaprog_main" begin
  @test !success(`$julia $build_app_jl --verbose`)  # no main.jl
end


text:
using Documenter

using MaxEntropyGraphs
#push!(LOAD_PATH, "../src/")

ci = get(ENV, "CI", "") == "true"

makedocs(sitename="MaxEntropyGraphs.jl",
         authors="Bart De Clerck",
         format = Documenter.HTML(prettyurls = ci),
         modules=[MaxEntropyGraphs],
         pages = [
            "Home" => "index.md",
            "Models" => "models.md",
            "Higher order metrics" => "derivedquantities.md",
            "GPU acceleration" => "GPU.md"
         ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
if ci
    @info "Deploying documentation to GitHub"
    deploydocs(
        repo = "github.com/B4rtDC/MaxEntropyGraphs.jl.git",
        devbranch = "main"
        )
end


text:
# Constructors:
###############
"""
    GBScalar{T}()

Create an unassigned GBScalar of type T.
"""
GBScalar{T}() where {T} = GBScalar{T}(libgb.GxB_Scalar_new(toGBType(T)))

"""
    GBScalar(v::T)

Create a GBScalar of type `T` and assign `v` to it.
"""
function GBScalar(v::T) where {T <: valid_union}
    x = GBScalar{T}()
    x[] = v
    return x
end

# Some Base and basic SparseArrays/LinearAlgebra functions:
###########################################################
Base.unsafe_convert(::Type{libgb.GxB_Scalar}, s::GBScalar) = s.p

function Base.copy(s::GBScalar{T}) where {T}
    return GBScalar{T}(libgb.GxB_Scalar_dup(s))
end

clear!(s::GBScalar) = libgb.GxB_Scalar_clear(s)

# Type dependent functions setindex and getindex:
for T ∈ valid_vec
    func = Symbol(:GxB_Scalar_setElement_, suffix(T))
    @eval begin
        function Base.setindex!(value::GBScalar{$T}, s::$T)
            libgb.$func(value, s)
        end
    end
    func = Symbol(:GxB_Scalar_extractElement_, suffix(T))
    @eval begin
        function Base.getindex(value::GBScalar{$T})
            libgb.$func(value)
        end
    end
end

Base.eltype(::Type{GBScalar{T}}) where{T} = T

function Base.show(io::IO, ::MIME"text/plain", s::GBScalar)
    gxbprint(io, s)
end

SparseArrays.nnz(v::GBScalar) = Int64(libgb.GrB_Scalar_nvals(v))


text:
using SafeTestsets

@safetestset "Constants Test" begin include("testsets/test_constants.jl") end
@safetestset "Dimensions Test" begin include("testsets/test_dimensions.jl") end
@safetestset "Body Test" begin include("testsets/test_bodies.jl") end
@safetestset "Two Body" begin include("testsets/test_twobody.jl") end


text:
using IterativeSolvers
using Test
using LinearMaps
using LinearAlgebra
using Random
using SparseArrays

#GMRES
@testset "GMRES" begin

srand(1234321)
n = 10

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n)
    b = rand(T, n)
    F = lu(A)
    tol = √eps(real(T))

    # Test optimality condition: residual should be non-increasing
    x, history = gmres(A, b, log = true, restart = 3, maxiter = 10, tol = tol);
    @test isa(history, ConvergenceHistory)
    @test all(diff(history[:resnorm]) .<= 0.0)

    # Left exact preconditioner
    x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, tol=tol, log=true)
    @test history.isconverged
    @test norm(F \ (A * x - b)) / norm(b) ≤ tol

    # Right exact preconditioner
    x, history = gmres(A, b, Pl=Identity(), Pr=F, maxiter=1, restart=1, tol=tol, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ tol
end

@testset "SparseMatrixCSC{$T}" for T in (Float64, ComplexF64)
    A = sprand(T, n, n, 0.5) + I
    b = rand(T, n)
    F = lu(A)
    tol = √eps(real(T))

    # Test optimality condition: residual should be non-increasing
    x, history = gmres(A, b, log = true, restart = 3, maxiter = 10);
    @test all(diff(history[:resnorm]) .<= 0.0)

    # Left exact preconditioner
    x, history = gmres(A, b, Pl=F, maxiter=1, restart=1, log=true)
    @test history.isconverged
    @test norm(F \ (A * x - b)) / norm(b) ≤ tol

    # Right exact preconditioner
    x, history = gmres(A, b, Pl = Identity(), Pr=F, maxiter=1, restart=1, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ tol
end

@testset "Linear operator defined as a function" begin
    A = LinearMap(cumsum!, 100; ismutating=true)
    b = rand(100)
    tol = 1e-5

    x = gmres(A, b; tol=tol, maxiter=2000)
    @test norm(A * x - b) / norm(b) ≤ tol
end
end


text:
function add_solver_xml!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
    # Add solution display
    if !anyints(m)
	m.xmlmodel = replace(m.xmlmodel, r"(?s)<wantint><![CDATA[yes]]></wantint>" => "<wantint><![CDATA[no]]></wantint>")
	# Lets hack this here. Maybe if we get fancy and offer other categories we should break this out a bit more
	m.xmlmodel = replace(m.xmlmodel, r"(?s)<category>milp</category>" => "<category>lp</category>")
    end
    # Add user options
    param_string = ""
    for key in keys(m.solver.params)
	param_string *= "$(key) $(m.solver.params[key])\n"
    end
    # Add user options
    m.xmlmodel = replace(m.xmlmodel, r"(?s)<param>.*</param>" => "<param><![CDATA[$(param_string)]]></param>")
end

function solution_status(m, args...)
    for s in args
	if match(Regex("SOLUTION STATUS\\W+?:\\W+?$s"), m.last_results) != nothing
	    return true
	end
    end
    return false
end

function parse_status!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
    if solution_status(m, "INTEGER_OPTIMAL", "OPTIMAL")
	m.status = OPTIMAL
    elseif solution_status(m, "PRIMAL_UNBOUNDED", "UNBOUNDED")
	m.status = UNBOUNDED
    elseif solution_status(m, "DUAL_INFEASIBLE")
	m.status = UNBNDORINF
    elseif solution_status(m, "DUAL_UNBOUNDED", "PRIMAL_INFEASIBLE", "INFEASIBLE")
	m.status = INFEASIBLE
    end
end

function parse_objective!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
    sci = match(r"PRIMAL\W+?OBJECTIVE\W+?:\W+?(-?[\d\.]+e[\+\-]\d+)", m.last_results).captures
    m.objVal = parse(Float64, sci[1])
end

function parse_solution!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
    for v in collect((m.match for m = eachmatch(r"V(\d+).+?(-?[\d\.]+e[\+\-]\d+)", m.last_results)))
	regmatch = match(r"V(\d+).+?(-?[\d\.]+e[\+\-]\d+)", v)
	m.solution[parse(Int, regmatch.captures[1])] = parse(Float64, regmatch.captures[2])
    end
end

function parse_duals!(::NEOSSolver{:MOSEK, :MPS}, m::MPSModel)
    m.duals = zeros(m.nrow)
    for v in collect((m.match for m = eachmatch(r"\d+\W+C(\d+).+?\n", m.last_results)))
	s = split(v)
	d_lb = parse(Float64, s[7])
	d_ub = -parse(Float64, s[8])
	m.duals[parse(Int, s[2][2:end])] = abs(d_lb) > abs(d_ub) ? d_lb : d_ub
    end

    m.reducedcosts = zeros(m.ncol)
    for v in collect((m.match for m = eachmatch(r"\d+\W+V(\d+).+?\n", m.last_results)))
	s = split(v)
	d_lb = parse(Float64, s[7])
	d_ub = -parse(Float64, s[8])
	m.reducedcosts[parse(Int, s[2][2:end])] = abs(d_lb) > abs(d_ub) ? d_lb : d_ub
    end

end


text:
const em = emulate

@test em(Int) != nothing
@test em(Integer) != nothing
@test em(Array) != nothing
@test em(Array{T,N} where T where N) != nothing
@test em(Tuple{Int8,Array{X,Y} where {X, Int8<:Y<:Int32}}) != nothing

@test em(Base.Bottom) == NewBottom
@test em(Any) == NewAny
# @test em(Tuple{Float16, Float64}) == NewDataType(NewTypeVar(:Tuple, BaseModule))


text:
module Soli

using ConstrainedRootSolvers: ReduceStepMethodND, SolutionToleranceND, find_peak
using DocStringExtensions: METHODLIST
using ..Emerald: VanGenuchten
using UnPack: @unpack

import ..Emerald: BrooksCorey

# export public functions
export soil_ψ_25

#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Sep-30: move this function out of BrooksCorey struct as an external method for the constructor
#
#######################################################################################################################################################################################################
"""
    BrooksCorey{FT}(name::String)

A constructor for BrooksCorey to create BrooksCorey type soil from VanGenuchten type
"""
BrooksCorey{FT}(vg::VanGenuchten{FT}) where {FT<:AbstractFloat} = (
    # generate data to fit
    _Θs   = range(vg.Θ_RES+FT(1e-2); stop=vg.Θ_SAT-FT(1e-2), length=30);
    _Ψ_vG = -1 .* soil_ψ_25.([vg], _Θs);

    # function to fit BrooksCorey parameters
    @inline _fit(x) = (
        _bc   = BrooksCorey{FT}(vg.TYPE, x[1], x[2], vg.Θ_SAT, vg.Θ_RES);
        _Ψ_BC = -1 .* soil_ψ_25.([_bc], _Θs);
        _diff = sum( (log.(_Ψ_BC) .- log.(_Ψ_vG)) .^ 2 );
        return -_diff
    );

    _st  = SolutionToleranceND{FT}([1e-3, 1e-6], 30);
    _ms  = ReduceStepMethodND{FT}(x_mins = FT[1e-3, 1e-6], x_maxs = FT[ 100, 1000], x_inis = [(2*vg.N-1) / (vg.N-1), 1 / (vg.Α)], Δ_inis = FT[0.1, 1e-3]);
    _sol = find_peak(_fit, _ms, _st);

    return BrooksCorey{FT}(vg.TYPE, _sol[1], _sol[2], vg.Θ_SAT, vg.Θ_RES);
)

#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Sep-30: create this function to work with two soil types using either VanGenuchten or BrooksCorey function
#
#######################################################################################################################################################################################################
"""
This function calculates soil metric potential from soil retention curve type and soil volumetric water potential. The supported methods are

$(METHODLIST)
"""
function soil_ψ_25 end

"""
    soil_ψ_25(vc::Union{BrooksCorey{FT}, VanGenuchten{FT}}, θ::FT) where {FT<:AbstractFloat}

Return the soil metric potential, given
- `vc` [`BrooksCorey`](@ref) or [`VanGenuchten`](@ref) type structure
- `θ` Soil volumetric water content (absolute value)

# Examples
```
soil_ψ_25(bc, FT(0.2))
soil_ψ_25(vg, FT(0.2))
```
"""
soil_ψ_25(vc::BrooksCorey{FT}, θ::FT) where {FT<:AbstractFloat} = (
    @unpack B, Ψ_SAT, Θ_RES, Θ_SAT = vc;

    # calculate effective θ
    _θ_e = min(1, max(0, (θ - Θ_RES) / (Θ_SAT - Θ_RES)));

    # if _θ_e >= 1, return 0
    if _θ_e >= 1
        return FT(0)
    end;

    return -Ψ_SAT / (_θ_e ^ B)
);

soil_ψ_25(vc::VanGenuchten{FT}, θ::FT) where {FT<:AbstractFloat} = (
    @unpack M, N, Α, Θ_RES, Θ_SAT = vc;

    # calculate effective θ
    _θ_e = min(1, max(0, (θ - Θ_RES) / (Θ_SAT - Θ_RES)));

    # if _θ_e >= 1, return 0
    if _θ_e >= 1
        return FT(0)
    end;

    return -1 * (_θ_e ^ (-1/M) - 1) ^ (1/N) / Α
);

end # module


text:
function simple(itr::Function, mgr::IterationManager=DefaultManager(1e-6, 10000))
    NS = 3000
    c = zeros(NS,4)

    # construct manager and state
    istate = DefaultState(c)  # initial state is the coefficients

    # Stages 2, 3, 4
    itr(mgr, istate) do c
        c_new = similar(c)

        # unpack coefficient vector from last time
        cadj, cnonadj, c1, c2 = c[:,1], c[:,2], c[:,3], c[:,4]

        c_new[:, 1] = rand()
        c_new[:, 2] = rand()
        c_new[:, 3] = rand()
        c_new[:, 4] = rand()

        # return new coefficients
        c_new
    end
end

function by_hand(maxiter=10000, tol=1e-6, verbose=true)
    NS = 3000
    c = zeros(NS,4)
    tol = 1e-6
    dist = 1.
    iter = 0
    elapsed = 0.0
    old_time = time()

    print_skip = div(maxiter, 5)

    while dist > tol && iter < maxiter
        # Stage 2: Iteration
        c_new = similar(c)

        # unpack coefficient vector from last time
        cadj, cnonadj, c1, c2 = c[:,1], c[:,2], c[:,3], c[:,4]

        c_new[:, 1] = rand()
        c_new[:, 2] = rand()
        c_new[:, 3] = rand()
        c_new[:, 4] = rand()

        # Stage 3: Between iteration processing
        dist = maxabs(c - c_new)
        iter += 1
        new_time = time()
        elapsed += new_time - old_time

        if verbose && iter % print_skip == 0
            println("Iteration: $iter\t dist: $(round(dist, 4))\t elapsed: $(elapsed)")
        end

        c = copy(c_new)
        old_time = new_time
    end
end


text:
function compute_new_ud_a!(grid, q_new, q, q_tendencies, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            a_predict = q[:a, k, i] + params[:Δt][1] * q_tendencies[:a, k, i]
            q_new[:a, k, i] = bound(a_predict, params[:a_bounds])
        end
    end
end

function compute_new_ud_w!(grid, q_new, q, q_tendencies, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    # Solve for updraft velocity
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            a_new_k = q_new[:a, k, i]
            ρ_k = aux[:ρ_0, k]
            w_i = q[:w, k, i]
            a_k = q[:a, k, i]
            ρa_k = ρ_k * a_k
            ρa_new_k = ρ_k * a_new_k
            ρaw_k = ρa_k * w_i
            w_predict =
                ρaw_k / ρa_new_k +
                params[:Δt][1] / ρa_new_k * q_tendencies[:w, k, i]
            q_new[:w, k, i] = bound(w_predict, params[:w_bounds])
        end
    end
end

function compute_new_ud_scalars!(grid, q_new, q, q_tendencies, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            a_k = q[:a, k, i]
            a_k_new = q_new[:a, k, i]
            ρ_k = aux[:ρ_0, k]
            ρa_k = ρ_k * a_k
            ρa_new_k = ρ_k * a_k_new
            θ_liq_predict =
                (
                    ρa_k * q[:θ_liq, k, i] +
                    params[:Δt][1] * q_tendencies[:θ_liq, k, i]
                ) / ρa_new_k
            q_tot_predict =
                (
                    ρa_k * q[:q_tot, k, i] +
                    params[:Δt][1] * q_tendencies[:q_tot, k, i]
                ) / ρa_new_k
            q_new[:θ_liq, k, i] = θ_liq_predict
            q_new[:q_tot, k, i] = q_tot_predict
        end
    end
end

#####
##### Update functions
#####

function assign_new_to_values!(grid, q_new, q, aux)
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems(grid), i in ud, name in (:w, :q_tot, :θ_liq)
        q_new[name, k, i] = q[name, k, i]
    end
end

function assign_values_to_new!(grid, q, q_new, aux)
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems(grid), i in ud
        @inbounds for name in (:a, :w, :q_tot, :θ_liq)
            q[name, k, i] = q_new[name, k, i]
        end
    end
    @inbounds for k in over_elems(grid)
        q[:tke, k, en] = q_new[:tke, k, en]
        q[:a, k, en] = 1 - sum([q_new[:a, k, i] for i in ud])
        q[:θ_liq, k, gm] = q_new[:θ_liq, k, gm]
        q[:q_tot, k, gm] = q_new[:q_tot, k, gm]
    end
end


text:
#kernels_pure.jl

function hpl2(A, b)
    F = LinAlg.generic_lufact!(A) #Pure Julia implementation, NO LAPACK OR BLAS calls
    F \ b
end