[Project Euler] 小数の法則
Project Euler Problem 210 やってたら、とても嫌な事態に遭遇。
何度検算しても「答え違う」と言われる。どんだけ debug してもとれないので仕方ないから他の人の答え拾ってきて計算結果比べてみたら、r という問題の大きさを調整する parameter が 1 から 10^8 のときまで答えが一致した挙句、実際の出題の r = 10^9 で初めて答えが不一致になるという嫌がらせとしか思えない事実が発覚。もっと正確には 1..5000 までは全て一致して、それ以降は 10^4, 10^5, ..., 10^8 で一致。
さて、こういうときは大抵精度の問題です。double が float になってたり、C/C++ だと割り算の前に cast か点忘れてたり (後者は例えば r が整数型のとき 1.0/r を 1/r と書いてたりする…)、Haskell だと Integer が Int になってたり。小さい入力でうまく動いて満足してたら大きい値で死にます。これぞ (浮動) 小数の法則。
(以下、ネタバレ注意)
最終的に、sqrt(n) より小さい最大の整数を返す int_lt_sqrt() が long double 使ってなかったから、って事に落ち着きました。
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <cassert>
using namespace std;
// Returns the greatest integer that is *strictly* less than sqrt (n).
// Note this is smaller than floor(sqrt(n)) if n is a square.
long int_lt_sqrt (long n)
{
long ret = floor (sqrt (n)); // 犯人
if (ret * ret == n)
--ret;
return ret;
}
int
main (int argc, char *argv[])
{
long r = atoi (argv[1]);
double radius = sqrt (2) * r/8;
long radius_sq = r*r / 32;
assert ((r*r) % 32 == 0);
long minj = ceil (- radius);
long maxj = -r/8-1; // floor (r/8 + radius);
// 円からできるだけ大きな正方形をくり出した跡に残る
// 欠片四つのうち、一つの格子点を数える。
// 円は境界を含まない open disk として扱う。
long gyoza = 0;
for (long j = minj; j <= maxj; ++j)
{
long s = int_lt_sqrt (radius_sq - j*j);
gyoza += 2*s + 1;
}
cout << gyoza << endl;
long disk = 4*gyoza + (r/4 + 1)*(r/4 + 1) - 4;
cout << (3*r*r/2 - r/4 + 1 + disk) << endl;
}で、link 先の code は何かサラッと ceil(sqrt()) とか書いてて、なんで大丈夫なのか謎。仮に精度の問題が何らかの理由でクリアされてたとしても、何かこれじゃ円の境界部の格子点も数えてしまいそうな気がするけど、計算結果は合ってる。何か納得いかね。
ちなみに link 先の人と同じく最初は Haskell でやってたけど、Haskell では long double が使えない*1から答えが間違ってる上に、どうしても 48 秒もかかる。上の C++ 版が 2 秒以下で走る事を考えると実に 20 倍以上の差が…。どれか忘れたけど、他にも浮動小数使うとやたら遅いことがあった気がするので、多分 fromInteger とかが入ってると上手く unbox できる箇所を推論出来ないに違いない。でも ghc-core 読めないので確認できませんでしたとさ。
ていうかこの問題が正答率低いのって、大勢が debug めんどくさくなって放り出すからだったりして。
-- Default integer type. Int or Integer.
type DInt = Int
-- Assume throughout that r `mod` 8 == 0.
s r = 2*r^2 + 2*r + 1
-- The open disk with diameter OC.
oc_disk !r = filling_square + 4 * moon r
where filling_square = ((r `div` 4) + 1)^2 - 4 -- (-4) removes corners
-- One of the moon shapes that remain after removing the maximal
-- square from the disk.
moon :: DInt -> DInt
moon r = rec 0 minj maxj
where !radius = sqrt(2) * fromIntegral r8
!r8 = r `div` 8
!radius_sq = (r^2) `div` 32
minj = ceiling (fromIntegral r8 - radius)
maxj = -1
rec acc j max | j > max = acc
| otherwise = rec (acc + maxi - mini + 1) (j+1) max
where maxi = r8 + s
mini = r8 - s
!s = int_lt_sqrt (radius_sq - (j-r8)^2)
-- The greatest integer that is *strictly* less than sqrt(x). Note
-- this is 1 less than the floor if x is a square.
int_lt_sqrt :: DInt -> DInt
int_lt_sqrt n | n < 0 = error $ "int_lt_sqrt: negative value " ++ show n
| s * s == n = s - 1
| otherwise = s
where s = round . sqrt $ fromIntegral n
-- Points O and C.
oc = 2
line_oc r = r + 1
p210 :: DInt -> DInt
p210 r = s r - belt_perp_oc r + oc_disk r - line_oc r + oc
