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 <- [0..n-1], mod (n*a) (a*a+b*b) == 0, mod (n*b) (a*a+b*b) == 0] sumDiv = sum.map f.divisors where f (a,b) | b == 0 = 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(/=0)$ [m * div lim (m*(a^2+b^2))| m <- [1..]] realSum lim = sum' [n*div lim n | n <- [1..lim]] sum' = foldl' (+) 0 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の効率が良くないわけですよ。
では効率化しよう。
といきたいところだが、もう寝る時間なので今日はここまで。
[追記]
少しだけ手を加えた。
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]
public class P153{ static long lim = 100000000; public static long gcd(long a,long b) { if (b==0)return a; else return gcd(b,a%b); } public static void main(String[] args){ long sum = 0; 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==0? a: gcd(b,a%b); } int main(){ long long sum = 0; 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 == 0 ? 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 = 0; 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 == 0 ? 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 = 0; 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秒未満になった．
