久しぶりに、プログラミングの記事です。
問題
ab(a＋b)(a－b)＝13c2
を満たす自然数(a,b,c)の組を見つけてください。
というものでした。
最初に数学的考証もせずに思うがままに組んだものがこちら。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
unsigned long long a,b,c,c2,m,n,count;
m = 13;
for (a=2,count=0; count<100; a++)
for (b=1; b<a; b++) {
n = a*b*(a+b)*(a-b);
if ( a*(a+b)*(a-b) > n || a*b*(a+b) > n ) {
fprintf(stderr,"*** overflow *** a=%llu,b=%llu\n",a,b);
a++;
b=0;
continue;
}
if ( n%m == 0 ) {
c2 = n/m;
c = sqrtl(c2)+0.5;
if ( c*c == c2 ) {
printf("a=%llu,b=%llu,c=%llu\n",a,b,c);
count++;
}
}
}
return EXIT_SUCCESS;
}
解は100個も出せば十分だろう。
エラー判定も一応しておくか。
と、解が見つかればいいなというだけのものでした。
続いて、13のところが、14の場合や、15の場合も見つけて欲しいとのことで、ああ、そこ変わるんだなと、少しプログラムを改変します。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char *argv[])
{
unsigned long long a,b,c,c2,m,n,count;
m = 13;
if ( argc == 2 ) {
m = atoi(argv[1]);
}
for (a=2,count=0; count<100; a++)
for (b=1; b<a; b++) {
n = a*b*(a+b)*(a-b);
if ( a*(a+b)*(a-b) > n || a*b*(a+b) > n ) {
fprintf(stderr,"*** overflow *** a=%llu,b=%llu\n",a,b);
a++;
b=0;
continue;
}
if ( n%m == 0 ) {
c2 = n/m;
c = sqrtl(c2)+0.5;
if ( c*c == c2 ) {
printf("a=%llu,b=%llu,c=%llu\n",a,b,c);
count++;
}
}
}
return EXIT_SUCCESS;
}
パラメータに値を渡すように変更します。
どうやら、合同数について、自然数(a,b,c)の組を求めたいんだなということに気が付きます。
100以下の合同数は、
5, 6, 7, 13, 14, 15, 21, 22, 23, 29, 30, 31, 34, 37, 38, 39, 41, 46, 47, 53, 55, 61, 62, 65, 69, 70, 71, 77, 78, 79, 85, 86, 87, 93, 94, 95
と36個あるので、これらについて調べてみることにします。
その中でも、合同数が素数の場合、
5, 7, 13, 23, 29, 31, 37, 41, 47, 53, 61, 71, 79
の後半は、64ビットの範疇で解が見つけられるか、怪しいところでもあります。
このままのプログラムでは、ちょっとだけ問題があるなと思い、プログラムを改変しますが、それについては後で説明します。
例えば、23の場合
>abc13 23 2>nul
a=24336,b=17689,c=72306780
a=42025,b=6647,c=144613560
a=48672,b=35378,c=289227120
a=73008,b=53067,c=650761020
a=84050,b=13294,c=578454240
a=97344,b=70756,c=1156908480
a=121680,b=88445,c=1807669500
a=126075,b=19941,c=1301522040
a=146016,b=106134,c=2603044080
a=168100,b=26588,c=2313816960
a=170352,b=123823,c=3543032220
a=210125,b=33235,c=3615339000
...
まず、プログラム的にcについて小さい順で求めているというわけではないこと。
また、cの2乗が64ビットの範疇に収まっていたとしても、取りこぼしが出てしまうという点です。
上記の中で一番大きいc、
c＝3615339000＜4294967296＝232
と32ビットで収まっていますが、2乗して23倍してしまうと、64ビットに収まらないことが解るかと思われます。
つまり、今までのプログラムではこういった解を取りこぼしていたということになります。
多倍長演算を組み込んだわけではなく、もっとゴリッゴリと場合分けをしました。
さて、100以下の合同数mについて自然数(a,b,c)を1個ずつ求めてみる。
m=5,a=5,b=4,c=6
m=6,a=2,b=1,c=1
m=7,a=16,b=9,c=60
m=13,a=325,b=36,c=9690
m=14,a=8,b=1,c=6
m=15,a=5,b=3,c=4
m=21,a=4,b=3,c=2
m=22,a=50,b=49,c=105
m=23,a=24336,b=17689,c=72306780
m=29,a=4901,b=4900,c=90090
m=30,a=5,b=1,c=2
m=31,a=1600,b=81,c=103320
m=34,a=9,b=8,c=6
m=37,a=777925,b=1764,c=4737551070
m=38,a=1250,b=289,c=118575
m=39,a=13,b=12,c=10
m=41,a=25,b=16,c=60
m=46,a=72,b=49,c=462
m=47,a=14561856,b=2289169,c=12111037689240
m=53,a=1873180325,b=1158313156,c=297855654284978790
m=55,a=125,b=44,c=1170
m=61,a=12079525,b=10227204,c=9147755349330
m=62,a=39200,b=22801,c=121068780
m=65,a=13,b=5,c=12
m=69,a=192,b=169,c=1976
m=70,a=9,b=5,c=6
m=71,a=3600,b=121,c=281820
m=77,a=2816,b=2809,c=63600
m=78,a=26,b=1,c=15
m=79,a=169000000,b=166952241,c=495683115837000
m=85,a=85,b=36,c=462
m=86,a=338,b=49,c=4641
m=87,a=17956,b=169,c=3353350
m=93,a=1444,b=75,c=49210
m=94,a=14112,b=529,c=3974124
m=95,a=1445,b=76,c=49062
ただし、このデータはa,b,cのいずれか一つでも最小の解を列挙していとは限らないことに注意してください。
もし手計算でやるならば、自明な解があるものであれば、どうにか見つけることは出来るのではなかろうか。
自明な解とは、a＝m, b＝m, a＋b＝m, a－b＝mといった条件で解が見つかるような場合でしょうか。
mが素数だと、見つけるのが大変かなとは思う。
さて、私のプログラマ遍歴でも晒してみると、
BASIC（自宅ではN88、学校ではF）、FORTRAN（77、90）、C（当時はTurbo-Cのアカデミアパックを使っていました）、Lisp、CASL、6502を大学時代にやって、CASLによる画像処理といった論文を書いたけど、特にどこかで見られるようにはなってないです。
社会人になって、8086、Java、PHP、Javascriptをやって、
実際に業務においてのプログラマと言えるような期間は1年もなくて、DOSのテストプログラマだったころに、8086やCでテストプログラムを組んだくらいで、あくまでもテスト用のプログラムなので、世に出ることはありません。
さて、私はプログラミングが出来るのか？と問われると、出来ない人からみれば、出来るんだろうけれども、今現在プログラマと呼ばれる人たちからは、どう見えるのだろうか。
おそらく、古いプログラマという認識になるのではなかろうか。
ではでは