Problem 210

210 Obtuse Angled Triangles
Problem 210 - Project Euler
格子点を数えるタイプの問題。
はじめ、問題をよく読まずに、正方形でなく、円の中の格子点を数えるのかと思った。
あとで、実は正方形だと気がついたが…
結局、円の中の格子点を数えることになった。

-- r is divisible 8
circleGrid :: Integer -> Integer
circleGrid r = 4 * sum (map grid [1..r'])
    where r' = floor $ sqrt 2 * fromIntegral (div r 8)
          grid h = ceiling.sqrt.fromIntegral $ 2 * (div r 8)^2 - h^2

main = print $ div (3*n*n) 2 + circleGrid n - 2 * (div n 8 - 1)
    where n = 10^9

探索範囲が広いから、遅い、遅い。
どうも、遅い。どうしたものか。
ルートをとるのが、マズイのかと思い、

circleGrid' :: Integer -> Integer
circleGrid' r = 4 * loop r' 0 0
    where r' = floor $ sqrt 2 * fromIntegral (div r 8)
          r2 = div (r*r) 32
          loop !h !x !c | x > r'         = c
                        | x^2 + h^2 < r2 = loop h (x+1) (c+h)
                        | otherwise      = loop (h-1) x c

こんなコードも書いてみたが、逆に遅くなった。
遅くなった理由は良く分からないが、この問題は手続き型のほうが良さそう。

追記

というわけで、Cで実装してみた。

#include <stdio.h>
#include <math.h>
#define N 1000000000LL
int main(){
  long long count, r, r2, x;
  count = 0;
  r = (long long) (sqrt(2)*(N/8));
  r2 = (N*N) / 32;
  for(x=1;x<=r;x++) count += (long long) ceil(sqrt (r2 - x*x));
  printf("%lld\n",4 * count + (3*N*N)/2 - N/4 + 2 );
}

はじめは、-lmオプションつけなくて、エラーがでて、途方に暮れていた。
勝手にリンクしてくれれば、良いのに。
もちろん、Haskellよりも、速くなりました。多分10倍ぐらい。