Skip to content

Instantly share code, notes, and snippets.

@carlobaldassi
Created June 28, 2012 15:50
Show Gist options
  • Save carlobaldassi/3012133 to your computer and use it in GitHub Desktop.
Save carlobaldassi/3012133 to your computer and use it in GitHub Desktop.
fft_powers tests
load("fft_powers.jl")
@assert fft_powers_approx([2,3,5], 37) == 40
@assert nextprod([2,3,5], 37) == 40
function bench_nextprod(ps, x)
tm = 1000
t = Array(Float64, tm)
for i = 1:tm
t[i] = @elapsed nextprod(ps, x)
end
println(" nextprod : $(mean(t)) ± $(std(t)/sqrt(tm))")
end
function bench_approx(ps, x)
tm = 1000
t = Array(Float64, tm);
for i = 1:tm
t[i] = @elapsed fft_powers_approx(ps, x)
end
println(" fft_powers_approx: $(mean(t)) ± $(std(t)/sqrt(tm))")
end
function bench(ps, x)
println("benchmark ps=$ps x=$x")
bench_approx(ps, x)
bench_nextprod(ps, x)
end
x = randi(200_000)
for i = 3:7
ps = PRIMES[1:i]
bench(ps, x)
end
load("sparse.jl", "glpk.jl")
let ps_opts=GLPSimplexParam(), opts=GLPIntoptParam(), lp=GLPProb()
opts["msg_lev"] = GLP_MSG_ERR
opts["tol_int"] = 1e-9
ps_opts["msg_lev"] = GLP_MSG_ERR
opts["br_tech"] = GLP_BR_FFV
global fft_powers_approx
function fft_powers_approx(ps::Vector, x::Real)
glp_erase_prob(lp)
glp_set_obj_dir(lp, GLP_MIN)
n = length(ps)
f0 = -log(x)
f = log(ps)
glp_set_obj_coef(lp, 0, f0)
glp_add_cols(lp, n)
for c = 1 : n
glp_set_obj_coef(lp, c, f[c])
glp_set_col_kind(lp, c, GLP_IV)
glp_set_col_bnds(lp, c, GLP_LO, 0., 0.)
end
glp_add_rows(lp, 1)
glp_set_mat_row(lp, 1, n, Int32[1:n], f)
glp_set_row_bnds(lp, 1, GLP_LO, -f0, 0.)
ret = glp_simplex(lp, ps_opts)
if ret != 0
error("fft_powers_approx failed")
end
ret = glp_intopt(lp, opts)
if ret != 0
error("fft_powers_approx failed")
end
sol = [ glp_mip_col_val(lp, c) for c = 1:n ]
return int(prod(ps .^ sol))
end
end
load("fft_powers.jl")
ps = PRIMES[1:4]
for i=1:1000_000
x = nextprod(ps, i)
@assert x >= i
y = fft_powers_approx(ps, i)
@assert y >= i
if (x != y)
println("diff: $i $x $y $(y-x)")
end
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment