search
LoginSignup
1

posted at

約数の総和の逆関数を作った

何言ってるのかよくわからないので例を挙げます。

たとえば、自然数 $60$ の約数を列挙すると、 $1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60$ です。この総和は $168$ ですね。
この関係を表す関数に約数関数 https://ja.wikipedia.org/wiki/%E7%B4%84%E6%95%B0%E9%96%A2%E6%95%B0 というものがあって、 $σ_1(60)=168$ と書いたり、あるいは単に $σ(60)=168$ と書いたりします。
この逆関数 $σ^{-1}$ は複数の値をとります。 $σ^{-1}(168)=60, 78, 92, 123, 143, 167$ など。

この値を全て列挙してみようという試みです。

python
from sympy import sieve
import copy


def sum_primes(n: int, primes: list):
    mydict = {}
    for p in primes:
        p_now = 1
        p_sum = 1
        mylist = []
        while True:
            p_now *= p
            p_sum += p_now
            if p_sum > n:
                break
            mylist.append(p_sum)
        mydict[p] = mylist
    return mydict


def sum_divisor_inverse(n: int):
    primes = list(sieve.primerange(n))
    spdict = sum_primes(n, primes)
    ans, _ = sdi(n, spdict)
    for s in ans:
        t = ""
        p = 1
        for ss in s:
            t = " * " + t
            if t == " * ":
                t = ""
            t = str(ss[0])+"^"+str(ss[1]) + t
            p *= ss[0]**ss[1]
        print(t, "=", p)
    return ans


def sdi(n: int, spdict: dict):
    if n == 1:
        return [[]], True
    f = False
    mylist = []
    divs = divisors(n)
    spdict_loop = copy.deepcopy(spdict)
    for k, v in spdict_loop.items():
        for i in range(len(v)):
            if v[i] in divs:
                spdict_cp = copy.deepcopy(spdict)
                del spdict_cp[k]
                s, flag = sdi(n // v[i], spdict_cp)
                if flag:
                    f = True
                    for ss in s:
                        ss.append([k, i+1, v[i]])
                        mylist.append(ss)
        del spdict[k]
    return mylist, f


def divisors(n: int):
    # https://qiita.com/LorseKudos/items/9eb560494862c8b4eb56 を参考にした。O(√N)
    lower_divisors, upper_divisors = [], []
    i = 1
    while i * i <= n:
        if n % i == 0:
            lower_divisors.append(i)
            if i != n // i:
                upper_divisors.append(n // i)
        i += 1
    return lower_divisors + upper_divisors[::-1]

$n=100000$ でも数秒かかるので、あんまり $n$ が大きいと使い物にならなそう。

sum_divisor_inverse(1680)を実行すると、以下のように print されます。たとえば、 $540 = 2^2・3^3・5^1$ の正の約数の和は $1680$ です。

output
2^1 * 3^1 * 139^1 = 834
2^1 * 3^3 * 13^1 = 702
2^2 * 3^1 * 59^1 = 708
2^2 * 3^3 * 5^1 = 540
2^2 * 7^1 * 29^1 = 812
2^2 * 11^1 * 19^1 = 836
2^2 * 239^1 = 956
2^3 * 7^1 * 13^1 = 728
3^1 * 13^1 * 29^1 = 1131
3^1 * 419^1 = 1257
3^3 * 41^1 = 1107
5^1 * 13^1 * 19^1 = 1235
11^1 * 139^1 = 1529
19^1 * 83^1 = 1577

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
What you can do with signing up
1