diff --git a/rough.jl b/rough.jl
new file mode 100644
index 00000000..df19984f
--- /dev/null
+++ b/rough.jl
@@ -0,0 +1,6 @@
+using Dice
+
+a = uniform(DistFix{5, -6}, 0.01, 1.0)
+pr(a)
+
+DistFix{5, -6}([false, true, true, true, true])
\ No newline at end of file
diff --git a/src/dist/dist.jl b/src/dist/dist.jl
index 91acdd9c..ea564bc4 100644
--- a/src/dist/dist.jl
+++ b/src/dist/dist.jl
@@ -25,6 +25,7 @@ include("misc.jl")
 include("number/uint.jl")
 include("number/int.jl")
 include("number/fix.jl")
+include("number/bitblast.jl")
 include("char.jl")
 include("string.jl")
 include("enum.jl")
diff --git a/src/dist/number/bitblast.jl b/src/dist/number/bitblast.jl
new file mode 100644
index 00000000..df56db13
--- /dev/null
+++ b/src/dist/number/bitblast.jl
@@ -0,0 +1,674 @@
+using Distributions
+using SymPy
+
+export bitblast, bitblast_linear, bitblast_sample, bitblast_exponential, bitblast_exact, unit_exponential, exponential, laplace, unit_gamma, shift_point_gamma, n_unit_exponentials, geometric, general_gamma
+
+#################################
+# primitives for distributions
+#################################
+
+"""
+    uniform(::Type{DistFix{W, F}}, start, stop; rtol=√eps(), kwargs...)
+
+Returns a uniform distribution of type `DistFix{W, F}` in the range [start, stop)
+"""
+function uniform(::Type{DistFix{W, F}}, start, stop; rtol=√eps(), kwargs...) where {W,F}
+    @assert -(2^(W-F-1)) <= start < stop <= 2^(W-F-1)
+    start = start > 0 ? start*(1+rtol) : start*(1-rtol)
+    stop =  stop  > 0 ? stop*(1-rtol)  : stop*(1+rtol) 
+    startint = floor(Int, start*2^float(F))
+    stopint = ceil(Int, (stop*2^float(F)))
+    mantissa = uniform(DistInt{W}, startint, stopint; kwargs...)
+    DistFix{W, F}(mantissa)
+ end
+ 
+
+function triangle(t::Type{DistFix{W, F}}, b::Int) where {W,F}
+    DistFix{W, F}(triangle(DistInt{W}, b))
+end
+
+###################################################
+# sound bit blasting for mixed gamma distributions
+###################################################
+
+# bit blasted exponential distribution on unit interval
+# reverse=true refers to LSB to MSB order of flips 
+function unit_exponential(t::Type{DistFix{W, F}}, beta::Float64; reverse=false) where W where F
+    if !reverse
+    DistFix{W, F}(vcat([false for i in 1:W-F], [flip(exp(beta/2^i)/(1+exp(beta/2^i))) for i in 1:F]))
+    else
+    DistFix{W, F}(vcat([false for i in 1:W-F], reverse([flip(exp(beta/2^i)/(1+exp(beta/2^i))) for i in F:-1:1])))
+    end
+end
+  
+function bitblast(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+    numpieces::Int, start::Float64, stop::Float64) where {W,F}
+
+    # count bits and pieces
+    @assert -(2^(W-F-1)) <= start < stop <= 2^(W-F-1)
+    f_range_bits = log2((stop - start)*2^float(F))
+    @assert isinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."
+    @assert ispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
+    intervals_per_piece = (2^Int(f_range_bits))/numpieces
+    bits_per_piece = Int(log2(intervals_per_piece))
+    # @show bits_per_piece
+
+    dist = truncated(dist, start, stop)
+
+    total_prob = 0
+    piece_probs = Vector(undef, numpieces) # prob of each piece
+    border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
+    linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
+
+    for i=1:numpieces
+    firstinter = start + (i-1)*intervals_per_piece/2^float(F) 
+    lastinter = start + (i)*intervals_per_piece/2^float(F)
+
+    piece_probs[i] = (cdf(dist, lastinter) - cdf(dist, firstinter))
+    total_prob += piece_probs[i]
+
+    border_probs[i] = [cdf(dist, firstinter + 1/2^float(F) ) - cdf(dist, firstinter), 
+            cdf(dist, lastinter) - cdf(dist, lastinter - 1/2^float(F) )]
+    @show bits_per_piece
+    linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2 * 2^(bits_per_piece)
+    end
+
+    PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
+    piecechoice = discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
+    slope_flips = Vector(undef, numpieces)
+
+    isdecreasing = Vector(undef, numpieces)
+    for i=numpieces:-1:1
+    iszero(linear_piece_probs[i]) && iszero(piece_probs[i]) && continue
+    a = border_probs[i][1]/linear_piece_probs[i]
+    isdecreasing[i] = a > 1/2^bits_per_piece
+    if isdecreasing[i]
+    slope_flips[i] = flip(2-a*2^bits_per_piece)
+    else
+    slope_flips[i] = flip(a*2^bits_per_piece)
+    end  
+    end
+
+    unif = uniform(DistFix{W,F},  bits_per_piece)
+    tria = triangle(DistFix{W,F}, bits_per_piece)
+    z = nothing
+    for i=1:numpieces
+    iszero(linear_piece_probs[i]) && continue
+    firstinterval = DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^float(F))
+    lastinterval = DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^float(F))
+    linear_dist = 
+    if isdecreasing[i]
+    (ifelse(slope_flips[i], 
+        (firstinterval + unif), 
+        (lastinterval - tria)))
+    else
+    firstinterval + ifelse(slope_flips[i], unif, tria)
+    end
+    z = if isnothing(z)
+    linear_dist
+    else
+    ifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)  
+    end
+    end
+    return z
+end
+
+
+function bitblast2(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+    numpieces::Int, start::Float64, stop::Float64, strategy=:linear) where {W,F}
+
+    @assert ispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
+    intervals_per_piece = (2^Int(f_range_bits))/numpieces
+    bits_per_piece = Int(log2(intervals_per_piece))
+
+    dist = truncated(dist, start, stop)
+
+    total_prob = 0
+    piece_probs = Vector(undef, numpieces) # prob of each piece
+    border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
+    linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
+
+    for i=1:numpieces
+    firstinter = start + (i-1)*intervals_per_piece/2^F 
+    lastinter = start + (i)*intervals_per_piece/2^F 
+
+    piece_probs[i] = (cdf(dist, lastinter) - cdf(dist, firstinter))
+    total_prob += piece_probs[i]
+
+    border_probs[i] = [cdf(dist, firstinter + 1/2^F ) - cdf(dist, firstinter), 
+            cdf(dist, lastinter) - cdf(dist, lastinter - 1/2^F )]
+    linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2 * 2^(bits_per_piece)
+    end
+
+    PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
+    piecechoice = discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
+    slope_flips = Vector(undef, numpieces)
+
+    isdecreasing = Vector(undef, numpieces)
+    for i=numpieces:-1:1
+    iszero(linear_piece_probs[i]) && iszero(piece_probs[i]) && continue
+    a = border_probs[i][1]/linear_piece_probs[i]
+    isdecreasing[i] = a > 1/2^bits_per_piece
+    if isdecreasing[i]
+    slope_flips[i] = flip(2-a*2^bits_per_piece)
+    else
+    slope_flips[i] = flip(a*2^bits_per_piece)
+    end  
+    end
+
+    unif = uniform(DistFix{W,F},  bits_per_piece)
+    tria = triangle(DistFix{W,F}, bits_per_piece)
+    z = nothing
+    for i=1:numpieces
+    iszero(linear_piece_probs[i]) && continue
+    firstinterval = DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^F)
+    lastinterval = DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^F)
+    linear_dist = 
+    if isdecreasing[i]
+    (ifelse(slope_flips[i], 
+        (firstinterval + unif), 
+        (lastinterval - tria)))
+    else
+    firstinterval + ifelse(slope_flips[i], unif, tria)
+    end
+    z = if isnothing(z)
+    linear_dist
+    else
+    ifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)  
+    end
+    end
+    return z
+end
+
+
+function bitblast(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+    start::Float64, stop::Float64, blast_strategy=:exact; kwargs...) where {W,F}
+    if blast_strategy == :linear
+    bitblast_linear(DistFix{W,F}, dist, start, stop; kwargs...)
+    else
+    error("Unknown bitblasting strategy: $strategy")
+    end
+end
+
+"Approximate a continuous distribution on an interval using a bit-blasted linear density"
+function bitblast_linear(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+    start::Float64, stop::Float64;
+    slope_flip = nothing, unif = nothing, tria = nothing) where {W,F}
+
+    dist = truncated(dist, start, stop)
+    firstprob = cdf(dist, start + 1/2^F ) - cdf(dist, start)
+    lastprob = cdf(dist, stop) - cdf(dist, stop - 1/2^F)
+    avgprob = (firstprob + lastprob)/2
+
+    @assert !iszero(firstprob) || !iszero(lastprob) "No probability mass found at the given boundaries."
+
+    @assert ispow2((stop - start)*2^F) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two."
+    num_bits = Int(log2((stop-start)*2^F))
+    isnothing(slope_flip) && (slope_flip = flip())
+    isnothing(unif) && (unif = uniform(DistFix{W,F}, num_bits))
+    isnothing(tria) && (tria = triangle(DistFix{W,F}, num_bits))
+
+    firstinterval = DistFix{W,F}(start)
+    lastinterval = DistFix{W,F}(start + (2^num_bits-1)/2^F)
+    if firstprob > avgprob # the slope is decreasing
+    slope_flip = flip_prob!(slope_flip, 2-firstprob/avgprob)
+    ifelse(slope_flip, (firstinterval + unif), (lastinterval - tria))
+    else # the slope is increasing
+    slope_flip = flip_prob!(slope_flip, firstprob/avgprob)
+    firstinterval + ifelse(slope_flip, unif, tria)
+    end  
+end
+
+# bit blasted exponential distribution on arbitrary interval
+# reverse=true refers to LSB to MSB order of flips
+# TODO: make the following function use unit_exponential
+function exponential(t::Type{DistFix{W, F}}, beta::Float64, start::Float64, stop::Float64; reverse=false) where W where F   
+range = stop - start
+@assert ispow2(range)
+
+new_beta = beta*range
+bits = Int(log2(range)) + F
+
+if !reverse
+bit_vector = vcat([false for i in 1:W - bits], [flip(exp(new_beta/2^i)/(1+exp(new_beta/2^i))) for i in 1:bits])
+else 
+bit_vector = vcat([false for i in 1:W - bits], reverse([flip(exp(new_beta/2^i)/(1+exp(new_beta/2^i))) for i in bits:-1:1]))
+end
+DistFix{W, F}(bit_vector) + DistFix{W, F}(start)
+end
+
+function beta(d::ContinuousUnivariateDistribution, start::Float64, stop::Float64, interval_sz::Float64)
+prob_start = cdf(d, start + interval_sz) - cdf(d, start)
+if prob_start == 0.0
+prob_start = eps(0.0)
+end
+prob_end = cdf(d, stop) - cdf(d, stop - interval_sz)
+result = (log(prob_end) - log(prob_start)) / (stop - start - interval_sz)
+
+if result ≈ Inf
+@show prob_start
+@show prob_end
+@show start, stop, interval_sz
+end
+result
+end
+
+function bitblast_exponential(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+numpieces::Int, start::Float64, stop::Float64, strategy=:exponential) where {W,F}
+
+# basic checks
+@assert start >= -(2^(W - F - 1))
+@assert stop <= (2^(W - F - 1))
+@assert start < stop
+a = Int(log2((stop - start)*2^F))
+@assert a isa Int 
+@assert ispow2(numpieces) "Number of pieces must be a power of two"
+piece_bits = Int(log2(numpieces))
+if piece_bits == 0
+piece_bits = 1
+end
+@assert typeof(piece_bits) == Int
+
+# preliminaries
+d = truncated(dist, start, stop)
+whole_bits = a
+point = F
+interval_sz = (2^whole_bits/numpieces)
+bits = Int(log2(interval_sz))
+areas = Vector(undef, numpieces)
+total_area = 0
+
+beta_vec = Vector(undef, numpieces)
+start_pts = Vector(undef, numpieces)
+stop_pts = Vector(undef, numpieces)
+
+
+
+# Figuring out end points
+for i=1:numpieces
+p1 = start + (i-1)*interval_sz/2^point 
+p3 = start + (i)*interval_sz/2^point 
+
+beta_vec[i] = beta(d, p1, p3, 2.0^(-F))
+# if (beta_vec[i] in [Inf, -Inf]) | isnan(beta_vec[i])
+#     beta_vec[i] = 0.0
+# end
+
+areas[i] = (cdf.(d, p3) - cdf.(d, p1))
+# @show p1, p2, p3, p4, areas[i]
+start_pts[i] = p1
+stop_pts[i] = p3
+
+total_area += areas[i]
+end
+
+
+# @show beta_vec
+
+rel_prob = areas/total_area
+
+b = discrete(DistUInt{piece_bits}, rel_prob)
+
+ans = DistFix{W, F}((2^(W-1)-1)/2^F)
+
+for i=numpieces:-1:1
+ans = ifelse( prob_equals(b, DistUInt{piece_bits}(i-1)), 
+  exponential(DistFix{W, F}, beta_vec[i], start_pts[i], stop_pts[i]),
+  ans)  
+end
+return ans
+end
+
+function continuous(t::Type{DistFix{W, F}}, d::ContinuousUnivariateDistribution, pieces::Int, start::Float64, stop::Float64, exp::Bool=false) where {W, F}
+c = if exp
+continuous_exp(DistFix{W, F}, d, pieces, start, stop)
+else 
+continuous_linear(DistFix{W, F}, d, pieces, start, stop)
+end
+return c
+end
+
+# https://en.wikipedia.org/wiki/Laplace_distribution
+function laplace(t::Type{DistFix{W, F}}, mean::Float64, scale::Float64, start::Float64, stop::Float64) where {W, F}
+@assert scale > 0
+
+beta1 = -1/scale
+e1 = exponential(DistFix{W, F}, beta1, mean, stop)
+
+beta2 = 1/scale
+e2 = exponential(DistFix{W, F}, beta2, start, mean)
+
+ifelse(flip(0.5), e1, e2)
+end
+
+#Helper function that returns exponentials 
+
+function shift_point_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64) where {W, F}
+DFiP = DistFix{W, F}
+if alpha == 0
+unit_exponential(DFiP, beta)
+else
+x1 = shift_point_gamma(DFiP, alpha - 1, beta)
+x2 = uniform(DFiP, 0.0, 1.0)
+observe(ifelse(flip(1/(1 + 2.0^(-F))), x2 < x1, true))
+x1
+end
+end
+
+#https://www.wolframalpha.com/input?i=sum+%28a*epsilon%29%5E2+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1
+function sum_qgp(β::Float64, ϵ::Float64)
+ans = (1/ϵ - 1)^2 * exp(β*ϵ*(2 + 1/ϵ))
+ans += (1/ϵ^2)*exp(β)
+ans += (2/ϵ - 2/ϵ^2 + 1)*exp(β * (1 + ϵ))
+ans -= exp(β*ϵ)
+ans -= exp(2*β*ϵ)
+ans *= ϵ^2
+ans /= (exp(β * ϵ) - 1)^3
+ans
+end
+
+#https://www.wolframalpha.com/input?i=sum+%28a*epsilon%29+*+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1
+function sum_agp(β::Float64, ϵ::Float64)
+ans = (1/ϵ - 1)*exp(β*(1 + ϵ))
+ans -= exp(β)/ϵ
+ans += exp(β*ϵ)
+ans *= ϵ
+ans /= (exp(β*ϵ) - 1)^2
+ans
+end
+
+#https://www.wolframalpha.com/input?i=sum+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1
+function sum_gp(β::Float64, ϵ::Float64)
+ans = (exp(β) - 1) / (exp(β*ϵ) - 1)
+ans
+end
+
+function sum_pgp(β::Float64, ϵ::Float64, p::Int)
+if p == 0
+sum_gp(β, ϵ)
+elseif p == 1
+sum_agp(β, ϵ)
+elseif p == 2
+sum_qgp(β, ϵ)
+else
+sum = 0
+for i = 0:ϵ:1-ϵ
+sum += i^p * exp(β*i)
+end
+sum
+end
+end
+
+function n_unit_exponentials(::Type{DistFix{W, F}}, betas::Vector{Float64}) where {W, F}
+DFiP = DistFix{W, F}
+l = length(betas)
+ans = Vector(undef, l)
+for i in 1:l
+ans[i] = Vector(undef, W)
+end
+for i in 1:W-F
+for j in 1:l
+ans[j][i] = false
+end
+end
+for i in F:-1:1
+for j in 1:l
+ans[j][i + W - F] = flip(exp(betas[j]/2^i)/(1+exp(betas[j]/2^i)))
+end
+end
+[DFiP(i) for i in ans] 
+end
+
+function exponential_for_gamma(α::Int, β::Float64)::Vector{Float64}
+if α == 0
+[]
+elseif α == 1
+[β, β, 0.0]
+else
+v = []
+for i in 1:α
+v = vcat(vcat([β], zeros(i-1)), v)
+end
+
+vcat(vcat(exponential_for_gamma(α-1, β), [0.0]), v)
+end
+end
+
+function gamma_constants(α::Int, β::Float64, ϵ::Float64)
+@syms varint
+@syms v2
+if α == 0
+[]
+else
+c1 = Float64(sympy.Poly(integrate(varint^α*exp(β*varint), (varint, 0, 1)), varint).coeffs().evalf()[1])
+c2 = [Float64(i) for i in sympy.Poly(simplify(v2*integrate(varint^(α-1)*exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2).coeffs()]
+p1 = 0
+for i in eachindex(c2)
+p1 += sum_pgp(β, ϵ, length(c2) + 1 - i) * c2[i]
+end
+p1 /= c1
+
+c2 = [Float64(i) for i in sympy.Poly(simplify(v2*integrate(varint^(α-1) * (varint - v2) *exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2).coeffs()]
+# @show c2
+# @show sympy.Poly(simplify(v2*integrate(varint^(α-1) * (varint - v2) *exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2)
+p2 = Vector(undef, α)
+for i in eachindex(c2)
+p2[i] = sum_pgp(β, ϵ, length(c2) - i) * c2[i]
+end
+# @show p2
+
+vcat([p1], p2, gamma_constants(α-1, β, ϵ))
+end
+end
+
+
+
+
+function unit_gamma(t::Type{DistFix{W, F}}, alpha::Int, beta::Float64; vec_arg=[], constants = [], discrete_bdd=[], coinflips = []) where {W, F}
+DFiP = DistFix{W, F}
+if alpha == 0
+unit_exponential(DFiP, beta)
+elseif alpha == 1
+
+t = (exp(beta*2.0^(-F))*(beta*2.0^(-F) - 1) + 1)*(1 - exp(beta)) / ((1 - exp(beta*2.0^(-F)))*(exp(beta) * (beta - 1) + 1))
+
+if (length(coinflips) != 0)
+coinflip = coinflips[1]
+else
+coinflip = flip(t)
+end
+
+
+if (length(vec_arg) != 0)
+(Y, Z, U) = vec_arg
+else
+(Y, Z, U) = n_unit_exponentials(DFiP, [beta, beta, 0.0])
+end
+observe(U < Y)
+
+
+
+final = ifelse(coinflip, Z, Y)
+final
+else 
+α = alpha
+β = beta
+if (length(vec_arg) != 0)
+vec_expo = vec_arg
+
+else
+discrete_bdd = Vector(undef, α)
+constants = gamma_constants(alpha, beta, 1/2^F)
+
+t = (exp(beta*2.0^(-F))*(beta*2.0^(-F) - 1) + 1)*(1 - exp(beta)) / ((1 - exp(beta*2.0^(-F)))*(exp(beta) * (beta - 1) + 1))
+f = flip(t)
+
+count = 0
+for i in α:-1:1
+  # @show constants
+  l = discrete(DistUInt{max(Int(ceil(log(i))), 1)}, normalize(constants[count + 2:count+i+1]))
+  count = count+i+1
+
+  discrete_bdd[α - i + 1] = l
+end
+vec_expo = n_unit_exponentials(DFiP, exponential_for_gamma(α, β))
+
+end
+
+seq = Int(α*(α^2 + 5)/6)
+# @show constant_flips
+x1 = unit_gamma(DFiP, alpha-1, beta, vec_arg=vec_expo[1:seq], constants=constants[α + 2:length(constants)], discrete_bdd=discrete_bdd[2:α])
+x2 = vec_expo[seq + 1]
+observe(x2 < x1)
+
+discrete_dist_vec = Vector(undef, α)
+count = seq+2
+for i in 1:α 
+x = vec_expo[count]
+count+=1
+for j in 1:α - i
+  observe(vec_expo[count] < x)
+  count+=1
+end
+discrete_dist_vec[i] = x
+end
+
+# l = discrete(DistUInt{Int(ceil(log(α)))}, normalize(constants[2:α+1]))
+l = discrete_bdd[1]
+t = DFiP(0.0)
+for i in 1:α
+t = ifelse(prob_equals(l, DistUInt{Int(ceil(log(α)))}(i-1)), discrete_dist_vec[i], t)
+end
+
+ifelse(flip(constants[1]), x1, t)
+end
+
+end
+
+function normalize(v)
+l = sum(v)
+[i/l for i in v]    
+end
+# # TODO: Write tests for the following function
+# function unit_concave(t::Type{DistFix{W, F}}, beta::Float64) where {W, F}
+#     @assert beta <= 0
+#     DFiP = DistFix{W, F}
+#     Y = uniform(DFiP, 0.0, 1.0)
+#     X = unit_exponential(DFiP, beta)
+#     observe((X < Y)| prob_equals(X, Y))
+#     Y
+# end
+
+######################################################################################################
+# bit blasting geometric distribution: https://en.wikipedia.org/wiki/Geometric_distribution
+######################################################################################################
+
+# Creates a geometric distribution over the integers in the range [0, stop-1]
+function geometric(::Type{DistFix{W, F}}, success::Float64, stop::Int) where {W, F}
+#TODO: use convert
+@assert ispow2(stop)
+bits = Int(log2(stop))
+@assert W - F > bits
+
+convert(DistFix{W, F}, DistFix{W, 0}(unit_exponential(DistFix{bits+1, bits}, log(1 - success)*2^bits).mantissa))
+end
+
+"""
+general_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64, ll::Float64, ul::Float64)
+
+Returns bitblast distribution for the density function x^alpha e^(beta*x)
+"""
+function general_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64, ll::Float64, ul::Float64) where {W, F}
+@assert ispow2(ul - ll)
+multiply = Int(log2(ul - ll))
+start = DistFix{W, F}(ul)
+new_type = DistFix{W, F + multiply} 
+
+DistFix{W, F}(unit_gamma(new_type, alpha, beta).mantissa.number.bits) + start
+end
+
+"""
+bitblast_sample(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+numpieces::Int, start::Float64, stop::Float64, offset::Float64, width::Float64)
+
+A modified version of the function bitblast that works with the assumption of lower bits being sampled.
+"""
+function bitblast_sample(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
+numpieces::Int, start::Float64, stop::Float64, offset::Float64, width::Float64) where {W,F}
+
+# count bits and pieces
+@assert -(2^(W-F-1)) <= start < stop <= 2^(W-F-1)
+f_range_bits = log2((stop - start)*2^F)
+@assert isinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."
+@assert ispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
+intervals_per_piece = (2^Int(f_range_bits))/numpieces
+bits_per_piece = Int(log2(intervals_per_piece))
+
+dist = truncated(dist, start, stop)
+
+total_prob = 0
+piece_probs = Vector(undef, numpieces) # prob of each piece
+border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
+linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
+
+for i=1:numpieces
+firstinter = start + (i-1)*intervals_per_piece/2^F 
+lastinter = start + (i)*intervals_per_piece/2^F 
+
+piece_probs[i] = 0
+# Warning: A potential source of high runtime
+for j=1:intervals_per_piece
+@show firstinter + offset + width + (j-1)/2^F
+@show firstinter + offset + (j-1)/2^F
+piece_probs[i] += cdf(dist, firstinter + offset + width + (j-1)/2^F) - cdf(dist, firstinter + offset + (j-1)/2^F)
+end
+total_prob += piece_probs[i]
+
+border_probs[i] = [cdf(dist, firstinter + offset+width) - cdf(dist, firstinter + offset), 
+  cdf(dist, lastinter -1/2^F + offset+width) - cdf(dist, lastinter - 1/2^F + offset)]
+linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2 * 2^(bits_per_piece)
+end
+
+
+@show piece_probs
+@show total_prob
+
+PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
+piecechoice = discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
+slope_flips = Vector(undef, numpieces)
+
+isdecreasing = Vector(undef, numpieces)
+for i=numpieces:-1:1
+iszero(linear_piece_probs[i]) && iszero(piece_probs[i]) && continue
+a = border_probs[i][1]/linear_piece_probs[i]
+isdecreasing[i] = a > 1/2^bits_per_piece
+if isdecreasing[i]
+slope_flips[i] = flip(2-a*2^bits_per_piece)
+else
+slope_flips[i] = flip(a*2^bits_per_piece)
+end  
+end
+
+unif = uniform(DistFix{W,F},  bits_per_piece)
+tria = triangle(DistFix{W,F}, bits_per_piece)
+z = nothing
+for i=1:numpieces
+iszero(linear_piece_probs[i]) && continue
+firstinterval = DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^F)
+lastinterval = DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^F)
+linear_dist = 
+if isdecreasing[i]
+  (ifelse(slope_flips[i], 
+      (firstinterval + unif), 
+      (lastinterval - tria)))
+else
+  firstinterval + ifelse(slope_flips[i], unif, tria)
+end
+z = if isnothing(z)
+  linear_dist
+else
+  ifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)  
+end
+end
+return z
+end
\ No newline at end of file
diff --git a/src/dist/number/fix.jl b/src/dist/number/fix.jl
index 372031f4..d7ce0632 100644
--- a/src/dist/number/fix.jl
+++ b/src/dist/number/fix.jl
@@ -1,13 +1,12 @@
 using Distributions
-using SymPy
 
-export DistFix, bitblast, bitblast_linear, bitblast_sample, bitblast_exponential, bitblast_exact, unit_exponential, exponential, laplace, unit_gamma, shift_point_gamma, n_unit_exponentials, geometric, general_gamma
+export DistFix
 
 ##################################
 # types, structs, and constructors
 ##################################
 
-"A signed fixed point number with `W` bits of which `F` are after the binary point"
+"A signed fixed point number with `W` bits in mantissa and exponent as 2^{-F}"
 struct DistFix{W, F} <: Dist{Int}
     mantissa::DistInt{W}
     function DistFix{W, F}(b) where {W,F}
@@ -29,21 +28,6 @@ function uniform(::Type{DistFix{W, F}}, n = W) where {W,F}
     DistFix{W, F}(uniform(DistInt{W}, n))
 end
 
-function uniform(t::Type{DistFix{W, F}}, start, stop; rtol=√eps(), kwargs...) where {W,F}
-    @assert -(2^(W-F-1)) <= start < stop <= 2^(W-F-1)
-    start = start > 0 ? start*(1+rtol) : start*(1-rtol)
-    stop =  stop  > 0 ? stop*(1-rtol)  : stop*(1+rtol) 
-    startint = floor(Int, start*2^F)
-    stopint = ceil(Int, (stop*2^F))
-    mantissa = uniform(DistInt{W}, startint, stopint; kwargs...)
-    DistFix{W, F}(mantissa)
- end
- 
-
-function triangle(t::Type{DistFix{W, F}}, b::Int) where {W,F}
-    DistFix{W, F}(triangle(DistInt{W}, b))
-end
-
 function Base.convert(::Type{DistFix{W2, F2}}, x::DistFix{W1, F1}) where {W1,W2,F1,F2}
     if (F1 == F2)
         return DistFix{W2, F2}(convert(DistInt{W2}, x.mantissa))
@@ -115,653 +99,3 @@ function Base.:(/)(x::DistFix{W, F}, y::DistFix{W, F}) where {W, F}
     z_mantissa = convert(DistInt{W}, z_mantissa_big)
     DistFix{W, F}(z_mantissa)
 end
-
-#################################
-# bit blasting continuous distributions
-#################################
-  
-function bitblast(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-                  numpieces::Int, start::Float64, stop::Float64) where {W,F}
-
-    # count bits and pieces
-    @assert -(2^(W-F-1)) <= start < stop <= 2^(W-F-1)
-    f_range_bits = log2((stop - start)*2^float(F))
-    @assert isinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."
-    @assert ispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
-    intervals_per_piece = (2^Int(f_range_bits))/numpieces
-    bits_per_piece = Int(log2(intervals_per_piece))
-    # @show bits_per_piece
-
-    dist = truncated(dist, start, stop)
-
-    total_prob = 0
-    piece_probs = Vector(undef, numpieces) # prob of each piece
-    border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
-    linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
-
-    for i=1:numpieces
-        firstinter = start + (i-1)*intervals_per_piece/2^float(F) 
-        lastinter = start + (i)*intervals_per_piece/2^float(F)
-
-        piece_probs[i] = (cdf(dist, lastinter) - cdf(dist, firstinter))
-        total_prob += piece_probs[i]
-
-        border_probs[i] = [cdf(dist, firstinter + 1/2^float(F) ) - cdf(dist, firstinter), 
-                        cdf(dist, lastinter) - cdf(dist, lastinter - 1/2^float(F) )]
-        @show bits_per_piece
-        linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2 * 2^(bits_per_piece)
-    end
-
-    PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
-    piecechoice = discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
-    slope_flips = Vector(undef, numpieces)
-
-    isdecreasing = Vector(undef, numpieces)
-    for i=numpieces:-1:1
-        iszero(linear_piece_probs[i]) && iszero(piece_probs[i]) && continue
-        a = border_probs[i][1]/linear_piece_probs[i]
-        isdecreasing[i] = a > 1/2^bits_per_piece
-        if isdecreasing[i]
-            slope_flips[i] = flip(2-a*2^bits_per_piece)
-        else
-            slope_flips[i] = flip(a*2^bits_per_piece)
-        end  
-    end
-
-    unif = uniform(DistFix{W,F},  bits_per_piece)
-    tria = triangle(DistFix{W,F}, bits_per_piece)
-    z = nothing
-    for i=1:numpieces
-        iszero(linear_piece_probs[i]) && continue
-        firstinterval = DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^float(F))
-        lastinterval = DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^float(F))
-        linear_dist = 
-            if isdecreasing[i]
-                (ifelse(slope_flips[i], 
-                    (firstinterval + unif), 
-                    (lastinterval - tria)))
-            else
-                firstinterval + ifelse(slope_flips[i], unif, tria)
-            end
-        z = if isnothing(z)
-                linear_dist
-            else
-                ifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)  
-            end
-    end
-    return z
-end
-
-
-function bitblast2(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-                  numpieces::Int, start::Float64, stop::Float64, strategy=:linear) where {W,F}
-
-    @assert ispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
-    intervals_per_piece = (2^Int(f_range_bits))/numpieces
-    bits_per_piece = Int(log2(intervals_per_piece))
-
-    dist = truncated(dist, start, stop)
-
-    total_prob = 0
-    piece_probs = Vector(undef, numpieces) # prob of each piece
-    border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
-    linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
-
-    for i=1:numpieces
-        firstinter = start + (i-1)*intervals_per_piece/2^F 
-        lastinter = start + (i)*intervals_per_piece/2^F 
-
-        piece_probs[i] = (cdf(dist, lastinter) - cdf(dist, firstinter))
-        total_prob += piece_probs[i]
-
-        border_probs[i] = [cdf(dist, firstinter + 1/2^F ) - cdf(dist, firstinter), 
-                        cdf(dist, lastinter) - cdf(dist, lastinter - 1/2^F )]
-        linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2 * 2^(bits_per_piece)
-    end
-
-    PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
-    piecechoice = discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
-    slope_flips = Vector(undef, numpieces)
-
-    isdecreasing = Vector(undef, numpieces)
-    for i=numpieces:-1:1
-        iszero(linear_piece_probs[i]) && iszero(piece_probs[i]) && continue
-        a = border_probs[i][1]/linear_piece_probs[i]
-        isdecreasing[i] = a > 1/2^bits_per_piece
-        if isdecreasing[i]
-            slope_flips[i] = flip(2-a*2^bits_per_piece)
-        else
-            slope_flips[i] = flip(a*2^bits_per_piece)
-        end  
-    end
-
-    unif = uniform(DistFix{W,F},  bits_per_piece)
-    tria = triangle(DistFix{W,F}, bits_per_piece)
-    z = nothing
-    for i=1:numpieces
-        iszero(linear_piece_probs[i]) && continue
-        firstinterval = DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^F)
-        lastinterval = DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^F)
-        linear_dist = 
-            if isdecreasing[i]
-                (ifelse(slope_flips[i], 
-                    (firstinterval + unif), 
-                    (lastinterval - tria)))
-            else
-                firstinterval + ifelse(slope_flips[i], unif, tria)
-            end
-        z = if isnothing(z)
-                linear_dist
-            else
-                ifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)  
-            end
-    end
-    return z
-end
-
-
-function bitblast(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-                  start::Float64, stop::Float64, blast_strategy=:exact; kwargs...) where {W,F}
-    if blast_strategy == :linear
-        bitblast_linear(DistFix{W,F}, dist, start, stop; kwargs...)
-    else
-        error("Unknown bitblasting strategy: $strategy")
-    end
-end
-
-"Approximate a continuous distribution on an interval using a bit-blasted linear density"
-function bitblast_linear(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-                  start::Float64, stop::Float64;
-                  slope_flip = nothing, unif = nothing, tria = nothing) where {W,F}
-
-    dist = truncated(dist, start, stop)
-    firstprob = cdf(dist, start + 1/2^F ) - cdf(dist, start)
-    lastprob = cdf(dist, stop) - cdf(dist, stop - 1/2^F)
-    avgprob = (firstprob + lastprob)/2
-
-    @assert !iszero(firstprob) || !iszero(lastprob) "No probability mass found at the given boundaries."
-
-    @assert ispow2((stop - start)*2^F) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two."
-    num_bits = Int(log2((stop-start)*2^F))
-    isnothing(slope_flip) && (slope_flip = flip())
-    isnothing(unif) && (unif = uniform(DistFix{W,F}, num_bits))
-    isnothing(tria) && (tria = triangle(DistFix{W,F}, num_bits))
-
-    firstinterval = DistFix{W,F}(start)
-    lastinterval = DistFix{W,F}(start + (2^num_bits-1)/2^F)
-    if firstprob > avgprob # the slope is decreasing
-        slope_flip = flip_prob!(slope_flip, 2-firstprob/avgprob)
-        ifelse(slope_flip, (firstinterval + unif), (lastinterval - tria))
-    else # the slope is increasing
-        slope_flip = flip_prob!(slope_flip, firstprob/avgprob)
-        firstinterval + ifelse(slope_flip, unif, tria)
-    end  
-end
-
-###################################################
-# sound bit blasting for mixed gamma distributions
-###################################################
-
-# bit blasted exponential distribution on unit interval
-# reverse=true refers to LSB to MSB order of flips 
-function unit_exponential(t::Type{DistFix{W, F}}, beta::Float64; reverse=false) where W where F
-    if !reverse
-        DistFix{W, F}(vcat([false for i in 1:W-F], [flip(exp(beta/2^i)/(1+exp(beta/2^i))) for i in 1:F]))
-    else
-        DistFix{W, F}(vcat([false for i in 1:W-F], reverse([flip(exp(beta/2^i)/(1+exp(beta/2^i))) for i in F:-1:1])))
-    end
-end
-
-# bit blasted exponential distribution on arbitrary interval
-# reverse=true refers to LSB to MSB order of flips
-# TODO: make the following function use unit_exponential
-function exponential(t::Type{DistFix{W, F}}, beta::Float64, start::Float64, stop::Float64; reverse=false) where W where F   
-    range = stop - start
-    @assert ispow2(range)
-
-    new_beta = beta*range
-    bits = Int(log2(range)) + F
-    
-    if !reverse
-        bit_vector = vcat([false for i in 1:W - bits], [flip(exp(new_beta/2^i)/(1+exp(new_beta/2^i))) for i in 1:bits])
-    else 
-        bit_vector = vcat([false for i in 1:W - bits], reverse([flip(exp(new_beta/2^i)/(1+exp(new_beta/2^i))) for i in bits:-1:1]))
-    end
-    DistFix{W, F}(bit_vector) + DistFix{W, F}(start)
-end
-
-function beta(d::ContinuousUnivariateDistribution, start::Float64, stop::Float64, interval_sz::Float64)
-    prob_start = cdf(d, start + interval_sz) - cdf(d, start)
-    if prob_start == 0.0
-        prob_start = eps(0.0)
-    end
-    prob_end = cdf(d, stop) - cdf(d, stop - interval_sz)
-    result = (log(prob_end) - log(prob_start)) / (stop - start - interval_sz)
-
-    if result ≈ Inf
-        @show prob_start
-        @show prob_end
-        @show start, stop, interval_sz
-    end
-    result
-end
-
-function bitblast_exponential(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-    numpieces::Int, start::Float64, stop::Float64, strategy=:exponential) where {W,F}
-
-    # basic checks
-    @assert start >= -(2^(W - F - 1))
-    @assert stop <= (2^(W - F - 1))
-    @assert start < stop
-    a = Int(log2((stop - start)*2^F))
-    @assert a isa Int 
-    @assert ispow2(numpieces) "Number of pieces must be a power of two"
-    piece_bits = Int(log2(numpieces))
-    if piece_bits == 0
-        piece_bits = 1
-    end
-    @assert typeof(piece_bits) == Int
-
-    # preliminaries
-    d = truncated(dist, start, stop)
-    whole_bits = a
-    point = F
-    interval_sz = (2^whole_bits/numpieces)
-    bits = Int(log2(interval_sz))
-    areas = Vector(undef, numpieces)
-    total_area = 0
-
-    beta_vec = Vector(undef, numpieces)
-    start_pts = Vector(undef, numpieces)
-    stop_pts = Vector(undef, numpieces)
-
-    
-
-    # Figuring out end points
-    for i=1:numpieces
-        p1 = start + (i-1)*interval_sz/2^point 
-        p3 = start + (i)*interval_sz/2^point 
-
-        beta_vec[i] = beta(d, p1, p3, 2.0^(-F))
-        # if (beta_vec[i] in [Inf, -Inf]) | isnan(beta_vec[i])
-        #     beta_vec[i] = 0.0
-        # end
-
-        areas[i] = (cdf.(d, p3) - cdf.(d, p1))
-        # @show p1, p2, p3, p4, areas[i]
-        start_pts[i] = p1
-        stop_pts[i] = p3
-
-        total_area += areas[i]
-    end
-
-
-    # @show beta_vec
-
-    rel_prob = areas/total_area
-
-    b = discrete(DistUInt{piece_bits}, rel_prob)
-
-    ans = DistFix{W, F}((2^(W-1)-1)/2^F)
-
-    for i=numpieces:-1:1
-        ans = ifelse( prob_equals(b, DistUInt{piece_bits}(i-1)), 
-                exponential(DistFix{W, F}, beta_vec[i], start_pts[i], stop_pts[i]),
-                ans)  
-    end
-    return ans
-end
-
-function continuous(t::Type{DistFix{W, F}}, d::ContinuousUnivariateDistribution, pieces::Int, start::Float64, stop::Float64, exp::Bool=false) where {W, F}
-    c = if exp
-        continuous_exp(DistFix{W, F}, d, pieces, start, stop)
-    else 
-        continuous_linear(DistFix{W, F}, d, pieces, start, stop)
-    end
-    return c
-end
-
-# https://en.wikipedia.org/wiki/Laplace_distribution
-function laplace(t::Type{DistFix{W, F}}, mean::Float64, scale::Float64, start::Float64, stop::Float64) where {W, F}
-    @assert scale > 0
-
-    beta1 = -1/scale
-    e1 = exponential(DistFix{W, F}, beta1, mean, stop)
-
-    beta2 = 1/scale
-    e2 = exponential(DistFix{W, F}, beta2, start, mean)
-
-    ifelse(flip(0.5), e1, e2)
-end
-
-#Helper function that returns exponentials 
-
-function shift_point_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64) where {W, F}
-    DFiP = DistFix{W, F}
-    if alpha == 0
-        unit_exponential(DFiP, beta)
-    else
-        x1 = shift_point_gamma(DFiP, alpha - 1, beta)
-        x2 = uniform(DFiP, 0.0, 1.0)
-        observe(ifelse(flip(1/(1 + 2.0^(-F))), x2 < x1, true))
-        x1
-    end
-end
-
-#https://www.wolframalpha.com/input?i=sum+%28a*epsilon%29%5E2+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1
-function sum_qgp(β::Float64, ϵ::Float64)
-    ans = (1/ϵ - 1)^2 * exp(β*ϵ*(2 + 1/ϵ))
-    ans += (1/ϵ^2)*exp(β)
-    ans += (2/ϵ - 2/ϵ^2 + 1)*exp(β * (1 + ϵ))
-    ans -= exp(β*ϵ)
-    ans -= exp(2*β*ϵ)
-    ans *= ϵ^2
-    ans /= (exp(β * ϵ) - 1)^3
-    ans
-end
-
-#https://www.wolframalpha.com/input?i=sum+%28a*epsilon%29+*+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1
-function sum_agp(β::Float64, ϵ::Float64)
-    ans = (1/ϵ - 1)*exp(β*(1 + ϵ))
-    ans -= exp(β)/ϵ
-    ans += exp(β*ϵ)
-    ans *= ϵ
-    ans /= (exp(β*ϵ) - 1)^2
-    ans
-end
-
-#https://www.wolframalpha.com/input?i=sum+e%5E%28beta+*+epsilon+*+a%29+from+a%3D0+to+a%3D2%5Eb-1
-function sum_gp(β::Float64, ϵ::Float64)
-    ans = (exp(β) - 1) / (exp(β*ϵ) - 1)
-    ans
-end
-
-function sum_pgp(β::Float64, ϵ::Float64, p::Int)
-    if p == 0
-        sum_gp(β, ϵ)
-    elseif p == 1
-        sum_agp(β, ϵ)
-    elseif p == 2
-        sum_qgp(β, ϵ)
-    else
-        sum = 0
-        for i = 0:ϵ:1-ϵ
-            sum += i^p * exp(β*i)
-        end
-        sum
-    end
-end
-
-function n_unit_exponentials(::Type{DistFix{W, F}}, betas::Vector{Float64}) where {W, F}
-    DFiP = DistFix{W, F}
-    l = length(betas)
-    ans = Vector(undef, l)
-    for i in 1:l
-        ans[i] = Vector(undef, W)
-    end
-    for i in 1:W-F
-        for j in 1:l
-            ans[j][i] = false
-        end
-    end
-    for i in F:-1:1
-        for j in 1:l
-            ans[j][i + W - F] = flip(exp(betas[j]/2^i)/(1+exp(betas[j]/2^i)))
-        end
-    end
-    [DFiP(i) for i in ans] 
-end
-
-function exponential_for_gamma(α::Int, β::Float64)::Vector{Float64}
-    if α == 0
-        []
-    elseif α == 1
-        [β, β, 0.0]
-    else
-        v = []
-        for i in 1:α
-            v = vcat(vcat([β], zeros(i-1)), v)
-        end
-
-        vcat(vcat(exponential_for_gamma(α-1, β), [0.0]), v)
-    end
-end
-
-function gamma_constants(α::Int, β::Float64, ϵ::Float64)
-    @syms varint
-    @syms v2
-    if α == 0
-        []
-    else
-        c1 = Float64(sympy.Poly(integrate(varint^α*exp(β*varint), (varint, 0, 1)), varint).coeffs().evalf()[1])
-        c2 = [Float64(i) for i in sympy.Poly(simplify(v2*integrate(varint^(α-1)*exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2).coeffs()]
-        p1 = 0
-        for i in eachindex(c2)
-            p1 += sum_pgp(β, ϵ, length(c2) + 1 - i) * c2[i]
-        end
-        p1 /= c1
-
-        c2 = [Float64(i) for i in sympy.Poly(simplify(v2*integrate(varint^(α-1) * (varint - v2) *exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2).coeffs()]
-        # @show c2
-        # @show sympy.Poly(simplify(v2*integrate(varint^(α-1) * (varint - v2) *exp(β*varint), (varint, v2, v2 + ϵ))/exp(β*v2)), v2)
-        p2 = Vector(undef, α)
-        for i in eachindex(c2)
-            p2[i] = sum_pgp(β, ϵ, length(c2) - i) * c2[i]
-        end
-        # @show p2
-        
-        vcat([p1], p2, gamma_constants(α-1, β, ϵ))
-    end
-end
-
-
-
-
-function unit_gamma(t::Type{DistFix{W, F}}, alpha::Int, beta::Float64; vec_arg=[], constants = [], discrete_bdd=[], coinflips = []) where {W, F}
-    DFiP = DistFix{W, F}
-    if alpha == 0
-        unit_exponential(DFiP, beta)
-    elseif alpha == 1
-        
-        t = (exp(beta*2.0^(-F))*(beta*2.0^(-F) - 1) + 1)*(1 - exp(beta)) / ((1 - exp(beta*2.0^(-F)))*(exp(beta) * (beta - 1) + 1))
-        
-        if (length(coinflips) != 0)
-            coinflip = coinflips[1]
-        else
-            coinflip = flip(t)
-        end
-        
-
-        if (length(vec_arg) != 0)
-            (Y, Z, U) = vec_arg
-        else
-            (Y, Z, U) = n_unit_exponentials(DFiP, [beta, beta, 0.0])
-        end
-        observe(U < Y)
-
-        
-        
-        final = ifelse(coinflip, Z, Y)
-        final
-    else 
-        α = alpha
-        β = beta
-        if (length(vec_arg) != 0)
-            vec_expo = vec_arg
-            
-        else
-            discrete_bdd = Vector(undef, α)
-            constants = gamma_constants(alpha, beta, 1/2^F)
-
-            t = (exp(beta*2.0^(-F))*(beta*2.0^(-F) - 1) + 1)*(1 - exp(beta)) / ((1 - exp(beta*2.0^(-F)))*(exp(beta) * (beta - 1) + 1))
-            f = flip(t)
-
-            count = 0
-            for i in α:-1:1
-                # @show constants
-                l = discrete(DistUInt{max(Int(ceil(log(i))), 1)}, normalize(constants[count + 2:count+i+1]))
-                count = count+i+1
-
-                discrete_bdd[α - i + 1] = l
-            end
-            vec_expo = n_unit_exponentials(DFiP, exponential_for_gamma(α, β))
-            
-        end
-
-        seq = Int(α*(α^2 + 5)/6)
-        # @show constant_flips
-        x1 = unit_gamma(DFiP, alpha-1, beta, vec_arg=vec_expo[1:seq], constants=constants[α + 2:length(constants)], discrete_bdd=discrete_bdd[2:α])
-        x2 = vec_expo[seq + 1]
-        observe(x2 < x1)
-
-        discrete_dist_vec = Vector(undef, α)
-        count = seq+2
-        for i in 1:α 
-            x = vec_expo[count]
-            count+=1
-            for j in 1:α - i
-                observe(vec_expo[count] < x)
-                count+=1
-            end
-            discrete_dist_vec[i] = x
-        end
-
-        # l = discrete(DistUInt{Int(ceil(log(α)))}, normalize(constants[2:α+1]))
-        l = discrete_bdd[1]
-        t = DFiP(0.0)
-        for i in 1:α
-            t = ifelse(prob_equals(l, DistUInt{Int(ceil(log(α)))}(i-1)), discrete_dist_vec[i], t)
-        end
-        
-        ifelse(flip(constants[1]), x1, t)
-    end
-
-end
-
-function normalize(v)
-    l = sum(v)
-    [i/l for i in v]    
-end
-# # TODO: Write tests for the following function
-# function unit_concave(t::Type{DistFix{W, F}}, beta::Float64) where {W, F}
-#     @assert beta <= 0
-#     DFiP = DistFix{W, F}
-#     Y = uniform(DFiP, 0.0, 1.0)
-#     X = unit_exponential(DFiP, beta)
-#     observe((X < Y)| prob_equals(X, Y))
-#     Y
-# end
-
-######################################################################################################
-# bit blasting geometric distribution: https://en.wikipedia.org/wiki/Geometric_distribution
-######################################################################################################
-
-# Creates a geometric distribution over the integers in the range [0, stop-1]
-function geometric(::Type{DistFix{W, F}}, success::Float64, stop::Int) where {W, F}
-    #TODO: use convert
-    @assert ispow2(stop)
-    bits = Int(log2(stop))
-    @assert W - F > bits
-
-    convert(DistFix{W, F}, DistFix{W, 0}(unit_exponential(DistFix{bits+1, bits}, log(1 - success)*2^bits).mantissa))
-end
-
-"""
-   general_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64, ll::Float64, ul::Float64)
-   
-Returns bitblast distribution for the density function x^alpha e^(beta*x)
-"""
-function general_gamma(::Type{DistFix{W, F}}, alpha::Int, beta::Float64, ll::Float64, ul::Float64) where {W, F}
-    @assert ispow2(ul - ll)
-    multiply = Int(log2(ul - ll))
-    start = DistFix{W, F}(ul)
-    new_type = DistFix{W, F + multiply} 
-
-    DistFix{W, F}(unit_gamma(new_type, alpha, beta).mantissa.number.bits) + start
-end
-
-"""
-    bitblast_sample(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-        numpieces::Int, start::Float64, stop::Float64, offset::Float64, width::Float64)
-
-A modified version of the function bitblast that works with the assumption of lower bits being sampled.
-"""
-function bitblast_sample(::Type{DistFix{W,F}}, dist::ContinuousUnivariateDistribution, 
-    numpieces::Int, start::Float64, stop::Float64, offset::Float64, width::Float64) where {W,F}
-
-    # count bits and pieces
-    @assert -(2^(W-F-1)) <= start < stop <= 2^(W-F-1)
-    f_range_bits = log2((stop - start)*2^F)
-    @assert isinteger(f_range_bits) "The number of $(1/2^F)-sized intervals between $start and $stop must be a power of two (not $f_range_bits)."
-    @assert ispow2(numpieces) "Number of pieces must be a power of two (not $numpieces)"
-    intervals_per_piece = (2^Int(f_range_bits))/numpieces
-    bits_per_piece = Int(log2(intervals_per_piece))
-
-    dist = truncated(dist, start, stop)
-
-    total_prob = 0
-    piece_probs = Vector(undef, numpieces) # prob of each piece
-    border_probs = Vector(undef, numpieces) # prob of first and last interval in each piece
-    linear_piece_probs = Vector(undef, numpieces) # prob of each piece if it were linear between end points
-
-    for i=1:numpieces
-        firstinter = start + (i-1)*intervals_per_piece/2^F 
-        lastinter = start + (i)*intervals_per_piece/2^F 
-
-        piece_probs[i] = 0
-        # Warning: A potential source of high runtime
-        for j=1:intervals_per_piece
-            @show firstinter + offset + width + (j-1)/2^F
-            @show firstinter + offset + (j-1)/2^F
-            piece_probs[i] += cdf(dist, firstinter + offset + width + (j-1)/2^F) - cdf(dist, firstinter + offset + (j-1)/2^F)
-        end
-        total_prob += piece_probs[i]
-
-        border_probs[i] = [cdf(dist, firstinter + offset+width) - cdf(dist, firstinter + offset), 
-                cdf(dist, lastinter -1/2^F + offset+width) - cdf(dist, lastinter - 1/2^F + offset)]
-        linear_piece_probs[i] = (border_probs[i][1] + border_probs[i][2])/2 * 2^(bits_per_piece)
-    end
-
-
-    @show piece_probs
-    @show total_prob
-
-    PieceChoice = DistUInt{max(1,Int(log2(numpieces)))}
-    piecechoice = discrete(PieceChoice, piece_probs ./ total_prob) # selector variable for pieces
-    slope_flips = Vector(undef, numpieces)
-
-    isdecreasing = Vector(undef, numpieces)
-    for i=numpieces:-1:1
-        iszero(linear_piece_probs[i]) && iszero(piece_probs[i]) && continue
-        a = border_probs[i][1]/linear_piece_probs[i]
-        isdecreasing[i] = a > 1/2^bits_per_piece
-        if isdecreasing[i]
-            slope_flips[i] = flip(2-a*2^bits_per_piece)
-        else
-            slope_flips[i] = flip(a*2^bits_per_piece)
-        end  
-    end
-
-    unif = uniform(DistFix{W,F},  bits_per_piece)
-    tria = triangle(DistFix{W,F}, bits_per_piece)
-    z = nothing
-    for i=1:numpieces
-        iszero(linear_piece_probs[i]) && continue
-        firstinterval = DistFix{W,F}(start + (i-1)*2^bits_per_piece/2^F)
-        lastinterval = DistFix{W,F}(start + (i*2^bits_per_piece-1)/2^F)
-        linear_dist = 
-            if isdecreasing[i]
-                (ifelse(slope_flips[i], 
-                    (firstinterval + unif), 
-                    (lastinterval - tria)))
-            else
-                firstinterval + ifelse(slope_flips[i], unif, tria)
-            end
-        z = if isnothing(z)
-                linear_dist
-            else
-                ifelse(prob_equals(piecechoice, PieceChoice(i-1)), linear_dist, z)  
-            end
-    end
-    return z
-end