text:
module TestFundamentals
using AlphaVantage
using Test

TEST_SLEEP_TIME =  parse(Float64, get(ENV, "TEST_SLEEP_TIME", "15"))
MAX_TESTS = parse(Int64, get(ENV, "MAX_TESTS", "1"))

@testset "Fundamentals" begin
    symbol = "IBM"

    @testset "Overview" begin
        data = company_overview(symbol, datatype="json")
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 59
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Income Statement" begin
        data = income_statement(symbol, datatype="json")
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Balance Sheet" begin
        data = balance_sheet(symbol, datatype="json")
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Cash Flow" begin
        data = cash_flow(symbol, datatype="json")
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Earnings" begin
        data = earnings(symbol, datatype="json")
        @test typeof(data) === Dict{String, Any}
        @test length(data) === 3
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "Listing Status" begin
        @testset "Default" begin
            data = listing_status()
            @test length(data) == 2
            sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
        end

        @testset "Delisted and Date" begin
            data = listing_status(state = "delisted", date = "2020-12-17")
            @test length(data) == 2
            sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
        end
    end

    @testset "Earnings Calendar" begin
        data = earnings_calendar(3)
        @test length(data) == 2
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end

    @testset "IPO Calendar" begin
        data = ipo_calendar()
        @test length(data) == 2
        sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
    end
end

end # module


text:
add_verbose_scope(:QuadraticExt)
add_assert_scope(:QuadraticExt)

##############################################################################
#
#  Sieves for primes and squarefree numbers
#
##############################################################################

function squarefree_up_to(n::Int; coprime_to::Array{fmpz,1}=fmpz[])

  list= trues(n)
  for x in coprime_to
    t=x
    while t<= n
      @inbounds list[Int(t)]=false
      t+=x
    end
  end
  i=2
  b=Base.sqrt(n)
  while i<=b
    @inbounds if list[i]
      j=i^2
      @inbounds if !list[j]
        i+=1
        continue
      else
        @inbounds list[j]=false
        t=2*j
        while t<= n
          @inbounds list[t]=false
          t+=j
        end
      end
    end
    i+=1
  end
  @inbounds return Int[i for i=1:n if list[i]]

end

function SqfSet(n::fmpz; coprime_to::Int=-1)
  return (i for i=-n:n if i!= -1 && i!=0 && issquarefree(i))
end

function primes_up_to(n::Int)

  list= trues(n)
  list[1]=false
  i=2
  s=4
  while s<=n
    list[s]=false
    s+=2
  end
  i=3
  b=sqrt(n)
  while i<=b
    if list[i]
      j=3*i
      s=2*i
      while j<= n
        list[j]=false
        j+=s
      end
    end
    i+=1
  end
  return Int[i for i=1:n if list[i]]

end

###########################################################################
#
#  Some useful functions
#
###########################################################################

#This function gets a quotient of the ray class group and the action on
# the ray class group
# In output, we get the quotient group as a ZpnGModule

function _action_on_quo(mq::GrpAbFinGenMap, act::Array{GrpAbFinGenMap,1})

  q=mq.header.codomain
  S,mS=snf(q)
  n=Int(S.snf[end])
  R=ResidueField(FlintZZ, n, cached=false)
  W=MatrixSpace(R, ngens(S), ngens(S), false)
  quo_action=Array{nmod_mat,1}(length(act))
  for s=1:length(act)
    quo_action[s]=W(mS.map*act[i].map*mS.imap)
  end
  return ZpnGModule(S, quo_action)

end

function _min_wild(D::Dict{NfOrdIdl, Int})

  res=1
  primes_done=fmpz[]
  for (p,v) in D
    s=minimum(p)
    if s in primes_done
      continue
    end
    push!(primes_done, s)
    res*=Int(s)^v
  end
  return res

end

function totally_positive_generators(mr::MapRayClassGrp, a::Int, wild::Bool=false)

  if isdefined(mr, :tame_mult_grp)
    tmg=mr.tame_mult_grp
    for (p,v) in tmg
      mr.tame_mult_grp[p]=(make_positive(v[1],a),v[2],v[3])
      @hassert :RayFacElem 1 istotally_positive(mr.tame_mult_grp[p][1])
    end
  end
  if wild && isdefined(mr, :wild_mult_grp)
    wld=mr.wild_mult_grp
    for (p,v) in wld
      x=v[1]
      for i=1:length(x)
        x[i]=make_positive(x[i],a)
        @hassert :RayFacElem 1 iscoprime(ideal(parent(x[i]), x[i]), ideal(parent(x[i]), a))
      end
      mr.wild_mult_grp[p]=(x,v[2],v[3])
    end
  end

end

function make_positive(x::NfOrdElem, a::Int)

  els=conjugates_arb(x)
  m=1
  for i=1:length(els)
    if isreal(els[i])
      y=BigFloat(midpoint(real(els[i]/a)))
      if y>0
        continue
      else
        m=max(m,1-ceil(Int,y))
      end
    end
  end
  @hassert :RayFacElem 1 iscoprime(ideal(parent(x),x+m*a), ideal(parent(x), a))
  @hassert :RayFacElem 1 istotally_positive(x+m*a)
  return x+m*a

end

function small_generating_set(Aut::Array{NfToNfMor, 1})
  K=Aut[1].header.domain
  a=gen(K)
  Identity = Aut[1]
  for i in 1:length(Aut)
    Au = Aut[i]
    if Au.prim_img == a
      Identity = Aut[i]
      break
    end
  end
  return  Hecke.small_generating_set(Aut, *, Identity)
end

function _are_there_subs(G::GrpAbFinGen,gtype::Array{Int,1})

  H=DiagonalGroup(gtype)
  H, _=snf(H)
  G1, _=snf(G)
  if length(G1.snf)<length(H.snf)
    return false
  end
  for i=0:length(H.snf)-1
    if !divisible(G1.snf[end-i],H.snf[end-i])
      return false
    end
  end
  return true
end

function issquarefree(n::Union{Int,fmpz})
  if iszero(n)
    throw(error("Argument must be non-zero"))
  end
  return isone(abs(n)) || maximum(values(factor(n).fac)) == 1
end

function absolute_automorphism_group(C::ClassField)
  L = number_field(C)
  K = base_field(C)
  autK = automorphisms(K)
  id = find_identity(autK, *)
  autK_gen = small_generating_set(autK, *, id)
  return absolute_automorphism_group(C, autK_gen)
end

function absolute_automorphism_group(C::ClassField, aut_gen_of_base_field::Array{NfToNfMor, 1})
  L = number_field(C)
  aut_L_rel = rel_auto(C)::Vector{NfRel_nsToNfRel_nsMor}
  rel_extend = NfRel_nsToNfRel_nsMor[ Hecke.extend_aut(C, a) for a in aut_gen_of_base_field ]
  rel_gens = vcat(aut_L_rel, rel_extend)
  return rel_gens::Array{NfRel_nsToNfRel_nsMor, 1}
end

function automorphism_groupQQ(C::ClassField)
  return rel_auto(C)
end

##########################################################################################################
#
#  Functions for conductors
#
##########################################################################################################

function tame_conductors_degree_2(O::NfOrd, bound::fmpz)

  K=nf(O)
  d=degree(O)
  b1=Int(root(bound,d))
  ram_primes=collect(keys(factor(O.disc).fac))
  sort!(ram_primes)
  filter!(x -> x!=2 ,ram_primes)
  list=squarefree_up_to(b1, coprime_to=vcat(ram_primes,2))

  extra_list=Tuple{Int, fmpz}[(1,fmpz(1))]
  for q in ram_primes
    tr=prime_decomposition_type(O,Int(q))
    e=tr[1][2]
    nq=fmpz(q)^(divexact(d,e))
    if nq> bound
      break
    end
    l=length(extra_list)
    for i=1:l
      n=extra_list[i][2]*nq
      if n> bound
        continue
      end
      push!(extra_list, (extra_list[i][1]*q, n))
    end
  end

  final_list=Tuple{Int,fmpz}[]
  l=length(list)
  for (el,norm) in extra_list
    for i=1:l
      if fmpz(list[i])^d*norm>bound
        continue
      end
      push!(final_list, (list[i]*el, fmpz(list[i])^d*norm))
    end
  end
  return final_list

end

function squarefree_for_conductors(O::NfOrd, n::Int, deg::Int ; coprime_to::Array{fmpz,1}=fmpz[])

  sqf= trues(n)
  primes= trues(n)

  #remove primes that can be wildly ramified or
  #that are ramified in the base field
  for x in coprime_to
    el=Int(x)
    t=el
    while t<= n
      @inbounds sqf[t]=false
      @inbounds primes[t]=false
      t+=el
    end
  end

  #sieving procedure

  if !(2 in coprime_to)
    dt=prime_decomposition_type(O,2)
    if gcd(2^dt[1][1]-1, deg)==1
      j=2
      while j<=n
        @inbounds sqf[j]=false
        @inbounds primes[j]=false
        j+=2
      end
    else
      i=2
      s=4
      while s<=n
        @inbounds primes[s]=false
        s+=2
      end
      s=4
      while s<=n
        @inbounds sqf[s]=false
        s+=4
      end
    end
  end
  i=3
  b=Base.sqrt(n)
  while i<=b
    if primes[i]
      dt=prime_decomposition_type(O,i)
      if gcd(deg,i^dt[1][1]-1)==1
        @inbounds primes[i]=false
        @inbounds sqf[i]=false
        j=i
        while j<= n
         @inbounds primes[j]=false
         @inbounds sqf[j]=false
         j+=i
        end
      else
        j=i
        while j<= n
          @inbounds primes[j]=false
          j+=i
        end
        j=i^2
        t=2*j
        while j<= n
          @inbounds sqf[j]=false
          j+=t
        end
      end
    end
    i+=2
  end
  while i<=n
    if primes[i]
      dt=prime_decomposition_type(O,i)
      if gcd(deg,i^dt[1][1]-1)==1
        @inbounds sqf[i]=false
        j=i
        while j<= n
         @inbounds sqf[j]=false
         j+=i
        end
      end
    end
    i+=2
  end

  if degree(O)==1
    i=2
    while i<=length(sqf)
      @inbounds sqf[i]=false
      i+=4
    end

  end

  return Int[i for i=1:length(sqf) if sqf[i]]

end

function conductors_tame(O::NfOrd, n::Int, bound::fmpz)

  if n==2
    return tame_conductors_degree_2(O,bound)
  end
  #
  #  First, conductors coprime to the ramified primes and to the
  #  degree of the extension we are searching for.
  #
  d=degree(O)
  K=nf(O)
  wild_ram=collect(keys(factor(fmpz(n)).fac))
  ram_primes=collect(keys(factor(O.disc).fac))
  filter!(x -> !divisible(fmpz(n),x), ram_primes)
  coprime_to=cat(1,ram_primes, wild_ram)
  sort!(ram_primes)
  m=minimum(wild_ram)
  k=divexact(n,m)
  b1=Int(root(fmpz(bound),Int(degree(O)*(m-1)*k)))

  list= squarefree_for_conductors(O, b1, n, coprime_to=coprime_to)

  extra_list=Tuple{Int, Int}[(1,1)]
  for q in ram_primes
    tr=prime_decomposition_type(O,Int(q))
    f=tr[1][1]
    nq=Int(q)^f
    if gcd(nq-1,n)==1
      continue
    end
    if nq> bound
      break
    end
    l=length(extra_list)
    for i=1:l
      no=extra_list[i][2]*nq
      if no> bound
        continue
      end
      push!(extra_list, (extra_list[i][1]*q, no))
    end
  end

  final_list=Tuple{Int,fmpz}[]
  l=length(list)
  for (el,norm) in extra_list
    for i=1:l
      if (list[i]^d)*norm>bound
        continue
      end
      push!(final_list, (list[i]*el, (list[i]^(Int((m-1)*k)*d))*norm))
    end
  end

  return final_list
end

function conductors(O::NfOrd, n::Int, bound::fmpz, tame::Bool=false)

  K=nf(O)
  d=degree(O)
  wild_ram=collect(keys(factor(fmpz(n)).fac))

  #
  # First, conductors for tamely ramified extensions
  #

  list=conductors_tame(O,n,bound)

  if tame
    return Tuple{Int, Dict{NfOrdIdl, Int}}[(x[1], Dict{NfOrdIdl, Int}()) for x in list]
  end
  #
  # now, we have to multiply the obtained conductors by proper powers of wildly ramified ideals.
  #
  wild_list=Tuple{Int, Dict{NfOrdIdl, Int}, fmpz}[(1, Dict{NfOrdIdl, Int}(),1)]
  for q in wild_ram
    lp=prime_decomposition(O,Int(q))
    fq=divexact(d,lp[1][2]*length(lp))
    l=length(wild_list)
    sq=fmpz(q)^(divexact(d,lp[1][2])) #norm of the squarefree part of the integer q
    #=
      we have to use the conductor discriminant formula to understand the maximal possible exponent of q.
      Let ap be the exponent of p in the relative discriminant, let m be the conductor and h_(m,C) the cardinality of the
      quotient of ray class group by its subgroup C.
      Then
          ap= v_p(m)h_(m,C)- sum_{i=1:v_p(m)} h_(m/p^i, C)
      Since m is the conductor, h_(m/p^i, C)<= h_(m,C)/q.
      Consequently, we get
        v_p(m)<= (q*ap)/(h_(m,C)*(q-1))
      To find ap, it is enough to compute a logarithm.
    =#
    nisc= gcd(q^(fq)-1,n)
    if nisc!=1
      nbound=n+n*lp[1][2]*valuation(n,q)-1
    else
      nbound=n+n*lp[1][2]*valuation(n,q)-div(fmpz(n), q^(valuation(n,q)))
    end
    obound=flog(bound,sq)
    bound_max_ap= min(nbound, obound)  #bound on ap
    bound_max_exp=div(q*bound_max_ap, n*(q-1)) #bound on the exponent in the conductor

    if nisc != 1
      fnisc=minimum(keys(factor(nisc).fac))
      nq=sq^((fnisc-1)*(divexact(n, fnisc)))
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (q*wild_list[s][1], wild_list[s][2], nn))
      end
    end
    for i=2:bound_max_exp
      d1=Dict{NfOrdIdl, Int}()
      for j=1:length(lp)
        d1[lp[j][1]]=i
      end
      nq= sq^(i*(q-1)*divexact(n,q))
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        d2=merge(max, d1, wild_list[s][2])
        if nisc!=1
          push!(wild_list, (q*wild_list[s][1], d2, nn))
        else
          push!(wild_list, (wild_list[s][1], d2, nn))
        end
      end
    end
  end

  #the final list
  final_list=Tuple{Int, Dict{NfOrdIdl, Int}}[]
  for (el, nm) in list
    for (q,d,nm2) in wild_list
      if nm*nm2 > bound
        continue
      end
      push!(final_list, (el*q,d))
    end
  end

  return final_list

end

function conductors_tameQQ(O::NfOrd, n::Int, bound::fmpz)

  #
  #  First, conductors coprime to the ramified primes and to the
  #  degree of the extension we are searching for.
  #

  wild_ram=collect(keys(factor(fmpz(n)).fac))
  m=minimum(wild_ram)
  k=divexact(n,m)
  b1=Int(root(fmpz(bound),Int((m-1)*k)))

  return squarefree_for_conductors(O, b1, n, coprime_to=wild_ram)

end

function conductorsQQ(O::NfOrd, n::Int, bound::fmpz, tame::Bool=false)

  K=nf(O)
  d=degree(O)
  wild_ram=collect(keys(factor(fmpz(n)).fac))

  #
  # First, conductors for tamely ramified extensions
  #

  list=conductors_tameQQ(O,n,bound)

  if tame
    return list
  end
  #
  # now, we have to multiply the obtained conductors by proper powers of wildly ramified ideals.
  #
  wild_list=Tuple{Int, Int, fmpz}[(1, 1,1)]
  for q in wild_ram
    l=length(wild_list)
    #=
      we have to use the conductor discriminant formula to understand the maximal possible exponent of q.
      Let ap be the exponent of p in the relative discriminant, let m be the conductor and h_(m,C) the cardinality of the
      quotient of ray class group by its subgroup C.
      Then
          ap= v_p(m)h_(m,C)- sum_{i=1:v_p(m)} h_(m/p^i, C)
      Since m is the conductor, h_(m/p^i, C)<= h_(m,C)/q.
      Consequently, we get
        v_p(m)<= (q*ap)/(h_(m,C)*(q-1))
      To find ap, it is enough to compute a logarithm.
    =#
    nisc= gcd(q-1,n)
    if nisc!=1
      nbound=n+n*valuation(n,q)-1
    else
      nbound=n+n*valuation(n,q)-div(fmpz(n), q^(valuation(n,q)))
    end
    obound=flog(bound,q)
    nnbound=valuation_bound_discriminant(n,q)
    bound_max_ap= min(nbound, obound, nnbound)  #bound on ap

    bound_max_exp=div(q*bound_max_ap, n*(q-1)) #bound on the exponent in the conductor

    if nisc != 1
      fnisc=minimum(keys(factor(nisc).fac))
      nq=q^((fnisc-1)*(divexact(n, fnisc)))
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (q*wild_list[s][1], wild_list[s][2], nn))
      end
    end
    t=(q-1)*divexact(n,q)
    for i=2:bound_max_exp
      nq= q^(i*t)
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (wild_list[s][1], wild_list[s][2]*q^i, nn))
      end
    end
  end

  #the final list
  final_list=Int[]
  exps=Int((minimum(wild_ram)-1)*divexact(n,minimum(wild_ram)))
  for el in list
    for (q,d,nm2) in wild_list
      if el^exps*nm2 > bound
        continue
      end
      push!(final_list, (el*q*d))
    end
  end

  return final_list

end

###############################################################################
#
#  Abelian extensions
#
###############################################################################

function abelian_extensions(O::Union{FlintIntegerRing, FlintRationalField},
                            gtype::Vector{Int}, discriminant_bound::fmpz;
                            real::Bool = false,
                            tame::Bool = false)

  Qx, x = PolynomialRing(FlintQQ, "x")
  K, _ = NumberField(x - 1, "a")
  OK = maximal_order(K)
  l = abelian_extensions(OK, gtype, discriminant_bound,
                         real = real,
                         tame = tame)

  newlist = Vector{NfAbsNS}(length(l))
  for j in 1:length(l)
    newlist[j], _ = number_field([Qx([coeff(coeff(f, i), 0) for i in 0:length(f)]) for f in l[j].abs_pol])
  end
  return newlist
end

function abelian_extensions(O::NfOrd, gtype::Array{Int,1}, absolute_discriminant_bound::fmpz; real::Bool=false, tame::Bool=false, with_autos::Type{Val{T}} = Val{false}) where T

  K=nf(O)
  @assert degree(K)==1
  n=prod(gtype)

  expo=lcm(gtype)

  #Getting conductors
  l_conductors=conductorsQQ(O,n,absolute_discriminant_bound, tame)
  @vprint :QuadraticExt 1 "Number of conductors: $(length(l_conductors)) \n"
  if with_autos==Val{true}
    fields = Tuple{NfRel_ns, Array{NfRel_nsToNfRel_nsMor,1}}[]
  else
    fields = NfRel_ns[]
  end

  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :QuadraticExt 1 "Conductor: $k \n"
    @vprint :QuadraticExt 1 "Left: $(length(l_conductors) - i)\n"
    r,mr=Hecke.ray_class_groupQQ(O,k,!real,expo)
    if !Hecke._are_there_subs(r,gtype)
      continue
    end
    ls=subgroups(r,quotype=gtype, fun= (x, y) -> quo(x, y, false)[2])
    for s in ls
      C=ray_class_field(mr, s)
      if Hecke._is_conductor_minQQ(C,n) && Hecke.discriminant_conductorQQ(O,C,k,absolute_discriminant_bound,n)
        @vprint :QuadraticExt 1 "New Field \n"
        L=number_field(C)
        if with_autos==Val{true}
          push!(fields,(L, automorphism_groupQQ(C)))
        else
          push!(fields, L)
        end
      end
    end
  end

  return fields

end

function abelian_normal_extensions(O::NfOrd, gtype::Array{Int,1}, absolute_discriminant_bound::fmpz; ramified_at_infplc::Bool=true, tame::Bool=false, absolute_galois_group::Symbol = :all,  with_autos::Type{Val{T}} = Val{false}, autos::Array{NfToNfMor,1}=NfToNfMor[]) where T

  d=degree(O)
  if d==1
    return abelian_extensions(O, gtype, absolute_discriminant_bound, real=!ramified_at_infplc, tame=tame, with_autos=with_autos)
  end

  K=nf(O)
  n=prod(gtype)
  inf_plc=InfPlc[]

  if with_autos==Val{true}
    fields=Tuple{NfRel_ns,Vector{NfRel_nsToNfRel_nsMor{nf_elem}}}[]
  else
    fields=NfRel_ns[]
  end

  bound = abs(div(absolute_discriminant_bound,discriminant(O)^n))
  if bound==0
    return fields
  end

  if mod(n,2)==0 && ramified_at_infplc
    inf_plc=real_places(K)
  end

  expo=lcm(gtype)
  C,mC=class_group(O)
  cgrp= gcd(n,order(C))!=1
  if cgrp
    S = prime_ideals_up_to(O, max(100,12*clog(abs(discriminant(O)),10)^2))
  end
  allow_cache!(mC)

  #
  # Getting a small set of generators
  # for the automorphisms group
  #
  if isempty(autos)
    Aut=Hecke.automorphisms(K)
    @assert length(Aut)==degree(O) #The field is normal over Q
    gens = Hecke.small_generating_set(Aut)
  else
    gens = autos
  end

  #Getting conductors
  l_conductors=conductors(O,n,bound, tame)
  @vprint :QuadraticExt 1 "Number of conductors: $(length(l_conductors)) \n"

  #Now, the big loop
  for (i, k) in enumerate(l_conductors)
    @vprint :QuadraticExt 1 "Conductor: $k \n"
    @vprint :QuadraticExt 1 "Left: $(length(l_conductors) - i)\n"
    r,mr=ray_class_group_quo(O,expo,k[1], k[2],inf_plc)
    if !_are_there_subs(r,gtype)
      continue
    end
    if cgrp
        mr.prime_ideal_cache = S
    end
    act=_act_on_ray_class(mr,gens)
    ls=stable_subgroups(r,gtype,act, op=(x, y) -> quo(x, y, false)[2])
    a=_min_wild(k[2])*k[1]
    totally_positive_generators(mr,a)
    for s in ls
      @hassert :QuadraticExt 1 order(codomain(s))==n
      C=ray_class_field(mr, s)
      if Hecke._is_conductor_min_normal(C,a) && Hecke.discriminant_conductor(O,C,a,mr,bound,n)
        @vprint :QuadraticExt 1 "New Field \n"
        L = number_field(C)
        if absolute_galois_group != :all
          autabs = absolute_automorphism_group(C, gens)
          G = generic_group(autabs, *)
          if absolute_galois_group == :_16T7
            if isisomorphic_to_16T7(G)
              @vprint :QuadraticExt 1 "I found a field with Galois group 16T7 (C_2 x Q_8)\n"
              @vtime :QuadraticExt 1 push!(fields, L)
            end
          elseif absolute_galois_group == :_8T5
            if isisomorphic_to_8T5(G)
              @vprint :QuadraticExt 1 "I found a field with Galois group 8T5 (Q_8)\n"
              @vtime :QuadraticExt 1 push!(fields, L)
            end
          else
            error("Group not supported (yet)")
          end
        else
          if with_autos==Val{true}
            push!(fields, (L,absolute_automorphism_group(C,gens)))
          else
            push!(fields, L)
          end
        end
      end
    end
  end

  return fields

end

###############################################################################
#
#  Quadratic Extensions of Q
#
###############################################################################

function quadratic_extensions(bound::Int; tame::Bool=false, real::Bool=false, complex::Bool=false, u::UnitRange{Int}=-1:0)

  @assert !(real && complex)
  Qx,x=PolynomialRing(FlintQQ, "x")
  sqf=squarefree_up_to(bound);
  if real
    deleteat!(sqf,1)
  elseif complex
    sqf=Int[-i for i in sqf]
  else
    sqf= vcat(sqf[2:end], Int[-i for i in sqf])
  end
  if tame
    filter!( x -> mod(x,4)==1, sqf)
    return ( number_field(x^2-x+divexact(1-i,4), cached=false)[1] for i in sqf)
  end
  final_list=Int[]
  for i=1:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list,sqf[i])
      continue
    end
    if mod(sqf[i],4)==1
      @views push!(final_list,sqf[i])
    end
  end
  if u==-1:0
    return ( mod(i,4)!=1 ? number_field(x^2-i, cached=false)[1] : number_field(x^2-x+divexact(1-i,4), cached=false)[1] for i in final_list)
  else
    return ( mod(final_list[i],4)!=1 ? number_field(x^2-final_list[i], cached=false)[1] : number_field(x^2-x+divexact(1-final_list[i],4), cached=false)[1] for i in u)
  end

end

function quadratic_extensions(bound::Int; tame::Bool=false, real::Bool=false, complex::Bool=false, with_autos::Type{Val{T}}=Val{false}) where T

  @assert !(real && complex)
  Qx,x=PolynomialRing(FlintQQ, "x")
  sqf=squarefree_up_to(bound)
  if real
    deleteat!(sqf,1)
  elseif complex
    sqf=Int[-i for i in sqf]
  else
    sqf= vcat(sqf[2:end], Int[-i for i in sqf])
  end
  if tame
    filter!( x -> mod(x,4)==1, sqf)
    return ( number_field(x^2-x+divexact(1-i,4), cached=false)[1] for i in sqf)
  end
  final_list=Int[]
  for i=1:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list,sqf[i])
      continue
    end
    if mod(sqf[i],4)==1
      @views push!(final_list,sqf[i])
    end
  end
  return ( mod(i,4)!=1 ? number_field(x^2-i, cached=false)[1] : number_field(x^2-x+divexact(1-i,4), cached=false)[1] for i in final_list)

end

function _quad_exts_as_ab_exts(bound::Int)

  Qy,y=PolynomialRing(FlintQQ, "y")
  K,a=NumberField(y-1)
  Kx,x=PolynomialRing(K,"x")
  sqf=squarefree_up_to(bound)
  sqf= vcat(sqf[2:end], Int[-i for i in sqf])
  final_list=Int[]
  for i=1:length(sqf)
    if abs(sqf[i]*4)< bound
      @views push!(final_list,sqf[i])
      continue
    end
    if mod(sqf[i],4)==1
      @views push!(final_list,sqf[i])
    end
  end
  return ( _quad_ext_with_auto(Kx,i) for i in final_list)

end

function _quad_ext_with_auto(Kx,a::Int)
  x=gen(Kx)
  if a % 4 == 1
    L, lg = number_field([x^2-x+divexact(1-a,4)])
    b=lg[1]
    mL=NfRel_nsToNfRel_nsMor(L,L, [1-b])
  else
    L, lg = number_field([x^2-a])
    b = lg[1]
    mL=NfRel_nsToNfRel_nsMor(L,L, [-b])
  end
  return (L, NfRel_nsToNfRel_nsMor[mL])

end

###############################################################################
#
#  C2 x C2 extensions of Q
#
###############################################################################

function C22_extensions(bound::Int)

  Qx,x=PolynomialRing(FlintZZ, "x")
  K,_=NumberField(x-1)
  Kx,x=PolynomialRing(K,"x", cached=false)
  b1=ceil(Int,Base.sqrt(bound))
  n=2*b1+1
  pairs=_find_pairs(bound)
  poszero=b1+1
  return (_ext(Kx,x,i-poszero,j-poszero) for i=1:n for j=i+1:n if pairs[i,j])

end

function _C22_exts_abexts(bound::Int)
  Qx,x=PolynomialRing(FlintZZ, "x")
  K,_=NumberField(x-1)
  Kx,x=PolynomialRing(K,"x", cached=false)
  pairs=_find_pairs(bound)
  b1=ceil(Int,Base.sqrt(bound))
  n=2*b1+1
  poszero=b1+1
  return (_ext_with_autos(Kx,x,i-poszero,j-poszero) for i=1:n for j=i+1:n if pairs[i,j])

end

function _ext_with_autos(Ox,x,i,j)

  y1=mod(i,4)
  pol1=Ox(1)
  if y1!=1
    pol1=x^2-i
  else
    pol1=x^2-x+divexact(1-i,4)
  end
  y2=mod(j,4)
  pol2=Ox(1)
  if y2!=1
    pol2=x^2-j
  else
    pol2=x^2-x+divexact(1-j,4)
  end
  L, lg= number_field([pol1,pol2])
  if y1!=1
    a1=-lg[1]
  else
    a1=1-lg[1]
  end
  if y2!=1
    a2=-lg[2]
  else
    a2=1-lg[2]
  end
  mL1=NfRel_nsToNfRel_nsMor(L,L, [a1,lg[2]])
  mL2=NfRel_nsToNfRel_nsMor(L,L, [lg[1],a2])
  return (L, [mL1,mL2])

end

function _find_pairs(bound::Int)

  b1=ceil(Int,Base.sqrt(bound))
  n=2*b1+1
  pairs=trues(n,n)
  poszero=b1+1

  for i=1:n
    pairs[i,poszero]=false
    pairs[i,b1+2]=false
    pairs[poszero,i]=false
    pairs[b1+2,i]=false
  end

  #sieve for squarefree
  list= trues(b1)
  i=2
  b=Base.sqrt(b1)
  while i<=b
    if list[i]
      j=i^2
      if !list[j]
        i+=1
        continue
      else
        list[j]=false
        t=2*j
        while t <= b1
          list[t]=false
          t+=j
        end
      end
    end
    i+=1
  end
  #now translate it into the matrix
  for i=1:b1
    if !list[i]
      pairs[poszero+i,1:n]=false
      pairs[1:n,poszero+i]=false
      pairs[poszero-i,1:n]=false
      pairs[1:n,poszero-i]=false
    end
  end
  #check
#=
  for i=1:n
    for j=1:n
      if pairs[i,j]
        @assert issquarefree(i-poszero)
        @assert issquarefree(j-poszero)
      end
    end
  end
=#
  #removing diagonal
  for i=1:n
    pairs[i,i]=false
  end

  #counting extensions
  for i=1:n
    for j=i+1:n
      if pairs[i,j]
        pairs[j,i]=false
        third=divexact((i-poszero)*(j-poszero), gcd(i-poszero, j-poszero)^2)
        if abs(third)>b1
          pairs[i,j]=false
        else
          y=third+poszero
          pairs[y,i]=false
          pairs[y,j]=false
          pairs[i,y]=false
          pairs[j,y]=false
          first=i-poszero
          second=j-poszero
          d1=_discn(first)
          d2=_discn(second)
          g1=d1*d2
          if abs(g1)<b1
            continue
          else
            d3=_discn(third)
            g2=d1*d3
            if abs(g2)<b1
              continue
            else
              g3=d2*d3
              if abs(g3)<b1
                continue
              else
                g=gcd(g1,g2,g3)
                if abs(g)<b1
                  continue
                elseif check_disc(first, second, third, bound)
                  pairs[i,j]=false
                end
              end
            end
          end
        end
      end
    end
  end
  return pairs

end

function check_disc(a::Int,b::Int,c::Int,bound::Int)

  if mod(a,4)!=2 && mod(b,4)!=2
    return true
  end
  if mod(a,4)==1 || mod(b,4)==1 || mod(c,4)==1
    return true
  end
  if 64*a*b*c>=bound
    return true
  else
    return false
  end

end

#given integers i,j, this function returns the extension Q(sqrt(i))Q(sqrt(j))
function _ext(Ox,x,i,j)

  y=mod(i,4)
  pol1=Ox(1)
  if y!=1
    pol1=x^2-i
  else
    pol1=x^2-x+divexact(1-i,4)
  end
  y=mod(j,4)
  pol2=Ox(1)
  if y!=1
    pol2=x^2-j
  else
    pol2=x^2-x+divexact(1-j,4)
  end
  return number_field([pol1,pol2])

end

# Given a squarefree integer n, this function computes the absolute value
# of the discriminant of the extension Q(sqrt(n))
function _discn(n::Int)

  x=mod(n,4)
  if x!=1
    return 4*n
  else
    return n
  end

end

function C22_extensions_tame_real(bound::Int)

  Qx,x=PolynomialRing(FlintZZ, "x")
  K,_=NumberField(x-1)
  Kx,x=PolynomialRing(K,"x", cached=false)

  b1=floor(Int,Base.sqrt(bound))
  n=b1
  pairs=trues(b1,b1)

  #sieve for squarefree number congruent to 1 mod 4
  i=2
  k=4
  while k<=b1
    pairs[1:i,i]=false
    pairs[1:k,k]=false
    pairs[i,i+1:n]=false
    pairs[k,k+1:n]=false
    i+=4
    k+=4
  end
  i=3
  b=Base.sqrt(b1)
  while i<=b
    if pairs[1,i]
      j=i^2
      if !pairs[1,j]
        i+=2
        continue
      else
        pairs[1:j,j]=false
        pairs[j,1:j]=false
        t=2*j
        while t <= b1
          pairs[t,t:n]=false
          pairs[1:t,t]=false
          t+=j
        end
      end
    end
    i+=2
  end
  k=3
  while k<=b1
    pairs[1:k,k]=false
    pairs[k,k:n]=false
    k+=4
  end
  pairs[1,1:n]=false

  #counting extensions
  for i=5:b1
    for j=i+1:b1
      if pairs[i,j]
        k=divexact(i*j, gcd(i, j)^2)
        if k>b1
          pairs[i,j]=false
        else
          if k>i
            pairs[i,k]=false
          else
            pairs[k,i]=false
          end
          if k>j
            pairs[j,k]=false
          else
            pairs[k,j]=false
          end
          if i*j*k<bound
            continue
          else
            pairs[i,j]=false
          end
        end
      end
    end
  end

  return (_ext(Kx,x,i,j) for i=5:b1 for j=i+1:b1 if pairs[i,j])

end

###############################################################################
#
#  Dycyclic Dic3
#
###############################################################################

#K is a C4 field
function Dic3_extensions(absolute_bound::fmpz, K::AnticNumberField)

  O=maximal_order(K)
  D=abs(discriminant(O))

  C,mC=class_group(O)
  allow_cache!(mC)
  cgrp=false
  if gcd(C.snf[end],3)!=1
    cgrp=true
    S = prime_ideals_up_to(O, max(100,100*clog(D,10)^2))
  end
  gens=Hecke.automorphisms(K)
  gens=small_generating_set(gens)

  #Getting conductors
  bound = div(absolute_bound, D^3)
  l_conductors=conductors(O,3, bound)
  @vprint :QuadraticExt "Number of conductors: $(length(l_conductors)) \n"

  res = []

  #Now, the big loop
  for k in l_conductors
    r,mr=ray_class_group_quo(O,3,k[1], k[2])
    if cgrp
      mr.prime_ideal_cache = S
    end
    act=_act_on_ray_class(mr,gens)
    ls=stable_subgroups(r,[3],act, op=(x, y) -> quo(x, y, false)[2])
    a=_min_wild(k[2])*k[1]
    totally_positive_generators(mr,a)
    for s in ls
      if _trivial_action(s,act,3)
        continue
      end
      C=ray_class_field(mr, s)
      if Hecke._is_conductor_min_normal(C,a) && Hecke.discriminant_conductor(O,C,a,mr,bound,3)
        @vprint :QuadraticExt 1 "New Field"
        L=number_field(C)
        SS=Hecke.simple_extension(L)[1]
        F=absolute_field(SS)[1]
        push!(res, (F.pol, prod([ fmpz(p)^e for (p, e) in C.absolute_discriminant])))
      end
    end
  end
  return res
end

###############################################################################
#
#  Dihedral group D5 with equation of degree 5
#
###############################################################################

function conductorsD5(O::NfOrd, bound_non_normal::fmpz)

  D=abs(discriminant(O))
  ram_primes=collect(keys(factor(O.disc).fac))
  coprime_to=cat(1,ram_primes, fmpz(5))
  sort!(ram_primes)
  b=root(bound_non_normal,2)
  b1=root(div(b,D),2)
  #
  # First, conductors for tamely ramified extensions
  #

  sqf_list= squarefree_for_conductors(O, Int(b1), 5, coprime_to=coprime_to)
  l=length(sqf_list)
  #
  # now, we have to multiply the obtained conductors by proper powers of wildly ramified ideals.
  #
  lp=prime_decomposition(O,5)
  final_list=Tuple{Int, Dict{NfOrdIdl, Int}}[]
  if 5 in ram_primes
    bound_max_exp=flog(b1,5)
    final_list=[(x, Dict{NfOrdIdl, Int}()) for x in sqf_list]
    for i=2:bound_max_exp
      isodd(i) && continue
      d1=Dict{NfOrdIdl, Int}()
      for j=1:length(lp)
        d1[lp[j][1]]=i
      end
      min_id=div(i,2)
      for s=1:l
        if sqf_list[s]*5^min_id<b1
          push!(final_list,(sqf_list[s],d1))
        end
      end
    end
  else
    bound_max_exp=flog(b1,5)
    final_list=[(x, Dict{NfOrdIdl, Int}()) for x in sqf_list]
    for i=2:bound_max_exp
      d1=Dict{NfOrdIdl, Int}()
      for j=1:length(lp)
        d1[lp[j][1]]=i
      end
      for s=1:l
        if sqf_list[s]*5^i<b1
          push!(final_list,(sqf_list[s],d1))
        end
      end
    end
  end
  return final_list

end

#The input are the absolute bound for the non-normal extension of degree 5 and the list of the quadratic fields
function D5_extensions(absolute_bound::fmpz, quad_fields)

  len=length(quad_fields)
  z = []
  for K in quad_fields
    len-=1

     @vprint :QuadraticExt 1 "Field: $K\n"
     @vprint :QuadraticExt 1 "Remaining Fields: $(len)\n"
    append!(z, single_D5_extensions(absolute_bound, K))
  end
  return z
end

function D5_extensions(absolute_bound::fmpz, quad_fields, f::IOStream)

  len=length(quad_fields)
  for K in quad_fields
    len-=1

    @vprint :QuadraticExt 1  "Field: $K\n"
    @vprint :QuadraticExt 1 "Remaining Fields: $(len)\n"
    for g in single_D5_extensions(absolute_bound, K)
      Base.write(f, "($g)\n" )
    end
  end
  return len
end

function single_D5_extensions(absolute_bound::fmpz, K::AnticNumberField)

  O=maximal_order(K)
  D=abs(discriminant(O))

  C,mC=class_group(O)
  allow_cache!(mC)
  cgrp=false
  if gcd(C.snf[end],5)!=1
    cgrp=true
    S = prime_ideals_up_to(O, max(100,100*clog(D,10)^2))
  end
  gens=Hecke.automorphisms(K)
  a=gen(K)
  if gens[1](a)==a
    deleteat!(gens,1)
  else
    deleteat!(gens,2)
  end

  z = []

  #Getting conductors
  l_conductors=conductorsD5(O,absolute_bound)
#  @vprint :QuadraticExt "Number of conductors: $(length(l_conductors)) \n"

  #Now, the big loop
  for k in l_conductors
    r,mr=ray_class_group_quo(O,5,k[1], k[2])
    if cgrp
      mr.prime_ideal_cache = S
    end
    act=_act_on_ray_class(mr,gens)
    ls=stable_subgroups(r,[5],act, op=(x, y) -> quo(x, y, false)[2])
    a=_min_wild(k[2])*k[1]
    for s in ls
      if _trivial_action(s,act,5)
        continue
      end
      C=ray_class_field(mr, s)
      if Hecke._is_conductor_min_normal(C,a)
        @vprint :QuadraticExt 1 "New Field"
        L=number_field(C)
        auto=Hecke.extend_aut(C, gens[1])
        pol=_quintic_ext(auto)
        push!(z, (pol,D^2*minimum(mr.modulus_fin)^4))
      end
    end
  end
  return z
end

function _quintic_ext(auto)#::NfRel_nsToNfRel_nsMor)

  #first, find a generator of the simple extension
  L=domain(auto)
  S,mS=simple_extension(L)
  x=mS(gen(S))

  #Find primitive element
  pr_el=x+auto(x)

  #Take minimal polynomial; I need to embed the element in the absolute extension
  pol=absolute_minpoly(pr_el)
  if degree(pol)==15
    return pol
  else
    pr_el=x*(auto(x))
    return absolute_minpoly(pr_el)
  end

end

###############################################################################
#
#  Dihedral group Dn with n odd
#
###############################################################################

function conductorsDn(O::NfOrd, n::Int, bound::fmpz, tame::Bool=false)

  K=nf(O)
  d=degree(O)
  wild_ram=collect(keys(factor(fmpz(n)).fac))
  ram_primes=collect(keys(factor(O.disc).fac))
  coprime_to=cat(1,ram_primes, wild_ram)
  sort!(ram_primes)
  m=minimum(wild_ram)
  k=divexact(n,m)
  b1=Int(root(fmpz(bound),Int(2*(minimum(wild_ram)-1)*k)))

  #
  # First, conductors for tamely ramified extensions
  #

  list= squarefree_for_conductors(O, b1, n, coprime_to=coprime_to)

  if tame
    error("Not yet implemented")
  end

  #
  # now, we have to multiply the obtained conductors by proper powers of wildly ramified ideals.
  #
  wild_list=Tuple{Int, Dict{NfOrdIdl, Int}, Int}[(1, Dict{NfOrdIdl, Int}(),1)]
  for q in wild_ram
    lp=prime_decomposition(O,q)
    l=length(wild_list)
    #=
      we have to use the conductor discriminant formula to understand the maximal possible exponent of q.
      Let ap be the exponent of p in the relative discriminant, let m be the conductor and h_(m,C) the cardinality of the
      quotient of ray class group by its subgroup C.
      Then
          ap= v_p(m)h_(m,C)- sum_{i=1:v_p(m)} h_(m/p^i, C)
      Since m is the conductor, h_(m/p^i, C)<= h_(m,C)/q.
      Consequently, we get
        v_p(m)<= (q*ap)/(h_(m,C)*(q-1))
      To find ap, it is enough to compute a logarithm.
    =#
    sq=fmpz(q)^(divexact(2,lp[1][2]))
    fq=divexact(2,length(lp)*lp[1][2])
    bound_max_ap=clog(bound,sq) #bound on ap
    bound_max_exp=divexact(q*bound_max_ap, n*(q-1)) #bound on the exponent in the conductor
    nisc= gcd(q^fq-1,n)!=1
    if !(q in ram_primes) && nisc
      for s=1:l
        nn=sq*wild_list[s][3]
        if nn>bound
          continue
        end
        push!(wild_list, (q*wild_list[s][1], wild_list[s][2], nn))
      end
    end
    for i=2:bound_max_exp
      if q in ram_primes && i % 2 ==1
        continue
      end
      d1=Dict{NfOrdIdl, Int}()
      for j=1:length(lp)
        d1[lp[j][1]]=i
      end
      nq= sq^i
      for s=1:l
        nn=nq*wild_list[s][3]
        if nn>bound
          continue
        end
        d2=merge(max, d1, wild_list[s][2])
        if nisc
          push!(wild_list, (q*wild_list[s][1], d2, nn))
        else
          push!(wild_list, (wild_list[s][1], d2, nn))
        end
      end
    end
  end

  #the final list
  final_list=Set(Tuple{Int, Dict{NfOrdIdl, Int}}[])
  for el in list
    for (q,d1,nm2) in wild_list
      if el^2*nm2 > bound
        continue
      end
      push!(final_list, (el*q,d1))
    end
  end

  return final_list

end

function Dn_extensions(n::Int, absolute_bound::fmpz; totally_real::Bool=false, complex::Bool=false, tame::Bool=false)

  if mod(n,2)==0
    error("Not yet implemented")
  end
  @assert !(totally_real && complex)

  bound_quadratic= Int(root(absolute_bound, n))
  list_quad=quadratic_extensions(bound_quadratic, tame=tame, real=totally_real, complex=complex)
  return Dn_extensions(n,absolute_bound, list_quad, tame=tame)

end

function Dn_extensions(n::Int, absolute_bound::fmpz, list_quad ; tame::Bool=false)
  @assert absolute_bound>0
  len=length(list_quad)
  fieldslist=Tuple{NfRel_ns,  Array{NfRel_nsToNfRel_nsMor{nf_elem},1},fmpz, Array{fmpz,1}}[]

  for K in list_quad
    len-=1

    @vprint :QuadraticExt 1 "Field: $K\n"
    @vprint :QuadraticExt 1 "Fields left: $len\n"
    O=maximal_order(K)
    D=abs(discriminant(O))
    if D^n>absolute_bound
      continue
    end
    bound = div(absolute_bound, abs(D)^n)

    C,mC=class_group(O)
    allow_cache!(mC)
    cgrp=false
    if gcd(C.snf[end],n)!=1
      cgrp=true
      S = prime_ideals_up_to(O, max(30,12*clog(D,10)^2))
    end

    gens=Hecke.automorphisms(K)
    a=gen(K)
    if gens[1](a)==a
      deleteat!(gens,1)
    else
      deleteat!(gens,2)
    end

    #Getting conductors
    l_conductors=conductorsDn(O,n,bound, tame)

    @vprint :QuadraticExt "Number of conductors: $(length(l_conductors)) \n"

    #Now, the big loop
    for k in l_conductors
      r,mr=ray_class_group_quo(O,n,k[1], k[2])
      if !_are_there_subs(r,[n])
        continue
      end
      if cgrp
        mr.prime_ideal_cache = S
      end
      act=_act_on_ray_class(mr,gens)
      ls=stable_subgroups(r,[n],act, op=(x, y) -> quo(x, y, false)[2])
      a=_min_wild(k[2])*k[1]
      for s in ls
        if _trivial_action(s,act,n)
          continue
        end
        C=ray_class_field(mr, s)
        if Hecke._is_conductor_min_normal(C,a) && Hecke.discriminant_conductor(O,C,a,mr,bound,n)
          @vprint :QuadraticExt 1 "\n New Field!\n"
          L=number_field(C)
          ram_primes=Set(collect(keys(factor(a).fac)))
          for p in keys(factor(O.disc).fac)
            push!(ram_primes, p)
          end
          aut=absolute_automorphism_group(C,gens)
          push!(fieldslist,(L,aut, evaluate(FacElem(C.absolute_discriminant)), collect(ram_primes)))
        end
      end
    end

  end

  return fieldslist

end

function _trivial_action(s::GrpAbFinGenMap, act::Array{GrpAbFinGenMap,1}, n::Int)

  @assert length(act)==1
  S=codomain(s)
  T,mT=snf(S)
  @assert T.snf[1]==n
  @assert ngens(T)==1
  new_el=mT\(s(act[1](s\mT(T[1]))))
  if new_el==(n-1)*T[1]
    return false
  else
    return true
  end

end

###############################################################################
#
#  C3 x D5 extensions
#
###############################################################################

function C3xD5_extensions(non_normal_bound::fmpz)

  bound_quadratic= Int(root(non_normal_bound, 6))
  list_quad=quadratic_extensions(bound_quadratic, complex=true)
  fieldslist=Tuple{AnticNumberField, Array{fmpz,1}}[]

  for K in list_quad

    @vprint :QuadraticExt 1 "Field: $K\n"
    O=maximal_order(K)
    D=abs(discriminant(O))

    bound5=root(non_normal_bound,3)
    bound = non_normal_bound^2

    C,mC=class_group(O)
    allow_cache!(mC)
    cgrp=false
    if gcd(C.snf[end],15)!=1
      cgrp=true
      S = prime_ideals_up_to(O, max(100,12*clog(D,10)^2))
    end
    gens=Hecke.automorphisms(K)
    a=gen(K)
    if gens[1](a)==a
      deleteat!(gens,1)
    else
      deleteat!(gens,2)
    end

    #Getting conductors
    l_conductors=conductors(O,15,bound)
    @vprint :QuadraticExt "Number of conductors: $(length(l_conductors)) \n"

    #Now, the big loop
    for (i, k) in enumerate(l_conductors)
      r,mr=ray_class_group_quo(O,15,k[1], k[2])
      if !_are_there_subs(r,[15])
        continue
      end
      if cgrp
        mr.prime_ideal_cache = S
      end
      act=_act_on_ray_class(mr,gens)
      ls=stable_subgroups(r,[15],act, op=(x, y) -> quo(x, y, false)[2])
      a=_min_wild(k[2])*k[1]
      for s in ls
        if !_right_actionD5C3(s,act)
          continue
        end
        C=ray_class_field(mr, s)
        if Hecke._is_conductor_min_normal(C,a) && Hecke.discriminant_conductor(O,C,a,mr,bound,15)
          @vprint :QuadraticExt "New Field \n"
          #Before computing the field, I check if the discriminant of the $D_5$ extension is compatible
          s1=codomain(s)
          q1,mq1=quo(s1,5, false)
          C1=ray_class_field(mr, Hecke._compose(s, mq1))
          cond=conductor(C1)[1]
          condint=minimum(cond)
          if condint^4*D^2>bound5
            @vprint :QuadraticExt "Too large (without computing it) :( \n"
            continue
          end
          L=number_field(C)
          ram_primes=Set(collect(keys(factor(a).fac)))
          for p in keys(factor(O.disc).fac)
            push!(ram_primes, p)
          end
          auto=extend_aut(C,gens[1])
          T,mT=simple_extension(L)
          x=mT(gen(T))
          if (auto*auto)(x)!=x
            auto=auto^3
          end
          pol=absolute_minpoly(x+auto(x))
          if degree(pol)!=15
            pol=absolute_minpoly(x*(auto(x)))
          end
          K,a=NumberField(pol, cached= false)
          if _is_discriminant_lower(K, collect(ram_primes), non_normal_bound)
            @vprint :QuadraticExt "The field is: $pol \n"
            push!(fieldslist, (number_field(pol)[1], collect(ram_primes)))
          end
        end
      end
    end

  end

  return fieldslist

end

function _right_actionD5C3(s::GrpAbFinGenMap, act::Array{GrpAbFinGenMap,1})

  @assert length(act)==1
  S=codomain(s)
  T,mT=snf(S)
  @assert T.snf[1]==15
  @assert ngens(T)==1
  new_el=mT\(s(act[1](s\mT(T[1]))))
  if 3*new_el==12*T[1] && 5*new_el==5*T[1]
    return true
  else
    return false
  end
end

###############################################################################
#
#  S3 x C5 extensions
#
###############################################################################
function S3xC5_extensions(non_normal_bound::fmpz)

  bound_quadratic= Int(root(non_normal_bound,5))
  list_quad=quadratic_extensions(bound_quadratic, complex=true)
  return S3xC5_extensions(non_normal_bound,collect(list_quad))
end

function S3xC5_extensions(non_normal_bound::fmpz, list_quad)

  fieldslist=Tuple{AnticNumberField, Array{fmpz,1}}[]

  for K in list_quad

    @vprint :QuadraticExt 1 "Field: $K\n"
    O=maximal_order(K)
    D=abs(discriminant(O))
    bound = non_normal_bound^2

    bound3=root(non_normal_bound, 5)
    C,mC=class_group(O)
    allow_cache!(mC)
    cgrp=false
    if gcd(C.snf[end],15)!=1
      cgrp=true
      S = prime_ideals_up_to(O, max(100,12*clog(D,10)^2))
    end
    gens=Hecke.automorphisms(K)
    a=gen(K)
    if gens[1](a)==a
      deleteat!(gens,1)
    else
      deleteat!(gens,2)
    end

    #Getting conductors
    l_conductors=conductors(O,15,bound)
    @vprint :QuadraticExt "Number of conductors: $(length(l_conductors)) \n"

    #Now, the big loop
    for (i, k) in enumerate(l_conductors)
      r,mr=ray_class_group_quo(O,15,k[1], k[2])
      if !_are_there_subs(r,[15])
        continue
      end
      if cgrp
        mr.prime_ideal_cache = S
      end
      act=_act_on_ray_class(mr,gens)
      ls=stable_subgroups(r,[15],act, op=(x, y) -> quo(x, y, false)[2])
      a=_min_wild(k[2])*k[1]
      for s in ls
        if !_right_actionC5S3(s,act)
          continue
        end
        C=ray_class_field(mr, s)
        if Hecke._is_conductor_min_normal(C,a) && Hecke.discriminant_conductor(O,C,a,mr,bound,15)
          @vprint :QuadraticExt 1  "\n New Field!\n"
          #Before computing the field, I check if the discriminant of the $S_3$ extension is compatible
          s1=codomain(s)
          q1,mq1=quo(s1,3, false)
          C1=ray_class_field(mr, _compose(s, mq1))
          cond=conductor(C1)[1]
          condint=minimum(cond)
          if condint^2*D>bound3
            @vprint :QuadraticExt "Too large :( \n"
            continue
          end
          L=number_field(C)
          ram_primes=Set(collect(keys(factor(a).fac)))
          for p in keys(factor(O.disc).fac)
            push!(ram_primes, p)
          end
          auto=extend_aut(C,gens[1])
          T,mT=simple_extension(L)
          x=mT(gen(T))
          if (auto*auto)(x)!=x
            auto=auto*auto*auto*auto*auto
          end
          pol=absolute_minpoly(x+auto(x))
          if degree(pol)!=15
            pol=absolute_minpoly(x*(auto(x)))
          end
          K=number_field(pol, cached=false)[1]
          if _is_discriminant_lower(K, collect(ram_primes), non_normal_bound)
            push!(fieldslist, (K, collect(ram_primes)))
            @vprint :QuadraticExt "New candidate! \n"
            @vprint :QuadraticExt "$(pol) \n"
          else
            @vprint :QuadraticExt "Too large :( \n"
          end
        end
      end
    end
  end

  return fieldslist

end

function _right_actionC5S3(s::GrpAbFinGenMap, act::Array{GrpAbFinGenMap,1})

  @assert length(act)==1
  S=codomain(s)
  T,mT=snf(S)
  @assert T.snf[1]==15
  @assert ngens(T)==1
  new_el=mT\(s(act[1](s\mT(T[1]))))
  if 3*new_el==3*T[1] && 5*new_el==10*T[1]
    return true
  else
    return false
  end
end

###############################################################################
#
#  Semidirect product of C9 and C4
#
###############################################################################

function C9semiC4(absolute_bound::fmpz; real::Bool=false)

  Qx,x=PolynomialRing(FlintQQ, "x")
  K,a=NumberField(x-1,"a")
  O=maximal_order(K)

  l=Hecke.abelian_normal_extensions(O,[4], root(absolute_bound, 9), ramified_at_infplc=!real)
  return C9semiC4(absolute_bound, l)

end

function C9semiC4(absolute_bound::fmpz, l)

  field=1
  for L in l
    S=Hecke.simple_extension(L)[1]
    K=Hecke.absolute_field(S)[1]
    K=simplify(K, canonical=true)[1]
    O=maximal_order(K)
    D=abs(discriminant(O))
    if D^9>absolute_bound
      continue
    end
    bo = ceil(Rational{BigInt}(absolute_bound//D^9))
    bound = FlintZZ(fmpq(bo))

    C,mC=class_group(O)
    allow_cache!(mC)
    cgrp=false
    if gcd(C.snf[end],3)!=1
      cgrp=true
      S = prime_ideals_up_to(O, max(100,12*clog(D,10)^2))
    end
    Aut=Hecke.automorphisms(K)
    @assert length(Aut)==degree(O) #The field is normal over Q
    gens = Hecke.small_generating_set(Aut)

    #Getting conductors
    l_conductors=conductors(O,9,bound, false)
    @vprint :QuadraticExt 1 "Conductors: $(length(l_conductors))\n"
    #Now, the big loop
    for (i, k) in enumerate(l_conductors)
      @vprint :QuadraticExt 1 "Doing $i\n"
      r,mr=ray_class_group_quo(O,9,k[1], k[2])
      if !_are_there_subs(r,[9])
        continue
      end
      if cgrp
        mr.prime_ideal_cache = S
      end
      @vprint :QuadraticExt 1 "Computing the action\n"
      act=_act_on_ray_class(mr,gens)
      @vprint :QuadraticExt 1 "Computing subgroups\n"
      ls=stable_subgroups(r,[9],act, op=(x, y) -> quo(x, y, false)[2])
      a=_min_wild(k[2])*k[1]
      for s in ls
        if _trivial_action(s,act,9)
          continue
        end
        C=ray_class_field(mr, s)
        if Hecke._is_conductor_min_normal(C,a) && Hecke.discriminant_conductor(O,C,a,mr,bound,9) && evaluate(FacElem(C.absolute_discriminant)) <= absolute_bound
          absolute_bound=evaluate(FacElem(C.absolute_discriminant))
          @vprint :QuadraticExt 1 "New Field with discriminant $absolute_bound"
          field=C
        end
      end
    end
  end
  #if typeof(field)==ClassField
  #  field=number_field(C)
  #end

  return field

end

###############################################################################
#
#  Sieving discriminants
#
###############################################################################

function _discriminant_bound(autos, ram_primes::Array{fmpz,1}, bound::fmpz)

  K=_to_non_normal(autos)
  @vprint :QuadraticExt 1 "Doing $(K)"
  #now, compute the discriminant of K. Since we know the ramified primes,
  #we know the primes dividing the discriminant and this is easier than computing the maximal order.
  return _is_discriminant_lower(K,ram_primes,bound)

end

function _is_discriminant_lower(K::AnticNumberField, ram_primes::Array{fmpz,1}, bound::fmpz)

  disc=fmpz(1)
  O=EquationOrder(K)
  for p in ram_primes
    OO=pmaximal_overorder(O,p)
    disc*=p^(valuation(discriminant(OO),p))
    if disc>bound
      return false
    end
  end
  return true

end

###############################################################################
#
#  From normal extension to the non-normal one using the trace/norm function
#
###############################################################################

function _to_non_normal(autos)#::Vector{NfRel_nsToNfRel_nsMor})

  #first, find a generator of the simple extension
  L=domain(autos[1])
  S,mS=simple_extension(L)
  x=mS(gen(S))

  #find an element of order 2 in the automorphisms
  Id=find_identity(autos,*)
  i=1
  while autos[i]==Id || autos[i]*autos[i]!=Id
    i+=1
  end

  #Find primitive element
  pr_el=x+autos[i](x)

  #Take minimal polynomial; I need to embed the element in the absolute extension
  pol=absolute_minpoly(pr_el)
  if degree(pol)==15
    return number_field(pol)[1]
  else
    pr_el=x*(autos[i](x))
    return number_field(absolute_minpoly(pr_el))[1]
  end

end

################################################################################
#
#  Valuation bounds for discriminants
#
################################################################################

function valuation_bound_discriminant(n::Int, p::Union{Integer, fmpz})
  # First compute the p-adic expansion of n
  S = Vector{typeof(p)}()
	q = typeof(p)(n)
  q, r = divrem(q, p)
  push!(S, r)
  while q >= p
    q, r = divrem(q, p)
    push!(S, r)
  end

	if !iszero(q)
		push!(S, q)
	end

	@assert sum(S[i + 1] * p^i for i in 0:length(S)-1) == n

	b = zero(typeof(p))

	for i in 0:length(S) - 1
		b = b + S[i + 1] * (i + 1) * p^i
		if !iszero(S[i + 1])
			b = b - 1
		end
	end

  return b
end

###############################################################################
#
#  From relative to absolute
#
###############################################################################

function _from_relative_to_abs(L::Tuple{NfRel_ns, Array{NfRel_nsToNfRel_nsMor{T},1}}) where T

  S,mS=simple_extension(L[1])
  K,mK=absolute_field(S, false)

  #First, we compute the maximal order of the absolute field.
  #We start from the maximal orders of the relative extension and of the base field.
  #FALSE: Since the computation of the relative maximal order is slow, I prefer to bring to the absolute field the elements
  # generating the equation order.
  O=maximal_order(L[1].base_ring)
  OL=EquationOrder(L[1])
  B=pseudo_basis(OL)

  @vprint :QuadraticExt 1 "Maximal Orders computed\n"
  #Then we consider the product basis
  basisL=Array{NfRel_nsElem, 1}(2*degree(L[1]))
  for i=1:degree(L[1])
    _assure_weakly_normal_presentation(B[i][2].num)
    basisL[2*i-1]=divexact(B[i][2].num.gen_one* B[i][1], B[i][2].den)
    basisL[2*i]=divexact(L[1].base_ring(B[i][2].num.gen_two)* B[i][1], B[i][2].den)
  end
  basisO=[L[1](O.basis_nf[i]) for i=1:degree(O)]

  nbasisL=[mK(mS\(el)) for el in basisL]
  nbasisO=[mK(mS\(el)) for el in basisO]

  cbasis=[x*y for x in nbasisL for y in nbasisO]
  append!(cbasis, [gen(K)^i for i=0:degree(K)-1])
  @vprint :QuadraticExt 1 "Product basis computed\n"
  #Now, we compute the maximal order. Then we simplify.
  O1=MaximalOrder(_order(K, cbasis))
  O1.ismaximal=1
  _set_maximal_order_of_nf(K,O1)
  Ks, mKs= simplify(K)

  #Now, we have to construct the maximal order of this field.
  #I am computing the preimages of mKs by hand, by inverting the matrix.
  M = zero_matrix(FlintZZ, degree(Ks), degree(Ks))
  prim_img=mKs(gen(Ks))
  M1=inv(basis_mat([prim_img^i for i=0:degree(Ks)-1]))
  basisO2=Array{nf_elem, 1}(degree(Ks))
  M=zero_matrix(FlintZZ, 1, degree(Ks))
  for i=1:length(basisO2)
    elem_to_mat_row!(M, 1, denominator(O1.basis_nf[i]), O1.basis_nf[i])
    mul!(M, M, M1.num)
    basisO2[i]=elem_from_mat_row(Ks, M, 1, M1.den*denominator(O1.basis_nf[i]))
  end
  O2=Order(Ks, basisO2, false)
  O2.ismaximal=1
  _set_maximal_order_of_nf(Ks,O2)
  @vprint :QuadraticExt 1 "MaximalOrder Computed. Now Automorphisms\n"

  #Now, the automorphisms.
  autos=Array{NfToNfMor, 1}(length(L[2]))
  el=mS(mK\(mKs(gen(Ks))))
  for i=1:length(L[2])
    x=mK(mS\(L[2][i](el)))
    elem_to_mat_row!(M, 1, denominator(x), x)
    mul!(M, M, M1.num)
    y=elem_from_mat_row(Ks, M, 1, M1.den*denominator(x))
    @assert Ks.pol(y)==0
    autos[i]=NfToNfMor(Ks,Ks,y)
  end

  @vprint :QuadraticExt 1 "Finished\n"
  return Ks, autos

end

###############################################################################
#
#  Read-Write
#
###############################################################################

function _write_fields(list::Array{Tuple{AnticNumberField, fmpz},1}, filename::String)
  f=open(filename, "a")
  for L in list
    x=([coeff(L[1].pol, i) for i=0:degree(L[1].pol)], L[2])
    Base.write(f, "$x\n")
  end
  close(f)
end

function _read_fields(filename::String)
  f=open(filename, "r")
  Qx,x=PolynomialRing(FlintQQ,"x")
  pols=Tuple{fmpq_poly, fmpz}[]
  for s in eachline(f)
    a=eval(parse(s))
	  push!(pols,(Qx(a[1]), a[2]))
	end
	close(f)
	return pols
end

ext: jl


text:
function gain(u,r,K,ind)
    return r[ind] .* u[ind] .* (K[ind] .- u[ind])
end

function holling2(u,r,kappa,A,ind)
    return r[ind] .* u[ind] .* (A' * u)[ind] ./ (kappa[ind] .+ (A' * u)[ind])
end

function dependentLoss(u,r,kappa,A,ind)
    return u[ind] .* sum((A[ind,:]' .* u) ./ (kappa .+ sum(A .* u, dims = 1)[:]), dims = 1)[:]
end

function webDynamic!(du,u,p,t)
    r,K,A,kappa,m,ind = p
    du[ind[1]] = gain(u,r,K,ind[1]) .- m[ind[1]] .* dependentLoss(u,r,kappa,A,ind[1])
    du[ind[2]] = holling2(u,r,kappa,A,ind[2]) .- m[ind[2]] .* dependentLoss(u,r,kappa,A,ind[2]) .- 0.01 .* u[ind[2]]
    du[ind[3]] = holling2(u,r,kappa,A,ind[3]) .- m[ind[3]] .* u[ind[3]]
end


text:
# le problème de Josephus

# listes circulaires doublement chaînées

mutable struct L
    valeur
    suivant
    precedent
end

# liste réduite à un élément
function make(v)
    r = L(v, nothing, nothing)
    r.suivant = r
    r.precedent = r
    r
end

# insertion après un élément donnée
function inserer_apres(l, v)
    e = make(v)
    e.suivant = l.suivant
    l.suivant = e
    e.suivant.precedent = e
    e.precedent = l
end

# suppression d'un élément donné
function supprimer(l)
    l.precedent.suivant = l.suivant
    l.suivant.precedent = l.precedent
end

# affichage
function afficher(l)
    p = l
    print(p.valeur)
    p = p.suivant
    while p != l
	print(p.valeur)
	p = p.suivant
    end
    print("\n")
end

# problème de Josephus

# construction de la liste circulaire 1,2,...,n;
# l'élément renvoyé est celui contenant 1
function cercle(n)
    l = make(1)
    i = n
    while i >= 2
	inserer_apres(l, i)
        i = i - 1
    end
    l
end

function josephus(n, p)
    # c est le joueur courant, 1 au départ
    c = cercle(n)

    # tant qu'il reste plus d'un joueur
    while c != c.suivant
	# on élimine un joueur
	for i = 1 : p-1
	    c = c.suivant
	end
	supprimer(c)
	c = c.suivant
    end
    c.valeur
end

println(josephus(7, 5)) # 6
println(josephus(5, 5)) # 2
println(josephus(5, 17)) # 4
println(josephus(13, 2)) # 11

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
abstract type AbelianExt end

mutable struct KummerExt <: AbelianExt
  zeta::nf_elem
  n::Int
  gen::Vector{FacElem{nf_elem, AnticNumberField}}

  AutG::GrpAbFinGen
  frob_cache::Dict{NfOrdIdl, GrpAbFinGenElem}
  frob_gens::Tuple{Vector{NfOrdIdl}, Vector{GrpAbFinGenElem}}
  gen_mod_nth_power::Vector{FacElem{nf_elem, AnticNumberField}}
  eval_mod_nth::Vector{nf_elem}

  function KummerExt()
    return new()
  end
end

function Base.show(io::IO, K::KummerExt)
  if isdefined(K.AutG, :snf)
    print(io, "KummerExt with structure $(K.AutG.snf)")
  else
    print(io, "KummerExt with structure $([K.AutG.rels[i, i] for i=1:ngens(K.AutG)])")
  end
end

function kummer_extension(n::Int, gen::Vector{FacElem{nf_elem, AnticNumberField}})
  K = KummerExt()
  k = base_ring(gen[1])
  L = maximal_order(k)
  zeta, o = torsion_units_gen_order(L)
  @assert o % n == 0

  K.zeta = k(zeta)^div(o, n)
  K.n = n
  K.gen = gen
  K.AutG = GrpAbFinGen(fmpz[n for i=gen])
  K.frob_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
  return K
end

function kummer_extension(exps::Array{Int, 1}, gens::Vector{FacElem{nf_elem, AnticNumberField}})
  K = KummerExt()
  k = base_ring(gens[1])
  L = maximal_order(k)
  zeta, o = torsion_units_gen_order(L)
  n = lcm(exps)
  @assert o % n == 0

  K.zeta = k(zeta)^div(o, n)
  K.n = n
  K.gen = gens
  K.AutG = DiagonalGroup(exps)
  K.frob_cache = Dict{NfOrdIdl, GrpAbFinGenElem}()
  return K
end

function kummer_extension(n::Int, gen::Array{nf_elem, 1})
  g = [FacElem(x) for x=gen]
  return kummer_extension(n, g)
end

###############################################################################
#
#  Base Field
#
###############################################################################

function base_field(K::KummerExt)
  return base_ring(K.gen[1])
end

###############################################################################
#
#  Degree
#
###############################################################################

function degree(K::KummerExt)
  return Int(order(K.AutG))
end

###############################################################################
#
#  From Kummer Extension to Number Field
#
###############################################################################

function number_field(K::KummerExt)
  k = base_field(K)
  kt, t = PolynomialRing(k, "t", cached = false)
  pols = Array{typeof(t), 1}(undef, length(K.gen))
  for i = 1:length(pols)
    p = kt()
    setcoeff!(p, Int(order(K.AutG[i])), k(1))
    setcoeff!(p, 0, -evaluate(K.gen[i]))
    pols[i] = p
  end
  return number_field(pols, check = false, cached = false)
end

###############################################################################
#
#  Computation of Frobenius automorphisms
#
###############################################################################

mutable struct NfToFqMor_easy <: Map{AnticNumberField, FqNmodFiniteField, HeckeMap, NfToFqMor_easy}
  header::MapHeader
  Fq::FqNmodFiniteField
  s::fq_nmod
  t::gfp_poly
  function NfToFqMor_easy(a::Map, k::AnticNumberField)
    r = new()
    r.Fq = codomain(a)
    r.header = MapHeader(k, r.Fq)
    r.s = r.Fq()
    r.t = PolynomialRing(GF(UInt(characteristic(r.Fq)), cached=false), cached=false)[1]()
    return r
  end
end

function image(mF::NfToFqMor_easy, a::FacElem{nf_elem, AnticNumberField}, quo::Int = 0)
  Fq = mF.Fq
  q = one(Fq)
  t = mF.t
  s = mF.s
  for (k, v) = a.fac
    vv = v
    if quo != 0
      vv = v %quo
      if vv < 0
        vv += quo
      end
    end
    @assert vv < order(Fq)  #please complain if this is triggered
    if !iszero(vv)
      if denominator(k) % characteristic(Fq) == 0
        throw(BadPrime(characteristic(Fq)))
      end
      _nf_to_fq!(s, k, Fq, t)
      if iszero(s)
        throw(BadPrime(1))
      end
      if vv < 0
        ccall((:fq_nmod_inv, :libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), s, s, Fq)
        vv = -vv
      end
      ccall((:fq_nmod_pow_ui, :libflint), Nothing, (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}), s, s, vv, Fq)
      mul!(q, q, s)
    end
  end
  return q
end

function image(mF::NfToFqMor_easy, a::nf_elem, n_quo::Int = 0)
  Fq = mF.Fq
  q = Fq()
  if denominator(a) % characteristic(Fq) == 0
    throw(BadPrime(characteristic(Fq)))
  end
  _nf_to_fq!(q, a, Fq, mF.t)
  return q
end

# the Frobenius at p in K:
#K is an extension of k, p a prime in k,
#returns a vector in (Z/nZ)^r representing the Frob
function can_frobenius(p::NfOrdIdl, K::KummerExt)
  @assert norm(p) % K.n == 1
  if haskey(K.frob_cache, p)
    return K.frob_cache[p]
  end
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end

  if nbits(minimum(p)) > 64
    error("Oops")
  end

  F, mF = ResidueFieldSmall(Zk, p)::Tuple{FqNmodFiniteField,NfOrdToFqNmodMor}
  #_mF = extend_easy(mF, number_field(Zk))
  mF = NfToFqMor_easy(mF, number_field(Zk))
  z_p = image(mF, K.zeta)^(K.n-1)

  # K = k(sqrt[n_i](gen[i]) for i=1:length(gen)), an automorphism will be
  # K[i] -> zeta^divexact(n, n_i) * ? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  aut = Array{fmpz, 1}(undef, length(K.gen))
  for j = 1:length(K.gen)
    ord_genj = Int(order(K.AutG[j]))
    ex = div(norm(p)-1, ord_genj)
    if isdefined(K, :gen_mod_nth_power)
      mu = image(mF, K.gen_mod_nth_power[j])^ex
    else
      mu = image(mF, K.gen[j], K.n)^ex  # can throw bad prime!
    end
    i = 0
    z_pj = z_p^divexact(K.n, ord_genj)
    while !isone(mu)
      i += 1
      @assert i <= K.n
      mul!(mu, mu, z_pj)
    end
    aut[j] = fmpz(i)
  end
  z = K.AutG(aut)
  K.frob_cache[p] = z
  return z
end

function can_frobenius1(p::NfOrdIdl, K::KummerExt)
  @assert norm(p) % K.n == 1
  if haskey(K.frob_cache, p)
    return K.frob_cache[p]
  end
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end

  if nbits(minimum(p)) > 64
    error("Oops")
  end

  F, mF = ResidueFieldSmall(Zk, p)
  mF = NfToFqMor_easy(mF, number_field(Zk))
  #mF = extend_easy(mF, number_field(Zk))
  z_p = mF(K.zeta)^(K.n-1)

  # K = k(sqrt[n_i](gen[i]) for i=1:length(gen)), an automorphism will be
  # K[i] -> zeta^divexact(n, n_i) * ? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p
  ex = div(norm(p)-1, K.n)
  gens = Vector{fq_nmod}(undef, length(K.gen))
  for i = 1:length(K.gen)
    if isdefined(K, :eval_mod_nth)
      gens[i] = image(mF, K.eval_mod_nth[i])^ex
      if iszero(gens[i])
        throw(BadPrime(p))
      end
    elseif isdefined(K, :gen_mod_nth_power)
      gens[i] = image(mF, K.gen_mod_nth_power[i])^ex
    else
      gens[i] = image(mF, K.gen[i], K.n)^ex  # can throw bad prime!
    end
  end

  aut = Array{fmpz, 1}(undef, length(K.gen))
  for j = 1:length(K.gen)
    mu = gens[j]
    i = 0
    while !isone(mu)
      i += 1
      @assert i <= K.n
      mul!(mu, mu, z_p)
    end
    aut[j] = fmpz(i)
  end
  z = K.AutG(aut)
  K.frob_cache[p] = z
  return z
end

#In this function, we are computing the image of $sqrt[n](g) under the Frobenius automorphism of p
function can_frobenius(p::NfOrdIdl, K::KummerExt, g::FacElem{nf_elem})
  Zk = order(p)
  if index(Zk) % minimum(p) == 0
    #index divisors and residue class fields don't agree
    # ex: x^2-10, rcf of 29*Zk, 7. 239 is tricky...
    throw(BadPrime(p))
  end

  if nbits(minimum(p)) > 64
    error("Oops")
  end

  F, mF = ResidueFieldSmall(Zk, p)
  mF = extend_easy(mF, nf(Zk))

  #K = sqrt[n](gen), an automorphism will be
  # K[i] -> zeta^? K[i]
  # Frob(sqrt[n](a), p) = sqrt[n](a)^N(p) (mod p) = zeta^r sqrt[n](a)
  # sqrt[n](a)^N(p) = a^(N(p)-1 / n) = zeta^r mod p

  z_p = inv(mF(K.zeta))
  @assert norm(p) % K.n == 1
  ex = div(norm(p)-1, K.n)

  mu = image(mF, g, K.n)^ex  # can throw bad prime!
  i = 0
  while true
    if isone(mu)
      break
    end
    i += 1
    @assert i <= K.n
    mul!(mu, mu, z_p)
  end
  return i
end


text:
abstract type AbstractPenalization end
abstract type AbstractPenalizer end
# Penalization API #
# ------------------------------------------------------------
penalty(ph::AbstractProgressiveHedging) = penalty(ph, ph.penalization)
initialize_penalty!(ph::AbstractProgressiveHedging) = initialize_penalty!(ph, ph.penalization)
update_penalty!(ph::AbstractProgressiveHedging) = update_penalty!(ph, ph.penalization)
# ------------------------------------------------------------
# Attributes #
# ------------------------------------------------------------
"""
    RawPenalizationParameter

An optimizer attribute used for raw parameters of the penalizer. Defers to `RawParameter`.
"""
struct RawPenalizationParameter <: PenalizationParameter
    name::Any
end
include("common.jl")
include("fixed.jl")
include("adaptive.jl")


text:
function _println_header(io::IO, B::LogBinner{N,T}) where {N, T}
    println("LogBinner{$(N),$(T)}")
end

function _println_body(io::IO, B::LogBinner{N,T}) where {N, T}
    n = B.count[1]
    print("| Count: ", n)
    if n > 0
        print("\n| Mean: ", round.(mean(B), digits=5))
        print("\n| StdError: ", round.(std_error(B), digits=5))
    end
end

Base.show(io::IO, B::LogBinner{N,T}) where {N, T} = begin
    _println_header(io, B)
    _println_body(io, B)
    nothing
end
Base.show(io::IO, m::MIME"text/plain", B::LogBinner{N,T}) where {N, T} = print(io, B)

Base.summary(io::IO, B::LogBinner{N,T}) where {N, T} = _println_header(io, B)
Base.summary(B::LogBinner{N,T}) where {N, T} = summary(stdout, B)


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

module RelationalDatasets

include("request.jl")
export load

include("convert.jl")

end # module


text:
module BLG_SQUID
    # Use README as the docstring of the module:
    @doc read(joinpath(dirname(@__DIR__), "README.md"), String) BLG_SQUID

    using Requires

    using CairoMakie, VegaLite
    using SharedArrays, Distributed
    using Quantica, StaticArrays, Parameters, LinearAlgebra, StatsBase
    using Baselet, Arpack, CSV, DataFrames, Dates
    using Optim, ProgressMeter, Interpolations, Random, LsqFit

    using PhysicalConstants.CODATA2018: ustrip, @u_str, h, ħ, k_B, m_e, e, μ_B

    export Params
    export nanoribbonS, nanoribbonSA, nanoribbonSZ, Params,  modelS, rectangle_weaklink,
        rectangle_randombounds_sc, ldosonlattice_averaged_sc, ldosonlattice
    export icϕ_exactdiag, fraunhofer_abs_exact, maxicϕ_exactdiag, fraunhofer_abs_exact_adaptive

    export Params_eff_quad, supercurrent_ef, supercurrent2_1dchains, spectrumvsphase_ef,
        fraunhofer_eff, plotfraunhof_vs_phase

    include("model.jl")
    include("nanoribbon.jl")
    include("bounded_sys.jl")
    include("ldos.jl")
    include("spectrum.jl")
    include("weaklink.jl")
    include("save.jl")
    include("adaptive_maxfinder.jl")
    include("transport.jl")
    include("effective_model.jl")

end


text:
# # Problems and the system

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
In specific problems that we wish to solve with immersed layers, there may
be other data and operators that we would like to cache. We do this with
an *extra cache*, which the user can define, along with a problem type associated
with this cache. The basic cache and the extra cache are generated and associated
together in a *system*.

There are a few basic ingredients to do this:
* **Create a problem type**, using the macro `@ilmproblem`. This mostly just serves as a means of dispatching correctly, but also holds the grid and bodies so they can be passed along when the caches are constructed.
* **Create an extra cache type**, making it a subtype of `AbstractExtraILMCache`. This can hold pretty much anything you want it to.
* **Extend the function `prob_cache(prob,base_cache)`** to serve as a constructor for your extra cache, when your problem type is passed in.

Optionally, you can also extend the function `solve` in order to perform the steps
of your algorithm. However, generically, you can just use pass in the system structure,
which holds the basic ILM cache and your extra cache, into any function.
=#

#=
## Example of problem and system use
We will demonstrate the use of problems and systems with the example
given in [A Dirichlet Poisson problem](@ref). Here, we will
assemble the various additional data structures and operators used to
solve this problem into an extra cache. We will also create a problem
type called `DirichletPoissonProblem`, which we make a subtype of
`AbstractScalarILMProblem`.
=#

using ImmersedLayers
#!jl using Plots
using UnPack

#=
### Create your problem type
For this job, we have a macro
=#
@ilmproblem DirichletPoisson scalar

#=
This generates a type `DirichletPoissonProblem` for dispatch, and an associated
constructor that we'll use later.
=#

#=
### Create your extra cache
Here, we'd like this extra cache to hold the Schur complement and the
filtering matrices, as well as some cache variables.
=#
struct DirichletPoissonCache{SMT,CMT,FRT,ST,FT} <: ImmersedLayers.AbstractExtraILMCache
   S :: SMT
   C :: CMT
   forcing_cache :: FRT
   fb :: ST
   fstar :: FT
end

#=
### Extend the `prob_cache` function
We need this to construct our extra cache. This will get called when the system
is constructed. We can make use of the data available
in the `base_cache` (which is already constructed before this is called)
in order to create the extra cache. (In the example here, we aren't quite ready to
use the `forcing_cache`, so we just keep it blank for now. We will use it later.)
=#
function ImmersedLayers.prob_cache(prob::DirichletPoissonProblem,base_cache::BasicILMCache)
    S = create_RTLinvR(base_cache)
    C = create_surface_filter(base_cache)
    forcing_cache = nothing
    fb = zeros_surface(base_cache)
    fstar = zeros_grid(base_cache)
    DirichletPoissonCache(S,C,forcing_cache,fb,fstar)
end

#=
### Extend the `solve` function
Here, we actually do the work of the algorithm, making use of all of the
operators and data structures that we have cached for efficiency.
The example below takes in some surface Dirichlet data `fbplus`,
and returns the solutions `f` and `s` (filtered).
=#
function ImmersedLayers.solve(fbplus,prob::DirichletPoissonProblem,sys::ILMSystem)
    @unpack extra_cache, base_cache = sys
    @unpack S, C, fb, fstar = extra_cache

    f = zeros_grid(base_cache)
    s = zeros_surface(base_cache)

    surface_divergence!(fstar,fbplus,base_cache)
    fb .= 0.5*fbplus

    inverse_laplacian!(fstar,base_cache)

    interpolate!(s,fstar,base_cache)
    s .= fb - s
    s .= -(S\s);

    regularize!(f,s,base_cache)
    inverse_laplacian!(f,base_cache)
    f .+= fstar;

    s .= C^6*s

    return f, s
end

#=
### Set up the grid, shape, and cache
We do this just as we did in [Immersed layer caches](@ref), but
now we don't create a cache, since it will be done internally.
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)
RadC = 1.0
Δs = 1.4*cellsize(g)
body = Circle(RadC,Δs)

#=
### Do the work
We do this in three steps:
- Create the problem instance. This is where we use the constructor generated by `@ilmproblem`
- Call `construct_system` to create the caches, assembled into a system
- Call `solve` to solve the problem.
Also, note that pretty much any function that accepts `base_cache`
as an argument also accepts `sys`.
=#
prob = DirichletPoissonProblem(g,body,scaling=GridScaling)
sys = construct_system(prob)

pts = points(sys)
f, s = solve(pts.u,prob,sys)

#!jl plot(f,sys)
#-
#!jl plot(s,sys)

#=
## More advanced use, with keyword arguments
In the previous example, we supplied the surface data as an argument directly to
the `solve` function. We also didn't have any right-hand side forcing. But in
a more general case we might have more complicated boundary conditions or forcing
information. This will be particularly true when we get to time varying problems.

So in this next example, we will demonstrate how the problem constructor can
be used to supply more sophisticated information to the `solve` function.
In particular, we will use three keywords:
* The `bc` keyword to supply boundary condition information.
* The `forcing` keyword to supply forcing.
* The `phys_params` keyword to supply physical parameters.

In general, the philosophy is that we wish to enable the user to provide a function that
supplies the required boundary or forcing data. There is a lot
of freedom in this approach, as we will demonstrate.
=#
#=
### Boundary conditions
Here, we wish to construct functions that return the value of `f`
on the outside and inside of the surface. These functions will
get called in the `solve` function, and they will be passed
the base cache, the physical parameters. So we can make use
of both of these to return the boundary data, e.g., here
we return the `x` coordinate of the surface points on the outside
of the surface, and zeros on the inside.
=#
get_fbplus(base_cache,phys_params) = points(base_cache).u
get_fbminus(base_cache,phys_params) = zeros_surface(base_cache)

#=
Using the `bc` keyword, we can pass these along to the problem specification
by various means. To make use of some routines that automate the
application of these functions, we must pass them in a `Dict`, with the
keys `exterior` and `interior`, respectively.
=#
bcdict = Dict("exterior"=>get_fbplus,"interior"=>get_fbminus)

#=
### Forcing
Forcing can come in several forms. It can be
* **area forcing**, in which the forcing is distributed over a whole region, possibly confined
* **line forcing**, in which the forcing is distributed along a curve (but singular in the orthogonal direction)
* **point forcing**, in which the forcing is singular at one or more points

In all cases, the user needs to supply some information about the **shape** of the
forcing region and some information about the instantaneous strength of the forcing,
in the form of a **model function**. These get bundled together into a **forcing model**.

For confined area forcing and for line forcing, the **shape** is supplied as a `Body` or `BodyList`,
similar to bodies (serving as the mask for confined-area forcing and as the curve
for line forcing). For point forcing, the **shape** is the set of discrete points
at which we want to apply forcing. We can either supply these points' coordinates directly
or provide them via a **position function** that returns the point coordinates. The advantage
of the latter is that it allows the points to be changed later without
regenerating the system.

In the example shown here, we apply point forcing at two points: one at (-1.0,1.5)
and the other at (1.0,1.5), with respective strengths -1 and 1.
We use a **position function** to provide the positions.
This function must have the signature as shown here, accepting
a general `state` argument, a time argument `t`, the cache for this forcing,
and physical parameters. Any of these can be used to help determine
the point positions. In this case, we just simply provide the coordinates
directly. However, it is straightforward to pass these in via the `state`
argument (or perhaps to update the positions by some velocity).

=#
function my_point_positions(state,t,fr::PointRegionCache,phys_params)
    x = [-1.0,1.0]
    y = [1.5,1.5]
    return x, y
end

#=
The **model function** returns the strength of the forcing. In this point
forcing example, it must return the strength of each point. It must be
set up as an in-place function, accepting the vector of strengths and
mutating (i.e. changing) this vector. Beyond that, it has the same
signature as the position function.
=#
function my_point_strengths!(σ,state,t,fr::PointRegionCache,phys_params)
    σ .= [-1.0,1.0]
end

#=
Now we bundle these together into the forcing model. We also tell
it what kind of regularized DDF we want to use. We use the $M_4'$ DDF here.
`pfm` will get passed into the problem setup via the `forcing` keyword.
=#
pfm = PointForcingModel(my_point_positions,my_point_strengths!;ddftype=CartesianGrids.M4prime);

#=
Now, we need to make sure that the forcing data is used. When
we generate the extra cache, we use the function `ForcingModelAndRegion`
to assemble the forcing's cache.
=#
function ImmersedLayers.prob_cache(prob::DirichletPoissonProblem,base_cache::BasicILMCache)
    @unpack phys_params, forcing = prob
    S = create_RTLinvR(base_cache)
    C = create_surface_filter(base_cache)
    forcing_cache = ForcingModelAndRegion(forcing,base_cache)
    fb = zeros_surface(base_cache)
    fstar = zeros_grid(base_cache)
    DirichletPoissonCache(S,C,forcing_cache,fb,fstar)
end

#=
We now redefine the `solve` function and use the boundary condition and forcing
information. For the forcing, we make use of the function `apply_forcing!`
which takes the forcing and automatically calculates the right-hand side field
=#
function ImmersedLayers.solve(prob::DirichletPoissonProblem,sys::ILMSystem)
    @unpack bc, forcing, phys_params, extra_cache, base_cache = sys
    @unpack gdata_cache = base_cache
    @unpack S, C, forcing_cache, fb, fstar = extra_cache

    f = zeros_grid(base_cache)
    s = zeros_surface(base_cache)

    ## apply_forcing! evaluates the forcing field on the grid and put
    ## the result in the `gdata_cache`.
    fill!(gdata_cache,0.0)
    apply_forcing!(gdata_cache,f,0.0,forcing_cache,phys_params)

    ## Get the prescribed jump in boundary data across the interface using
    ## the functions we supplied via the `Dict`.
    prescribed_surface_jump!(fb,sys)

    ## Evaluate the double-layer term and add it to the right-hand side
    surface_divergence!(fstar,fb,base_cache)
    fstar .+= gdata_cache

    ## Intermediate solution
    inverse_laplacian!(fstar,base_cache)

    ## Get the prescribed average of boundary data on the interface using
    ## the functions we supplied via the `Dict`.
    prescribed_surface_average!(fb,sys)

    ## Correction
    interpolate!(s,fstar,base_cache)
    s .= fb - s
    s .= -(S\s);

    regularize!(f,s,base_cache)
    inverse_laplacian!(f,base_cache)
    f .+= fstar;

    return f, C^6*s
end

#=
### Set up the problem
Now we specify the problem, create the system, and solve it, as before,
but now supplying the boundary condition and forcing information with the `bc`
and `forcing` keywords:
=#
prob = DirichletPoissonProblem(g,body,scaling=GridScaling,bc=bcdict,forcing=pfm)
sys = construct_system(prob)
f, s = solve(prob,sys)
#!jl plot(f,sys)

#=
So we get the additional features from the sources. Now, suppose we wish to change the
the boundary conditions or source points? We can do it easily without regenerating the
cache and system, simply by redefining our bc and forcing model functions. For
example, to create an internal solution, with surface data equal to the $y$ coordinate,
and four sources inside.
=#
get_fbplus(base_cache,phys_params) = zeros_surface(base_cache)
get_fbminus(base_cache,phys_params) = points(base_cache).v

function my_point_positions(state,t,fr::PointRegionCache,phys_params)
    x = [-0.2,0.2,-0.2,0.2]
    y = [0.2,0.2,-0.2,-0.2]
    return x, y
end
function my_point_strengths!(σ,state,t,fr::PointRegionCache,phys_params)
    σ .= [-1.0,1.0,-1.0,1.0]
end

#=
We can solve immediately without having to reconstruct the system, so
it's very fast.
=#
@time f, s = solve(prob,sys)
#!jl plot(f,sys)

#md # ## Problem types and functions
#md #
#md # ```@docs
#md # AbstractScalarILMProblem
#md # AbstractVectorILMProblem
#md # @ilmproblem
#md # BasicScalarILMProblem
#md # BasicVectorILMProblem
#md # prob_cache
#md # ```

#md # ## System types and functions
#md #
#md # ```@docs
#md # ILMSystem
#md # construct_system
#md # update_system
#md # update_system!
#md # ```

#md # ## Forcing functions
#md # ```@docs
#md # AreaForcingModel
#md # LineForcingModel
#md # PointForcingModel
#md # AreaRegionCache
#md # LineRegionCache
#md # PointRegionCache
#md # ForcingModelAndRegion
#md # arcs(::LineRegionCache)
#md # mask(::AreaRegionCache)
#md # points(::PointRegionCache)
#md # apply_forcing!
#md # ```


text:
"""
    co_identity(x)
Identity function. Already defined in Julia as `identity`, specialized for scalars in the `comparison` layer.
"""
co_identity(x; param=nothing, dom_size=0, nvars=0) = identity(x)

"""
    co_abs_diff_val_param(x; param)
Return the absolute difference between `x` and `param`.
"""
co_abs_diff_val_param(x; param, dom_size=0, nvars=0) = abs(x - param)

"""
    co_val_minus_param(x; param)
Return the difference `x - param` if positive, `0.0` otherwise.
"""
co_val_minus_param(x; param, dom_size=0, nvars=0) = max(0.0, x - param)

"""
    co_param_minus_val(x; param)
Return the difference `param - x` if positive, `0.0` otherwise.
"""
co_param_minus_val(x; param, dom_size=0, nvars=0) = max(0.0, param - x)

"""
    co_euclidian_param(x; param, dom_size)
Compute an euclidian norm with domain size `dom_size`, weigthed by `param`, of a scalar.
"""
function co_euclidian_param(x; param, dom_size, nvars=0)
    return x == param ? 0.0 : (1.0 + abs(x - param) / dom_size)
end

"""
    co_euclidian(x; dom_size)
Compute an euclidian norm with domain size `dom_size` of a scalar.
"""
function co_euclidian(x; param=nothing, dom_size, nvars=0)
    return co_euclidian_param(x; param=0.0, dom_size=dom_size)
end

"""
    co_abs_diff_val_vars(x; nvars)
Return the absolute difference between `x` and the number of variables `nvars`.
"""
co_abs_diff_val_vars(x; param=nothing, dom_size=0, nvars) = abs(x - nvars)

"""
    co_val_minus_vars(x; nvars)
Return the difference `x - nvars` if positive, `0.0` otherwise, where `nvars` denotes the numbers of variables.
"""
co_val_minus_vars(x; param=nothing, dom_size=0, nvars) = co_val_minus_param(x; param=nvars)

"""
    co_vars_minus_val(x; nvars)
Return the difference `nvars - x` if positive, `0.0` otherwise, where `nvars` denotes the numbers of variables.
"""
co_vars_minus_val(x; param=nothing, dom_size=0, nvars) = co_param_minus_val(x; param=nvars)

"""
    comparison_layer(param = false)
Generate the layer of transformations functions of the ICN. Iff `param` value is set, also includes all the parametric comparison with that value. The operations are mutually exclusive, that is only one will be selected.
"""
function comparison_layer(param=false)
    comparisons = LittleDict{Symbol,Function}(
        :identity => co_identity,
        :euclidian => co_euclidian,
        :abs_diff_val_vars => co_abs_diff_val_vars,
        :val_minus_vars => co_val_minus_vars,
        :vars_minus_val => co_vars_minus_val,
    )

    if param
        comparisons_param = LittleDict{Symbol,Function}(
            :abs_diff_val_param => co_abs_diff_val_param,
            :val_minus_param => co_val_minus_param,
            :param_minus_val => co_param_minus_val,
            :euclidian_param => co_euclidian_param,
        )
        comparisons = LittleDict{Symbol,Function}(union(comparisons, comparisons_param))
    end

    return Layer(comparisons, true)
end


text:
function GLONASSAlmanac(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jdouble, arg15::TimeScale, arg2::jint, arg3::jint, arg4::jint, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return GLONASSAlmanac((jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble, TimeScale, jint, jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function GLONASSAlmanac(arg0::jint, arg1::jint, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jdouble, arg2::jint, arg3::jint, arg4::jint, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return GLONASSAlmanac((jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble, jint, jint, jint, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_date(obj::GLONASSAlmanac)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_delta_i(obj::GLONASSAlmanac)
    return jcall(obj, "getDeltaI", jdouble, ())
end

function get_delta_t(obj::GLONASSAlmanac)
    return jcall(obj, "getDeltaT", jdouble, ())
end

function get_delta_t_dot(obj::GLONASSAlmanac)
    return jcall(obj, "getDeltaTDot", jdouble, ())
end

function get_e(obj::GLONASSAlmanac)
    return jcall(obj, "getE", jdouble, ())
end

function get_frequency_channel(obj::GLONASSAlmanac)
    return jcall(obj, "getFrequencyChannel", jint, ())
end

function get_gamma_n(obj::GLONASSOrbitalElements)
    return jcall(obj, "getGammaN", jdouble, ())
end

function get_glo2utc(obj::GLONASSAlmanac)
    return jcall(obj, "getGlo2UTC", jdouble, ())
end

function get_glo_offset(obj::GLONASSAlmanac)
    return jcall(obj, "getGloOffset", jdouble, ())
end

function get_gps2_glo(obj::GLONASSAlmanac)
    return jcall(obj, "getGPS2Glo", jdouble, ())
end

function get_health(obj::GLONASSAlmanac)
    return jcall(obj, "getHealth", jint, ())
end

function get_iod(obj::GLONASSOrbitalElements)
    return jcall(obj, "getIOD", jint, ())
end

function get_lambda(obj::GLONASSAlmanac)
    return jcall(obj, "getLambda", jdouble, ())
end

function get_n4(obj::GLONASSAlmanac)
    return jcall(obj, "getN4", jint, ())
end

function get_na(obj::GLONASSAlmanac)
    return jcall(obj, "getNa", jint, ())
end

function get_pa(obj::GLONASSAlmanac)
    return jcall(obj, "getPa", jdouble, ())
end

function get_time(obj::GLONASSAlmanac)
    return jcall(obj, "getTime", jdouble, ())
end

function get_tn(obj::GLONASSOrbitalElements)
    return jcall(obj, "getTN", jdouble, ())
end

function get_x(obj::GLONASSOrbitalElements)
    return jcall(obj, "getX", jdouble, ())
end

function get_x_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getXDot", jdouble, ())
end

function get_x_dot_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getXDotDot", jdouble, ())
end

function get_y(obj::GLONASSOrbitalElements)
    return jcall(obj, "getY", jdouble, ())
end

function get_y_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getYDot", jdouble, ())
end

function get_y_dot_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getYDotDot", jdouble, ())
end

function get_z(obj::GLONASSOrbitalElements)
    return jcall(obj, "getZ", jdouble, ())
end

function get_z_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getZDot", jdouble, ())
end

function get_z_dot_dot(obj::GLONASSOrbitalElements)
    return jcall(obj, "getZDotDot", jdouble, ())
end


text:
using Statistics: std #We have some problem in the calculation of the std

function _discrete_fourier_transform(series, word_length, norm, inverse_sqrt_win_size, normalise)
    _length = length(series)
    output_length = ceil(Int, word_length/2)
    start = norm ? 1 : 0

    _std = 1

    if normalise
        s = std(series)  # some issue with std statement ??
        if s!=0
            _std = s
        end
    end

    dft = [ sum( [ [  series[n] * cos(2 * π * (n-1) * i / _length),
                     -series[n] * sin(2 * π * (n-1) * i / _length) ]
                                                   for n=1:_length ] ) # sum end
                                                   for i = start:(start + output_length-1) ] #end of dft

    dft = vcat(dft...)

    if normalise
        dft = dft * (inverse_sqrt_win_size/ _std)
    end
    return dft
end

# self._create_word(), self._add_to_pyramid(),  self.levels,  self._add_to_bag(), self.bigrams,
# self.save_words, self.words.append, self.window_size, self.word_length,

function transform(self, X, y=None)
    # check_is_fitted()
    # X = check_X(X, enforce_univariate=True)
    # X = tabularize(X, return_array=True)
    # bags = pd.DataFrame()
    # dim = []
    for i in range(X.shape[0])
        dfts = _mft(X[i, :])
        bag = Dict() # {} This is how they use the ditionary in the python
        last_word = -1
        repeat_words = 0
        words = []
        for window in range(dfts.shape[0])  # use size(dfts, 1)
            word = _create_word(dfts[window])
            append!(words, word)
            # repeat_word = (_add_to_pyramid(bag, word, last_word,      # adding aproproate if eles statement
            #                                     window -
            #                                     Int(repeat_words/2))
            #                if levels > 1 else
            #                _add_to_bag(bag, word, last_word))
            if repeat_word
                repeat_words += 1
            else
                last_word = word.word
                repeat_words = 0
            end
            if bigrams
                if ((window - window_size) >= 0) && (window > 0)
                    # bigram = words[window - window_size]\  # Adding the biagram finding its library
                    #    .create_bigram(word, word_length)
                    if levels > 1
                        bigram = (bigram, 0)
                    end
                    bag[bigram] = bag.get(bigram, 0) + 1
                end
            end
        if save_words
            words.append(words)
        end
    end
        dim.append(pd.Series(bag))
    bags[0] = dim
    end
end

function _mft(series, word_length, window_size, norm, inverse_sqrt_win_size, series_length )
    start_offset = norm ? 2 : 0
    lengthl = (word_length + word_length) ÷ 2

    phis = [  [cos(2 * π * (-((i * 2) + start_offset) / 2) / window_size),
               -sin(2 * π * (-((i * 2) + start_offset) / 2) / window_size) ]
                                          for i = 0:((lengthl÷2)-1) ] #end of dft
    phis = vcat(phis...)
    endl = max(1, series_length - window_size + 1)
    stds = _calc_incremental_mean_std(series, endl, window_size)
    transformed = zeros(Float64, endl, lengthl)
    mft_data = []  # Find way to do catching
    @inbounds for i = 0:endl-1
        if i > 0
            @inbounds for n = 1:2:lengthl
                real = mft_data[n] + series[i + window_size] - series[i]
                imag = mft_data[n + 1]
                mft_data[n] = real * phis[n] - imag * phis[n + 1]
                mft_data[n + 1] = real * phis[n + 1] + phis[n] * imag
            end
        else
            mft_data = _discrete_fourier_transform(series[1:window_size], word_length,
                                               norm, inverse_sqrt_win_size, false) #normalise=false
        end
        normalising_factor = ( (1 / ( stds[i + 1] > 0 ? stds[i + 1] : 1) ) * inverse_sqrt_win_size)
        transformed[i + 1, :] = mft_data * normalising_factor
    end
    return transformed
end

function _calc_incremental_mean_std(series, endl, window_size)
    means = zeros(endl)
    stds = zeros(endl)
    window = series[1:window_size]
    series_sum = sum(window)
    square_sum = sum(window.*window)

    r_window_length = 1 / window_size
    means[1] = series_sum * r_window_length
    buf = square_sum * r_window_length - means[1]^2
    stds[1] = buf > 0 ? sqrt(buf) : 0

    @inbounds for w =1:(endl-1)
        series_sum += series[w + window_size] - series[w]
        means[w + 1] = series_sum * r_window_length
        square_sum += series[w + window_size]^2 - series[w]^2
        buf = square_sum * r_window_length - means[w + 1]^2
        stds[w + 1] = buf > 0 ? sqrt(buf) : 0
    end
    return stds
end
# self.series_length  self.window_size, self.n_instances, self.alphabet_size,  self.word_length

function _mcb(X, window_size, alphabet_size, word_length, norm, inverse_sqrt_win_size, normalise)
    n_instances, series_length = size(X)
    num_windows_per_inst = ceil(series_length / window_size)
    dft = zeros(Float64, n_instances, num_windows_per_inst, word_length)
    @inbounds for k=1:n_instances
        @inbounds for i=1:num_windows_per_inst-1
            dft[k, i, :] = _discrete_fourier_transform(view(X, k, ((i-1)*window_size + 1):(i*window_size)),
                                                word_length, norm, inverse_sqrt_win_size, normalise)
        end
        dft[k, end, :] = _discrete_fourier_transform(view(X, k, (series_length-window_size + 1):series_length),
                                                word_length, norm, inverse_sqrt_win_size, normalise)
    end
    total_num_windows = n_instances * num_windows_per_inst
    breakpoints = zeros(Float64, word_length, alphabet_size)
    @inbounds for letter in range(word_length)                        # sortring algo, making a array and finding the aproproate ramge to place
        sort!(dft[((letter-1)*total_num_windows + 1):(letter*total_num_windows)])
        bin_index = 0
        target_bin_depth = total_num_windows / alphabet_size
        @inbounds for bp in range(alphabet_size - 1)
            bin_index += target_bin_depth
            breakpoints[letter][bp] = dft[int(bin_index + (letter-1)*total_num_windows)]  # Check indexing here
        end
        breakpoints[letter][alphabet_size] = maxintfloat(Float64)   # check this out in the system
    return breakpoints
end

# self.window_size, self.window_size
#=
function _mcb_dft(series, num_windows_per_inst, window_size, series_length,
                        word_length, norm, inverse_sqrt_win_size, normalise)
    # Splits individual time series into windows and returns the DFT for
    # each
    Z = zeros(Float64, num_windows_per_inst, window_size)
    @inbounds for i=1:num_windows_per_inst-1
        Z[i, :] = _discrete_fourier_transform(series[((i-1)*window_size + 1):(i*window_size)],
                                            word_length, norm, inverse_sqrt_win_size, normalise)
    end
    Z[end, :] = _discrete_fourier_transform(series[ (series_length-window_size + 1):series_length],
                                            word_length, norm, inverse_sqrt_win_size, normalise)
    return Z # find the way to have two indiaces for he same loops @inoans may help
end
=#
# self.word_length, self.alphabet_size, self.breakpoints

function _create_word(self, dft)
    word = _BitWord()                       # check out this function
    for i in range(self.word_length)
        for bp in range(self.alphabet_size)
            if dft[i] <= self.breakpoints[i][bp]
                word.push(bp)
                break
            end
        end
    end
    return word
end

# self.remove_repeat_words

function _add_to_bag(self, bag, word, last_word)
    if self.remove_repeat_words && word.word == last_word
        return false
    end
    bag[word.word] = bag.get(word.word, 0) + 1   # check out this get fucntion
    return true
end

# self.levels, self.remove_repeat_words, self.window_size,  self.series_length, self.level_weights

function _add_to_pyramid(self, bag, word, last_word, window_ind)
    if self.remove_repeat_words && word.word == last_word  # check if python `and` and `&&` from the julia are same
        return false
    end
    start = 0
    for i in range(self.levels)
        num_quadrants = pow(2, i)                # Check out this fucntion
        quadrant_size = self.series_length / num_quadrants
        pos = window_ind + Int((self.window_size / 2))
        quadrant = start + (pos / quadrant_size)
        bag[(word.word, quadrant)] = (bag.get((word.word, quadrant), 0)
                                      + self.level_weights[i])
        start += num_quadrants
    end
    return True
end

function BOSSfit(X,  y, min_window, alphabet_size, word_length, norm, inverse_sqrt_win_size,
                 normalise, series_length, max_win_len_prop, max_ensemble_size)

    max_window_searches = series_length ÷ 4
    max_window = series_length * max_win_len_prop
    win_inc = (max_window - min_window) ÷ max_window_searches
    inverse_sqrt_win_size = 1 / sqrt(window_size)
    for normalise in [true, false]
        for win_size = min_window:win_inc:(max_window + 1)
            boss = _mcb(X, win_size, alphabet_size, word_length, norm, inverse_sqrt_win_size, normalise)
            _individual_train_acc(outputfrom_boss, y, )
end

function _train_predict(self, train_num)
    test_bag = self.transformed_data[train_num]
    best_dist = sys.float_info.max
    nn = None

    for n, bag in enumerate(self.transformed_data)
        if n == train_num
            continue
        end
        dist = boss_distance(test_bag, bag, best_dist)
        if dist < best_dist
            best_dist = dist
            nn = self.class_vals[n]
        end
    end
    return nn
end

function _individual_train_acc(self, boss, y, train_size, lowest_acc)
    correct = 0
    required_correct = int(lowest_acc * train_size)

    for i in range(train_size)
        if correct + train_size - i < required_correct
            return -1
        end
        c = boss._train_predict(i)

        if c == y[i]
            correct += 1
        end
    end
    return correct / train_size
end

"""
julia> linespace(2,4,5, Array{Int})
5-element Array{Int64,1}:
 2
 2
 3
 3
 4
"""
linespace(start,stop,_length, _type) =  convert(_type, floor.(LinRange(start,stop,_length)))


text:
export SeisBin,
SeisGeometry,
SeisPatch,
SeisPatchProcess,
SeisProcess,
SeisProcessHeaders,
SeisSort,
SeisUnPatch,
SeisWindow,
SeisWindowHeaders
include("SeisBin.jl")
include("SeisGeometry.jl")
include("SeisPatch.jl")
include("SeisUnPatch.jl")
include("SeisPatchProcess.jl")
include("SeisProcess.jl")
include("SeisProcessHeaders.jl")
include("SeisSort.jl")
include("SeisWindow.jl")
include("SeisWindowHeaders.jl")