友愛数を列挙する

雑記帳

Mathematicaで友愛数を列挙するプログラム例として以下のようなものが見受けられる。

yakuwa[n_] := DivisorSigma[1, n] - n;
Do[If[(yakuwa[yakuwa[k]] == k) && (yakuwa[k] != k), Print[{k, yakuwa[k]}]], {k, 1, 1000}];

しかし、Doでループを回してPrintで書き出していくのはMathematica的に美しくないと思う。

Mathematicaなら関数型プログラミングとパターンマッチを用いるのが良いと思うので、私なら以下のように書く。

Cases[NestList[DivisorSigma[1, #] - # &, #, 2] & /@ Range[100000], {a_, b_, a_} /; a < b -> {a, b}]

実行速度もこちらのほうが1割程度速い。

この書き方から一般の社交数を求めるプログラムに拡張するのは簡単で、例えば5個組の社交数を探す場合は以下のようになる。

Cases[NestList[DivisorSigma[1, #] - # &, #, 5] & /@ Range[100000], {a_, b__, a_} /; a < Min[b] -> {a, b}]

ただし、私の方法はメモリを大量に使う。