在蔡炎龍老師的貼文上看到，今年 2025 是完美平方 (perfect square) 年
自然想起了一個問題：那其他完美平方年是多少呢？
先設定一個範圍，就考慮 20 世紀後吧 ( ≥ 1900 A.C)。
老實說我數學不好，要我去算出來可能不太行，但應該是可以寫程式去快速找出所有完美平方年，順手寫了一下，也把過程記錄在這個 blog 。
定義問題
基本上我們可以這樣定義完美平方年：
其實一但有了這個定義，要寫程式找出所有可能的解就很簡單了，這邊簡單用 Python 解一下：
from itertools import product
solutions = []
for x, y in product(range(19, 100), range(0, 100)):
if (x + y)**2 == (100*x+y):
solutions.append((x, y))
print(solutions) # [(20, 25), (30, 25), (98, 1)]碰！我們就可以知道 20 世紀後的完美平方年有 2025 、3025 跟 9801。
但我發現這個題目用 ORTools 寫也是個不錯的練習，就來試試看用它來解看看。
ORTools
這邊我會用到 ORTools 裡的 SAT Solver。
起初我想說簡單的寫法不就是像這樣：
from ortools.sat.python import cp_model
model = cp_model.CpModel()
var_x = model.new_int_var(19, 99, "var_x")
var_y = model.new_int_var(0, 99, "var_y")
model.add((100*var_x + var_y) == (var_x + var_y) * (var_x + var_y))結果馬上報錯：
原因在於 SAT solver 不容許 (x + y) * (x + y) 這樣的 expression。
這時候就可以用上一個技巧，引入額外的 proxy variable 去解決這個問題：
model = cp_model.CpModel()
var_x = model.new_int_var(19, 99, "var_x")
var_y = model.new_int_var(0, 99, "var_y")
var_sum = model.new_int_var(0, 198, "var_sum")
var_sq = model.new_int_var(1900, 9999, "var_sq")這邊我引入了 var_sum 跟 var_sq 這兩個 proxy variables ，把原本兩個變數的問題進行擴增，然後再加上額外的條件：
model.add(var_sum == (var_x + var_y))
model.add_multiplication_equality(var_sq, var_sum, var_sum)
model.add(var_sq == 100*var_x + var_y)var_sum == (var_x + var_y)去限制var_sum必須等於兩變數總和。model.add_multiplication_equality(...)去限制var_sq必須等於var_sum * var_sum。- 最後再讓
var_sq必須等於100*var_x + var_y。
簡單的來說，就是透過適當的限制條件在這些 proxy variable 上，讓這個擴增問題的解變成跟原問題相同，進而繞過 SAT Solver 的一些限制。
最後來檢查一下 model：
model.validate()
# ""如果看到空字串的話代表沒問題，但如果有 error message 就可能要再利用上面的 proxy variable 的技巧改寫去修正。
最後透過 solution callback 就可以拿到所有的解了：
# declare callback
class SolutionCollector(cp_model.CpSolverSolutionCallback):
def __init__(self, variables):
cp_model.CpSolverSolutionCallback.__init__(self)
self._variables = variables[:]
self._solutions = []
def on_solution_callback(self):
self._solutions.append(
tuple(
self.value(var)
for var in self._variables
)
)
@property
def solutions(self):
return self._solutions[:]
# create solver
solver = cp_model.CpSolver()
solver.parameters.enumerate_all_solutions = True
collector = SolutionCollector([var_x, var_y, var_sq])
# solve the model
status = solver.solve(model, solution_callback=collector)
# print solutions
if status in [cp_model.FEASIBLE, cp_model.OPTIMAL]:
print(collector.solutions)
# [(20, 25, 2025), (30, 25, 3025), (98, 1, 9801)]結語
新年剛過待業在家，剛好看到有趣的問題，而且也能藉此介紹一個在使用 ORTools 時很實用的 proxy variable 技巧。我自己的經驗是這在求解更複雜的 SAT 問題時，常常可以繞過 SAT solver 的限制。如果是想透過 ORTools 解最佳化問題的話，這是個很實用的技巧。
不過我想我們大家應該都看不到下一個完美平方年了 (3025) 🤣
Happy Python Programming!