「カット・アンド・スクエア」問題 別解

　CodeIQにて「カット・アンド・スクエア」問題というのを出題しました。

　問題文はこちら。解法はこちらです。

　以下では、解法pdfで言及した別解について解説します。

 

　上 n/2 桁の数を a、下 n/2 桁の数を b とおくと、題意は次のようになります。

10n/2 a + b = a2 + b2

　ここで、A = 2a - 10n/2，B = 2b - 1 と変数変換を行うと、題意の式は次のようになります。

A2 + B2 = 10n + 1

　整数を２平方和で表す問題に帰着されます。

 

　この問題を解くときのポイントは、「①10n + 1 を素因数分解する」「②素因数分解した素因数のそれぞれを２平方和で表す」の２点です。

 

n = 6 の場合

　まず「①10n + 1 を素因数分解する」ことを考えます。単純には試し割りのアルゴリズムが有効です。コードは詳細します。106 + 1 = 101 × 9901 となります。

　次に「②素因数分解した素因数のそれぞれを２平方和で表す」ことを考えます。例えば x2 + y2 = 101 となる整数組 (x, y) を見つけるということです。これも単純には、x の値をまず一つ固定し、y = √(101 - x2) が整数となるかチェックするという方法が考えられます。結果、102 + 12 = 101 と 992 + 102 = 9901 が得られます。

　①②ができれば、あとは単純です。ブラーマグプタの二平方恒等式を使います。

A2 + B2 = 101 × 9901

　　　 = (102 + 12) × (992 + 102)

　　　 = (10・99 - 1・10)2 + (10・10 + 1・99)2 = 9802 + 1992 ･･･ (1)

　　　 = (10・99 + 1・10)2 + (10・10 - 1・99)2 = 10002 + 12 ･･･ (2)

　あとは (A, B) から (a, b) への変換を行い、題意を満たすかをチェックすればOKです。結果、(a, b) = (990, 100) が正解と分かります。

 

n = 10 の場合

　n = 10 の場合も同様です。①素因数分解の結果は、1010 + 1 = 101 × 3541 × 27961 となります。

　②素因数のそれぞれを２平方和で表した結果は、102 + 12 = 101 と 542 + 252 = 3541 と 1442 + 852 = 27961 となります。

　ここから先の (A, B) の組を見つけるところがちょっとややこしいです。ブラーマグプタの二平方恒等式では、２つの２平方和の積が、２平方和として２通りに表せます。これが３つになると、計４通りの組み合わせが存在することになります。

　結果、(A, B) = (1, 100000), (64700 76249), (48320, 87551), (19801, 98020) の４つが解となります（A, B とも正の場合のみ記載）。n = 6 の場合と同様、(a, b) への変換を行い、題意を満たすかをチェックすればOKです。ただし A または B が負になる場合ももれなくチェックしましょう。

　(a, b) = (17650, 38125), (25840, 43776), (74160, 43776), (82350, 38125), (99010, 9901) の５つが正解です。

　以上の方針だと、n = 16 程度までなら、一瞬の実行時間で答えを得ることができます。計算量は O(√(10n+1の最大の素因数)) です。

 

さらに大きな n の場合

　n がさらに大きな値の場合の高速化を考えます。

　①素因数分解については、上では単純な試し割りのアルゴリズムを紹介しましたが、ポラード・ロー素因数分解法などの高速なアルゴリズムが有効です。

　②素数の２平方和への分解については、セレの方法などが存在します。H. Davenport 著「The Higher Arithmetic: An Introduction to the Theory of Numbers」の p.120~122 に詳しく書かれてあります。上記 Amazon のページでキーワード「Serret」でなか見検索すると確認できると思います。同ページのルジャンドルの方法とともに、参照コードをこちらに書きました。

　この方針で、n = 36 以下の全ての n に対する答えを ideone 環境で約 1.2 秒で得ることができました（Python）。コードはこちらです。

　10n + 1 が大きな素因数を持つ場合はこの方法でも高速に解けません。n = 38 の場合（1038 + 1 = 101 × 722817036322379041 × 1369778187490592461）は、上記コードでも数分待っても答えは得られませんでした。