Problem 153

http://projecteuler.net/index.php?section=problems&id=153

とりあえず書いたが、速く動かない。（10^8ってかなり大きい気がするんですが。）

import Number
import Data.List
-- ナーイブな実装
divisors n = [(a,b)| a <- [1..n], b <- [..n-1], mod (n*a) (a*a+b*b) == , mod (n*b) (a*a+b*b) == ]
sumDiv = sum.map f.divisors
where f (a,b) | b ==  = a
| otherwise = 2*a
sumDivs n = sum.map sumDiv $ [1..n]
gussianSum lim = sum' [appearSum a b lim | a <- [1..sqrt'$div lim 2], b <- [a..sqrt'$lim-a^2], gcd a b == 1]
where sqrt' = floor.sqrt.fromIntegral
appearSum a b lim | a == b = 2*a*appear a b lim
| otherwise = 2*(a+b)*appear a b lim
appear a b lim = sum'.takeWhile(/=)$ [m * div lim (m*(a^2+b^2))| m <- [1..]]
realSum lim = sum' [n*div lim n | n <- [1..lim]]
sum' = foldl' (+) 
main = print$ realSum (10^5) + gussianSum (10^5)

そこで、プロファイルしてみると

Mon Dec 22 02:24 2008 Time and Allocation Profiling Report  (Final)
p153.exe +RTS -p -h -RTS
total time  =        0.32 secs   (16 ticks @ 20 ms)
total alloc =  68,996,996 bytes  (excludes profiling overheads)
COST CENTRE                    MODULE               %time %alloc
appear                         Main                  56.3   65.6
sum'                           Main                  18.8   12.8
realSum                        Main                  12.5   15.1
gussianSum                     Main                  12.5    5.9
individual    inherited
COST CENTRE              MODULE                                               no.    entries  %time %alloc   %time %alloc
MAIN                     MAIN                                                   1           0   0.0    0.0   100.0  100.0
main                    Main                                                 174           1   0.0    0.0     0.0    0.0
CAF                     Main                                                 168           7   0.0    0.0   100.0  100.0
appear                 Main                                                 184           0   0.0    2.6     0.0    2.6
gussianSum             Main                                                 180           0   0.0    0.0     0.0    0.0
main                   Main                                                 175           0   0.0    0.0   100.0   97.4
gussianSum            Main                                                 178         225  12.5    5.9    81.3   75.0
appearSum            Main                                                 181       23880   0.0    0.6    62.5   66.9
appear              Main                                                 182       23880  56.3   63.0    62.5   66.3
sum'               Main                                                 183       23880   6.3    3.3     6.3    3.3
sum'                 Main                                                 179           1   6.3    2.2     6.3    2.2
realSum               Main                                                 176           1  12.5   15.1    18.8   22.4
sum'                 Main                                                 177           1   6.3    7.3     6.3    7.3
CAF                     GHC.ConsoleHandler                                   166           1   0.0    0.0     0.0    0.0
CAF                     GHC.Num                                              151           1   0.0    0.0     0.0    0.0
CAF                     GHC.Handle                                           113           2   0.0    0.0     0.0    0.0

つまり、関数appearの効率が良くないわけですよ。

では効率化しよう。

といきたいところだが、もう寝る時間なので今日はここまで。

1 [追記]
2 [追記2]
3 [追記3]
4 [追記4]

[追記]

少しだけ手を加えた。

appear a b lim = sum' [m * div lim (m*n)| m <- [1..div lim $ 2*n]]+ t (div lim $ 2*n) (div lim n)
where n = a^2+b^2
t m n = div (n*(n+1)-m*(m+1)) 2

あんまり変わらん。

速い計算機で回してもまだ、約90秒。

変な書き方すれば、もっと速くなるはずだけど…

それもいやだ。

全体は

import Data.List
gussianSum lim = sum [appearSum a b lim | a <- [1..sqrt'$div lim 2], b <- [a..sqrt'$lim-a^2], gcd a b == 1]
where sqrt' = floor.sqrt.fromIntegral
appearSum a b lim | a == b = 2*a*appear a b
| otherwise = 2*(a+b)*appear a b
appear a b = let n = a^2+b^2
in sum [m * div lim (m*n)| m <- [1..div lim $ 2*n]]+ t (div lim $ 2*n) (div lim n)
t m n = div (n*(n+1)-m*(m+1)) 2
realSum lim = sum [n*div lim n | n <- [1..lim]]
main = print$ realSum (10^8) + gussianSum (10^8)

[追記2]

javaで同じアルゴリズムを実装したら一分切った。

public class P153{
static long lim = 100000000;
public static long gcd(long a,long b) {
if (b==)return a;
else return gcd(b,a%b);
}
public static void main(String[] args){
long sum = ;
for(int n=1;n<=lim;sum+=n*(lim/n++));
for(int m=1;m<=lim/2;sum+=2*m*(lim/(2*m++)));
for(int a=1,u=(int)Math.sqrt(lim/2);a<=u;a++)
for(int b=a+1,v=(int)Math.sqrt(lim-a*a);b<=v;b++)
if(gcd(a,b)==1){
long n=a*a+b*b,s=lim/n,t=lim/(2*n);
sum+=2*(a+b)*(s*(s+1)-t*(t+1))/2;
for(int m=1;m<=t;sum+=2*(a+b)*m*(lim/(n*m++)));
}
System.out.println(sum);
}
}

やっぱりこういうのは手続き型向けかな。

[追記3]

c++で実装したらもっと速くなった。

#include <iostream>
#include <cmath>
using namespace std;
long lim = 100000000;
long gcd(long a,long b) {
return b==? a: gcd(b,a%b);
}
int main(){
long long sum = ;
for(int n=1;n<=lim;n++)sum+=n*(lim/n);
for(int m=1;m<=lim/2;m++)sum+=2*m*(lim/(2*m));
for(int a=1,u=(int)sqrt(lim/2);a<=u;a++)
for(int b=a+1,v=(int)sqrt(lim-a*a);b<=v;b++)
if(gcd(a,b)==1){
long n=a*a+b*b,s=lim/n,t=lim/(2*n);
sum+=2*(a+b)*(s*(s+1)-t*(t+1))/2;
for(int m=1;m<=t;sum+=2*(a+b)*m*(lim/(n*m++)));
}
cout << sum << endl;
}

とおもったら、出力が異なる？なぜ？c++のことはよくわからないからなー

[追記4]

和のとりかたを工夫して，計算量を減らした．

さらに，前の計算結果を利用するためにmapを使用．

#include <iostream>
#include <cmath>
#include <map>
using namespace std;
map<int, long long> a;
int gcd(int a,int b) { return b ==  ? a: gcd(b, a % b); }
inline long long sum (long long n) { return n * (n + 1) / 2; }
inline long long ap(int n) {
map<int, long long>::iterator it = a.find(n);
if(it != a.end()) return it->second;
long long c = ;
for (int k = 1, l = n / ((int)sqrt(n) + 1); k <= l; k ++) c += k * (n / k);
for (int k = 1, l = sqrt(n); k <= l; k++) c += k * (sum(n / k) - sum(n / (k + 1)));
return a[n] = c;
}
int main(){
const int u = 100000000;
long long sum = ap(u) + 2 * ap(u / 2);
for (int a = 1, ma = (int)sqrt(u / 2); a <= ma; a++)
for(int b = a + 1, mb = (int)sqrt(u - a * a); b <= mb; b++)
if(gcd(a, b) == 1) sum += 2 * (a + b) * ap(u / (a*a + b*b));
cout << sum << endl;
}

しかし，STLのmapは木で実装されている．hash_mapが欲しい．あった．

unordered_mapというらしい．よくわからないが，まだ，正式採用されていないっぽい．

#include <iostream>
#include <cmath>
#include <tr1/unordered_map>
using namespace std;
tr1::unordered_map<int, long long> a;
int gcd(int a,int b) { return b ==  ? a: gcd(b, a % b); }
inline long long sum (long long n) { return n * (n + 1) / 2; }
inline long long ap(int n) {
tr1::unordered_map<int, long long>::iterator it = a.find(n);
if(it != a.end()) return it->second;
long long c = ;
for (int k = 1, l = n / ((int)sqrt(n) + 1); k <= l; k ++) c += k * (n / k);
for (int k = 1, l = sqrt(n); k <= l; k++) c += k * (sum(n / k) - sum(n / (k + 1)));
return a[n] = c;
}
int main(){
const int u = 100000000;
long long sum = ap(u) + 2 * ap(u / 2);
for (int a = 1, ma = (int)sqrt(u / 2); a <= ma; a++)
for(int b = a + 1, mb = (int)sqrt(u - a * a); b <= mb; b++)
if(gcd(a, b) == 1) sum += 2 * (a + b) * ap(u / (a*a + b*b));
cout << sum << endl;
}

実行時間は5秒未満になった．

コンテンツ

[追記]

[追記2]

[追記3]

[追記4]