約数関数の整理

Eulerは完全数や友愛数を考えるにあたって、すべての約数の和を表す約数関数を定義し、駆使しました。(友愛数とは何かについては別の機会に紹介します。)
Maximaではは約数関数はdivsumという名前です。

ではグラフを書いてみましょう。

upper : 500$
ns : makelist(n, n, 1, upper)$
plot2d([[discrete, ns, map(divsum, ns)], 2*x], [x, 0, upper], [style, points, line])$

赤線上の点が完全数です。

約数関数は次のような性質を持ちます。

1. pが素数であることと、が成り立つことは同値である。
2. Nが完全数であることと、が成り立つことは同値である。
3. pが素数のとき、が成り立つ。
4. もしaとbが互いに素なら、が成り立つ。

Maximaのdivsum関数は、引数が整数の時には値を計算してくれますが、引数が式の場合にはそのまま(名詞形)で、上記のような性質を利用した整理をしてくれません。
では、整理する関数を作ってみましょう。

matchdeclare(xx, all)$
defrule(dsrule0, divsum(xx, 1), dssimp(divsum(xx, 1)))$

matchdeclare(pp, lambda([x], integerp(x) and primep(x)))$
matchdeclare(nn, lambda([x], (integerp(x) or featurep(x, integer)) and x > 0))$
defrule(dsrule1, divsum(pp^nn, 1), (pp^(nn + 1) - 1)/(pp - 1))$

factorall(expr) := block([e],
  e : factor(expr),
  if atom(expr) then e
  else if inpart(e, 0) = "*" then map(factorall, e)
  else if inpart(e, 0) = "^" then factorall(inpart(e, 1))^inpart(e, 2)
  else e)$

base(factored) :=
  if not atom(factored) and inpart(factored, 0) = "^"
  then inpart(factored, 1)
  else factored$

dssimp(expr) := block([arg],
  /* expr should be a form of "divsum(argument, 1)" by dsrule0. */
  arg : factorall(inpart(expr, 1)),
  if inpart(arg, 0) = "*" then block([bases],
    bases : maplist(base, arg),
    if apply("and", maplist(integerp, bases))
    then if (apply(ezgcd, bases))[1] = 1
         then apply1(map(divsum, arg), dsrule1)
         else divsum(arg)
    else if readonly("Are ", bases, " relatively prime?[yes, no]: ") = 'yes
    then apply1(map(divsum, arg), dsrule1)
    else divsum(arg)
    )
  else apply1(divsum(arg), dsrule1)
  )$

divsumsimp(expr) := applyb1(expr, dsrule0)$

使用した機能

• integerp
• featurep
• atom
• inpart
• ezgcd
• readonly
• applyb1

関数divsumsimpは、先に書いた約数関数の性質を用いて、引数式の中の関数divsumの整理を試みます。
こんな感じになります。

declare(a, integer, m, integer, n, integer)$
assume(m > 0, n > 0)$
divsumsimp(divsum(a*6^m*10^n));

Are [a,3,5,2] relatively prime?[yes, no]: yes

道具が一つ作れました。