numpyで全4桁パターンから10を作る計算式を列挙する

はじめに

4桁の数字を使って、四則計算を組み合わせて10を作る計算は、昔から暇つぶしのクイズなどで出題されている。
今回は、この出題の全パターンについて10を求める計算式を列挙するアルゴリズムを、pythonとくにnumpyを使用して実装することを考える。
(アルゴリズムについては先駆者がより効率的なものを紹介していると思う。本記事はnumpyでこういう風にも書けるよ、という紹介記事である。)

基本方針は以下。

• 全出題パターンを列挙する
• 各出題パターンについて、計算式の探索を行う
• 解が存在するパターンについて、計算式を生成する

それぞれ作っていく前に、前準備を行う。

(前準備)被りのないパターンの数の計算

※ちょっと高校数学な話。苦手な方は読み飛ばしても大丈夫です。

4桁から10を計算するパズルにおいて、順番が入れ替わっただけの4桁は同じ計算方法で解ける問題となる。これらは重複問題といえるので除去したい。

4桁からなる全パターンは$10^4$パターン存在するが、その中から被りのない問題は何パターンあるだろうか？
これは、重複組み合わせと呼ぶ数学の問題に帰着する。

以下で、そのパターン数を得る公式を導出してみよう。高校数学のおさらいである。

考え方

9個のNと、ピックしたい桁数だけの個数(例として2個)のSを並べて以下のような文字の列を作る。

N,N,N,N,...,N,S,S

この文字の列は、どんなNとSの並びであっても指定桁数(桁の数値は順不同)のパターンと1対1対応するように解読方法を定めることができる。

その解読方法とは、

• Sの左側にあるNの数をカウントして桁の数を求める。
• 表れた桁をSの数だけ集める。並びは順不同

である。たとえば、

S,N,N,N,S,...,N

としたとき、これは[0,3]を
選択したことに対応させる、と考える。

N,N,S,S,N,...,N

これは、[2,2]に対応する。

この表記ルールにおいて、NとSの並びを変えると[00]から[99]までの全ての桁を表現できる。また、

S,N,S,N,N,...,N
は[0,1]と並びによらず対応するルールのため、[01]と[10]はともにこの1つのパターンにまとめられる。

そのため、同じ4桁を並び替えただけのパターンをダブルカウントしなかったときの全パターンの個数は、先ほど定義したNとS(4個)の並び替えで可能なパターンの総数と同じだといえる。

さて、総数を求めてみよう。
Nは最大の番号である9個、Sは桁数だけである。その組み合わせは「9+桁数個の玉の中から桁数個の玉をSとして選ぶ」という手順で作れるので、単なる組み合わせ計算になる。

組み合わせの公式は、

$${ }_nC_r=\frac{n!}{(n-r)!r!}$$

であることを利用しよう。

2桁の場合: 11!/(9!2!)=55パターン

3桁の場合: 12!/(9!3!)=220パターン

4桁の場合: 13!/(9!4!)=715パターン

求まった。

検算

0,1,2,3,4の5つの数字を組み合わせて順不同の2桁を選ぶとき、パターンを列挙しよう。

[00][01][02][03][04][11][12][13][14][22][23][24][33][34][44]
合計、5+4+3+2+1=15通りある。
先ほどの計算で考えると、Nは最大の数である4、Sが2個なので、

(4+2)!/4!2!=15。

確かに、この計算で正しく求まっているようだ。

ということで、重複を含めると10000パターンあった四桁の数字の組み合わせは、実は715パターン見れば充分だということになる。

まあ、この計算は今回実は直接実装上必須ではないのだが、次の節で全パターン列挙する関数でパターン数に過不足がないかチェックできるので、知っていて損はない話だと思う。

(前準備2)組み合わせ、順列の列挙処理

4桁のパターンは715通りあることは分かったが、今度はこれを全パターン網羅する処理を考える。

順番入れ替えで同じパターンにならない4桁の数字の列挙は、以下のようにすれば求まる。

• 必ず左の桁以上になる数字を右にくっつけるように4桁生成する(同じ数字を使ってはいけない場合は、必ず左の桁+1以上)

これを網羅的に実現するコードは、愚直に書けば以下になる。

    def straightCombination(n=10, r=4, s=[], overlap=False):
        # nCrで列挙される全パターンのうち、sで指定した数字までを左詰め固定したリストを得る
        if len(s)<r:
            if len(s)==0:
                s0=0 if overlap else -1
            else:
                s0=s[-1]
            q=[]
            if not overlap:
                s0+=1
            for t in range(s0,n):
                # 再帰的に、入力した桁を左に詰めした条件でのパターンを得る
                q+=straightCombination(n, r, s+[t])
            # これで得られるのは、入力桁の条件であり得る全パターン
            return q
        else:
            # 欲しい桁数に到達したらリストを返す
            return [s]

この実装は桁数の数だけネストさせたforループを回すことで、欲しいパターンを片端から全探査する処理になっている。

一方、pythonはforループは速くないので、こういうネストの深いforで解決する実装はあまり望ましくない。

いろいろなパターンをforを使って個別計算する代わりに、numpyの行列演算はこういった処理に向いている。
行列演算を使った実装例は以下。

    def makeCombination(n, r, allowOverlap=True):
        # 最初の一桁は全数字が該当
        q=np.arange(n)[None]
        for _ in range(1, r):
            # 0,1,2,3...を0,0,0,...,1,1,1,...という風に拡張する
            q_ex=q.repeat(n,axis=-1)
            # 0,1,2,3,...,0,1,2,3,...という並びを作る
            p=np.tile(np.arange(n),q.shape[1])
            # qとpをつなげて桁を1つ増やしたパターン列挙にする
            qp=np.concatenate((q_ex,p[None]),axis=0)
            # 欲しいパターン条件は、1つ前の桁より最後の桁の方が大きい場合なので抜粋
            if allowOverlap:
                q=qp[:,(qp[-2]<=qp[-1])]
            else:
                q=qp[:,(qp[-2]<qp[-1])]
            # これをr回繰り返してr桁にする
        return q

パターンを列挙するコツは、
4桁から3桁を選ぶパターンを作る場合を例にすると以下である。

2桁目まで選んだとき、そのパターンが(0,1),(0,2),(0,3),(1,2),(1,3),(2,3)だったとして、
(0,1,0)から(0,1,3)までの4通り、(0,2,0)から(0,2,3)までの4通り、と単純に最終桁に0,1,2,3をあてはめてしまう。
得られるパターンは合計6*4=24通りになる。

この中で、桁が昇順になってるパターンを抜き出すと
(0,1,2),(0,1,3),(1,2,3)の3通りが抽出され、3桁までを選んだパターンが得られる。あとは、この作業を冒頭に戻って繰り返せば何桁でも生成できる、という風に作られたアルゴリズムである。

新しい桁の抽出ルールを以下のように「今までの桁に使われていない数字」にすれば、順列の列挙関数の出来上がりである。

    def makePermutation(n, r):
        q=np.arange(n)[None]
        for _ in range(1, r):
            q_ex=q.repeat(n,axis=-1)
            p=np.tile(np.arange(n),q.shape[1])
            qp=np.concatenate((q_ex,p[None]),axis=0)
            for i in range(qp.shape[0]-1):
                qp=qp[:,qp[i]!=qp[-1]]
            q=qp
        return q

全出題パターンの列挙

これは前準備2で定義したmakeCombinationを使うだけである。

numbers=makeCombination(10,4)

これで全715パターンが生成される。

解答の探索

四則演算のみが許されているので、計算に使う数字a,bが決まっているなら、計算式は以下の6パターンのいずれかが使える。

1. a+b
2. a-b
3. b-a
4. a/b
5. b/a
6. a*b

aとbは、最初4桁から2桁を選んで計算し、計算結果として1つの数字を得る。

流れ的には、

(1) 4桁から2桁を選ぶ
(2) 6種類の演算パターンのどれかを使って計算する。
(3) 選ばなかった2桁と計算で作った1つ、合計3つの数字から2つを選ぶ
(4) (2)と同じ処理
(5) 残り2つの数字を使って(2)をおこなう

という手順で、ある適当に四則演算を組み合わせたときの計算結果が得られる。
これを全パターン網羅すれば、指定した4桁の数字から10は作れるのかわかる。

なお、全パターン試すには桁数字の重複を考えなければ

$$ { }_4C_2 * 6 * { }
_3C_2 * 6 * { }_2C_2 * 6$$

通りの試行回数が必要である。

上記手法をpython, とくにnumpyを使って実装してみよう。

まずは、計算に使う2桁を選ぶ全パターンを生成する処理

    def permutePatterns(num):
        assert num>1
        # 引数に指定した桁数の数列に対して、
        # 四則演算に使う2桁を昇順に配置し、
        # その末尾に残りの桁を昇順に配置した際の
        # 全パターンを列挙する
        # それには、まず桁数を自由に並べてできる全パターンを列挙する
        p=makePermutation(num,num)
        # 最初の2桁は演算に使用する。その2桁が昇順になっているパターンのみ抽出
        p=p[:,p[0]<p[1]]
        # 抽出されたパターンのなかから、さらに
        # 残りの桁が昇順に並んでいるパターンを抽出
        # なお、最初に順列になっている全パターンから始めているので、残りの桁と前半2桁が重複する心配はない
        for i in range(2,num-1):
            p=p[:,p[i]<p[i+1]]
        return p

使う四則演算パターンを全列挙する関数

    def divide(a,b):
        # 0除算を迂回した割り算
        q=np.zeros(b.shape)
        q[b!=0]=a[b!=0]/b[b!=0]
        # 0割結果はnanを代入
        q[b==0]=np.nan
        return q
    
    def patternCalc(a,b,pat=None):
        # 指定した二項に対して、全6パターンの四則演算結果を返す
        ress=[
            a+b,
            a-b,
            b-a,
            a*b,
            divide(a,b),
            divide(b,a)
        ]
        if pat is None:
            # patを指定しない場合は、6通り全部返す
            return np.stack(ress)
        # patで指定した計算結果のみ返す
        return ress[pat]

これを使って、任意の桁の数字列から四則計算で10を作れるパターンを全探索する関数を作ると以下。

    def processCalc(mat):
        # 四則演算に使う二項を全パターン×全6通りの計算結果を得る
        # 1パターンとは、例として(3,4,5,6)のうち(3,4)を四則演算に使う数字に選び、かけ算を適用すると
        # (3*4,5,6)=(12,5,6)と3項が残る
        # これを全パターン列挙した場合のパターン数は
        # 入力桁数が4項なら4P2*6パターン=36通りである
        # matは第一成分が桁並び、第二成分が上記でいう4P2に相当するパターン、第三成分が探索に指定した4桁の数字のバリエーションに相当した3階のテンソルである
        # 前半2項を6通り計算して、6通りの1項の結果を得る。第一成分は項に対応させて1、第二成分が6通りの四則演算結果となる。
        a=patternCalc(mat[0],mat[1])[None]
        # 選ばなかった残りの項を、第二成分(四則演算の各結果)に沿って同じ値になるように拡張する
        ar=mat[2:,None].repeat(6,axis=1)
        # 四則演算で得られた1桁と、演算しなかった残りの桁を結合して、入力より1項少ない状態で出力する。
        res=np.concatenate((a,ar),axis=0)
        return res
    
    def nearby(v,e):
        return np.abs(v-e)<1e-5
    
    # 
    def getJudgeMatrix(n, dim=4):
        # 全数字の組み合わせパターンを列挙
        numbers=makeCombination(n,dim)
        mnum=numbers.shape[1]
        r=numbers
        for k in range(dim,1,-1):
            # 2桁を選ぶ組み合わせ列挙
            r=r[permutePatterns(k)]
            # 四則演算全パターン列挙
            r=processCalc(r)
            # 四則演算により、項は1つ減る
            # 項が1つ減った計算の組み合わせを、全パターン第二成分に押し込める
            r=r.reshape(k-1,-1)
            # パターン数は指数関数的に増えていく様子は、shapeを見れば確認できる
            print(r.shape)
        # 最後の1項になるとループは終了
        # 列成分がmnumずつ括ることで、行成分は項の数が1つになるまでの計算に使った桁の順番と使った四則計算式の種類の全パターンの数だけの次元数になる
        r=r.reshape(-1,mnum)
        # パターンの通し番号を用意する
        idxs=np.arange(r.shape[0])[:,None].repeat(mnum,axis=1)
        # 足して10にならなかったパターンに対応する番号は、最後のパターンの番号よりも大きくする
        # 計算結果は浮動小数点であり丸め誤差が載るので、10との距離が非常に近い値も、足して10になるとみなす
        idxs[nearby(r,10)==False]=r.shape[0]
        # まだ通し番号がついているパターンのなかで、最も若い番号を代表計算結果とする。番号が若い順のほうが加算減算優先のため簡単な計算式が出やすい。
        minidx=idxs.min(axis=0)
        # 逆に、最も大きい番号側は難しい式が出やすい。こっちも参考例としてピックアップする
        # maxが有効値になるように、足して10にならないパターンに対応する番号を-1に移し替える
        idxs[idxs==r.shape[0]]=-1
        maxidx=idxs.max(axis=0)
        # 数字の組み合わせごとに、足して10になるパターンの数をカウントする
        ans_cnt=np.where(idxs!=-1,1,0).sum(axis=0)
        # 解が存在する数字の対応する場所を列挙
        ans_num=np.where(minidx<r.shape[0])
        # 最大最小と総数のリストから、解が存在するパターンのみ抽出
        ans_min=minidx[ans_num]
        ans_max=maxidx[ans_num]
        ans_cnt=ans_cnt[ans_num]
        ans_idx=numbers[:,ans_num[0]]
        # 結合して並び替えて終了
        ansvi=np.stack((*ans_idx,ans_min,ans_max,ans_cnt))
        return ansvi.T

これは、例えばある特定の4桁(3,2,4,5)から10を作る計算の流れが以下だとして、それ全数字パターン、全計算経路で全探索している処理に相当する。

• (3,4,2,5)と並び替える
• 先頭2つを選んで、4-3=1とする。使った項を消して計算結果に置き換え(1,2,5)をつくる
• (1,5,2)と並び替えて先頭２つを1x5と計算し置き換え、(5,2)にする
• (5,2)を5x2=10として10を生成

こうすると、計算に使う数字の並びや掛け算や引き算を使う順序を変えると10にならなくなるパターンが無数に出てくる。上で定義した関数は、組み合わせが10になるならないに問わず全計算パターンをベクトル演算で並列的に網羅して、あとで10になる組み合わせを抽出する処理をすることで全解を抽出している。

計算式の構築

計算式を作る方法は先ほど行った例と手順を使えば、

• (3,4,2,5)と並び替える
• 先頭2つを選んで、"4-3"を実施したので、("4-3",2,5)をつくる
• ("4-3",5,2)と並び替えて先頭２つを"4-3"x5した。かけ算より引き算を優先する必要があるので括弧をつけて("(4-3)x5",2)とする
• (4-3)x5x2、が最終的な答えになる

ようは、10を作ったパターンを計算式にするには、各ステップでどの数字をどの計算に通したかがわかれば構築できる。先ほど全計算パターンを列挙したとあるので、そのパターンに通し番号をつけてあれば、後から計算式構築ができる、ということになる。

先ほどの関数後半の

        # パターンの通し番号を用意する
        idxs=np.arange(r.shape[0])[:,None].repeat(mnum,axis=1)

から

        # maxが有効値になるように、足して10にならないパターンに対応する番号を-1に移し替える
        idxs[idxs==r.shape[0]]=-1

は、足して10にならないパターンを捨てて10になるパターンの番号を抽出する処理に対応している。

あとはパターン番号に沿って、先ほどと同じ計算探索順で今度は文字列が出るような
関数を以下のように用意すれば、計算式が列挙される。

    def splitNumber(num,mds):
        res=np.zeros((len(mds)), dtype=int)
        for k, md in enumerate(mds):
            res[k]=num%md
            num=num//md
        return res
    
    # make "a" to "(a)", if necessary.
    def termToStr(term, lv):
        # termは数字、またはtupleが想定。数字の場合はndimをメンバにもつndarrayか、もしくは[]が使えないはずなので、以下条件で見分ける
        if hasattr(term, "ndim") or not hasattr(term, "__getitem__"):
            # 数字だった場合はいったんterm_lv=0のtupleにする
            term=("%d"%term,0)
        term_val, term_lv=term
        # termlvより次の演算のlvの方が低いとき、そのまま次の項をつなげると、項内の計算式末尾の項と次の項で先に計算されてしまう。そうならないように、ここまでの項の内容を()でくくる
        if term_lv>lv:
            return "("+term_val+")"
        return term_val
    
    
    def toPhrase(a, b, m):
        ci="+--*//"    # パターン番号に対応する四則演算の記号
        co="001001" # 二項を使用する順番
        lvsl="444222" # 演算子のレベル(左の項)
        lvsr="433211" # (右の項)
        lvl,lvr=[int(lvs[m]) for lvs in (lvsl,lvsr)]
        pi=0 if co[m]=="0" else 1
        qi=1-pi
        # 2項と演算記号を、必要に応じて括弧で括ってから並べる
        ab=a,b
        p=ab[pi]
        q=ab[qi]
        return termToStr(p,lvl)+ci[m]+termToStr(q,lvr), lvl+(lvl%2)
    
    def parse(v, m, d):
        vi=v
        # パターン番号から、対応する計算順序と用いた演算記号の組み合わせパターンの番号に分解する
        mods=[]
        for k in range(d, 1, -1):
            mods+=[permutePatterns(k).shape[1],6]
        mi=splitNumber(m,mods)
        # パターン番号は1度二項演算を行うごとに2個消費する。d桁の数字が1項の数字までまとまるまで、ループをまわす
        for k in range(d,1,-1):
            # k項の数字から2項を選ぶ組み合わせパターンリストを取得
            ps=permutePatterns(k)
            # 組み合わせパターン番号に対応する項の並びを取得
            vi=[vi[p] for p in ps[:,mi[(d-k)*2]]]
            # 先頭2項を計算に使った計算経過tupleを生成し、[0]の項の内容を置き換え
            vi[0]=toPhrase(*vi[:2],mi[(d-k)*2+1])
            # [1]は計算で消費したので、取り除く
            vi.pop(1)
        # 最後に残った1項が最終的な計算式になる
        # [0]に格納された計算経過tupleのうち、計算式を格納しているのは[0]なので、[0][0]を返す
        return vi[0][0]

ここで、括弧をつけるタイミングについて補足する。
演算式a+b,axbなどに新しい項cを計算する際、aとbの演算結果を保存するように記述すると

• a+b+c: 括弧不要
• axb+c: 括弧不要
• (a+b)xc: 括弧必要
• c-(a+b): 括弧必要
• c/(axb): 括弧必要

という例で示す形になる。演算子記述の優先度を整理すると

• lv0: 演算子による結合がなく、括弧つける必要なし
• lv1: 割り算の右辺。+-x/で結ばれた項は括弧をつけてから配置
• lv2: 掛け算もしくは割り算の左辺。+-で結ばれた項は括弧をつけてから配置
• lv3: 引き算の右辺。+-で結ばれた項は括弧をつけてから配置
• lv4: 足し算もしくは引き算の左辺

というレベルを定義すれば判断できる。

たとえばa+bはレベル4の項であり、xcはレベル2の演算である。そのままa+bxcと結合するとbxcが先に計算されてしまうので、括弧でくくる。

これを判断しているのがtoPhraseおよびtermToStr関数である。

実践

あとは、ここまで組み上げたパーツを直列実行するだけである。
以下は、4桁から10をつくることができる組み合わせを全探索し回答の数(同じ数字の重複は許容する)を出し、その中で足し算優先およびかけ算優先で計算式を構築した場合の解答例を標準出力に出すプログラムである。

    # 4桁で10をつくる計算式を探索する処理
    dim=4
    getTens=getJudgeMatrix(10, dim)
    # 探索が終わったので、各結果について10をもとめる計算式に変換
    for g in getTens:
        _num=g[:dim]
        _min,_max,_cnt=g[dim:]
        nArr=_num
        nStr="".join(["%d"%d for d in nArr])
        s=nStr+": "
        s+="%15s"%parse(_num, _min, dim)
        if _min!=_max:
            t="%15s"%parse(_num, _max, dim)
            s+=", "+t
        if _cnt>2:
            s+=", (+%4d ans.)." % (_cnt-2)
        print(s)

実行結果は長いので省略。全部見せてしまうと、我々の娯楽が一つ消えることになるので、実行にはそういう意味で注意である。

あとはお遊び

難しい出題パターン

このプログラムを動かすと、10を作るのは可能だが難しい4桁がわかる。
いくつか難しいパターンを例示しよう。

1158
1337
3478

いずれも、解は単独である上に分数を経由するという意味で、難しいパターンだといえる。分数を経由する演算、というだけでパターンはかなり絞られるので大ヒントを出してしまったかもしれないが。

トリッキーな解探し

このプログラムはmax側の番号の探索パターンに除算を割り当てているので、比較的トリッキーな式になりがちである。
個人的に面白いと思った解は以下。

• 2789: 素直な解=2+7+9-8, トリッキーな解: (8x9-2)/7
  • 敢えて掛け算で大きな値に飛ばしている点に好感が持てる
• 1589: 素直な解=(1+9)/5+8, トリッキーな解=8/(9/5-1)
  • 分数使用

ただ、眺めてみるとだいたいワンパターンな組み合わせが多く、トリッキーな答えにも飽きがきやすい。解を眺めていると自力で計算する気が起きなくなってくる点、やはりこの手の問題は計算機任せで自動で求めるべきではないのかもしれない。

おわりに

今回は4桁から10を求める問題の解の全列挙という問題を通して、実はnumpyで組み合わせ列挙関数をどう組むかに重きを置いて練習がてら組んだものだったが、以外と数学的側面と実装工夫的側面で面白かったので結構時間を使ってしまった。
numpyでの順列、組み合わせパターン全列挙処理の実装は他の問題にも役立ちそうなため、この手法は覚えておきたいと思う。

なお、このプログラムは4桁でない場合でも解答が求まるように組んであるが、5桁以上では組み合わせの総数が指数関数的に大きくなるため消費メモリが激しく増加する。そのため、この方法は桁数が大きい処理にはお勧めできない点、注意である。