求素数幂的公倍数 (2, 3, 5, 7...) > N 并最小化指数

Question

我这样做的原因是找到一个比 N 大的最接近的数，它是质数幂的公倍数 --- 以便能够使用 FFTW。

据我了解，这是一个 optimization/linear 编程问题。我设法用数学公式表示：

minimize(a log(2) + b log(3) + c log(5) + ...)

Constraints::
  a log(2) + b log(3) + c log(5) + ..... >= log (N)  (known parameter: N)
  {a, b, c ...} >= 0      (i.e. are positive)
  {a, b, c ...} % 1 == 0  (i.e. are integers)

其中 a、b、c 是指数并且是整数。

我想用数字实现同样的东西。最好使用 scipy.optimize.

编辑：我根据评论稍微修改了方程式。即使不存在实现，算法也会有所帮助。

Answer 1

根据@sascha，a conversation in Scipy-Dev mailing list scipy 尚未实现混合整数线性规划 (MILP a.k.a.MIP) 求解器。

The Python wiki contains a list of LP / MILP python packages. I ended up using pulp and its built-in solver. By modifying the code from a similar problem and examples found in a paper，解决方案如下所示。

#!/usr/bin/env python
"""Find the closest multiple of prime powers greater than N."""
from pulp import *
import numpy as np


N = 1000
bases = [2, 3, 5, 7]
debug = True

prob = LpProblem("FFTW Gridsize Problem", LpMinimize)

exp_max = np.ceil(np.log2(N))
exps = LpVariable.dicts('exponent', bases, 0, exp_max, cat=LpInteger)
log_N_new = LpVariable("log_grid_size", 0)
log_N_new = lpSum(exp * np.log(base) for exp, base in zip(exps.values(), bases))

prob += log_N_new  # Target to be minimized
# Subject to:
prob += log_N_new >= np.log(N), "T1"

if debug:
    print('bases =', bases)
    print('exponents =', exps)
    print('log_N_new =', log_N_new)
    prob.writeLP("OptimizationModel.lp")

prob.solve()

if debug:
    print("Status:", LpStatus[prob.status])
    for v in prob.variables():
        print(v.name, "=", v.varValue)

N_new = np.prod(
    [base**v.varValue for base, v in zip(bases, prob.variables())])
print(N_new)

求素数幂的公倍数 (2, 3, 5, 7...) > N 并最小化指数

Find common multiple of prime powers (2, 3, 5, 7...) > N and minimize the exponents

python

optimization

scipy

linear-programming

integer-programming