fft_powers tests

fft_benchmarks.jl

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

fft_powers.jl

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

fft_test.jl

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