text:
###############################################################################
#### AbstractAlgebra overrides ################################################

function AbstractAlgebra.factor(a::Generic.MPoly{nf_elem})
  return AbstractAlgebra.MPolyFactor.mfactor_char_zero(a)
end

function AbstractAlgebra.factor_squarefree(a::Generic.MPoly{nf_elem})
  return AbstractAlgebra.MPolyFactor.mfactor_squarefree_char_zero(a)
end

function AbstractAlgebra.factor(a::Generic.Poly{nf_elem})
  return Hecke.factor(a)
end

function AbstractAlgebra.MPolyFactor.mfactor_choose_eval_points!(
  alphas::Vector,
  A::E,
  mainvar::Int,
  minorvars::Vector{Int},
  size::Int) where E <: Generic.MPoly{nf_elem}

  n = length(minorvars)
  R = parent(A)
  K = base_ring(R)
  d = degree(K)

  size = clamp(size, 2, 100)

  for i in 1:n
    # choose integers for now
    zero!(alphas[i])
    Hecke.set_den!(alphas[i], one(ZZ))
    for j in 0:0
      Hecke.set_num_coeff!(alphas[i], j, rand_bits(ZZ, 2 + rand(0:size)))
    end
  end
end

function AbstractAlgebra.MPolyFactor.hlift_have_lcs(
  A::E,
  Auf::Vector{E},       # univariate factorization in mainvar
  lcs::Vector{E},       # lead coeffs of factors
  mainvar::Int,
  minorvars::Vector{Int},   # length n
  alphas::Vector
) where E <: Generic.MPoly{nf_elem}

  return hlift_have_lcs_generic(A, Auf, lcs, mainvar, minorvars, alphas)
#  return hlift_have_lcs_crt(A, Auf, lcs, mainvar, minorvars, alphas)
end

###############################################################################
#### implementation ###########################################################

#### the original generic method using nf_elem arithmetic #####################

function hlift_have_lcs_generic(
  A::E,                 # lead coeff should equal prod(lcs)
  Auf::Vector{E},       # univariate factorization in mainvar
  lcs::Vector{E},       # lead coeffs of factors
  mainvar::Int,
  minorvars::Vector{Int},   # length n
  alphas::Vector
) where E <: Generic.MPoly{nf_elem}

  R = parent(A)
  K = base_ring(R)
  n = length(minorvars)
  r = length(Auf)
  @assert n > 0
  @assert r > 1

  lc_evals = zeros(R, n + 1, r)
  for j in 1:r
    lc_evals[n + 1, j] = lcs[j]
    for i in n:-1:1
      lc_evals[i, j] = AbstractAlgebra.MPolyFactor.eval_one(lc_evals[i + 1, j],
                                                       minorvars[i], alphas[i])
    end
  end

  A_evals = zeros(R, n + 1)
  A_evals[n + 1] = A
  for i in n:-1:1
    A_evals[i] = AbstractAlgebra.MPolyFactor.eval_one(A_evals[i + 1],
                                                       minorvars[i], alphas[i])
  end

  fac = zeros(R, r)
  for j in 1:r
    @assert is_constant(lc_evals[1, j])
    fac[j] = Auf[j]*divexact(lc_evals[1, j],
                           AbstractAlgebra.MPolyFactor.get_lc(Auf[j], mainvar))
  end

  tdegs = degrees(A_evals[n + 1])
  liftdegs = [tdegs[minorvars[i]] for i in 1:n]

  for i in 2:n+1
    tfac = [AbstractAlgebra.MPolyFactor.set_lc(fac[j],
                                         mainvar, lc_evals[i, j]) for j in 1:r]

    ok, fac = AbstractAlgebra.MPolyFactor.hliftstep(tfac, mainvar,
                          minorvars[1:i-1], liftdegs, alphas, A_evals[i], true)
    if !ok
      return false, fac
    end
  end

  return true, fac
end

#### modular approach #########################################################

function nf_elem_split(a::nf_elem)
  den = denominator(a)
  num = Hecke.Globals.Zx(parent(defining_polynomial(parent(a)))(a)*den)
  return num, den
end

#=
  Theorem 6.2 in Martin Lee's thesis
  return (lcc::fmpz, bits::Int) such that for any monic divisor g of A:
    lcc*g is in Z[alpha][x1,x2,...]
    height(lcc*g) <= 2^bits
=#
function mfactor_martin_bounds(A::Generic.MPoly{nf_elem})
  R = parent(A)
  K = base_ring(R)
  Alen = length(A)
  @assert Alen > 0

  # mu in Z[alpha] is the defining polynomial
  mu = defining_polynomial(K)
  mu = Hecke.Globals.Zx(denominator(mu) * mu)

  # define f to be scale*A so that f in Z[alpha][X,x1,...,xn]
  # f will not be computed explicitly - we only need height(f) and leading_coefficient(f)
  leadf, scale1 = nf_elem_split(coeff(A, 1))
  heightf = height(leadf)
  scale = scale1
  for i in 2:Alen
    num, den = nf_elem_split(coeff(A, i))
    g = gcd(scale, den)
    den = divexact(den, g)
    heightf = max(heightf*den, height(num)*divexact(scale, g))
    scale *= den
  end
  leadf *= divexact(scale, scale1)

  # factors will be reconstructed with leading coeff df*D
  df = resultant(mu, leadf)
  D = resultant(mu, derivative(mu))

  degs = degrees(A)
  prod_degs = prod(fmpz(i) + 1 for i in degs)
  sum_degs = sum(fmpz(i) for i in degs)
  s = sum(i > 0 ? 1 : 0 for i in degs)  # count of present variables
  k = fmpz(degree(mu))
  bits = sum_degs + k + nbits(prod_degs)
  bits += cdiv(k*nbits(cdiv((k + 1)^7*heightf^2*height(mu)^8, leading_coefficient(mu)^2)) - s, 2)

  return (df*D, Int(bits))
end

function is_pairwise_coprime(v::Vector{E}) where E
  n = length(v)
  for i in 1:n-1, j in i+1:n
    if !isone(gcd(v[i],v[j]))
      return false
    end
  end
  return true
end

function hlift_have_lcs_crt(
  A::E,
  Auf::Vector{E},       # univariate factorization in mainvar
  lcs::Vector{E},       # lead coeffs of factors
  mainvar::Int,
  minorvars::Vector{Int},   # length n
  alphas::Vector
) where E <: Generic.MPoly{nf_elem}

  R = parent(A)
  K = base_ring(R)
  n = length(minorvars)
  r = length(Auf)
  @assert n > 0
  @assert r > 1

  nf_lcc, bits = mfactor_martin_bounds(A)

  m = fmpz(1)
  Af = zeros(R, r)

  for p in PrimesSet(Hecke.p_start, -1)

    local pA, pAuf, plcs, palphas

    if divisible(nf_lcc, p)
        continue
    end

    me = Hecke.modular_init(K, p)
    s = length(me.fld)

    try
      pA = modular_proj(fq_nmod, A, me)

      # make sure no leading coeff vanishes
      plcs = [modular_proj(fq_nmod, lc, me) for lc in lcs]
      plcs = [[plcs[j][i] for j in 1:r] for i in 1:s]
      ok = true
      for i in 1:s, j in 1:r
        ok = ok && !iszero(plcs[i][j])
      end
      if !ok
        continue
      end

      # make sure univariate factorizations remain pairwise prime
      pAuf = [modular_proj(fq_nmod, f, me) for f in Auf]
      pAuf = [[pAuf[j][i] for j in 1:r] for i in 1:s]
      ok = true
      for i in 1:s
        ok = ok && is_pairwise_coprime(pAuf[i])
      end
      if !ok
        continue
      end

      palphas = [deepcopy(Hecke.modular_proj(alphas[j], me)) for j in 1:n]
      palphas = [[palphas[j][i] for j in 1:n] for i in 1:s]

    catch ee
      if ee <: Hecke.BadPrime
        continue
      end
      rethrow(ee)
    end

    pAf = [eltype(pA)[] for j in 1:r]
    for i in 1:s
      ok, t = AbstractAlgebra.MPolyFactor.hlift_have_lcs(
                       pA[i], pAuf[i], plcs[i], mainvar, minorvars, palphas[i])
      if !ok
        return false
      end

      for j in 1:r
        push!(pAf[j], (nf_lcc*inv(leading_coefficient(t[j]))*t[j]))
      end
    end

    for j in 1:r
      t = modular_lift(pAf[j], me)
      Af[j] = Hecke.induce_crt(Af[j], m, t, fmpz(p), true)[1]
    end
    m *= p

    if nbits(m) > bits + 1
      break
    end
  end

  # scale the Af[j] to match the given lcs if possible
  for j in 1:r
    ok, c = divides(lcs[j], AbstractAlgebra.MPolyFactor.get_lc(Af[j], mainvar))
    if !ok || !is_constant(c)
      return false
    end
    Af[j] *= c
  end

  ok = (A == prod(f for f in Af))
  return ok, Af
end


text:
# Weave readme
using Pkg
cd("c:/git/TidyStanza/")
Pkg.activate("c:/git/TidyStanza/readme-env")

using Weave

weave("README.jmd", out_path = :pwd, doctype = "github")

if false
    tangle("README.jmd")
end

ext: jl


text:
function set_mesh!(vis)
    l,w,h = 0.4,0.4,1.0
    body = Sphere(Point3f0(0,0,0), 7l/16)
    setobject!(vis["robot"]["torso"]["body"], body,
                            MeshPhongMaterial(color=colorant"yellow"))
    axle = Cylinder(Point3f0(0,0,0), Point3f0(0,w/2,0), 0.03f0)
    setobject!(vis["robot"]["torso"]["Laxle"], axle,
                            MeshPhongMaterial(color=colorant"black"))
    setobject!(vis["robot"]["torso"]["Raxle"], axle,
                            MeshPhongMaterial(color=colorant"black"))
    settransform!(vis["robot"]["torso"]["Laxle"], Translation(0,+l/4,0))
    settransform!(vis["robot"]["torso"]["Raxle"], Translation(0,-3l/4,0))

    foot = HyperSphere(Point3f0(0,0,0f0), 0.05f0)
    Lleg = Cylinder(Point3f0(0,+l/2,0), Point3f0(0,+l/2,1), 0.03f0)
    Rleg = Cylinder(Point3f0(0,-l/2,0), Point3f0(0,-l/2,1), 0.03f0)
    setobject!(vis["robot"]["torso"]["Lleg"]["geom"], Lleg,
                        MeshPhongMaterial(color=colorant=colorant"blue"))
    setobject!(vis["robot"]["torso"]["Rleg"]["geom"], Rleg,
                        MeshPhongMaterial(color=colorant=colorant"green"))

    Lfoot = setobject!(vis["robot"]["torso"]["Lleg"]["geom"]["Lfoot"]["geom"],
                    foot, MeshPhongMaterial(color=colorant"firebrick"))
    setobject!(vis["robot"]["torso"]["Rleg"]["geom"]["Rfoot"]["geom"],
                    foot, MeshPhongMaterial(color=colorant"firebrick"))
    settransform!(vis["robot"]["torso"]["Lleg"]["geom"]["Lfoot"]["geom"],
                    Translation(0,+l/2,h))
    settransform!(vis["robot"]["torso"]["Rleg"]["geom"]["Rfoot"]["geom"],
                    Translation(0,-l/2,h))
return Lfoot
end

function animate_walk!(vis, x)
    xb,yb,zb = x[1], x[2], ones(length(x[1]))
    xl,yl,zl = x[3], x[4], x[5]
    xr,yr,zr = x[6], x[7], x[8]
    dt = 0.1
    anim = Animation()
    for i=1:length(xb)
        atframe(anim, i) do
            Llen = norm(SA[xl[i]-xb[i], zl[i]-zb[i]])
            Rlen = norm(SA[xr[i]-xb[i], zr[i]-zb[i]])
            θl = atan(xl[i]-xb[i], zl[i]-zb[i])
            θr = atan(xr[i]-xb[i], zr[i]-zb[i])
            settransform!(vis["robot"]["torso"]["Lleg"], LinearMap(RotY(θl)))
            settransform!(vis["robot"]["torso"]["Rleg"], LinearMap(RotY(θr)))
            settransform!(vis["robot"]["torso"], Translation(xb[i],yb[i],zb[i]))
            settransform!(vis["robot"]["torso"]["Lleg"]["geom"],
                                    LinearMap(Diagonal(SA[1,1,Llen])))
            settransform!(vis["robot"]["torso"]["Rleg"]["geom"],
                                    LinearMap(Diagonal(SA[1,1,Rlen])))
        end
    end
    setanimation!(vis, anim)
end

function animate_walking_trajectory(X)
    vis  = Visualizer()
    setprop!(vis["/Background"], "top_color", colorant"white")
    setprop!(vis["/Background"], "bottom_color", colorant"white")
    setprop!(vis["/Grid"], "visible",false)
    set_mesh!(vis)
    animate_walk!(vis, X)
    render(vis)
end

# X = move_to_position([0.0,0.0,1.0], goal_position=[5.0, 0.0, 1.0])
# animate_walking_trajectory(X)


text:
# system/z-test.jl

using Test

@test ir_lm_test_getindex_cumsum()
@test hcat_lm(:test)


text:
using Base64

function addcol!(df::DataFrame, k, val)
    col_val = Vector{Union{Missing, typeof(val)}}(undef, nrow(df))
    col_val .= missing
    col_val[end] = val
    setproperty!(df, k, col_val)
end

function setcol!(df::DataFrame, i, k, val)
    if string(k) ∈ names(df)
    	T = eltype(df[k])
        df[i, k] = _convert(T,val)
    else
        addcol!(df, k, val)
    end
end

Base.push!(df::DataFrame) = push!(df, [missing for i=1:ncol(df)])

_convert(T::Type{Union{Missing, String}}, val::JSON3.Object{Base.CodeUnits{UInt8,String},SubArray{UInt64,1,Array{UInt64,1},Tuple{UnitRange{Int64}},true}}) = begin
    ks = keys(val)
    if length(ks) == 1 && first(ks) == Symbol("py/b64")
        return convert(T, String(base64decode(val[first(ks)])))
    end
end
_convert(T::Type{<:Union{Missing, String}}, val::Number) = "$val"
_convert(T::Type, val) = convert(T, val)


text:
module Hygienic

Base.@kwdef struct ExpansionContext
    already_used::Set{Symbol} = Set{Symbol}()
    parent_context::Union{Nothing,ExpansionContext} = nothing
    name_map::Dict{Symbol,Symbol} = Dict{Symbol,Symbol}()
end
ExpansionContext(ctx::ExpansionContext) =
    ExpansionContext(; parent_context = ctx)

function get_parent_sym(ctx::ExpansionContext, s::Symbol)
    if ctx.parent_context === nothing
        push!(ctx.already_used, s)
    end
    get_sym(ctx.parent_context, s)
end

get_sym(ctx::ExpansionContext, s::Symbol) =
    haskey(ctx.name_map, s) ? ctx.name_map[s] : get_parent_sym(ctx, s)
    get_sym(::Nothing, s::Symbol) = s # the symbol is not aliased because it is not defined in the current expr

function set_sym!(ctx::ExpansionContext, s::Symbol)
    if s ∈ ctx.already_used
        error("$s was used before definition, use a different name to prevent confusion.")
    end

    if haskey(ctx.name_map, s)
        error("$s was already defined, note that branching is not supported.")
    end
    ctx.name_map[s] = gensym(s)
end

function sanitize_arglist!(ctx, ex::Expr)
    if ex.head == :tuple
       map!(
            s -> s isa Symbol ? set_sym!(ctx, s) : s,
            ex.args,
            ex.args,
       )
       return ex
    end

    error("Unvalid arglist type $(ex.head)")
end

isexpr(ex::Expr, heads) = any(h -> Meta.isexpr(ex, h), heads)
isexpr(ex::Expr, s::Symbol) = Meta.isexpr(ex, s)

function sanitize!(ctx, ex::Expr)::Expr
    if ex.head == :(=)
        value = ex.args[2]
        if ex.args[2] isa Expr
            sanitize!(ctx, ex.args[2])
        elseif value isa Symbol
            ex.args[2] = get_sym(ctx, ex.args[2])
        end

        assignee = ex.args[begin]
        if assignee isa Symbol
            ex.args[begin] = set_sym!(ctx, assignee)
        elseif Meta.isexpr(assignee, :tuple)
            sanitize_arglist!(ctx, assignee)
        elseif Meta.isexpr(assignee, :ref, 1)
            if assignee.args[1] isa Symbol
                assignee.args[1] = get_sym(ctx, assignee.args[1])
            end
        elseif Meta.isexpr(assignee, :$)
            nothing
        else
            dump(assignee)
            error("Unhandled assignee $(assignee)")
        end

        return ex
    elseif ex.head == :function
        func_ctx = ExpansionContext(ctx)
        sanitize_arglist!(func_ctx, ex.args[1])
        sanitize!(func_ctx, ex.args[2])
        return ex
    elseif ex.head == :(->)
        ctx = ExpansionContext(ctx)
        if ex.args[1] isa Expr
            sanitize_arglist!(ctx, ex.args[1])
        elseif ex.args[1] isa Symbol
            ex.args[1] = set_sym!(ctx, ex.args[1])
        end
        if ex.args[2] isa Expr
            sanitize!(ctx, ex.args[2])
        elseif ex.args[2] isa Symbol
            ex.args[2] = get_sym(ctx, ex.args[2])
        end
        return ex
    elseif isexpr(ex, (:ref,))
        if ex.args[1] isa Symbol
            ex.args[1] = get_sym(ctx, ex.args[1])
        elseif ex.args[1] isa Expr
            sanitize!(ctx, ex.args[1])
        end
        return ex
    elseif isexpr(ex, (:quote,))
        sanitize!(ctx, ex.args[1])
        return ex
    elseif isexpr(ex, :local)
        if ex.args[1] isa Expr
            sanitize!(ctx, ex.args[1])
        else
            map!(s -> set_sym!(ctx, s), ex.args, ex.args)
        end
        return ex
    elseif isexpr(ex, :do)
        sanitize!(ctx, ex.args[1])
        sanitize!(ctx, ex.args[2])
        return ex
    elseif isexpr(ex, (:block, :call, :if, :elseif, :&&, :||, :curly, :tuple, :vect))
        map!(
            a -> a isa Expr ? sanitize!(ctx, a) : a isa Symbol ? get_sym(ctx, a) : a,
            ex.args,
            ex.args,
        )
        return ex
    elseif isexpr(ex, :macrocall)
        # TODO: Support DSL macros
        map!(
            a -> a isa Expr ? sanitize!(ctx, a) : a isa Symbol ? get_sym(ctx, a) : a,
            ex.args,
            ex.args,
        )
        return ex
    elseif isexpr(ex, :.)
        child_ex = ex
        while isexpr(child_ex, :.) && child_ex.args[1] isa Expr
            child_ex = child_ex.args[1]
        end
        if child_ex.args[1] isa Symbol
            child_ex.args[1] = get_sym(ctx, child_ex.args[1])
        else
            error("Invalid selector $(child_ex.args[1])")
        end
        return ex
        return ex
    elseif isexpr(ex, :$)
        return ex
    else
        dump(ex)
        error("Unhandled expr type $(ex.head)")
    end
    error("unreachable")
end

"""
    @hygienize(ex::Expr)

A macro to transform a quote by replacing variable definitions with gensymed symbols.

```julia
julia> using Hygienic

julia> @hygienize quote
           x = 1
           x + y
       end
quote
    #= REPL[2]:2 =#
    var"##x#274" = 1
    #= REPL[2]:3 =#
    var"##x#274" + y
end
```
"""
macro hygienize(ex)
    @assert Meta.isexpr(ex, :quote) "Expected quote but got $(ex)"

    sanitize!(ExpansionContext(), ex)

    esc(ex)
end

let # For Precompilation?
    @hygienize quote
        x = 1
        x + y
    end
end

export @hygienize

end


text:
function generateposition()
    # Placeholder knights
    rank = ['♘', '♘', '♘', '♘', '♘', '♘', '♘', '♘']
    lrank = length(rank)

    # Check if a space is available
    isfree(x::Int) = rank[x] == '♘'

    # Place the King
    rank[indking = rand(2:lrank-1)] = '♔'

    # Place rooks
    rank[indrook = rand(filter(isfree, 1:lrank))] = '♖'
    if indrook > indking
        rank[rand(filter(isfree, 1:indking-1))] = '♖'
    else
        rank[rand(filter(isfree, indking+1:lrank))] = '♖'
    end

    # Place bishops
    rank[indbish = rand(filter(isfree, 1:8))] = '♗'
    pbish = filter(iseven(indbish) ? isodd : iseven, 1:lrank)
    rank[rand(filter(isfree, pbish))] = '♗'

    # Place queen
    rank[rand(filter(isfree, 1:lrank))] = '♕'
    return rank
end

@show generateposition()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
##########
# Tokens #
##########

module Kinds
@enum(Kind,
    EOF,
    IDENTIFIER,
    NUMBER,
    begin_keywords,
        DEF,
        EXTERN,
        IF,
        THEN,
        ELSE,
        FOR,
        IN,
        VAR,
    end_keywords,
    begin_operators,
        PLUS,
        MINUS,
        SLASH,
        STAR,
    end_operators,
    begin_comparisons,
        LESS,
        GREATER,
        EQUAL,
    end_comparisons,
    begin_delimiters,
        LPAR,
        RPAR,
        SEMICOLON,
        COMMA,
        LBRACE,
        RBRACE,
    end_delimiters,
    )
end

# Lookup-table for keywords
const KEYWORD_KINDS = Dict{String, Kinds.Kind}()
for kind in instances(Kinds.Kind)
    if Kinds.begin_keywords < kind < Kinds.end_keywords
        KEYWORD_KINDS[lowercase(string(kind))] = kind
    end
end

struct Token
    kind::Kinds.Kind
    val::String
end
Token(k::Kinds.Kind) = Token(k, "")

function Base.show(io::IO, t::Token)
    print(io, string(t.kind))
    if !isempty(t.val)
        print(io, ": ", t.val)
    end
end

#########
# Lexer #
#########

mutable struct Lexer{IO_t <: IO}
    io::IO_t
    last_char::Char
    buffer_scratch::IOBuffer
end

Lexer(io::IO) = Lexer(io, EOF_CHAR, IOBuffer())
Lexer(str::String) = Lexer(IOBuffer(str))

function Base.collect(l::Lexer)
    tokens = Token[]
    while true
        t = gettok(l)
        push!(tokens, t)
        t.kind == Kinds.EOF && break
    end
    return tokens
end

const EOF_CHAR = convert(Char,typemax(UInt32))

# Reading
readchar(io::IO) = eof(io) ? EOF_CHAR : read(io, Char)
function readchar(l::Lexer, dowrite = true)
    if l.last_char != EOF_CHAR
        c = l.last_char
        l.last_char = EOF_CHAR
    else
        c = readchar(l.io)
    end
    dowrite && write(l.buffer_scratch, c)
    return c
end
discardchar(l::Lexer) = readchar(l, false)

function readwhile(f, l::Lexer, dowrite=true)
    while true
        c = readchar(l, false)
        if !f(c)
            l.last_char = c
            break
        end
        dowrite && write(l.buffer_scratch, c)
    end
end
discardwhile(f, l::Lexer) = readwhile(f, l, false)

function gettok(l::Lexer)::Token
    t = _gettok(l)
    seek(l.buffer_scratch, 0)
    return t
end

function _gettok(l::Lexer)::Token
    discardwhile(isspace, l)
    if l.last_char == '#'
        discardwhile(x -> (x == 'r' || x == '\n'), l)
        discardchar(l) # discard the new line
    end
    c = readchar(l)
    if c == EOF_CHAR; return Token(Kinds.EOF)
    elseif isalpha(c); return lex_alpha(l)
    # number
    elseif (c == '.' || isdigit(c)); return lex_number(l)
    # delimiters
    elseif c == '('; return Token(Kinds.LPAR)
    elseif c == ')'; return Token(Kinds.RPAR)
    elseif c == '{'; return Token(Kinds.LBRACE)
    elseif c == '}'; return Token(Kinds.RBRACE)
    elseif c == ';'; return Token(Kinds.SEMICOLON)
    elseif c == ','; return Token(Kinds.COMMA)
    # comparisons
    elseif c == '<'; return Token(Kinds.LESS)
    elseif c == '>'; return Token(Kinds.GREATER)
    elseif c == '='; return Token(Kinds.EQUAL)
    # operators
    elseif c == '+'; return Token(Kinds.PLUS)
    elseif c == '-'; return Token(Kinds.MINUS)
    elseif c == '/'; return Token(Kinds.SLASH)
    elseif c == '*'; return Token(Kinds.STAR)
    else error("failed to handle character $c in lexer")
    end
end

function lex_alpha(l::Lexer)
    isvalididentifier(x) = isalpha(x) || isdigit(x)
    readwhile(isvalididentifier, l)
    str = String(take!(l.buffer_scratch))
    if haskey(KEYWORD_KINDS, str)
        return Token(KEYWORD_KINDS[str])
    else
        return Token(Kinds.IDENTIFIER, str)
    end
end

# TODO: Better lexing of numbers
function lex_number(l)
    isvalid_digit = x -> isdigit(x) || x == '.'
    readwhile(isvalid_digit, l)
    num = String(take!(l.buffer_scratch))
    return Token(Kinds.NUMBER, num)
end


text:
using BenchmarkTools
import Base: *, length, size, getindex

println("""
    1. A one-hot vector can be described by its length and the index of
    its hot bit.
""")

println("""
    2. A one-hot vector composite type is defined with the above two
    fields and an inner constructor method that throws an error when the
    arguments are in the wrong order.
""")

struct OneHot
    index::UInt
    length::UInt
    function OneHot(index, length)
        if index > length
            error("index should be smaller than length")
        end
        new(index, length)
    end
end

println("OneHot(6, 9) = $(OneHot(6, 9))\n")

println("""
    3. To make the innersum() function work, a length() and a getindex()
    methods are added along with a method for multiplying a matrix with
    a one-hot vector.
""")

*(a::AbstractMatrix, v::OneHot) = a[:, v.index]
length(v::OneHot) = v.length
getindex(v::OneHot, i::Integer) = i == v.index

function innersum(A, vs)
    t = zero(eltype(A)) # generic!
    for v in vs
        y = A * v
        for i = 1:length(vs[1])
            t += v[i] * y[i]
        end
    end
    return t
end

A = rand(3, 3)
vs = [rand(3) for _ = 1:10]
vs_onehot = [OneHot(rand(1:3), 3) for _ = 1:10]

println("innersum(A, vs) = $(innersum(A, vs))")
println("innersum(A, vs_onehot) = $(innersum(A, vs_onehot))\n")

println("""
    4. Comparison of the speed of innersum() with a OneHot and innersum()
    with a Vector{Float64}:
""")

print("@btime (innersum(A, vs)):")
@btime (innersum(A, vs))
print("@btime (innersum(A, vs_onehot)):")
@btime (innersum(A, vs_onehot))
println()

println("""
    5. A OneHotVector composite type is defined. It is identical to
    OneHot but is declared to be a subtype of AbstractVector{Bool}.
""")

struct OneHotVector <: AbstractVector{Bool}
    index::UInt
    length::UInt
    function OneHotVector(index, length)
        if index > length
            error("index should be smaller than length")
        end
        new(index, length)
    end
end

size(v::OneHotVector) = (v.length,)
getindex(v::OneHotVector, i::Integer) = i == v.index

println("OneHotVector(6, 9) = $(OneHotVector(6, 9))\n")

println("""
    6. Comparison of the speed of innersum() with a OneHot, with a
    OneHotVector and a Vector{Float64}.
""")

vs_onehotvector = [OneHotVector(rand(1:3), 3) for _ = 1:10]

print("@btime (innersum(A, vs)):")
@btime (innersum(A, vs))
print("@btime (innersum(A, vs_onehot)):")
@btime (innersum(A, vs_onehot))
print("@btime (innersum(A, vs_onehotvector)):")
@btime (innersum(A, vs_onehotvector))


text:
import DocSeeker: dynamicsearch

function firstN(matches, desired, N = 3)
  binds = map(x -> x[2].name, matches[1:N])
  for d in desired
    if !(d in binds)
      return false
    end
  end
  return true
end

@testset "dynamicsearch" begin
  @test firstN(dynamicsearch("sine"), ["sin", "sind", "asin"], 20)
  @test firstN(dynamicsearch("regular expression"), ["match", "eachmatch", "replace"], 20)

  @test dynamicsearch("Real")[1][2].name == "Real"
  @test length(dynamicsearch("Real")[1][2].text) > 0
  @test dynamicsearch("regex")[1][2].name == "Regex"

  let downloadsearch = dynamicsearch("download")
    dfound = 0
    for d in downloadsearch
      if d[2].name == "download" && d[2].mod == "Base"
        dfound += 1
      end
    end
    @test dfound == 1
  end
end


text:
using Documenter, PAPA

makedocs(;
    modules=[PAPA],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://qiplab.bbn.com/matthewware/PAPA.jl/blob/{commit}{path}#L{line}",
    sitename="PAPA.jl",
    authors="Luke C. Govia",
    assets=String[],
)


text:
using Revise, ApproxOperator

elements, nodes = importmsh("./msh/bar.msh")
nₚ = length(nodes[:x])
nₑ = length(elements["Domain"])

type = (SNode,:Quadratic1D,:□,:CubicSpline)
sp = RegularGrid(nodes[:x],nodes[:y],nodes[:z],n = 2,γ = 1)
elements["Domain"] = ReproducingKernel{type...,:Seg2}(elements["Domain"])
elements["NBC"] = ReproducingKernel{type...,:Poi1}(elements["NBC"])
elements["EBC"] = ReproducingKernel{type...,:Poi1}(elements["EBC"])
sp(elements["Domain"])
sp(elements["NBC"])
sp(elements["EBC"])
elements["DomainS"] = ReproducingKernel{type...,:Seg2}(elements["Domain"])
s = 0.25*ones(nₚ)
nodes[:s₁] = s
nodes[:s₂] = s
nodes[:s₃] = s

set𝓖!(elements["Domain"],:SegRK3,:∂1,:∂x,:∂y,:∂z)
set𝓖!(elements["DomainS"],:SegGI2,:∂1,:∂x,:∂y,:∂z)
set𝓖!(elements["NBC"],:PoiGI1,:∂1)
set𝓖!(elements["EBC"],:PoiGI1,:∂1,:∂x,:∂y,:∂z)

elements["EBC"] = ReproducingKernel{type...,:Seg2}(elements["Domain"],elements["EBC"],sharing=false)

set𝒏!(elements["Domain"])

set𝝭!(elements["Domain"])
set∇̃𝝭!(elements["DomainS"],elements["Domain"])
set∇𝝭!(elements["EBC"])
set𝝭!(elements["NBC"])

r = 3
prescribe!(elements["Domain"],:b,(x,y,z)->-r*(r-1)*x^abs(r-2))
prescribe!(elements["NBC"],:t,(x,y,z)->r*x^abs(r-1))
prescribe!(elements["EBC"],:g,(x,y,z)->x^r)

coefficient = (:k=>1.0,:α=>1e3)
ops = [Operator(:∫∇v∇udΩ,coefficient...),
       Operator(:∫vbdΩ,coefficient...),
       Operator(:∫vtdΓ,coefficient...),
       Operator(:∫∇𝑛vgdΓ,coefficient...),
       Operator(:∫vgdΓ,coefficient...),
       Operator(:H₁)]

k = zeros(nₚ,nₚ)
f = zeros(nₚ)

ops[1](elements["DomainS"],k)
ops[2](elements["Domain"],f)
ops[3](elements["NBC"],f)
ops[4](elements["EBC"],k,f)
ops[5](elements["EBC"],k,f)
d = k\f

push!(nodes,:d=>d)
set𝓖!(elements["Domain"],:SegGI5,:∂1,:∂x,:∂y,:∂z)
prescribe!(elements["Domain"],:u,(x,y,z)->x^r)
prescribe!(elements["Domain"],:∂u∂x,(x,y,z)->r*x^abs(r-1))
set∇𝝭!(elements["Domain"])
h1, l2 = ops[6](elements["Domain"])

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
## sample_method.jl : types to represent sample methods (MC and QMC)
#
# Representation of a sampling method.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

abstract type AbstractSampleMethod end

## MC ##
"""
    MC()

Return a Monte Carlo sample method.

See also: [`QMC`](@ref)
"""
struct MC <: AbstractSampleMethod end

## QMC ##
"""
    QMC()

Return a Quasi-Monte Carlo sample method.

See also: [`MC`](@ref)
"""
struct QMC <: AbstractSampleMethod end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Bison"
version = v"3.5.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.gnu.org/gnu/bison/bison-3.5.2.tar.xz", "24e273db9eb6da8bbb6f0648284d0724a5cbd6268a163db402f961350a4e50dd"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/bison-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# Disable windows for now
platforms = filter(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("bison", :bison_exe)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using Base.Math: @horner
using Base.MPFR: ROUNDING_MODE
#useful constants
const acc0 = [5.0e-15 , 5.0e-7 , 5.0e-4] #accuracy options
const big1 = [25.0 , 14.0, 10.0]
const e0 = [0.25e-3 , 0.25e-1 , 0.14]
const x0=[31.0 , 17.0 , 9.7]
const alog10 = log(10)
const rt2pin = 1.0/sqrt(2*pi)
const rtpi = sqrt(pi)
const exparg = -745.1
const stirling_coef = [1.996379051590076518221, -0.17971032528832887213e-2, 0.131292857963846713e-4, -0.2340875228178749e-6, 0.72291210671127e-8, -0.3280997607821e-9, 0.198750709010e-10, -0.15092141830e-11, 0.1375340084e-12, -0.145728923e-13, 0.17532367e-14, -0.2351465e-15, 0.346551e-16, -0.55471e-17, 0.9548e-18, -0.1748e-18, 0.332e-19, -0.58e-20]
const auxgam_coef = [-1.013609258009865776949, 0.784903531024782283535e-1, 0.67588668743258315530e-2, -0.12790434869623468120e-2, 0.462939838642739585e-4, 0.43381681744740352e-5, -0.5326872422618006e-6, 0.172233457410539e-7, 0.8300542107118e-9, -0.10553994239968e-9, 0.39415842851e-11, 0.362068537e-13, -0.107440229e-13, 0.5000413e-15, -0.62452e-17, -0.5185e-18, 0.347e-19, -0.9e-21]

#----------------COEFFICIENTS FOR TEMME EXPANSION------------------

const d00 = -.333333333333333E+00
const d0=[.833333333333333E-01 , -.148148148148148E-01 , .115740740740741E-02 , .352733686067019E-03 , -.178755144032922E-03 , .391926317852244E-04]
const d10 = -.185185185185185E-02
const d1=[-.347222222222222E-02 , .264550264550265E-02 , -.990226337448560E-03 , .205761316872428E-03]
const d20 = .413359788359788E-02
const d2=[-.268132716049383E-02 , .771604938271605E-03]
const d30 = .649434156378601E-03
const d3=[.229472093621399E-03 , -.469189494395256E-03]
const d40 = -.861888290916712E-03
const d4=[.784039221720067E-03]
const d50 = -.336798553366358E-03
const d5=[-.697281375836586E-04]
const d60 = .531307936463992E-03
const d6=[-.592166437353694E-03]
const d70 = .344367606892378E-03
const d80 = -.652623918595309E-03
#Source of logmxp1(x): https://github.com/JuliaStats/StatsFuns.jl/blob/master/src/basicfuns.jl
# The kernel of log1pmx
# Accuracy within ~2ulps for -0.227 < x < 0.315
function _log1pmx_ker(x::Float64)
    r = x/(x+2.0)
    t = r*r
    w = @horner(t,
                6.66666666666666667e-1, # 2/3
                4.00000000000000000e-1, # 2/5
                2.85714285714285714e-1, # 2/7
                2.22222222222222222e-1, # 2/9
                1.81818181818181818e-1, # 2/11
                1.53846153846153846e-1, # 2/13
                1.33333333333333333e-1, # 2/15
                1.17647058823529412e-1) # 2/17
    hxsq = 0.5*x*x
    r*(hxsq+w*t)-hxsq
end
"""
    log1pmx(x::Float64)
Return `log(1 + x) - x`
Use naive calculation or range reduction outside kernel range.  Accurate ~2ulps for all `x`.
"""
function log1pmx(x::Float64)
    if !(-0.7 < x < 0.9)
        return log1p(x) - x
    elseif x > 0.315
        u = (x-0.5)/1.5
        return _log1pmx_ker(u) - 9.45348918918356180e-2 - 0.5*u
    elseif x > -0.227
        return _log1pmx_ker(x)
    elseif x > -0.4
        u = (x+0.25)/0.75
        return _log1pmx_ker(u) - 3.76820724517809274e-2 + 0.25*u
    elseif x > -0.6
        u = (x+0.5)*2.0
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        u = (x+0.625)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    end
end
"""
    logmxp1(x::Float64)
Return `log(x) - x + 1` carefully evaluated.
"""
function logmxp1(x::Float64)
    if x <= 0.3
        return (log(x) + 1.0) - x
    elseif x <= 0.4
        u = (x-0.375)/0.375
        return _log1pmx_ker(u) - 3.55829253011726237e-1 + 0.625*u
    elseif x <= 0.6
        u = 2.0*(x-0.5)
        return _log1pmx_ker(u) - 1.93147180559945309e-1 + 0.5*u
    else
        return log1pmx(x - 1.0)
    end
end

"""
   rgamma1pm1(a)

   Computation of ``1/Gamma(a+1) - 1`` for `-0.5<=a<=1.5` : ``1/\\Gamma (a+1) - 1``
   Uses the relation `gamma(a+1) = a*gamma(a)`.
"""
function rgamma1pm1(a::Float64)
    t=a
    rangereduce = a > 0.5
    t = rangereduce ? a-1 : a #-0.5<= t <= 0.5
    if t == 0.0
        return 0.0
    elseif t < 0.0
        top = @horner(t , -.422784335098468E+00 , -.771330383816272E+00 , -.244757765222226E+00 , .118378989872749E+00 , .930357293360349E-03 , -.118290993445146E-01 , .223047661158249E-02 , .266505979058923E-03 , -.132674909766242E-03)
        bot = @horner(t , 1.0 , .273076135303957E+00 , .559398236957378E-01)
        w = top/bot
        return rangereduce ? t*w/a : a*(w+1)
    else
        top = @horner(t , .577215664901533E+00 , -.409078193005776E+00 , -.230975380857675E+00 , .597275330452234E-01 , .766968181649490E-02 , -.514889771323592E-02 , .589597428611429E-03)
        bot = @horner(t , 1.0 , .427569613095214E+00 , .158451672430138E+00 , .261132021441447E-01 , .423244297896961E-02)
        w = top/bot
        return rangereduce ? (t/a)*(w - 1.0) : a*w
    end
end
"""
    rgammax(a,x)

Evaluation of ``1/\\Gamma(a) e^{-x} x^{a}``.
"""
function rgammax(a::Float64,x::Float64)
    if x == 0.0
        return 0.0
    elseif a >= 20.0
        u =x/a
        if u == 0.0
            return 0.0
        end
        t = (1.0/a)^2
        t1 = (((0.75*t - 1.0)*t + 3.5)*t - 105.0)/(a*1260.0) #Using stirling Series : https://dlmf.nist.gov/5.11.1
        t1 = t1 + a*logmxp1(u)
        if t1 >= exparg
            return rt2pin*sqrt(a)*exp(t1)
        end
    else
        t = a*log(x) - x
        if t < exparg
            return 0.0
        elseif a >= 1.0
            return exp(t)/gamma(a)
        else
            return (a*exp(t))*(1.0 + rgamma1pm1(a))
        end
    end
end
"""
    auxgam(x)

Compute function `g` in ``1/\\Gamma(x+1) = 1+x*(x-1)*g(x)``, `-1 <= x <= 1`.
"""
function auxgam(x::Float64)
    if x < 0
        return -(1.0+(1+x)*(1+x)*auxgam(1+x))/(1.0-x)
    else
        t = 2*x - 1.0
        return chepolsum(t, auxgam_coef)
    end
end

"""
    loggamma1p(a)

Compute ``log(\\Gamma(1+a))`` for `-1 < a <= 1`.
"""
function loggamma1p(x::Float64)
    return -log1p(x*(x-1.0)*auxgam(x))
end

"""
    chepolsum(n,x,a)

Computes a series of Chebyshev Polynomials given by: `a[1]/2 + a[2]T1(x) + .... + a[n]T{n-1}(X)`.
"""
function chepolsum(x::Float64, a::Array{Float64,1})
    n = length(a)
    if n == 1
        return a[1]/2.0
    elseif n == 2
        return a[1]/2.0 + a[2]*x
    else
        tx = 2*x
        r = a[n]
        h = a[n-1] + r*tx
        for k = n-2:-1:2
            s=r
            r=h
            h=a[k]+r*tx-s
        end
        return a[1]/2.0 - r + h*x
    end
end

"""
    stirling(x)

Compute ``\\ln{\\Gamma(x)} - (x-0.5)*\\ln{x} + x - \\ln{(2\\pi)}/2`` using stirling series for asymptotic part
in https://dlmf.nist.gov/5.11.1
"""
function stirling(x::Float64)
    if x < floatmin(Float64)*1000.0
        return floatmax(Float64)/1000.0
    elseif x < 1
        return loggamma1p(x) - (x+0.5)*log(x) + x - log((2*pi))/2.0
    elseif x < 2
        return loggamma1p(x-1) - (x-0.5)*log(x) + x - log((2*pi))/2.0
    elseif x < 3
        return loggamma1p(x-2) - (x-0.5)*log(x) + x  - log((2*pi))/2.0 + log(x-1)
    elseif x < 12
        z=18.0/(x*x)-1.0
        return chepolsum(z, stirling_coef)/(12.0*x)
    else
        z = 1.0/(x*x)
        if x < 1000
            return @horner(z, 0.25721014990011306473e-1, 0.82475966166999631057e-1, -0.25328157302663562668e-2, 0.60992926669463371e-3, -0.33543297638406e-3, 0.250505279903e-3)/((0.30865217988013567769 + z)*x)
        else
            return (((-z/1680.0+1.0/1260.0)*z-1.0/360.0)*z+1.0/12.0)/x
        end
    end
end
@doc raw"""
    gammax(x)

```math
\operatorname{gammax}(x) = \begin{cases}e^{\operatorname{stirling}(x)}\quad\quad\quad \text{if} \quad x>0,\\
\frac{\Gamma(x)}{\sqrt{2 \pi}e^{-x + (x-0.5)\operatorname{log}(x)}},\quad \text{if} \quad x\leq 0.
\end{cases}
```
"""
function gammax(x::Float64)
    if x >= 3
        return exp(stirling(x))
    elseif x > 0
        return gamma(x)/(exp(-x+(x-0.5)*log(x))*sqrt(2*pi))
    else
        return floatmax(Float64)/1000.0
    end
end
"""
    lambdaeta(eta)

Compute the value of eta satisfying ``eta^{2}/2 = \\lambda-1-\\ln{\\lambda}``.
"""
function lambdaeta(eta::Float64)
    s = eta*eta*0.5
    if eta == 0.0
        la = 1
    elseif eta < -1.0
        r = exp(-1-s)
        la = @horner(r, 0.0, 1.0, 1.0, 1.5, 8.0/3.0, 125.0/24.0, 15.0/5.0)
    elseif eta < 1.0
        r = eta
        la = @horner(r, 1.0, 1.0, 1.0/3.0, 1.0/36.0, -1.0/270.0, 1.0/4320.0, 1.0/17010.0)
    else
        r = 11 + s
        l = log(r)
        la = r+l
        r = 1.0/r
        ak1 = 1.0
        ak2 = (2 - l)*0.5
        ak3 = (@horner(l,6,-9,2))/6.0
        ak4 = -(@horner(l,-12,36,-22,3))/12.0
        ak5 = (@horner(l,60,-300,350,-125,12))/60.0
        ak6 = -(@horner(l,-120,900,-1700,1125,-274,20))/120.0
        la = la + l*@horner(r,0.0, ak1, ak2, ak3, ak4, ak5, ak6)
    end
    r = 1
    if (eta>-3.5 && eta<-0.03) || (eta>0.03 && eta<40)
        r=1
        q=la
        while r>1.0e-8
            la = q*(s+log(q))/(q-1.0)
            r = abs(q/la-1)
            q=la
        end
    end
    return la
end
"""
Computing the first coefficient for the expansion :
```math
\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3}
```
Refer Eqn (3.12) in the paper
"""
function coeff1(eta::Float64)
    if abs(eta) < 1.0
        coeff1 = @horner(eta, -3.333333333438e-1, -2.070740359969e-1, -5.041806657154e-2, -4.923635739372e-3, -4.293658292782e-5) / @horner(eta, 1.000000000000e+0, 7.045554412463e-1, 2.118190062224e-1, 3.048648397436e-2, 1.605037988091e-3)
    else
        la = lambdaeta(eta)
        coeff1 = log(eta/(la-1.0))/eta
    end
    return coeff1
end
"""
Computing the second coefficient for the expansion :
```math
\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3}
```
Refer Eqn (3.12) in the paper
"""
function coeff2(eta::Float64)

    if eta < -5.0
        x=eta*eta
        lnmeta = log(-eta)
        coeff2 = (12.0 - x - 6.0*lnmeta*lnmeta)/(12.0*x*eta)
    elseif eta < -2.0
        coeff2 = @horner(eta, -1.72847633523e-2, -1.59372646475e-2, -4.64910887221e-3, -6.06834887760e-4, -6.14830384279e-6) / @horner(eta, 1.00000000000e+0, 7.64050615669e-1, 2.97143406325e-1, 5.79490176079e-2, 5.74558524851e-3)
    elseif eta < 2.0
        coeff2 = @horner(eta, -1.72839517431e-2, -1.46362417966e-2, -3.57406772616e-3, -3.91032032692e-4, 2.49634036069e-6) / @horner(eta, 1.00000000000e+0, 6.90560400696e-1, 2.49962384741e-1, 4.43843438769e-2, 4.24073217211e-3)
    elseif eta < 1000.0
        x = 1.0/eta
        coeff2 = @horner(x, 9.99944669480e-1, 1.04649839762e+2, 8.57204033806e+2, 7.31901559577e+2, 4.55174411671e+1) / @horner(x, 1.00000000000e+0, 1.04526456943e+2, 8.23313447808e+2, 3.11993802124e+3, 3.97003311219e+3)
    else
        coeff2 = -1.0/(12.0*eta)
    end
    return coeff2
end
"""
Computing the third and last coefficient for the expansion :
```math
\\epsilon (\\eta_{0},a) = \\epsilon_{1} (\\eta_{0},a)/a + \\epsilon_{2} (\\eta_{0},a)/a^{2} + \\epsilon_{3} (\\eta_{0},a)/a^{3}
```
Refer Eqn (3.12) in the paper
"""
function coeff3(eta::Float64)
    if eta < -8.0
        x=eta*eta
        y = log(-eta)/eta
        coeff3=(-30.0+eta*y*(6.0*x*y*y-12.0+x))/(12.0*eta*x*x)
    elseif eta < -4.0
        coeff3 = (@horner(eta, 4.95346498136e-2, 2.99521337141e-2, 6.88296911516e-3, 5.12634846317e-4, -2.01411722031e-5) / @horner(eta, 1.00000000000e+0, 7.59803615283e-1, 2.61547111595e-1, 4.64854522477e-2, 4.03751193496e-3))/(eta*eta)
    elseif eta < -2.0
        coeff3 = @horner(eta, 4.52313583942e-3, 1.20744920113e-3, -7.89724156582e-5, -5.04476066942e-5, -5.35770949796e-6) / @horner(eta, 1.00000000000e+0, 9.12203410349e-1, 4.05368773071e-1, 9.01638932349e-2, 9.48935714996e-3)
    elseif eta < 2.0
        coeff3 = @horner(eta, 4.39937562904e-3, 4.87225670639e-4, -1.28470657374e-4, 5.29110969589e-6, 1.57166771750e-7) / @horner(eta, 1.00000000000e+0, 7.94435257415e-1, 3.33094721709e-1, 7.03527806143e-2, 8.06110846078e-3)
    elseif eta < 10.0
        x=1.0/eta
        coeff3 = (@horner(x, -1.14811912320e-3, -1.12850923276e-1, 1.51623048511e+0, -2.18472031183e-1, 7.30002451555e-2) / @horner(x, 1.00000000000e+0, 1.42482206905e+1, 6.97360396285e+1, 2.18938950816e+2, 2.77067027185e+2))/(eta*eta)
    elseif eta < 100.0
        x=1.0/eta
        coeff3 = (@horner(x, -1.45727889667e-4, -2.90806748131e-1, -1.33085045450e+1, 1.99722374056e+2, -1.14311378756e+1) / @horner(x, 1.00000000000e+0, 1.39612587808e+2, 2.18901116348e+3, 7.11524019009e+3, 4.55746081453e+4))/(eta*eta)
    else
        eta3 = eta*eta*eta
        coeff3 = -log(eta)/(12.0*eta3)
    end
    return coeff3
end

"""
    gamma_inc_cf(a, x, ind)

Computes ``P(a,x)`` by continued fraction expansion given by :
```math
R(a,x) * \\frac{1}{1-\\frac{z}{a+1+\\frac{z}{a+2-\\frac{(a+1)z}{a+3+\\frac{2z}{a+4-\\frac{(a+2)z}{a+5+\\frac{3z}{a+6-\\dots}}}}}}}
```
Used when `1 <= a <= BIG` and `x < x0`.

External links: [DLMF](https://dlmf.nist.gov/8.9.2)

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_cf(a::Float64, x::Float64, ind::Integer)
    acc = acc0[ind + 1]
    tol = 4.0*acc
    a2nm1 = 1.0
    a2n = 1.0
    b2nm1 = x
    b2n = x + (1.0 - a)
    c = 1.0
    while true
       a2nm1 = x*a2n + c*a2nm1
       b2nm1 = x*b2n + c*b2nm1
       c = c + 1.0
       t = c - a
       a2n = a2nm1 + t*a2n
       b2n = b2nm1 + t*b2n
       a2nm1 = a2nm1/b2n
       b2nm1 = b2nm1/b2n
       a2n = a2n/b2n
       b2n = 1.0
       if abs(a2n - a2nm1/b2nm1) < tol*a2n
           break
       end
    end
    q = rgammax(a,x)*a2n
    return (1.0 - q, q)
end
"""
    gamma_inc_taylor(a, x, ind)

Compute ``P(a,x)`` using Taylor Series for P/R given by :
```math
R(a,x)/a * (1 + \\sum_{n=1}^{\\infty} x^{n}/((a+1)(a+2)...(a+n)))
```
Used when `1 <= a <= BIG` and `x <= max{a, ln 10}`.

External links: [DLMF](https://dlmf.nist.gov/8.11.2)

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_taylor(a::Float64, x::Float64, ind::Integer)
    acc = acc0[ind + 1]
    wk = zeros(30)
    flag = false
    apn = a + 1.0
    t = x/apn
    wk[1] = t
    loop=2
    for indx = 2:20
       apn = apn + 1.0
       t = t*(x/apn)
       if t <= 1.0e-3
           loop = indx
           flag = true
           break
       end
       wk[indx] = t
    end
    if !flag
        loop = 20
    end
    sm = t
    tol = 0.5*acc #tolerance
    while true
       apn = apn+1.0
       t = t*(x/apn)
       sm = sm + t
       if t <= tol
           break
       end
    end
    for j = loop-1:-1:1
       sm += wk[j]
    end
    p = (rgammax(a,x)/a)*(1.0 + sm)
    return (p, 1.0-p)
end
"""
    gamma_inc_asym(a, x, ind)

Compute ``P(a,x)`` using asymptotic expansion given by:
```math
R(a,x)/a * (1 + \\sum_{n=1}^{N-1}(a_{n}/x^{n} + \\Theta _{n}a_{n}/x^{n}))
```
where `R(a,x) = rgammax(a,x)`. Used when `1 <= a <= BIG` and `x >= x0`.

External links: [DLMF](https://dlmf.nist.gov/8.11.2)

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_asym(a::Float64, x::Float64, ind::Integer)
    wk = zeros(30)
    flag = false
    acc = acc0[ind + 1]
    amn = a-1.0
    t=amn/x
    wk[1]=t
    loop=2
    for indx = 2 : 20
       amn = amn-1.0
       t=t*(amn/x)
       if abs(t) <= 1.0e-3
           loop = indx
           flag = true
           break
       end
       wk[indx]=t
    end
    if !flag
        loop = 20
    end
    sm = t
    while true
       if abs(t) < acc
           break
       end
       amn=amn-1.0
       t=t*(amn/x)
       sm=sm+t
    end
    for j = loop-1:-1:1
       sm += wk[j]
    end
    q = (rgammax(a,x)/x)*(1.0 + sm)
    return (1.0 - q, q)
end
"""
    gamma_inc_taylor_x(a,x,ind)

Computes ``P(a,x)`` based on Taylor expansion of ``P(a,x)/x**a`` given by:
```math
J = -a * \\sum_{1}^{\\infty} (-x)^{n}/((a+n)n!)
```
and ``P(a,x)/x**a`` is given by :
```math
(1 - J)/ \\Gamma(a+1)
```
resulting from term-by-term integration of `gamma_inc(a,x,ind)`.
This is used when `a < 1` and `x < 1.1` - Refer Eqn (9) in the paper.

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_taylor_x(a::Float64, x::Float64, ind::Integer)
    acc = acc0[ind + 1]
    l=3.0
    c=x
    sm= x/(a + 3.0)
    tol = 3.0*acc/(a + 1.0)
    while true
       l=l+1.0
       c=-c*(x/l)
       t=c/(a+l)
       sm=sm+t
       if abs(t) <= tol
           break
       end
    end
    temp = a*x*((sm/6.0 - 0.5/(a + 2.0))*x + 1.0/(a + 1.0))
    z = a*log(x)
    h = rgamma1pm1(a)
    g = 1.0 + h
    if (x < 0.25 && z > -.13394) || a < x/2.59
       l = expm1(z)
       w = 1.0+l
       q = max((w*temp - l)*g - h, 0.0)
       return (1.0 - q, q)
    else
       w = exp(z)
       p = w*g*(1.0 - temp)
       return (p, 1.0 - p)
    end
end
"""
    gamma_inc_minimax(a,x,z)

Compute ``P(a,x)`` using minimax approximations given by :
```math
1/2 * erfc(\\sqrt{y}) - e^{-y}/\\sqrt{2\\pi*a}* T(a,\\lambda)
``` where
```math
T(a,\\lambda) = \\sum_{0}^{N} c_{k}(z)a^{-k}
```

This is a higher accuracy approximation of Temme expansion, which deals with the region near `a ≈ x` with `a` large.
Refer Appendix F in the paper for the extensive set of coefficients calculated using Brent's multiple precision arithmetic(set at 50 digits) in BRENT, R. P. A FORTRAN multiple-precision arithmetic package, ACM Trans. Math. Softw. 4(1978), 57-70 .

External links: [DLMF](https://dlmf.nist.gov/8.12.8)

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_minimax(a::Float64, x::Float64, z::Float64)
    l = x/a
    s = 1.0 - l
    y = -a*logmxp1(l)
    c = exp(-y)
    w = 0.5*erfcx(sqrt(y))

    if abs(s) <= 1.0e-3
        c0 = @horner(z , d00 , d0[1] , d0[2] , d0[3])
        c1 = @horner(z , d10 , d1[1] , d1[2] , d1[3])
        c2 = @horner(z , d20 , d2[1] , d2[2])
        c3 = @horner(z , d30 , d3[1] , d3[2])
        c4 = @horner(z , d40 , d4[1])
        c5 = @horner(z , d50 , d5[1])
        c6 = @horner(z , d60 , d6[1])

        t = @horner(z , c0 , c1 , c2 , c3 , c4 , c5 , c6 , d70 , d80)
        if l < 1.0
            p = c*(w - rt2pin*t/sqrt(a))
            return (p , 1.0 - p)
        else
            q = c*(w + rt2pin*t/sqrt(a))
            return (1.0 - q , q)
        end
    end
    #---USING THE MINIMAX APPROXIMATIONS---
    c0 = @horner(z , -.333333333333333E+00 , -.159840143443990E+00 , -.335378520024220E-01 , -.231272501940775E-02)/(@horner(z , 1.0 , .729520430331981E+00 , .238549219145773E+00, .376245718289389E-01 , .239521354917408E-02 , -.939001940478355E-05 , .633763414209504E-06) )
    c1 = @horner(z , -.185185185184291E-02 , -.491687131726920E-02 , -.587926036018402E-03 , -.398783924370770E-05)/(@horner(z , 1.0 , .780110511677243E+00 , .283344278023803E+00 , .506042559238939E-01 , .386325038602125E-02))
    c2 = @horner(z , .413359788442192E-02 , .669564126155663E-03)/(@horner(z , 1.0 , .810647620703045E+00 , .339173452092224E+00 , .682034997401259E-01 , .650837693041777E-02 , -.421924263980656E-03))
    c3 = @horner(z , .649434157619770E-03 , .810586158563431E-03)/(@horner(z , 1.0 , .894800593794972E+00, .406288930253881E+00 , .906610359762969E-01 , .905375887385478E-02 , -.632276587352120E-03))
    c4 = @horner(z , -.861888301199388E-03 , -.105014537920131E-03)/(@horner(z , 1.0 , .103151890792185E+01 , .591353097931237E+00 , .178295773562970E+00 , .322609381345173E-01))
    c5 = @horner(z , -.336806989710598E-03, -.435211415445014E-03)/(@horner(z , 1.0 , .108515217314415E+01 , .600380376956324E+00 , .178716720452422E+00))
    c6 = @horner(z , .531279816209452E-03 , -.182503596367782E-03)/(@horner(z , 1.0 , .770341682526774E+00 , .345608222411837E+00))
    c7 = @horner(z , .344430064306926E-03, .443219646726422E-03)/(@horner(z , 1.0 , .115029088777769E+01 , .821824741357866E+00))
    c8 = @horner(z , -.686013280418038E-03 , .878371203603888E-03)

    t = @horner(1.0/a , c0 , c1 , c2 , c3 , c4 , c5 , c6 , c7 , c8)
    if l < 1.0
        p = c*(w - rt2pin*t/sqrt(a))
        return (p , 1.0 - p)
    else
        q = c*(w + rt2pin*t/sqrt(a))
        return (1.0 - q , q)
    end
end
"""
    gamma_inc_temme(a, x, z)

Compute ``P(a,x)`` using Temme's expansion given by :
```math
1/2 * erfc(\\sqrt{y}) - e^{-y}/\\sqrt{2\\pi*a}* T(a,\\lambda)
``` where
```math
T(a,\\lambda) = \\sum_{0}^{N} c_{k}(z)a^{-k}
```
This mainly solves the problem near the region when `a ≈ x` with a large, and is a lower accuracy version of the minimax approximation.

External links: [DLMF](https://dlmf.nist.gov/8.12.8)

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_temme(a::Float64, x::Float64, z::Float64)
    l = x/a
    y = -a*logmxp1(x/a)
    c = exp(-y)
    w = 0.5*erfcx(sqrt(y))
    c0 = @horner(z , d00 , d0[1] , d0[2] , d0[3] , d0[4] , d0[5] , d0[6])
    c1 = @horner(z , d10 , d1[1] , d1[2] , d1[3] , d1[4])
    c2 = @horner(z , d20 , d2[1])
    t = @horner(1.0/a , c0 , c1 , c2)
    if l < 1.0
        p = c*(w - rt2pin*t/sqrt(a))
        return (p , 1.0 - p)
    else
        q = c*(w + rt2pin*t/sqrt(a))
        return (1.0 - q , q)
    end
end
"""
    gamma_inc_temme_1(a, x, z, ind)

Computes ``P(a,x)`` using simplified Temme expansion near ``y=0`` by :
```math
E(y) - (1 - y)/\\sqrt{2\\pi*a} * T(a,\\lambda)
```
where
```math
E(y) = 1/2 - (1 - y/3)*(\\sqrt(y/\\pi))
```
Used instead of it's previous function when ``\\sigma <= e_{0}/\\sqrt{a}``.

External links: [DLMF](https://dlmf.nist.gov/8.12.8)
"""
function gamma_inc_temme_1(a::Float64, x::Float64, z::Float64, ind::Integer)
    iop = ind + 1
    l = x/a
    y = -a * logmxp1(l)
    if a*eps()*eps() > 3.28e-3
        throw(DomainError((a,x,ind,"P(a,x) or Q(a,x) is computationally indeterminant in this case.")))
    end
    c = exp(-y)
    w = 0.5*erfcx(sqrt(y))
    u = 1.0/a
    if iop == 1
        c0 = @horner(z , d00 , d0[1] , d0[2] , d0[3])
        c1 = @horner(z , d10 , d1[1] , d1[2] , d1[3])
        c2 = @horner(z , d20 , d2[1] , d2[2])
        c3 = @horner(z , d30 , d3[1] , d3[2])
        c4 = @horner(z , d40 , d4[1])
        c5 = @horner(z , d50 , d5[1])
        c6 = @horner(z , d60 , d6[1])

        t = @horner(u , c0 , c1 , c2 , c3 , c4 , c5 , c6 , d70 , d80)

    elseif iop == 2
        c0 = @horner(d00 , d0[1] , d0[2])
        c1 = @horner(d10 , d1[1])
        t = @horner(u , c0 , c1 , d20)

    else
        t = @horner(z , d00 , d0[1])

    end
    if l < 1.0
        p = c*(w - rt2pin*t/sqrt(a))
        return (p , 1.0 - p)
    else
        q = c*(w + rt2pin*t/sqrt(a))
        return (1.0 - q , q)
    end
end
"""
    gamma_inc_fsum(a,x)

Compute using Finite Sums for ``Q(a,x)`` when `a >= 1 && 2a` is integer.
Used when `a <= x <= x0` and `a = n/2`.
Refer Eqn (14) in the paper.

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc)
"""
function gamma_inc_fsum(a::Float64, x::Float64)
    if isinteger(a)
        sm = exp(-x)
        t = sm
        N = 1
        c=0.0
        i = trunc(Int,a )
    else
        rtx = sqrt(x)
        sm = erfc(rtx)
        t = exp(-x)/(rtpi*rtx)
        N=0
        c=-0.5
        i = trunc(Int,a )
    end
    for lp = N : i-1
        if i == 0
            break
        end
        c += 1.0
        t = (x*t)/c
        sm += t
    end
    q = sm
    return (1.0 - q, q)

end
"""
    gamma_inc_inv_psmall(a,p)

Compute `x0` - initial approximation when `p` is small.
Here we invert the series in Eqn (2.20) in the paper and write the inversion problem as:
```math
x = r\\left[1 + a\\sum_{k=1}^{\\infty}\\frac{(-x)^{n}}{(a+n)n!}\\right]^{-1/a},
```
where ``r = (p\\Gamma(1+a))^{1/a}``
Inverting this relation we obtain ``x = r + \\sum_{k=2}^{\\infty}c_{k}r^{k}``.
"""
function gamma_inc_inv_psmall(a::Float64, p::Float64)
    logr = (1.0/a)*(log(p) + logabsgamma(a + 1.0)[1])
    r = exp(logr)
    ap1 = a + 1.0
    ap1² = (ap1)*ap1
    ap1³ = (ap1)*ap1²
    ap1⁴ = ap1²*ap1²
    ap2 = a+2.0
    ap2² = ap2*ap2
    ck1 = 1.0
    ck2 = 1.0/(1.0+a)
    ck3 = 0.5*(3*a+5)/(ap1²*(a+2))
    ck4 = (1.0/3.0)*(@horner(a,31,33,8))/(ap1³*ap2*(a+3))
    ck5 = (1.0/24.0)*(@horner(a, 2888, 5661, 3971, 1179, 125))/(ap1⁴*ap2²*(a+3)*(a+4))
    x0 = @horner(r, 0.0, ck1, ck2, ck3, ck4, ck5)
    return x0
end
"""
    gamma_inc_inv_qsmall(a,q)

Compute `x0` - initial approximation when `q` is small from ``e^{-x_{0}} x_{0}^{a} = q \\Gamma(a)``.
Asymptotic expansions Eqn (2.29) in the paper is used here and higher approximations are obtained using
```math
x \\sim x_{0} - L + b \\sum_{k=1}^{\\infty} d_{k}/x_{0}^{k}
```
where ``b = 1-a``, ``L = \\ln{x_0}``.
"""
function gamma_inc_inv_qsmall(a::Float64, q::Float64)
    b=1.0-a
    eta=sqrt(-2/a*log(q*gammax(a)*sqrt(2*pi)/sqrt(a)))
    x0 = a*lambdaeta(eta)
    l = log(x0)

    if a > 0.12 || x0 > 5
        r = 1.0/x0
        ck1 = l - 1.0
        ck2 = (@horner(l, @horner(b, 2, 3), @horner(b, -2, -2), 1))/2.0
        ck3 = (@horner(l, @horner(b, -12, -24, -11), @horner(b, 12, 24, 6), @horner(b, -6, -9), 2))/6.0
        ck4 = (@horner(l, @horner(b, 72, 162, 120, 25), @horner(b, -72, -168, -114, -12), @horner(b, 36, 84, 36), @horner(b, -12, -22), 3))/12.0
        x0 = x0 - l + b * r * @horner(r, ck1, ck2, ck3, ck4)
    else
        r = 1.0/x0
        l² = l*l
        ck1 = l - 1.0
        x0 = x0 - l + b * r * ck1
    end
    return x0
end
"""
    gamma_inc_inv_asmall(a,p,q,pcase)

    Compute `x0` - initial approximation when `a` is small.
    Here the solution `x` of ``P(a,x)=p`` satisfies ``x_{l} < x < x_{u}``
    where ``x_{l} = (p\\Gamma(a+1))^{1/a}`` and ``x_{u} = -\\log{(1 - p\\Gamma(a+1))}``, and is used as starting value for Newton iteration.
    """
function gamma_inc_inv_asmall(a::Float64, p::Float64, q::Float64, pcase::Bool)
    logp = (pcase) ? log(p) : log1p(-q)
    return exp((1.0/a)*(logp +loggamma1p(a)))
end
"""
    gamma_inc_inv_alarge(a,porq,s)

Compute `x0` - initial approximation when `a` is large.
The inversion problem is rewritten as :
```math
0.5 \\operatorname{erfc}(\\eta \\sqrt{a/2}) + R_{a}(\\eta) = q
```
For large values of `a` we can write: ``\\eta(q,a) = \\eta_{0}(q,a) + \\epsilon(\\eta_{0},a)``
and it is possible to expand:
```math
\\epsilon(\\eta_{0},a) = \\epsilon_{1}(\\eta_{0},a)/a + \\epsilon_{2}(\\eta_{0},a)/a^{2} + \\epsilon_{3}(\\eta_{0},a)/a^{3} + ...
```
which is calculated by coeff1, coeff2 and coeff3 functions below.
This returns a tuple `(x0,fp)`, where `fp` is computed since it's an approximation for the coefficient after inverting the original power series.
"""
function gamma_inc_inv_alarge(a::Float64, porq::Float64, s::Integer)
    r = erfcinv(2*porq)
    eta = s*r/sqrt(a*0.5)
    eta += (coeff1(eta) + (coeff2(eta) + coeff3(eta)/a)/a)/a
    x0 = a*lambdaeta(eta)
    fp = -sqrt(a/(2*pi))*exp(-0.5*a*eta*eta)/gammax(a)
    return (x0,fp)
end
# Reference : 'Computation of the incomplete gamma function ratios and their inverse' by Armido R DiDonato , Alfred H Morris.
# Published in Journal: ACM Transactions on Mathematical Software (TOMS)
# Volume 12 Issue 4, Dec. 1986 Pages 377-393
# doi>10.1145/22721.23109

"""
    gamma_inc(a,x,IND)

Returns a tuple ``(p, q)`` where ``p + q = 1``, and
``p=P(a,x)`` is the Incomplete gamma function ratio given by:
```math
P(a,x)=\\frac{1}{\\Gamma (a)} \\int_{0}^{x} e^{-t}t^{a-1} dt.
```
and ``q=Q(a,x)`` is the Incomplete gamma function ratio given by:
```math
Q(x,a)=\\frac{1}{\\Gamma (a)} \\int_{x}^{\\infty} e^{-t}t^{a-1} dt.
```

`IND ∈ [0,1,2]` sets accuracy: `IND=0` means 14 significant digits accuracy, `IND=1` means 6 significant digit, and `IND=2` means only 3 digit accuracy.

External links: [DLMF](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)

See also [`gamma(z)`](@ref SpecialFunctions.gamma), [`gamma_inc_inv(a,p,q)`](@ref SpecialFunctions.gamma_inc_inv)
"""
function gamma_inc(a::Float64,x::Float64,ind::Integer)
    iop = ind + 1
    acc = acc0[iop]
    if a<0.0 || x<0.0
        throw(DomainError((a,x,ind,"`a` and `x` must be greater than 0 ---- Domain : (0,inf)")))
    elseif a==0.0 && x==0.0
        throw(DomainError((a,x,ind,"`a` and `x` must be greater than 0 ---- Domain : (0,inf)")))
    elseif a*x==0.0
        if x<=a
            return (0.0,1.0)
        else
            return (1.0,0.0)
        end
    end

    if a >= 1.0
        if a >= big1[iop]
            l = x/a
            if l == 0.0
                return (0.0,1.0)
            end
            s = 1.0 - l
            z = -logmxp1(l)
            if z >= 700.0/a
                return (0.0,1.0)
            end
            y = a*z
            rta = sqrt(a)
            if abs(s) <= e0[iop]/rta
                z = sqrt(z + z)
                if l < 1.0
                    z=-z
                end
                return gamma_inc_temme_1(a, x, z, ind)
            end

            if abs(s) <= 0.4
                if abs(s) <= 2.0*eps() && a*eps()*eps() > 3.28e-3
                    throw(DomainError((a,x,ind,"P(a,x) or Q(a,x) is computationally indeterminant in this case.")))
                end
                c = exp(-y)
                w = 0.5*erfcx(sqrt(y))
                u = 1.0/a
                z = sqrt(z + z)
                if l < 1.0
                    z=-z
                end
                if iop == 1
                    return gamma_inc_minimax(a,x,z)
                elseif iop == 2
                    return gamma_inc_temme(a,x,z)
                else
                    t = @horner(z , d00 , d0[1] , d0[2] , d0[3])
                    return gamma_inc_temme_1(a, x, z, ind)
                end
            end
        elseif a > x || x >= x0[iop] || !isinteger(2*a)
            r = rgammax(a,x)
            if r == 0.0
                if x <= a
                    return (0.0,1.0)
                else
                    return (1.0,0.0)
                end
            end
            if x <= max(a,alog10)
                return gamma_inc_taylor(a, x, ind)
            elseif x < x0[iop]
                return gamma_inc_cf(a, x, ind)
            else
                return gamma_inc_asym(a, x, ind)
            end
        else
            return gamma_inc_fsum(a,x)

        end
    elseif a == 0.5
        if x >= 0.25
            q = erfc(sqrt(x))
            return ( 1.0 - q , q )
        end
        p = erf(sqrt(x))
        return ( p , 1.0 - p )
    elseif x < 1.1
        return gamma_inc_taylor_x(a, x, ind)
    end
    r = rgammax(a,x)
    if r == 0.0
        return (1.0, 0.0)
    else
        return gamma_inc_cf(a, x, ind)
    end

end

function gamma_inc(a::BigFloat,x::BigFloat,ind::Integer) #BigFloat version from GNU MPFR wrapped via ccall
    z = BigFloat()
    ccall((:mpfr_gamma_inc, :libmpfr), Int32 , (Ref{BigFloat} , Ref{BigFloat} , Ref{BigFloat} , Int32) , z , a , x , ROUNDING_MODE[])
    q = z/gamma(a)
    return (1.0 - q, q)
end
gamma_inc(a::Float32,x::Float32,ind::Integer) = ( Float32(gamma_inc(Float64(a),Float64(x),ind)[1]) , Float32(gamma_inc(Float64(a),Float64(x),ind)[2]) )
gamma_inc(a::Float16,x::Float16,ind::Integer) = ( Float16(gamma_inc(Float64(a),Float64(x),ind)[1]) , Float16(gamma_inc(Float64(a),Float64(x),ind)[2]) )
gamma_inc(a::Real,x::Real,ind::Integer) = (gamma_inc(float(a),float(x),ind))
gamma_inc(a::Integer,x::Integer,ind::Integer) = gamma_inc(Float64(a),Float64(x),ind)
gamma_inc(a::AbstractFloat,x::AbstractFloat,ind::Integer) = throw(MethodError(gamma_inc,(a,x,ind,"")))

#EFFICIENT AND ACCURATE ALGORITHMS FOR THECOMPUTATION AND INVERSION OF THE INCOMPLETEGAMMA FUNCTION RATIOS by Amparo Gil, Javier Segura, Nico M. Temme
#SIAM Journal on Scientific Computing 34(6) (2012), A2965-A2981
# arXiv:1306.1754

"""
    gamma_inc_inv(a,p,q)

Inverts the `gamma_inc(a,x)` function, by computing `x` given `a`,`p`,`q` in ``P(a,x)=p`` and ``Q(a,x)=q``.

External links: [DLMF](https://dlmf.nist.gov/8.2.4), [Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)

See also: [`gamma_inc(a,x,ind)`](@ref SpecialFunctions.gamma_inc).
"""
function gamma_inc_inv(a::Float64, p::Float64, q::Float64)
    if p < 0.5
        pcase = true
        porq = p
        s = -1
    else
        pcase = false
        porq = q
        s = 1
    end
    haseta = false

    logr = (1.0/a)*( log(p) + logabsgamma(a + 1.0)[1] )
    if logr < log(0.2*(1+a)) #small value of p
        x0 = gamma_inc_inv_psmall(a,p)
    elseif ((q < min(0.02, exp(-1.5*a)/gamma(a))) && (a < 10)) #small q
        x0 = gamma_inc_inv_qsmall(a, q)
    elseif abs(porq - 0.5) < 1.0e-05
        x0 = a - 1.0/3.0 + (8.0/405.0 + 184.0/25515.0/a) / a
    elseif abs(a-1.0) < 1.0e-4
        x0 = pcase ? -log1p(-p) : -log(q)
    elseif a < 1.0 # small value of a
        x0 = gamma_inc_inv_asmall(a, p, q, pcase)
    else    #large a
        haseta = true
        (x0,fp) = gamma_inc_inv_alarge(a,porq,s)
    end

    t = 1
    x = x0
    n = 1
    logabsgam = logabsgamma(a)[1]
    #Newton-like higher order iteration with upper limit as 15.
    while t > 1.0e-15 && n < 15
        x = x0
        x² = x*x
        if !haseta
            dlnr = (1.0-a)*log(x) + x + logabsgam
            if dlnr > log(floatmax(Float64)/1000.0)
                n=20
            else
                r = exp(dlnr)
            end
        else
            r = -(1/fp)*x
        end

        (px,qx) = gamma_inc(a,x,0)
        ck1 = (pcase) ? -r*(px-p) : r*(qx-q)
        ck2 = (x-a+1.0)/(2.0*x)
        ck3 = (@horner(x, @horner(a, 1, -3, 2), @horner(a, 4, -4), 2))/(6*x²)
        r = ck1
        if a > 0.1
            x0 = @horner(r, x, 1.0, ck2, ck3)
        else
            if a > 0.05
                x0 = @horner(r, x, 1.0, ck2)
            else
                x0 = x + r
            end
        end

        t=abs(x/x0 - 1.0)
        n+=1
        x=x0

    end
    return x
end

gamma_inc_inv(a::Float32, p::Float32, q::Float32) = Float32( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )
gamma_inc_inv(a::Float16, p::Float16, q::Float16) = Float16( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )
gamma_inc_inv(a::Real, p::Real, q::Real) = ( gamma_inc_inv(float(a), float(p), float(q)) )
gamma_inc_inv(a::Integer, p::Integer, q::Integer) = ( gamma_inc_inv(Float64(a), Float64(p), Float64(q)) )
gamma_inc_inv(a::AbstractFloat, p::AbstractFloat, q::AbstractFloat) = throw(MethodError(gamma_inc_inv,(a,p,q,"")))


text:
module LineEdit
    using Terminals

    import Terminals: raw!, width, height, cmove, Rect, Size, getX,
                      getY, clear_line, beep

    import Base: ensureroom, peek

    abstract TextInterface

    export run_interface, Prompt, ModalInterface, transition, reset_state, edit_insert

    immutable ModalInterface <: TextInterface
        modes
    end

    type MIState
        interface::ModalInterface
        current_mode
        aborted::Bool
        mode_state
    end

    type Mode <: TextInterface

    end

    type Prompt <: TextInterface
        prompt
        first_prompt
        prompt_color::ASCIIString
        keymap_func
        keymap_func_data
        input_color
        complete
        on_enter
        on_done
        hist
    end

    immutable InputAreaState
        num_rows::Int64
        curs_row::Int64
    end

    type PromptState
        terminal::TextTerminal
        p::Prompt
        input_buffer::IOBuffer
        ias::InputAreaState
        indent::Int
    end

    input_string(s::PromptState) = bytestring(pointer(s.input_buffer.data),s.input_buffer.size)

    abstract HistoryProvider
    abstract CompletionProvider

    type EmptyCompletionProvider <: CompletionProvider
    end

    type EmptyHistoryProvider <: HistoryProvider
    end

    reset_state(::EmptyHistoryProvider) = nothing

    completeLine(c::EmptyCompletionProvider,s) = []

    terminal(s::IO) = s
    terminal(s::PromptState) = s.terminal

    for f in [:terminal,:edit_insert,:on_enter,:add_history,:buffer,:edit_backspace,:(Base.isempty),
            :replace_line,:refreshMultiLine,:input_string,:completeLine,:edit_move_left,:edit_move_right,
            :update_display_buffer]
        @eval ($f)(s::MIState,args...) = $(f)(s.mode_state[s.current_mode],args...)
    end

    function common_prefix(completions)
        ret = ""
        i = nexti = 1
        cc,nexti = next(completions[1],1)
        while true
            for c in completions
                if i > length(c) || c[i] != cc
                    return ret
                end
            end
            ret = ret*string(cc)
            if i >= length(completions[1])
                return ret
            end
            i = nexti
            cc,nexti = next(completions[1],i)
        end
    end

    function completeLine(s::PromptState)
        (completions,partial) = completeLine(s.p.complete,s)
        if length(completions) == 0
            beep(LineEdit.terminal(s))
        elseif length(completions) == 1
            # Replace word by completion
            prev_pos = position(s.input_buffer)
            seek(s.input_buffer,prev_pos-sizeof(partial))
            edit_replace(s,position(s.input_buffer),prev_pos,completions[1])
        else
            p = common_prefix(completions)
            if length(p) > 0 && p != partial
                # All possible completions share the same prefix, so we might as
                # well complete that
                prev_pos = position(s.input_buffer)
                seek(s.input_buffer,prev_pos-sizeof(partial))
                edit_replace(s,position(s.input_buffer),prev_pos,p)
            else
                # Show available completions
                colmax = maximum(map(length,completions))
                num_cols = div(width(LineEdit.terminal(s)),colmax+2)
                entries_per_col = div(length(completions),num_cols)+1
                println(LineEdit.terminal(s))
                for row = 1:entries_per_col
                    for col = 0:num_cols
                        idx = row + col*entries_per_col
                        if idx <= length(completions)
                            cmove_col(LineEdit.terminal(s),(colmax+2)*col)
                            print(LineEdit.terminal(s),completions[idx])
                        end
                    end
                    println(LineEdit.terminal(s))
                end
            end
        end
    end

    clear_input_area(s) = (clear_input_area(s.terminal,s.ias); s.ias = InputAreaState(0,0))
    function clear_input_area(terminal,state::InputAreaState)
        #println(s.curs_row)
        #println(s.num_rows)

        # Go to the last line
        if state.curs_row < state.num_rows
            cmove_down(terminal,state.num_rows-state.curs_row)
        end

        # Clear lines one by one going up
        for j=0:(state.num_rows - 2)
            clear_line(terminal)
            cmove_up(terminal)
        end

        # Clear top line
        clear_line(terminal)
    end

    prompt_string(s::PromptState) = s.p.prompt
    prompt_string(s::String) = s

    refreshMultiLine(s::PromptState) = s.ias = refreshMultiLine(s.terminal,buffer(s),s.ias,s,indent=s.indent)

    function refreshMultiLine(terminal,buf,state::InputAreaState,prompt = "";indent = 0)
        cols = width(terminal)

        clear_input_area(terminal,state)

        curs_row = -1 #relative to prompt
        curs_col = -1 #absolute
        curs_pos = -1 # 1 - based column position of the cursor
        cur_row = 0
        buf_pos = position(buf)
        line_pos = buf_pos
        # Write out the prompt string
        write_prompt(terminal,prompt)
        prompt = prompt_string(prompt)

        seek(buf,0)

        llength = 0

        l=""

        plength = length(prompt)
        pslength = length(prompt.data)
        # Now go through the buffer line by line
        while cur_row == 0 || (!isempty(l) && l[end] == '\n')
            l = readline(buf)
            cur_row += 1
            # We need to deal with UTF8 characters. Since the IOBuffer is a bytearray, we just count bytes
            llength = length(l)
            slength = length(l.data)
            if cur_row == 1 #First line
                if line_pos < slength
                    num_chars = length(l[1:line_pos])
                    curs_row = div(plength+num_chars-1,cols)+1
                    curs_pos = (plength+num_chars-1)%cols+1
                end
                cur_row += div(plength+llength-1,cols)
                line_pos -= slength
                write(terminal,l)
            else
                # We expect to be line after the last valid output line (due to
                # the '\n' at the end of the previous line)
                if curs_row == -1
                    if line_pos < slength
                        num_chars = length(l[1:line_pos])
                        curs_row = cur_row+div(indent+num_chars-1,cols)
                        curs_pos = (indent+num_chars-1)%cols+1
                    end
                    line_pos -= slength #'\n' gets an extra pos
                    cur_row += div(llength+indent-1,cols)
                    cmove_col(terminal,indent+1)
                    write(terminal,l)
                    # There's an issue if the last character we wrote was at the very right end of the screen. In that case we need to
                    # emit a new line and move the cursor there.
                    if curs_pos == cols
                        write(terminal,"\n")
                        cmove_col(terminal,1)
                        curs_row+=1
                        curs_pos=0
                        cur_row+=1
                    end
                else
                    cur_row += div(llength+indent-1,cols)
                    cmove_col(terminal,indent+1)
                    write(terminal,l)
                end

            end

        end

        seek(buf,buf_pos)

        # If we are at the end of the buffer, we need to put the cursor one past the
        # last character we have written

        if curs_row == -1
            curs_pos = ((cur_row == 1 ? plength : indent)+llength-1)%cols+1
            curs_row = cur_row
        end

        # Same issue as above. TODO: We should figure out
        # how to refactor this to avoid duplcating functionality.
        if curs_pos == cols
            write(terminal,"\n")
            cmove_col(terminal,1)
            curs_row+=1
            curs_pos=0
            cur_row+=1
        end

        # Let's move the cursor to the right position
        # The line first
        n = cur_row-curs_row
        if n>0
            cmove_up(terminal,n)
        end

        #columns are 1 based
        cmove_col(terminal,curs_pos+1)

        flush(terminal)

        # Updated cur_row,curs_row
        return InputAreaState(cur_row,curs_row)
    end

    # Edit functionality

    char_move_left(s::PromptState) = char_move_left(s.input_buffer)
    function char_move_left(buf::IOBuffer)
        while position(buf)>0
            seek(buf,position(buf)-1)
            c = peek(buf)
            if ((c&0x80) == 0) || ((c&0xc0) == 0xc0)
                break
            end
        end
    end

    function edit_move_left(s::PromptState)
        if position(s.input_buffer)>0
            #move to the next UTF8 character to the left
            char_move_left(s.input_buffer)
            refresh_line(s)
        end
    end

    char_move_right(s) = char_move_right(buffer(s))
    function char_move_right(buf::IOBuffer)
        while position(buf) != buf.size
            seek(buf,position(buf)+1)
            if position(buf)==buf.size
                break
            end
            c = peek(buf)
            if ((c&0x80) == 0) || ((c&0xc0) == 0xc0)
                break
            end
        end
    end

    const non_word_chars = " \t\n\"\\'`@\$><=:;|&{}()[].,+-*/?%^~"

    function char_move_word_right(s)
        while !eof(s.input_buffer) && !in(read(s.input_buffer,Char),non_word_chars)
        end
    end

    function char_move_word_left(s)
        while position(s.input_buffer) > 0
            char_move_left(s)
            c = peek(s.input_buffer)
            if c < 0x80 && in(char(c),non_word_chars)
                read(s.input_buffer,Uint8)
                break
            end
        end
    end

    function edit_move_right(s)
        if position(s.input_buffer)!=s.input_buffer.size
            #move to the next UTF8 character to the right
            char_move_right(s)
            refresh_line(s)
        end
    end

    ## Move line up/down
    # Querying the terminal is expensive, memory access is cheap
    # so to find the current column, we find the offset for the start
    # of the line.
    function edit_move_up(s)
        buf = buffer(s)
        npos = rsearch(buf.data,'\n',position(buf))
        if npos == 0 #we're in the first line
            return false
        end
        # We're interested in character count, not byte count
        offset = length(bytestring(buf.data[(npos+1):(position(buf))]))-1
        npos2 = rsearch(buf.data,'\n',npos-1)
        seek(buf,npos2+1)
        for _ = 1:offset
            pos = position(buf)
            if read(buf,Char) == '\n'
                seek(buf,pos)
                break
            end
        end
        refresh_line(s)
        return true
    end

    function edit_move_down(s)
        buf = buffer(s)
        npos = rsearch(buf.data[1:buf.size],'\n',position(buf))
        # We're interested in character count, not byte count
        offset = length(bytestring(buf.data[(npos+1):(position(buf))]))-1
        npos2 = search(buf.data[1:buf.size],'\n',position(buf)+1)
        if npos2 == 0 #we're in the last line
            return false
        end
        seek(buf,npos2+1)
        for _ = 1:offset
            pos = position(buf)
            if eof(buf) || read(buf,Char) == '\n'
                seek(buf,pos)
                break
            end
        end
        refresh_line(s)
        return true
    end

    function charlen(ch::Char)
        if ch < 0x80
            return 1
        elseif ch < 0x800
            return 2
        elseif ch < 0x10000
            return 3
        elseif ch < 0x110000
            return 4
        end
        error("Corrupt UTF8")
    end

    function edit_replace(s,from,to,str)
        room = length(str.data)-(to-from)
        ensureroom(s.input_buffer,s.input_buffer.size + room)
        ccall(:memmove, Void, (Ptr{Void},Ptr{Void},Int), pointer(s.input_buffer.data,to+room+1),pointer(s.input_buffer.data,to+1),s.input_buffer.size-to)
        s.input_buffer.size += room
        seek(s.input_buffer,from)
        write(s.input_buffer,str)
    end

    function edit_insert(s::PromptState,c)
        str = string(c)
        edit_insert(s.input_buffer,str)
        if !('\n' in str) && eof(s.input_buffer) &&
            ((position(s.input_buffer) + length(s.p.prompt) + sizeof(str) - 1) < width(LineEdit.terminal(s)))
            #Avoid full update
            write(LineEdit.terminal(s),str)
        else
            refresh_line(s)
        end
    end

    # TODO: Don't use memmove
    function edit_insert(buf::IOBuffer,c)
        if eof(buf)
            write(buf,c)
        else
            s = string(c)
            ensureroom(buf,buf.size-position(buf)+sizeof(s))
            oldpos = position(buf)
            ccall(:memmove, Void, (Ptr{Void},Ptr{Void},Int), pointer(buf.data,position(buf)+1+sizeof(s)), pointer(buf.data,position(buf)+1),
                buf.size-position(buf))
            buf.size += sizeof(s)
            write(buf,c)
        end
    end

    function edit_backspace(s::PromptState)
        if edit_backspace(s.input_buffer)
            refresh_line(s)
        else
            beep(LineEdit.terminal(s))
        end

    end
    function edit_backspace(buf::IOBuffer)
        if position(buf) > 0 && buf.size>0
            oldpos = position(buf)
            char_move_left(buf)
            ccall(:memmove, Void, (Ptr{Void},Ptr{Void},Int), pointer(buf.data,position(buf)+1), pointer(buf.data,oldpos+1),
                buf.size-oldpos)
            buf.size -= oldpos-position(buf)
            return true
        else
            return false
        end
    end

    function edit_delete(s)
        buf = buffer(s)
        if buf.size>0 && position(buf) < buf.size
            oldpos = position(buf)
            char_move_right(s)
            ccall(:memmove, Void, (Ptr{Void},Ptr{Void},Int), pointer(buf.data,oldpos+1), pointer(buf.data,position(buf)+1),
                buf.size-position(buf))
            buf.size -= position(buf)-oldpos
            seek(buf,oldpos)
            refresh_line(s)
        else
            beep(LineEdit.terminal(s))
        end
    end

    function replace_line(s::PromptState,l::IOBuffer)
        s.input_buffer = l
    end

    function replace_line(s::PromptState,l)
        s.input_buffer.ptr = 1
        s.input_buffer.size = 0
        write(s.input_buffer,l)
    end

    history_prev(::EmptyHistoryProvider) = ("",false)
    history_next(::EmptyHistoryProvider) = ("",false)
    history_search(::EmptyHistoryProvider,args...) = false
    add_history(::EmptyHistoryProvider,s) = nothing
    add_history(s::PromptState) = add_history(mode(s).hist,s)

    function history_prev(s,hist)
        (l,ok) = history_prev(mode(s).hist)
        if ok
            replace_line(s,l)
            refresh_line(s)
        else
            beep(LineEdit.terminal(s))
        end
    end
    function history_next(s,hist)
        (l,ok) = history_next(mode(s).hist)
        if ok
            replace_line(s,l)
            refresh_line(s)
        else
            beep(LineEdit.terminal(s))
        end
    end

    refresh_line(s) = refreshMultiLine(s)

    default_completion_cb(::IOBuffer) = []
    default_enter_cb(_) = true

    write_prompt(terminal,s::PromptState) = write_prompt(terminal,s,s.p.prompt)
    function write_prompt(terminal,s::PromptState,prompt)
        @assert terminal == LineEdit.terminal(s)
        write(terminal,s.p.prompt_color)
        write(terminal,prompt)
        write(terminal,Base.text_colors[:normal])
        write(terminal,s.p.input_color)
    end
    write_prompt(terminal,s::ASCIIString) = write(terminal,s)

    function normalize_key(key)
        if isa(key,Char)
            return string(key)
        elseif isa(key,Integer)
            return string(char(key))
        elseif isa(key,String)
            if in('\0',key)
                error("Matching \\0 not currently supported.")
            end
            buf = IOBuffer()
            i = start(key)
            while !done(key,i)
                (c,i) = next(key,i)
                if c == '*'
                    write(buf,'\0')
                elseif c == '^'
                    (c,i) = next(key,i)
                    write(buf,uppercase(c)-64)
                elseif c == '\\'
                    (c,i) == next(key,i)
                    if c == 'C'
                        (c,i) == next(key,i)
                        @assert c == '-'
                        (c,i) == next(key,i)
                        write(buf,uppercase(c)-64)
                    elseif c == 'M'
                        (c,i) == next(key,i)
                        @assert c == '-'
                        (c,i) == next(key,i)
                        write(buf,'\e')
                        write(buf,c)
                    end
                else
                    write(buf,c)
                end
            end
            return takebuf_string(buf)
        end
    end

    # Turn an Dict{Any,Any} into a Dict{'Char',Any}
    # For now we use \0 to represent unknown chars so that they are sorted before everything else
    # If we ever actually want to mach \0 in input, this will have to be
    # reworked
    function normalize_keymap(keymap)
        ret = Dict{Char,Any}()
        for key in keys(keymap)
            newkey = normalize_key(key)
            current = ret
            i = start(newkey)
            while !done(newkey,i)
                (c,i) = next(newkey,i)
                if haskey(current,c)
                    if !isa(current[c],Dict)
                        println(ret)
                        error("Conflicting Definitions for keyseq "*escape_string(newkey)*" within one keymap")
                    end
                elseif done(newkey,i)
                    if isa(keymap[key],String)
                        current[c] = normalize_key(keymap[key])
                    else
                        current[c] = keymap[key]
                    end
                    break
                else
                    current[c] = Dict{Char,Any}()
                end
                current = current[c]
            end
        end
        ret
    end

    keymap_gen_body(keymaps,body::Expr,level) = body
    keymap_gen_body(keymaps,body::Function,level) = keymap_gen_body(keymaps,:($(body)(s)))
    keymap_gen_body(keymaps,body::Char,level) = keymap_gen_body(keymaps,keymaps[body])
    keymap_gen_body(keymaps,body::Nothing,level) = nothing
    function keymap_gen_body(keymaps,body::String,level)
        if length(body) == 1
            return keymap_gen_body(keymaps,body[1],level)
        end
        current = keymaps
        for c in body
            if haskey(current,c)
                if isa(current[c],Dict)
                    current = current[c]
                else
                    return keymap_gen_body(keymaps,current[c],level)
                end
            elseif haskey(current,'\0')
                return keymap_gen_body(keymaps,current['\0'],level)
            else
                error("No match for redirected key $body")
            end
        end
        error("No exact match for redirected key $body")
    end

    keymap_gen_body(a,b) = keymap_gen_body(a,b,1)
    function keymap_gen_body(dict,subdict::Dict,level)
        block = Expr(:block)
        bc = symbol("c"*string(level))
        push!(block.args,:($bc=read(LineEdit.terminal(s),Char)))

        if haskey(subdict,'\0')
            last_if = keymap_gen_body(dict,subdict['\0'],level+1)
        else
            last_if = nothing
        end

        for c in keys(subdict)
            if c == '\0'
                continue
            end
            cblock = Expr(:if,:($bc==$c))
            push!(cblock.args,keymap_gen_body(dict,subdict[c],level+1))
            if isa(cblock,Expr)
                push!(cblock.args,last_if)
            end
            last_if = cblock
        end

        push!(block.args,last_if)
        return block
    end

    export @keymap

    # deep merge where target has higher precedence
    function keymap_merge!(target::Dict,source::Dict)
        for k in keys(source)
            if !haskey(target,k)
                target[k] = source[k]
            elseif isa(target[k],Dict)
                keymap_merge!(target[k],source[k])
            else
                # Ignore, target has higher precedence
            end
        end
    end

    fixup_keymaps!(d,l,s,sk) = nothing
    function fixup_keymaps!(dict::Dict, level, s, subkeymap)
        if level > 1
            for d in dict
                fixup_keymaps!(d[2],level-1,s,subkeymap)
            end
        else
            if haskey(dict,s)
                if isa(dict[s],Dict) && isa(subkeymap,Dict)
                    keymap_merge!(dict[s],subkeymap)
                end
            else
                dict[s] = deepcopy(subkeymap)
            end
        end
    end

    function add_specialisations(dict,subdict,level)
        default_branch = subdict['\0']
        if isa(default_branch,Dict)
            for s in keys(default_branch)
                if s == '\0'
                    add_specialisations(dict,default_branch,level+1)
                end
                fixup_keymaps!(dict,level,s,default_branch[s])
            end
        end
    end

    fix_conflicts!(x) = fix_conflicts!(x,1)
    fix_conflicts!(others,level) = nothing
    function fix_conflicts!(dict::Dict,level)
        # needs to be done first for every branch
        if haskey(dict,'\0')
            add_specialisations(dict,dict,level)
        end
        for d in dict
            if d[1] == '\0'
                continue
            end
            fix_conflicts!(d[2],level+1)
        end
    end

    function keymap_prepare(keymaps)
        if isa(keymaps,Dict)
            keymaps = [keymaps]
        end
        push!(keymaps,{"*"=>:(error("Unrecognized input"))})
        @assert isa(keymaps,Array) && eltype(keymaps) <: Dict
        keymaps = map(normalize_keymap,keymaps)
        map(fix_conflicts!,keymaps)
        keymaps
    end

    function keymap_unify(keymaps)
        if length(keymaps) == 1
            return keymaps[1]
        else
            ret = Dict{Char,Any}()
            for keymap in keymaps
                keymap_merge!(ret,keymap)
            end
            fix_conflicts!(ret)
            return ret
        end
    end

    macro keymap(func, keymaps)
        dict = keymap_unify(keymap_prepare(isa(keymaps,Expr)?eval(keymaps):keymaps))
        body = keymap_gen_body(dict,dict)
        esc(quote
            function $(func)(s,data)
                $body
                return :ok
            end
        end)
    end

    const escape_defaults = {
        # Ignore other escape sequences by default
        "\e*" => nothing,
        "\e[*" => nothing,
        # Also ignore extended escape sequences
        # TODO: Support tanges of characters
        "\e[1**" => nothing,
        "\e[2**" => nothing,
        "\e[3**" => nothing,
        "\e[4**" => nothing,
        "\e[5**" => nothing,
        "\e[6**" => nothing,
        "\e[1~" => "\e[H",
        "\e[4~" => "\e[F",
        "\e[7~" => "\e[H",
        "\e[8~" => "\e[F",
        "\eOH"  => "\e[H",
        "\eOF"  => "\e[F",
    }

    function write_response_buffer(s::PromptState,data)
        offset = s.input_buffer.ptr
        ptr = data.respose_buffer.ptr
        seek(data.respose_buffer,0)
        write(s.input_buffer,readall(data.respose_buffer))
        s.input_buffer.ptr = offset+ptr-2
        data.respose_buffer.ptr = ptr
        refresh_line(s)
    end

    type SearchState
        terminal
        histprompt
        #rsearch (true) or ssearch (false)
        backward::Bool
        query_buffer::IOBuffer
        respose_buffer::IOBuffer
        ias::InputAreaState
        #The prompt whose input will be replaced by the matched history
        parent
        SearchState(a,b,c,d,e) = new(a,b,c,d,e,InputAreaState(0,0))
    end

    terminal(s::SearchState) = s.terminal

    function update_display_buffer(s::SearchState,data)
        history_search(data.histprompt.hp,data.query_buffer,data.respose_buffer,data.backward,false) || beep(LineEdit.terminal(s))
        refresh_line(s)
    end

    function history_next_result(s::MIState,data::SearchState)
        #truncate(data.query_buffer,s.input_buffer.size - data.respose_buffer.size)
        history_search(data.histprompt.hp,data.query_buffer,data.respose_buffer,data.backward,true) || beep(LineEdit.terminal(s))
        refresh_line(data)
    end

    function history_set_backward(s::SearchState,backward)
        s.backward = backward
    end

    function refreshMultiLine(s::SearchState)
        buf = IOBuffer()
        write(buf,pointer(s.query_buffer.data),s.query_buffer.ptr-1)
        write(buf,"': ")
        offset = buf.ptr
        ptr = s.respose_buffer.ptr
        seek(s.respose_buffer,0)
        write(buf,readall(s.respose_buffer))
        buf.ptr = offset+ptr-1
        s.respose_buffer.ptr = ptr
        refreshMultiLine(s.terminal,buf,s.ias,s.backward ? "(reverse-i-search)`" : "(i-search)`")
    end

    function reset_state(s::SearchState)
        if s.query_buffer.size != 0
            s.query_buffer.size = 0
            s.query_buffer.ptr = 1
        end
        if s.respose_buffer.size != 0
            s.respose_buffer.size = 0
            s.query_buffer.ptr = 1
        end
        reset_state(s.histprompt.hp)
    end

    type HistoryPrompt <: TextInterface
        hp::HistoryProvider
        keymap_func::Function
        HistoryPrompt(hp) = new(hp)
    end

    init_state(terminal,p::HistoryPrompt) = SearchState(terminal,p,true,IOBuffer(),IOBuffer())

    state(s::MIState,p) = s.mode_state[p]
    state(s::PromptState,p) = (@assert s.p == p; s)
    mode(s::MIState) = s.current_mode
    mode(s::PromptState) = s.p
    mode(s::SearchState) = @assert false

    function accept_result(s,p)
        parent = state(s,p).parent
        replace_line(state(s,parent),state(s,p).respose_buffer)
        transition(s,parent)
    end

    function setup_search_keymap(hp)
        p = HistoryPrompt(hp)
        pkeymap = {
            "^R"    => :( LineEdit.history_set_backward(data,true); LineEdit.history_next_result(s,data) ),
            "^S"    => :( LineEdit.history_set_backward(data,false); LineEdit.history_next_result(s,data) ),
            "\r"    => s->accept_result(s,p),
            "\t"    => nothing, #TODO: Maybe allow tab completion in R-Search?

            # Backspace/^H
            '\b'    => :(LineEdit.edit_backspace(data.query_buffer)?LineEdit.update_display_buffer(s,data):beep(LineEdit.terminal(s))),
            127     => '\b',
            "^C"    => s->transition(s,state(s,p).parent),
            "^D"    => s->transition(s,state(s,p).parent),
            # ^A
            1       => s->(accept_result(s,p); move_line_start(s)),
            # ^E
            5       => s->(accept_result(s,p); move_line_end(s)),
            # Try to catch all Home/End keys
            "\e[H"  => s->(accept_result(s,p); move_input_start(s)),
            "\e[F"  => s->(accept_result(s,p); move_input_end(s)),
            "*"     => :(LineEdit.edit_insert(data.query_buffer,c1);LineEdit.update_display_buffer(s,data))
        }
        @eval @LineEdit.keymap keymap_func $([pkeymap, escape_defaults])
        p.keymap_func = keymap_func
        keymap = {
            "^R"    => s->( state(s,p).parent = mode(s); state(s,p).backward = true; transition(s,p) ),
            "^S"    => s->( state(s,p).parent = mode(s); state(s,p).backward = false; transition(s,p) ),
        }
        (p,keymap)
    end

    keymap(state,p::HistoryPrompt) = p.keymap_func
    keymap_data(state,::HistoryPrompt) = state

    Base.isempty(s::PromptState) = s.input_buffer.size == 0

    on_enter(s::PromptState) = s.p.on_enter(s)

    move_input_start(s) = (seek(buffer(s),0); refresh_line(s))
    move_input_end(s) = (seekend(buffer(s)); refresh_line(s))
    function move_line_start(s)
        buf = buffer(s)
        curpos = position(buf)
        if curpos == 0
            return
        end
        if buf.data[curpos] == '\n'
            move_input_start(s)
        else
            seek(buf,rsearch(buf.data,'\n',curpos-1))
        end
        refresh_line(s)
    end
    function move_line_end(s)
        buf = buffer(s)
        curpos = position(buf)
        if eof(buf)
            return
        end
        c = read(buf,Char)
        if c == '\n'
            move_input_end(s)
            return
        end
        seek(buf,curpos)
        pos = search(buffer(s).data,'\n',curpos+1)
        if pos == 0
            move_input_end(s)
            return
        end
        seek(buf,pos-1)
        refresh_line(s)
    end

    function commit_line(s)
        LineEdit.move_input_end(s)
        println(LineEdit.terminal(s))
        LineEdit.add_history(s)
        LineEdit.state(s,LineEdit.mode(s)).ias =
            LineEdit.InputAreaState(0,0)
    end

    const default_keymap =
    {
        # Tab
        '\t' => s->begin
            buf = buffer(s)
            # Yes, we are ignoring the possiblity
            # the we could be in the middle of a multi-byte
            # sequence, here but that's ok, since any
            # whitespace we're interested in is only one byte
            if position(buf) != 0
                c = buf.data[position(buf)]
            else
                c = '\n'
            end
            if c == ' ' || c == '\n' || c == '\t'
                edit_insert(s," "^4)
                return
            end
            LineEdit.completeLine(s)
            LineEdit.refresh_line(s)
        end,
        # Enter
        '\r' => quote
            if LineEdit.on_enter(s)
                LineEdit.commit_line(s)
                return :done
            else
                LineEdit.edit_insert(s,'\n')
            end
        end,
        '\n' => '\r',
        # Backspace/^H
        '\b' => edit_backspace,
        127 => '\b',
        # ^D
        4 => quote
            if LineEdit.buffer(s).size > 0
                LineEdit.edit_delete(s)
            else
                println(LineEdit.terminal(s))
                return :abort
            end
        end,
        # ^B
        2 => edit_move_left,
        # ^F
        6 => edit_move_right,
        # Meta Enter
        "\e\r" => :(LineEdit.edit_insert(s,'\n')),
        # Simply insert it into the buffer by default
        "*" => :( LineEdit.edit_insert(s,c1) ),
        # ^U
        21 => :( truncate(LineEdit.buffer(s),0); LineEdit.refresh_line(s) ),
        # ^K
        11 => :( truncate(LineEdit.buffer(s),position(LineEdit.buffer(s))); LineEdit.refresh_line(s) ),
        # ^A
        1 => move_line_start,
        # ^E
        5 => move_line_end,
        # Try to catch all Home/End keys
        "\e[H"  => move_input_start,
        "\e[F"  => move_input_end,
        # ^L
        12 => :( Terminals.clear(LineEdit.terminal(s)); LineEdit.refresh_line(s) ),
        # ^W (#edit_delte_prev_word(s))
        23 => :( error("Unimplemented") ),
        # ^C
        "^C" => s->begin
            move_input_end(s);
            LineEdit.refresh_line(s);
            print(LineEdit.terminal(s), "^C\n\n");
            transition(s,:reset);
            LineEdit.refresh_line(s)
        end,
        # Right Arrow
        "\e[C" => edit_move_right,
        # Left Arrow
        "\e[D" => edit_move_left,
        # Up Arrow
        "\e[A" => edit_move_up,
        # Down Arrow
        "\e[B" => edit_move_down,
        # Bracketed Paste Mode
        "\e[200~" => s->begin
            ps = state(s,mode(s))
            input = readuntil(ps.terminal,"\e[201~")[1:(end-6)]
            input = replace(input,'\r','\n')
            if position(buffer(s)) == 0
                indent = Base.indentation(input)[1]
                input = Base.unindent(input[(indent+1):end],indent)
            end
            edit_insert(s,input)
        end,
    }

    function history_keymap(hist)
        return {
            # ^P
            16 => :( LineEdit.history_prev(s,$hist) ),
            # ^N
            14 => :( LineEdit.history_next(s,$hist) ),
            # Up Arrow
            "\e[A" => :( LineEdit.edit_move_up(s) || LineEdit.history_prev(s,$hist) ),
            # Down Arrow
            "\e[B" => :( LineEdit.edit_move_down(s) || LineEdit.history_next(s,$hist) )
        }
    end

    function deactivate(p::Union(Prompt,HistoryPrompt),s::Union(SearchState,PromptState))
        clear_input_area(s)
        s
    end

    function activate(p::Union(Prompt,HistoryPrompt),s::Union(SearchState,PromptState))
        s.ias = InputAreaState(0,0)
        refresh_line(s)
    end

    function activate(p::Union(Prompt,HistoryPrompt),s::MIState)
        @assert p == s.current_mode
        activate(p,s.mode_state[s.current_mode])
    end
    activate(m::ModalInterface,s::MIState) = activate(s.current_mode,s)

    function transition(s::MIState,mode)
        if mode == :abort
            s.aborted = true
            return
        end
        if mode == :reset
            reset_state(s)
            return
        end
        s.mode_state[s.current_mode] = deactivate(s.current_mode,s.mode_state[s.current_mode])
        s.current_mode = mode
        activate(mode,s.mode_state[mode])
    end

    function reset_state(s::PromptState)
        if s.input_buffer.size != 0
            s.input_buffer.size = 0
            s.input_buffer.ptr = 1
        end
        s.ias = InputAreaState(0,0)
    end

    function reset_state(s::MIState)
        for (mode,state) in s.mode_state
            reset_state(state)
        end
    end

    @LineEdit.keymap default_keymap_func [LineEdit.default_keymap,LineEdit.escape_defaults]

    function Prompt(prompt;
        first_prompt = prompt,
        prompt_color="",
        keymap_func = default_keymap_func,
        keymap_func_data = nothing,
        input_color="",
        complete=EmptyCompletionProvider(),
        on_enter=default_enter_cb,on_done=()->nothing,hist=EmptyHistoryProvider())
        Prompt(prompt,first_prompt,prompt_color,keymap_func,keymap_func_data,input_color,complete,on_enter,on_done,hist)
    end

    function run_interface(::Prompt)

    end

    init_state(terminal,prompt::Prompt) = PromptState(terminal,prompt,IOBuffer(),InputAreaState(1,1),length(prompt.prompt))

    function init_state(terminal,m::ModalInterface)
        s = MIState(m,m.modes[1],false,Dict{Any,Any}())
        for mode in m.modes
            s.mode_state[mode] = init_state(terminal,mode)
        end
        s
    end

    function run_interface(terminal,m::ModalInterface)
        s = init_state(terminal,m)
        while !s.aborted
            p = s.current_mode
            buf,ok = prompt!(terminal,m,s)
            s.mode_state[s.current_mode].p.on_done(s,buf,ok)
        end
    end

    buffer(s::PromptState) = s.input_buffer
    buffer(s::SearchState) = s.query_buffer

    keymap(s::PromptState,prompt::Prompt) = prompt.keymap_func
    keymap_data(s::PromptState,prompt::Prompt) = prompt.keymap_func_data
    keymap(ms::MIState,m::ModalInterface) = keymap(ms.mode_state[ms.current_mode],ms.current_mode)
    keymap_data(ms::MIState,m::ModalInterface) = keymap_data(ms.mode_state[ms.current_mode],ms.current_mode)

    function prompt!(terminal,prompt,s=init_state(terminal,prompt))
        raw!(terminal,true)
        isa(terminal,UnixTerminal) && Terminals.Unix.enable_bracketed_paste(terminal)
        try
            start_reading(terminal)
            activate(prompt,s)
            while true
                state = keymap(s,prompt)(s,keymap_data(s,prompt))
                if state == :abort
                    stop_reading(terminal)
                    return (buffer(s),false)
                elseif state == :done
                    stop_reading(terminal)
                    return (buffer(s),true)
                else
                    @assert state == :ok
                end
            end
        finally
            raw!(terminal,false) && Terminals.Unix.disable_bracketed_paste(terminal)
        end
    end
end