約数和配列を用いて友愛数を求める
友愛数とは、約数の和がお互いの数自身になる数の組を言います。例えば220と284は友愛数です。220は1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110を約数として持ちその和は284となり、284は約数として1, 2, 4, 71, 142となりその和は220となり友愛数であることが分かります。
プログラムは完全数 - 約数和配列で使った約数和配列を利用しています。ただし、完全数はその数自身で判定が可能だったので配列を繰り返し使用することで大きな値まで探索が可能でしたが、友愛数は二つの数が関係することからこの方法を使用することが出来ません。大きな範囲を探索するには更なるプログラムの工夫が必要です。
#include <stdio.h>
#include <stdlib.h>
#define DIVISORSUMLIST_SIZE_MAX (100000000)
long long divisorsumlist[DIVISORSUMLIST_SIZE_MAX];
int main( int ac, char *av[] )
{
long long divisorsumlistsize, i, j, base, n, n2, tmp, min_n, max_n;
/* コマンドラインから探索範囲を決定する */
if( ac < 3 ) {
fprintf( stderr, "usage : yuuai min_n max_n\n" );
return( 1 );
}
min_n = strtoll( av[1], NULL, 10 );
max_n = strtoll( av[2], NULL, 10 );
if( min_n < 1 || max_n < 1 ) {
fprintf( stderr, "bad parameter\n" );
return( 1 );
}
/* DIVISORSUMLIST_SIZE_MAX / 2 ごとの整数区間を調べる */
/* 約数和配列自体は最大 DIVISORSUMLIST_SIZE_MAX まで求める(ペア数字が配列内に収まりやすくするため) */
for( base = min_n; base <= max_n; base += DIVISORSUMLIST_SIZE_MAX / 2 ) {
/* 整数区間配列の大きさを決める */
if( max_n - base + 1 > DIVISORSUMLIST_SIZE_MAX / 2 )
divisorsumlistsize = DIVISORSUMLIST_SIZE_MAX / 2;
else
divisorsumlistsize = max_n - base + 1;
/* 配列を初期化する */
for( i = 0; i < divisorsumlistsize * 2; i++ )
divisorsumlist[i] = 0;
/* 約数の和を求める */
for( i = 1; i <= ( base + divisorsumlistsize * 2 ) / 2; i++ ) {
j = ( base + i - 1 ) / i * i;
if( j < i * 2 )
j = i * 2;
for( ; j < base + divisorsumlistsize * 2; j += i )
divisorsumlist[j - base] += i;
}
/* 約数の和が互いの数自身になれば友愛数である */
for( i = 0; i < divisorsumlistsize; i++ ) {
n = base + i;
n2 = divisorsumlist[i];
if( n2 <= n )
continue;
if( n2 < base + divisorsumlistsize * 2 ) {
if( divisorsumlist[n2 - base] == n )
printf( "%lld %lld\n", n, n2 );
}
else {
/* 配列の範囲外の場合は個別計算する */
tmp = 0;
for( j = 1; j <= n2 / j; j++ ) {
if( n2 % j == 0 ) {
tmp += j;
if( n2 / j != j && j != 1 )
tmp += ( n2 / j );
}
}
if( n == tmp )
printf( "%lld %lld\n", n, n2 );
}
}
}
return( 0 );
}
両方ともに1千万以下の友愛数は全部で100組見つかりました。これまでに見つかった友愛数は全て偶数同士もしくは奇数同士らしいのですが、実際1千万以下の範囲で見つかった100組も全て偶数もしくは奇数同士でした。偶数と奇数の組み合わせの友愛数があるかどうかはまだ分かっていないそうです。
(2012/2/16追記)
メモリが潤沢に使用できるようになりましたので、配列を大きくして1億以下まで求めるようにしました。1億以下同士の友愛数は231組ありました。
(2020/4/25追記)
ソースコードを少し整理しました。
(2021/10/5追記)
配列を更に大きくして10億以下の範囲で求めました。10億以下同士の友愛数は564組あり、最大の組は(957871508 998051692)でした。
10億以下ではメモリを4GB(4bytes x 10億)使用しましたが、100億以下となると8バイト整数を使用するため80GBものメモリが必要になってしまい、そうした環境を調達するのは難しそうです。将来、メモリが更に安価になったら挑戦してみたいです。
ちなみに素数を求める - エラトステネスのふるいを繰り返し使うのように約数和配列を繰り返し使う方法は、友愛数の組の両方が同じ約数和配列に収まってる場合のみしか求められないため不適です。大抵の友愛数は近しい範囲の値に収まっているのですが、例外があり得ないとは言い切れません。
(2021/11/10追記)
クラウドサービスを使えば巨大なメモリを搭載した環境も手軽に調達できることに気づいたので、AWS にて m5.8xlarge(メモリ128GB)のインスタンスを作成して100億以下の範囲で友愛数を探索しました。1時間あたり1.984 USDとなかなか痺れる単価ですが、せいぜい数時間の利用ですから実際にはそれほど厳しいわけではありません。
100億以下同士の友愛数は1391組あり、最大の組は(9958985128 9994662872)でした。
(2026/2/5追記)
当たり前だと思って記載してなかったら自分でも混乱してしまったので明記しておきます。プログラムは 32bit 整数(int)を用いていますので、100億以下の範囲ではオーバーフローしてしまいます。一時的に int を 64bit 整数(long long)に書き換えて計算を実行しています。
(2026/3/13追記)
int を long long に変更し、ソースコードを整理するとともに、約数和配列を繰り返し使用するように修正しました。
約数和配列を繰り返し使用する場合、ペアの数字が約数和配列に無い場合に対応出来ないという問題があります。これは素直に約数和をその場で計算することで対処しました。次節「約数和配列を用いずに友愛数を求める」のコードを合わせたようなものですね。計算時間は不利になりますが、巨大なメモリがなくても広範囲に探索が行えるようになります。
なお、約数和配列の範囲を計算範囲そのままにするとペアの数字が約数和配列に収まらないケースが多発するため、約数和配列は計算範囲の2倍用意するようにしました。実装としては計算範囲の方を半分に区切るようにしています。
(2026/3/27追記)
単純な方法で約数和を求める処理を√nまでしか計算しないようにして高速化しました。
(2026/3/27追記その2)
約数和配列を繰り返し使用することでメモリが潤沢に無くても広範囲に友愛数を探索できるようになりました。ただし、約数和配列の範囲外になった場合に単純な方法で約数和を計算している処理が非常に重いため、実際のところとしては探索範囲をそこまで広げることが出来ていません。今回は10億以下までを探索しましたので、一覧を更新しました。合わせて組の数も掲載しています。100億以下の掲載は一旦取り下げています。
約数和配列の範囲外になった場合に単純な方法を使用するメリットの一つとして、ペアが両方とも約数和配列に収まっていなくても探索可能なことです。例えば1億以下で探索した場合、従来は見つからなかった(95791430 115187002)という組を発見することが出来るようになりました。
約数和配列を用いずに友愛数を求める
(2020/9/9追記)
順序が逆になってしまいましたが、約数和配列を用いず、常に約数を求めながら友愛数を求める版を作りました。速度的には非常に遅くなってしまいますが、探索範囲を配列サイズに制限されないというメリットがあります。
(2021/10/4追記)
約数を調べる際、n/2までしか調べなくていいことをうっかりしてたので修正しました。
(2026/3/27追記)
int を long long に変更したのがこちらには反映されていませんでした。作業漏れです。
合わせて、単純な方法で約数和を求める処理を√nまでしか計算しないようにして高速化しました。
#include <stdio.h>
#include <stdlib.h>
int main( int ac, char *av[] )
{
long long i, n, n2, tmp, min_n, max_n;
/* コマンドラインから探索範囲を決定する */
if( ac < 3 ) {
fprintf( stderr, "usage : yuuai min_n max_n\n" );
return( 1 );
}
min_n = strtoll( av[1], NULL, 10 );
max_n = strtoll( av[2], NULL, 10 );
if( min_n < 1 || max_n < 1 ) {
fprintf( stderr, "bad parameter\n" );
return( 1 );
}
/* 探索範囲の数を調べる */
for( n = min_n; n <= max_n; n++ ) {
n2 = 0;
for( i = 1; i <= n / i; i++ ) {
if( n % i == 0 ) {
n2 += i;
if( n / i != i && i != 1 )
n2 += ( n / i );
}
}
if( n2 <= n )
continue;
tmp = 0;
for( i = 1; i <= n2 / 2; i++ ) {
if( n2 % i == 0 )
tmp += i;
}
if( n == tmp )
printf( "%lld %lld\n", n, n2 );
}
return( 0 );
}
友愛数の一覧
小さい数字が10万以下の友愛数は以下の通りです。
(220 284) (1184 1210) (2620 2924) (5020 5564) (6232 6368) (10744 10856) (12285 14595) (17296 18416) (63020 76084) (66928 66992) (67095 71145) (69615 87633) (79750 88730)
10万以上の友愛数は数が多いので別ファイルとして置いておきます。
友愛数の組の数一覧
| 範囲 | 友愛数の組の数 |
|---|---|
| 小さい数字が1000以下 | 1個 |
| 小さい数字が1万以下 | 5個 |
| 小さい数字が10万以下 | 13個 |
| 小さい数字が100万以下 | 42個 |
| 小さい数字が1000万以下 | 108個 |
| 小さい数字が1億以下 | 236個 |
| 小さい数字が10億以下 | 586個 |