任意の整数に於ける約数ペアの合計が最小になる組み合わせ

なんだかコメ欄で凄い揉めていて,なんで揉めているのかさっぱり分からないけど.

Getting multiples of a number with least sum in O(1) – StackOverflow

最初,タイトルから誤解していて,単調増大な数列とその逆な数列のペアの合計値が最小になる組み合わせかと思って,楽勝じゃないかと思ったんだけど,任意の整数における約数となると「そりゃあ無理だ……」という話.

例えば,100という整数の場合,

import sympy

for i in sympy.divisors(100):
    print(i, '\tx\t', 100//i)

1 x 100
2 x 50
4 x 25
5 x 20
10 x 10
20 x 5
25 x 4
50 x 2
100 x 1

12の場合,

for i in sympy.divisors(12):
    print(i, '\tx\t', 12//i)

1 x 12
2 x 6
3 x 4
4 x 3

6 x 2
12 x 1

これら約数ペアの合計が最小になるのは,大体真ん中になるのはみれば分かる.
任意の整数が平方数であれば,平方根を求めれば丁度真ん中になる.
つまり,平方数であれば,という条件があれば,O(1)になる.

任意の整数を考えるなら,普通に上から順に計算していけばO(n)だけど,
平方数であれば,平方根を求めれば丁度真ん中である事が分かっているので,
平方数じゃなくても,平方根を求めれば,大体真ん中である事が分かる.
大体真ん中,から,探索していけば,O(sqrt(n))になる.

import math


m = 100

for i in range(int(math.sqrt(m)), m):
    if m % i == 0:
        print(i, int(m / i))
        break

10 10

m = 12

for i in range(int(m**(1/2)), m):
    if m % i == 0:
        print(i, int(m / i))
        break

3 4

%%timeit
for i in range(math.ceil(math.sqrt(m)), 0, -1):
    if m % i == 0:
        #print(i, int(m / i))
        break

1.3 µs ± 14.6 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%%timeit
for i in range(int(m**(1/2)), m):
    if m % i == 0:
        #print(i, int(m / i))
        break

1.16 µs ± 19.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

そもそも,任意の整数の約数について,O(1)で分かる方法があるとすれば,それはO(1)で素因数分解が解けるって話で,この問題の本質は,日本語では余り聞かないけど「整数分解(integer factorization)」の困難性と同じ.

という話で揉めているのかなあ.よく分からないけど.

数理的には,「O(sqrt(n))」という話で終わりだけど,じゃあモジュール使いまくってでも,Pythonで最も高速に処理をするには?という話に置き換えれば.

素数なら → 1 x m

平方数なら → sqrt(m) x sqrt(m)

それ以外なら → 約数の合計が最小なものを返す

と処理するのが速いんじゃないだろうか.

from sympy import divisors, isprime


m = 1001
sqrt_m = int(m**(1/2))

if isprime(m):
    print(f'1 x {m}')
elif m % sqrt_m:
    a = divisors(m)
    print(f'{a[len(a)//2 - 1:len(a)//2  + 1]}')
else:
    print(f'{sqrt_m} x {sqrt_m}')

m = 12 100 541 900 1681 1987 10000について比較.

1.3 µs ± 140 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

27.4 µs ± 2.87 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

1.17 µs ± 28 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

2.1 µs ± 128 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

63.1 µs ± 2.11 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

3.15 µs ± 343 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

1.22 µs ± 10.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

2 µs ± 172 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

1.27 µs ± 40.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

2.96 µs ± 201 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

261 µs ± 23.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

3.01 µs ± 190 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

1.22 µs ± 22.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

1.98 µs ± 163 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

ちなみに,

%%timeit

if isprime(m):
    a = 1
    #print(f'1 x {m}')
else:
    for i in range(int(m**(1/2)), m):
        if m % i == 0:
            break
    #print(f'{a[len(a)//2 - 1:len(a)//2  + 1]}')

3.03 µs ± 319 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

広告
カテゴリー: 未分類 パーマリンク

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中