23の倍数で各桁の和が23の数を11¹²まで数える

再帰関数＋lru_cache

[P294: Sum of digits - experience #23] (https://projecteuler.net/problem=294) は表題のように簡単そうな問題ですが、数がとんでもなく大きいのでギブアップでした。

まず普通に再帰関数で書くとこんな感じですが時間がかかって問題外。

def S(nd,m,ds): #桁 10**nd, m: 23で割った余り, ds: 桁の和
  if nd == 0:
    return 1 if (m % 23 == 0) and ds == 23 else 0
  ret = 0
  for d in range(min(23-ds+1,10)):
    ret += S(nd-1,(m*10+d)%23,ds+d)
  return ret

ND = 42
print((f"ND = {ND}, answer = {S(ND,0,0)})

そこでお約束のlru_cacheを設定すると劇的に改善して11⁴までは何とか1分以内で終了。でも先は遠い、、、

from functools import lru_cache

@lru_cache(maxsize=None)

動的計画法（1条件）

この手の問題は動的計画法かなと思いましたが、あまり得意ではないので、まず条件の１つの各桁の和が23のみを書きました。

import numpy as np
n, D = 9, 23
dp = np.array([0] * (n+1)*(D+1)).reshape(n+1,D+1)
dp[0, 0] = 1

for i in range(n):
  for j in range(D+1):
    for k in range(min(D-j+1,10)):
        dp[i+1,j+k] += dp[i,j]

print(dp[n,D])

動的計画法（2条件）

この問題では以下の2つの条件があるので、

各桁の和が23
23で割り切れる

DP tableを(23+1 x 23+1)の2次元mapにして、X軸で23で割った余り、Y軸で桁の数の合計を、各桁ごとにカウントするようにしました。計算には一つ前の桁の値しか使わないので、2つのmapを交互に使うようにしました。(使った後で0でクリアしてます)
でもこれも1分で解けるのはk=11⁴までで、今回はちょっとお手上げでした。

import numpy as np

def clearMap(map):
  for i in range(len(map[0])):
    for j in range(len(map)):
      map[i,j] = 0  # clear next table for reuse 

def dpMap(n, D, M):
  BF, D1, D1 = 2, D+1, D+1   # 桁、桁の和、余り、
  dp = np.array([0] * BF*D1*D1).reshape(BF, D1, D1)
  dp[0, 0, 0] = 1
  for i in range(n):  # 桁
    clearMap(dp[(i+1)%2]) # 次のmapを0でクリア
    for d in range(D1):  # 桁の和
      for m in range(D1):  # 余り
        for k in range(min(D-d+1,10)):  # i+1番目の桁の数 (d+k<=D)
          dp[(i+1)%2,d+k,(m*10+k)%D] += dp[i%2,d,m]
          dp[(i+1)%2,d+k,(m*10+k)%D] %= M
  return dp[(i+1)%2,D,0]

n = 42
ans = dpMap(n, 23, 10**9)
print(f"n = {n}, ans = {ans}")