４．高速化のための別のアルゴリズム

前記の

013, 019, 036, 056, 059, 079, 236, 247, 678, 789

について解を持つかどうか調べてみよう。
1., 3. で説明したアルゴリズムは結局 10ⁿ のオーダーだから、
探査範囲を１桁拡げるのに10倍の時間がかかる。
これではすぐに行き詰まってしまう。
何か他の方法はないか？

ここでは、

Ｉ．ヴァルディ、『Mathematica計算の愉しみ』、（株）トッパン、1991
（原題：『Computational Recreations in Mathematica』, Ilan Vardi)

で紹介されているアルゴリズムについて説明する。

2n桁、2n−１桁の数の平方根は、上位ｎ桁によって支配される
下ｎ桁、ｎ−１桁は、四捨五入によって１の位に影響を与える程度である
そこで以下のようなアルゴリズムを考える。

◎  ３種類の数からなるｎ桁の数をａとする
    ｃ＝［ａに０をｎ−１個くっつけた数］とする（ｃ＝ａ×10^(n−1)）
    ｄ＝isqrt(c)
    (d＋1)^2 が３種類の数かどうか調べる
      このとき、下位ｎ桁のみ調べればよい
      ∵  上位ｎ桁は３種類の数となることがわかっているので
          四捨五入の影響も考えて１桁追加した
          下位ｎ桁のみを調べればよい
    ３種類なら、見つかったことになる

    ｃ＝［ａに０をｎ個くっつけた数］とする（ｃ＝ａ×10^n）
    ｄ＝isqrt(c)
    (d＋1)^2 が３種類の数かどうか調べる
      このとき、下位ｎ＋１桁のみ調べればよい
    ３種類なら、見つかったことになる

このアルゴリズムだと、ａに１桁追加すると、探査範囲が２桁増える。
しかも、ａを１桁増やすことによるａの範囲は、元の範囲の３倍である。
（特に３種類のうち１つが０の場合は、２倍となる）。
前のアルゴリズムでの増加が10倍だったことを考えると、
これはかなりの進歩である。

さて、◎をいかにして生成するか？

以下のような数を順次生成していきたい。

• 0,1,2 のとき、
  1 , 2 , 10 , 11 , 12 , 20 , 21 , 22 , 100 , 101 , ……
• 1,2,3 のとき、
  1 , 2 , 3 , 11 , 12 , 13 , 21 , 22 , 23 , 31 , 32 , 33 , 111 , ……

アルゴリズムは以下のとおりである。

まず対象となる３文字を T$(0)、T$(1)、T$(2)に小さい順に格納しておく。
３種類の数字からなる数を S$ とする。桁数をｉとする。

○  ｉ桁目に着目する

    ｉ桁目がT$(0)なら
        ｉ桁目をT$(1)に置き換えて
        ｉ−１桁から１桁まではそのまま

    ｉ桁目がT$(1)なら
        ｉ桁目をT$(2)に置き換えて
        ｉ−１桁から１桁まではそのまま

    ｉ桁目がT$(2)なら
        ｉ桁目をT$(0)に置き換えて
        ｉ＝ｉ−１として、○のサブルーチンへ

        ｉが０となった場合は
            最上位桁に（最左端に）T$(0)を付加する
            T$(0) が０の時は、T$(1)を付加する

プログラムは以下のとおり。

 10   ' nnuma.ub non zero
 20   ' print=print+"nnum2.txt"
 30   dim T$(3)
 40   T$(0)=cutlspc(str(4))
 50   T$(1)=cutlspc(str(5))
 60   T$(2)=cutlspc(str(6))
 70   S$=T$(0)
 80   gosub 90:gosub 180:goto 80
 90   U$=S$:S$="":I=len(U$)
100   if I=0 then S$=T$(0)+S$:return
110   X$=mid(U$,I,1)
120   if X$=T$(0) then S$=T$(1)+S$:goto 150
130   if X$=T$(1) then S$=T$(2)+S$:goto 150
140   if X$=T$(2) then S$=T$(0)+S$:dec I:gosub 100:return
150   for J=I-1 to 1 step -1:S$=mid(U$,J,1)+S$:next J:return
160   end
170   '
180   A=val(S$):I=len(S$)
190   Wa=A*10^(I-1):B=(isqrt(Wa)+1)^2:gosub 230
200   Wa=Wa*10:B=(isqrt(Wa)+1)^2:gosub 230
210   return
220   '
230   if B@10=0 then return
240   C$=cutlspc(str(B)):D=len(C$)
250   for J=D to D-I step -1:E$=mid(C$,J,1)
260     for K=0 to 2
270       if T$(K)=E$ then cancel for:goto 300
280     next K
290     cancel for:return
300   next J
310   print isqrt(B),B
320   return

160行目のENDは、いらない。
バグがあった時の暴走止めでデバッグ時のなごりである。まあ、あっても悪くはない。
そんなことより、140行目の「gosub 100」というのがあまり美しくない。
90行目から150行目までが、上記３種類の数字を生成するサブルーチンなのだが、
140行目の「gosub 100」では、その途中の100行目に飛び込んでいる。
サブルーチンのラベルを目指して飛び込んでいるのではなく、100という行番号を目指して飛び込んでいる。
つまりgotoと同じである。

プログラム自体は正しく動作するので、とりあえずこれでよしとする。

実行結果は以下のとおり。
b=a², a ≤ 10¹² の範囲で、次のような解が見つかった。

数字	a	b (=a2)
036	1740871344	3030633036360366336
059	70778174997	5009550055905955950009
	771395165003	595050500590005595990009
236	14945311394	223362332663626223236
	179619783606	32263266662666266363236
247	88002411582	7744424444247727742724

013 , 019 , 056 , 079 , 678 , 789 に解はあるのか？

('00/09/23 追記) 10¹⁵まで探査が完了。
('03/07/13 追記) 10²⁰まで探査が完了。
('04/05/16 追記) 12 May 2004, Milos Tatarevic からメールがあり、「個別解」の０を含まないパターンについて、
10²⁵まで探査が完了したとのこと。結果はこちらに記す。