ＰＥ１２ 三角数の約数の数を求める処理が遅い

Question

###追記
解決したのは数ヶ月前でしたが別言語で改めて実装してみたところ
@swordoneさんの回答に近い実装で今のところ最速となりました。
当時はまだ理解が足りていなかったのでBAを変更させていただきました。

要点としては・・・
１．約数の個数は素因数分解で求める
２．１つ前で求めた解をキャッシュして使う　※@swordoneさんコードのlow,high変数

約数の個数と素因数分解の参考情報
約数の個数の公式と平方数の性質
素因数分解と約数の個数と総和の求め方を説明!|数学勉強法

結果が出るまで時間が掛かり過ぎる

参考にした例題元　リンクURL
結果が出るまで10分以上掛かり何度も中断しているので改善点などアドバイス頂けたら幸いです。
※ソースコードは結果が即出るように変数を設定しています。

ソースコード

	/**
	 * 
	 * 問題 https://projecteuler.net/problem=12
	 * ※以下Google翻訳
	 * 三角数のシーケンスは、自然数を加算することによって生成されます。
	 * したがって、7 番目の三角形の数は1 + 2 + 3 + 4 + 5 + 6 + 7 = 28となります。
	 * 最初の10の項は次のようになります。1、3、6、10、15、21、28、36、45、55、...
	 * 最初の7つの三角数の要素をリストアップしましょう。
	 * 
	 * 1：1 
	 * 3：1,3 
	 * 6：1,2,3,6 
	 * 10：1,2,5,10 
	 * 15：1,3,5,15 
	 * 21：1,3,7,21 
	 * 28：1,2,4,7,14,28
	 * 
	 * 28が5つの除数を持つ最初の三角形の数であることが分かります。
	 * 最初の三角形の数には何百もの約数がありますか？
	 * 
	 * 原文：What is the value of the first triangle number to have over five hundred divisors?
	 * 
	 * ---以下推測---
	 * 
	 * 28: 1,2,4,7,14,28 
	 *       2,4,7,14,28 ===>>> 5
	 *       
	 * over five hundred divisors?      
	 * 約数が500以上の三角数を求める?
	 * 
	 */
	public static void main(String[] args) {

		long start = System.currentTimeMillis();

//		func(5, 1); 		//0ms		//Answer:28
//		func(300, 1); 		//2245ms	//Answer:2162160
		func(500,12375);	//205ms		//Answer:76576500//即答用
//		func(500,1);		//未達

		long end = System.currentTimeMillis();
		System.out.println((end - start) + "ms");
		System.out.println();

	}

	/**
	 * num個以上の約数を持つ最初の三角数を求める
	 * @param num : 約数の数
	 * @param cnt : 探索開始初期値1
	 */
	static void func(int num, int cnt) {

		//探索開始初期値
		int counter = cnt;
		//		即答用		
		//		int counter = 12375;
		//		ログ c:12375 t:76576500 d:575
		//		Answer:76576500

		while (true) {

			//三角数
			int triangle = 1;

			for (int i = 2; i <= counter; i++) {
				triangle += i;
			}

			//奇数はスキップ
			if (triangle % 2 == 1) {
				counter++;
				continue;
			}

			//約数
			int divisors = 0;

			for (int i = 2; i <= triangle; i++) {
				if (triangle % i == 0) {
					divisors++;
				}
			}

			//ログ
//			System.out.println("Log  c:" + counter + " t:" + triangle + " d:" + divisors);

			//出力
			if (divisors >= num) {
				System.out.println("Answer:" + triangle);
				return;
			}
			counter++;
		}
	}

出力結果（※即答用に変数を指定した場合）
Answer:76576500
204ms

---

試したこと

			//奇数はスキップ
			if(triangle % 2 == 1) {
				counter++;
				continue;
			}

• 大きい解は偶数のみなので暫定的に奇数はスキップするようにしたところ300個探索で約2倍速度upしました。
• counter探索開始初期値を変更することで解を即答出来るようにしました。

わからないこと

• 例題なので答えがすぐに表示される解法があるのではないか
• ボトルネックになっているのはどこなのか

初歩的な総当たりループで探索し初期値を手動で変更などしているので本来の求め方では無いと思います。
指数倍的に時間が掛るのはループに原因があるのかなと思いますがアドバイスやヒントなど頂けましたら幸いです。

Comments

[退会済みユーザー]

等差数列の和は　（最小値＋最大値）＊（最大値／２）でだせるからループなんて必要ないんだけどね

[opyon]

ありがとうございます。n(n+1)/2の公式を取り入れ作り直すとこまで出来たので、追加でその公式も考慮してみたいと思います。

[opyon]

と思ったのですが、n(n+1)/2の公式が三角数ですでに等差数列の和ですよね？ 約数の個数を求めるのに（例えば6なら1,2,3,6の計4個） 三角数＝等差数列の和なので敢えて使うところが無いのかなと思っています。 こちらの勘違いだったらご指摘ください。

[退会済みユーザー]

The sequence of triangle numbers is generated by adding the natural numbers. So the 7th triangle number would be 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28. The first ten terms would be:

[退会済みユーザー]

この部分

[退会済みユーザー]

三角数を求めてる場所がループで求めてるからループする必要ないってこと

[opyon]

はいそうです。その代りにn(n+1)/2を公式に取り入れて作り直したのでご確認ください。コメントした時にはコード掲載してなかったので説明が分かりづらかったと思います。すみません。

Answer 1

約数の個数の性質として、「互いに素なaとbについて、a*bの約数の個数は、aの約数の個数×bの約数の個数に等しい」ということがあります。

三角数はn(n+1)/2となりますが、nとn+1は互いに素なので、「nとn+1のうち、奇数の方」と「nとn+1のうち、偶数を2で割ったもの」で分けて約数の個数を計算して、それをかければ、一気に計算量を減らせます。

    static int dividors(int num) {
        int divisors = 1;

        for (int i = 2; i <= num; i++) {
            if (num % i == 0) {
                divisors++;
            }
        }
        return divisors;
    }
    
    static void func(int num, int cnt) {
        for(int n = cnt; ;++n) {
            int x, y;
            if(n % 2 == 1) {
                x = n;
                y = (n + 1) / 2;
            } else {
                x = n / 2;
                y = n + 1;
            }
            if(dividors(x) * dividors(y) > 500) {
                System.out.println("Answer:" + ((long)x * y));
                return;
            }
        }
    }

paiza.ioで実行した結果

Comments

[opyon]

ありがとうございます。
公式を利用した計算を考えてみます。

[opyon]

ベストアンサーとても迷いましたが@Stars1024さんに送りました。
一人にしか送れないのがとても惜しいです。
奇数偶数で分けたり、平方根までの計算だったり、数学やアルゴリズムに慣れている方のアドバイスはとても参考になります。
ありがとうございました。

Answer 2

参考にした問題のページにアクセスして、実際に問題を見たところ,「約数の個数が500を超える最小の三角数を求める」と書いてあったのでそれを解くコードを作成しました。

約数の個数を求めるのに工夫したところ

例えば、nが2で割り切れたら、n / 2もnの約数です。(約数の性質)
つまり、nの平方根の整数部分までの約数の個数を求めそれを2倍すればいいのです。
ここで、注意すべき点はnが平方数の場合です。~~この場合はさらに1加えればいいでしょう。
~~
この場合は平方数の部分が重複するので1減らす必要があります。

コードはまだ作成していないけれど、他に思いついた方法について

約数の個数を求める方法なら素因数分解して求めるのが一番 はやいと思います。
ただ、まだ実装できていないので、まだそこは考えます。
素因数分解された形に表示するだけならできますが、それぞれの因数の個数を求めるとなると
おそらく配列に素数を順に格納させ、それぞれの個数を数え上げる必要があると思います。
その時の処理速度が気になりますが....

とりあえず、求めたい個数が500を超える最小の三角数を求めることはできたので
その部分のコードを載せておきます。

import java.util.*;

public class Main {
    public static void main(String[] args){
        
        Scanner sc = new Scanner(System.in);
        int i = 1;
        while(number_divisors(i * (i + 1) / 2) < 500){
            i++;
        }
        System.out.println(i + "番目の三角数:" + i * (i + 1) / 2);        
    }
    public static int number_divisors(int n)
    {
        int sum = 2;
        for(int i = 2; i <= Math.sqrt(n); i++){
            if(n % i == 0){
                sum += 2;
            }
        }
        if(is_square(n)){
            sum--;
        }
        return sum;
    }
    public static boolean is_square(int n)
    {
        return n == (int)Math.pow((int)Math.sqrt(n),2);
    }
}

<実行結果>
12375番目の三角数:76576500

Comments

[退会済みユーザー]

素因数分解をして約数の個数を求める方法を実装し、もしそちらの方の処理速度が速かったら、そちらを載せます。

[opyon]

ありがとうございます。
公式を利用した計算を考えてみます。
コードはチラ見だけして丸パクリはしないように注意します。
どうしても詰まったらカンニングさせてもらいますね。

[opyon]

@Stars1024さんの平方根まで求め2倍にする方法と、
@maisumakunさんの奇数偶数で分ける方法と、
本質的には同じような意味なのでしょうか？
数学に疎いもので理解が足りずすみません。

[opyon]

よく考えたら奇数偶数で分けてもそれぞれの約数の個数を計算しているので、
試しに平方根まで求め2倍にする方法を実装してみました。
爆速です！
コード修正したものを後ほど差し替えておきます。

[退会済みユーザー]

>@Stars1024さんの平方根まで求め2倍にする方法と、
@maisumakunさんの奇数偶数で分ける方法と、
本質的には同じような意味なのでしょうか？

私が書いた「平方根まで求める」は約数の性質を使っており、
maisumakunさんがおっしゃった方法は約数の個数の性質を使っています。
maisumakunさんの回答にもありますが、「互いに素である２つの自然数m,nに対して
(mの約数の個数) * (nの約数の個数) = (m*nの約数の個数)」が成り立ちます。
ｎ,n + 1は互いに素であるのですが、
(i) nが偶数 n + 1が奇数
(ii)nが奇数 n + 1が偶数
の2つの場合があります。そのため、「偶数奇数で分ける」必要があります。

結論を言うと、厳密には異なると思います。前述の通り「約数の性質」と「約数の個数の性質」だからです。

[opyon]

なるほどそうですね、数学的に厳密な理解は出来ていませんが、実装してみることで得たことは多いです。
「約数の性質」と「約数の個数の性質」の合わせ技で10分以上でも答えが出なかったものがわずか7msで出るようになりました。
ありがとうございました。

[opyon]

今更なのですが
func(500)などの大きい数値なら問題なく答えが出るのですが、
func(5)などの小さい数値で求めると正しい結果が得られないことに気づきました。
@Stars1024さんの「平方根まで求める」部分のコードをそのまま組み込んでも同じなので調べているところです。
ちなみに約数の個数を全ループで探す元の方法に戻すと大丈夫でした。

[退会済みユーザー]

func(5) = 28であっていますよね?　確認したところ私が書いたコードでは 28と出力されました。

[opyon]

@Stars1024さんの全コードをまとめて使うと大丈夫だと思います。
「平方根まで求める」部分のみ利用させて頂くと結果が違うのでこちらの実装方法と使い方が間違っているのだと思います。

[opyon]

＠Stars1024さんの全コードで試してみました。n=1,2,3で試すと全て1になります。
n=4で6の結果が出ますので、3までの処理が必要かもしれません。

「平方根まで求める」機能をなんとか自分でも理解して実装してみたいと思います。

[x_x]

> この場合はさらに1加えればいい
引くのでは？ 平方根を二重にカウントしているのですよね？

[opyon]

（元が訂正されたので他の方が見る時に誤解の無いように削除しておきます）

[opyon]

すみませんリンクありましたので貼っておきます。
https://mathtrain.jp/numberofd

[x_x]

number_divisors()のことを言っているのですが……
number_divisors(3) // 2 ……あってる
number_divisors(4) // 5 ……3が正しい
number_divisors(5) // 2 ……あってる
……
number_divisors(12) // 6 ……あってる
number_divisors(16) // 7 ……5が正しい

[退会済みユーザー]

すみません。x_xさんのおっしゃる通り平方数の場合は約数がその平方根になるとき重複するので
１減らす必要があります。
最初の平方根までやる部分で
for(int i = 2; i <= Math.sqrt(n); i++)
と i<=Math.sqrt(n)のように 「<=」不等号に=がついてるためです。見落としていました。
x_xさん、ご指摘ありがとうございます。

[退会済みユーザー]

その個所の訂正はしました。

[opyon]

両者様ともありがとうございます。

[opyon]

検索で見に来た方の為にBAを変更させていただきました。

Answer 3

三角数は1からnまでの総和なのでn*(n+1)/2で求められます。
1番目:12/2=1
2番目:23/2=3
3番目:3*4/2=6
という具合に。

そして約数の個数は、素因数分解できれば、各素因数の個数+1をすべて掛けることで求められます。
また、隣り合った自然数は互いに素なので、その積の約数の個数は各自然数の約数の個数の積になります(maisumakunさんの回答通り)。

これらの事実を利用すれば、各自然数の約数の個数を素因数分解しながら求めることで、高速に求められるかと思います。

private static int numOfDivisor(int n) {
    if (n < 2) return 1;
    int answer = 1;

    // 偶数だった場合は2で割る(ビットシフト)
    if ((n & 1)  == 0)
        n >>= 1;

    // nが√Integer.MAX_VALUE以上の素数だったりすると、この継続条件は良くないが、
    // 簡単なのでこちらで進める
    for (int i = 2; i * i <= n; i++) {
        int count = 0;
        while (n % i == 0) {
            n /= i;
            count++;
        }
        answer *= count + 1;
    }
    // forを抜けて残ったnが1でないときは、そのnは素数なので1+1=2を掛ける
    if (n != 1)
        answer *= 2;
    return answer;
}

で、500個以上の三角数を探すのはこういうコード

public static void main(String[] args) {
    int n = 500;
    System.out.println(n + "個以上の約数を持つ最初の三角数は" + func(n));
}

public static int func(int num) {
    // 2番目の三角数2*3/2=3に対応する、low:2/2=1の約数の個数1、high:3の約数の個数2
    int low = 1;
    int high = 2;
    int n = 2;
    while (low * high < num) {
        n++;
        low = high;
        high = numOfDivisor(n + 1);
    }
    return n * (n + 1) / 2;
}

ideoneでの実行結果
0.09秒で解決

Comments

[opyon]

ありがとうございます。
すでに7msで結果出ているのでと思ったのですがコードを読んでみると知らない技法など使っていてとても参考になります。

[opyon]

かなり前の回答ですが・・・
別言語のC++で改めて速度を追求してみたところ、
こちらの回答に近い実装で最速となりました。(平均で約１ミリ秒)

検索で見に来た方の為にBAを変更させていただきました。

Answer 4

python での回答例を web で検索し、速度を比較してみました。

１. https://qiita.com/hoxo_m/items/ad17361b7276f1c6527b
２. http://freak-da.hatenablog.com/entry/20100729/p1
３. https://qiita.com/yukkyo/items/29cd0d4c87ba29819f40
４. https://stackoverflow.com/questions/17743851

どれも 一瞬で終わるレベルですが、 このなかでは 2 番目のものが一番 速そうです。
アルゴリズムを研究して java に移植するとよいと思います。
上のものは "Project Euler 12 python" で google 検索してみつけたものですが、 "Project Euler 12 java" で検索すれば、 java での解をいろいろ見つけることができると思います。

py12_1.py

# See https://qiita.com/hoxo_m/items/ad17361b7276f1c6527b
import math

# limit = 5
limit = 500

def compute_factors(num):
  factors = [1]
  if (num == 1): return factors
  for i in range(2, int(math.floor(math.sqrt(num))) + 1):
    if(num % i == 0):
      factors.append(i)
  factors_rev = list(factors)
  factors_rev.reverse()
  for i in factors_rev:
    fac = num // i
    if(fac not in factors):
      factors.append(fac)
  return factors

triangular_nums = [1]
factors = []
n = 2
while True:
  triangular_num = triangular_nums[n - 1 -1] + n
  triangular_nums.append(triangular_num)
  factors = compute_factors(triangular_num)
  if len(factors) >= limit:
    break
  n += 1

result = triangular_nums[-1]

print(result)
print(result == 76576500)
print(n)
print(len(factors), factors[:10])

py12_2.py

# http://freak-da.hatenablog.com/entry/20100729/p1

import math
import sys

fact = {1:1, 2:2} # 約数の個数をキャッシュする辞書

def num_factors(n):
    if n in fact: return fact[n]

    r = int(math.sqrt(n))

    num = 0
    if r * (r + 1) == n:
        num = num_factors(r) * num_factors(r + 1)
    else:
        i = 1
        while i < r:
            if n % i == 0:
                num += 1

            i += 1

        num *= 2
        if r * r == n: num += 1

    fact[n] = num
    return num

def main():
    n = 2
    while True:
        if n % 2 == 0:
            fn = num_factors(n / 2) * num_factors(n + 1)
        else:
            fn = num_factors(n) * num_factors((n + 1) / 2)

        if fn >= 500:
            print("%d = (%d * %d) / 2, it has %d factors" % (n * (n + 1) / 2, n, n + 1, fn))
            return

        fact[n * (n + 1) / 2] = fn

        n += 1

if __name__ == '__main__':
    main()

py12_3.py

# See  Phttps://qiita.com/yukkyo/items/29cd0d4c87ba29819f40

# problem 12
import math
import collections

# 数字を素因数分解する
def trial_division_sqrt(n):
    prime_count = collections.Counter()
    for i in range(2, int(math.sqrt(n)) + 2):
        while n % i == 0:
            n /= i
            prime_count[i] += 1
    if n > 1:
        prime_count[int(n)] += 1
    return prime_count

# 約数の個数を求める
def count_divisor(num):
    ans_count = 1
    count_one = 0
    # 素因数分解して各素数の個数のリストを求める
    prime_count = trial_division_sqrt(num)

    # 各リスト内の個数をカウント
    for v in prime_count.values():
        if v == 1:
            count_one += 1
        else:
            # 同じ素数が複数個含まれている場合
            ans_count *= v + 1
    return ans_count * pow(2, count_one)

ans_num = 0
count = 0
divisors = 500

while(True):
    count += 1
    ans_num += count
    if count_divisor(ans_num) > divisors:
        break
ans_num

py12_4.py

# See https://stackoverflow.com/questions/17743851

import math

def get_factors(n):
    return sum(2 for i in range(1, round(math.sqrt(n)+1)) if not n % i)

def generate_triangles(limit):
    l = 1
    while l <= limit:
        yield sum(range(l + 1))
        l += 1

def test_triangles():
    triangles = generate_triangles(100000)
    for i in triangles:
        if get_factors(i) > 499:
            return i

print(test_triangles())

Comments

[opyon]

ありがとうございます。
検索でここへ辿り着いた方へ参考になると思います。
自分は他言語覚えるほど余裕はないですが他言語覚える時は参考にさせていただきます。

[opyon]

2番が速いとのことでしたのでリンク先の解説読んでみたところ
主に①「三角数の公式」②「公式nの偶数判定」③「公式n+1の再利用」を実装してるようです。
自分もアドバイス頂いた①と②は実装出来ましたが③はまだだったのでやってみます。

[opyon]

2番目を参考に③「公式n+1の再利用」を実装してみましたが、逆に遅くなりました。
7ms⇢136ms
追加したのはmapだけですので全体のコードは省略しますが読み書きでオーバーヘッド？が発生してるのかなと思います。

static HashMap<Integer, Integer> map = new HashMap<>();

単純に1,3,6,10毎に約数をmap.putして再利用時にmap.getするだけの3行追加しただけです。
再利用の理屈はよくわかったのですが再利用するためにアクセスするより単純計算させておいたほうが速いようです。

また別の機会で改善策として利用できそうなのでとても参考になりました。
ありごとうございました。

Answer 5

Javaなど不要。そう、シェル芸ならね！

$ seq 100000|awk '$0=$1+s,s=$0'|gfactor|awk '{delete a;for(i=2;i<=NF;i++){a[$i]++};s=1;for(i in a){s*=a[i]+1};if(s>=500){print $1,s"個";exit}}'
76576500: 576個

なお実行時間

$ time (seq 100000|awk '$0=$1+s,s=$0'|gfactor|awk '{delete a;for(i=2;i<=NF;i++){a[$i]++};s=1;for(i in a){s*=a[i]+1};if(s>=500){print $1,s"個";exit}}')
76576500: 576個

real	0m0.080s
user	0m0.106s
sys	0m0.016s

Comments

[opyon]

ありがとうございます。1行？で解を求められるなんて凄いですね。

Answer 6

頂いたアドバイスで作り直したところ、かなり改善できました。
ありがとうございました。
※追加で平方根の処理を実装したところ更に改善できました。

---

	public static void main(String[] args) {
		long start = System.currentTimeMillis();

		func(500);

		long end = System.currentTimeMillis();
		System.out.println((end - start) + "ms");
	}

	static void func(int num) {

		int n = 0; //自然数n
		int t = 0; //n番目の三角数（総和）

		//約数の個数の性質「互いに素なaとbについて、
		//a*bの約数の個数は、aの約数の個数×bの約数の個数に等しい」
		//nとn+1は互いに素なのでnかn+1のどちらかが偶数になる

		int x = 0; //(偶数奇数で分けて偶数の方を2で割る)
		int y = 0; //(偶数奇数で分けて偶数の方を2で割る)
		int d = 0; //分けて計算した約数の個数を掛けた総数

		while (d < num) {

			//自然数n
			n++;

			//三角数の公式
			//n(n+1)/2

			//n       1 2 3  4  5  6  7  8  9 10
			//三角数  1 3 6 10 15 21 28 36 45 55
			//約数　  1 2 4  4  4  4  6  9  6  4

			//偶数奇数判定
			if (n % 2 == 0) {
				x = n / 2;
				y = n + 1;
			} else {
				x = n;
				y = (n + 1) / 2;
			}

			//個別に約数の個数を求め掛け合わせて総数を求める
			d = getD(x) * getD(y);
		}

		//ログ出力用
		t = n * (n + 1) / 2;

		System.out.print("last n:" + (n));
		System.out.print(" t:" + t);
		System.out.print(" d:" + d);
		System.out.println();
		System.out.println("Ans:" + t);
	}

	//約数の個数を求める
	static int getD(final int t) {

		int d = 0;

		final int rootT = (int) Math.sqrt(t);

		for (int i = 1; i <= rootT; i++) {
			if (t % i == 0) {
				d++;
			}
		}

		d *= 2;

		if (t == rootT * rootT) {
			d--;
		}

		return d;
	}

出力結果

last n:12375 t:76576500 d:576
Ans:76576500
0ms