３辺ａ、ｂ、ｃが有理数となる直角３角形の面積ｇとして表される数を合同数(congruum)と呼ぶ。
ａ、ｂ、ｃが自然数の場合、ａ、ｂが直角を挟む辺とすると、直角３角形の面積ｓは、
ｓ＝ａｂ／２
ｋを、ｓ＝k2ｇ、（ｋは１以上の自然数、ｇは２乗因子を含まない）と定義し、両辺をk2で割ると、
ｇ＝(a/k)(b/k)／２
よって、合同数の定義は、
「３辺ａ、ｂ、ｃが自然数となる直角３角形の面積ｓから２乗因子を取り除いた数として表される自然数ｇ」
と云いかえることができる。
３辺が互いに素な自然数となる直角３角形の３辺ａ、ｂ、ｃは以下の式で表される。
ａ＝m2−n2
ｂ＝2mn
ｃ＝m2＋n2
（ｍ、ｎは任意の自然数）
この直角３角形の面積 ｓ は、ｓ＝mn(m2−n2)となる。
ｓの２乗因子をｋをすると、ｓ＝k 2ｇ＝mn(m2−n2) 。
例えば、
ｍ＝２、ｎ＝１
に対し、
ａ＝m2−n2＝22−12＝３
ｂ＝2mn＝2・2・1＝４
となり、
ｓ＝mn(m2−n2)＝2・1（22−12）＝６
であるから、６は合同数である。
ここで、関数 T(n) を「任意の自然数ｎに対し、ｎの２乗因子を取り除いたもの」とすると、
合同数の定義式は、
ｇ＝T(m)T(n)T(m+n)T(m-n)
となる。
ある値以下のｍ、ｎについてｇ＜1000となるケースを探してみる。メインループは、
for m=1 to 10^6
for n=1 to m-1
g=t(m)*t(n)*t(m+n)*t(m-n)
if g<1000 then print g,m,n
next n
next m
となる。
２乗因子を取り除くためのアルゴリズムは以下のとおり。
10 fnSub(N) 20 local P,W=1 30 P=prmdiv(N):if P=0 then return(N) 40 N=N\P:if N@P=0 then N=N\P else W=W*P 50 if N=1 then return(W) else goto 30
まず、最小素因数で１回割る。
２回目も割れるときは２乗因子、割れない場合は２乗因子ではないのでｗに値をセットする。
ｎが１になったときのｗの値が、ｎから２乗因子を取り除いた値である。
これを前のメインループに入れて実行してもよいが、もう少しメインループの方の効率化を考える。
これは明らかであろう。
ｍ、ｎともに偶数のときは、上のとおり。
ｍ、ｎともに奇数のとき、
ｍ＋ｎ
ｍ−ｎ
は共に偶数であり、共通因子２を持つ。
よって、Ｍ＝(ｍ＋ｎ)/2、Ｎ＝(ｍ−ｎ)/2 とおくと、Ｍ＜ｍ、Ｎ＜ｎ。
従って、ｍ、ｎを調べる前に、既にＭ、Ｎが調べられている。
上記２つの条件を満たすｍ、ｎについて、
ｍ、ｎ、ｍ＋ｎ、ｍ−ｎ
の内、任意の２つは、互いに素となる（これも容易に確かめられる）。
この時、ｍ、ｎ、ｍ＋ｎ、ｍ−ｎは、
ｍ、ｎのどちらかは偶数
残りの３個は奇数
となる。このときの解を原始的な解と呼ぶ。
また、原始的な解 (m, n) から (m+n, m-n) を作ると、上記により、これも解になる。
（ただし、奇数同士の組となる）。これを同値な解と呼ぶことにする。
よって、配列 T(i)、ｉ＝1,2,…… にｉから２乗因子を取り除いたものを格納しておき、
ｍ、ｎの代わりに、T(m)、T(n)、T(m+n)、T(m-n) を使えばよい。
ｇ＝T(m)T(n)T(m+n)T(m-n)
は２乗因子を持たない。
また、ｇの４つの因子 T(m)、T(n)、T(m+n)、T(m-n) について、１つでも1000を越えると、ｇ＞1000となる。
よって、1000を越えた時点で次のｍ、ｎに移ってよい。
以上、1.〜4.で処理速度は格段に向上する。原始的な解を求めるためのプログラムは以下のとおり。
10 ' cong_1.ub 20 K%=1000 30 L%=10000:dim T%(2*L%+1) 40 for I%=1 to 2*L%+1:T%(I%)=fnSub(I%):next I% 50 ' 60 for M%=2 to L%:if T%(M%)>=K% then 160 70 for N%=1 to M%-1:if T%(N%)>=K% then 150 80 if even(M%+N%) then 150 90 if gcd(M%,N%)>1 then 150 100 if T%(M%+N%)>=K% then 150 110 if T%(M%-N%)>=K% then 150 120 G=T%(M%)*T%(N%)*T%(M%+N%)*T%(M%-N%) 130 if G>=K% then 150 140 print G,M%,N% 150 next N% 160 next M% 170 end 180 ' 190 fnSub(N) 200 local P,W=1 210 P=prmdiv(N):if P=0 then return(N) 220 N=N\P:if N@P=0 then N=N\P else W=W*P 230 if N=1 then return(W) else goto 210
M%、N% を integer を超す範囲までループで廻したい場合は、まず、M%、N% を 多倍長型に変更する。
それから、UBASIC では配列として割り当てられる個数に制限があるので、
ある値以上については、T%( ) に格納するのではなく、直接 fnSub(N) を呼び出す形とする。
また、この m、n 及び、k2g＝mn(m2−n2) となる k に対して、
X=mg/n
Y=kg2/n2
とおくと、点 (X, Y) は、楕円曲線 Y2=X3-g2X 上の有理点となる。
かつ、楕円曲線上の有理点 (X, Y), X, Y∈Ｑ は、(x/z2, y/z3),
x, y, z∈Ｚ と表すことができる。
解となる M%、N% が見つかった時、この x, y, z も合わせて計算することにする。
『数論における未解決問題集』によると、1000以下の合同数は全部で361個ある。
100万以下のｍ、ｎについて調べたところ、そのうち、198個について解が求められた。
実行結果は以下のとおり。
解が見つからなかったのは、以下の163個。
47, 53, 61, 79, 101, 103, 118, 127, 134, 142, 157, 166, 167, 173, 181, 183, 191, 197, 199, 206, 213, 223, 229, 237, 247, 262, 263, 269, 271, 277, 278, 287, 293, 302, 303, 309, 311, 317, 326, 327, 334, 335, 349, 358, 359, 365, 367, 373, 382, 383, 389, 397, 398, 407, 413, 415, 421, 431, 437, 439, 446, 447, 453, 454, 461, 463, 469, 471, 478, 479, 485, 487, 493, 501, 502, 503, 511, 519, 526, 533, 541, 542, 543, 557, 559, 566, 573, 583, 589, 591, 597, 599, 607, 613, 614, 623, 631, 638, 647, 653, 655, 661, 662, 677, 678, 679, 685, 687, 695, 701, 703, 717, 718, 719, 727, 733, 742, 743, 757, 758, 766, 767, 773, 781, 789, 797, 807, 815, 821, 822, 823, 829, 831, 838, 839, 853, 862, 863, 871, 877, 878, 886, 887, 893, 894, 902, 911, 917, 919, 926, 933, 941, 942, 958, 965, 967, 974, 982, 983, 989, 991, 997, 998
| 前 | この章の目次 | 次 |
|---|