yukicoder : No.413 +5,000,000pts

これは自分の環境にはつらい…。
http://yukicoder.me/problems/no/413

問題

正整数dに対しt^2+t≦dとなる最大の整数tを求める関数として、C++11で以下の関数calcを作成した。
1≦d≦10^18の範囲で、下記関数が不正解を返すdを10^5個列挙せよ。

typedef long long ll;

ll calc(ll d) {
  return (ll)((-1 + sqrt(1 + 4*d)) / 2.0);
}

解法

普通doubleは仮数部が52bitなので、(1+4*d)が52bit以上となるような整数を指定すると精度の問題でsqrtが返す値が真の平方根とずれる。
そのためには1+4*dが2^58前後となるdを適当に10^5個列挙すればよい。

なお、上記calc関数はgcc 32bit環境ではsqrtが誤差を生じない点に注意。
gccのデフォルトでは32bit環境では小数の途中計算を80bitのレジスタで行うため、仮数部が64bitありlong longの値の精度を維持出来てしまう。
(自分は32bit OS利用のため、これに気づくまで一向にcalcが不正な値を返す例が見つからず大幅に時間を食った。)

32bit環境では以下の対策で小数の計算を64bitで行うことで回避できる。一番下は試していない。

• g++のオプションに-mfpmath=sseを付ける(64bit OSだとデフォルトがこちらなので気にしなくてよい)
• 一旦値をdouble型の変数に代入し、確実に仮数部53bitによる誤差を生じさせる(以下のコードではvolatileを用いて、コンパイラによる最適化を回避している。)
• なぜか途中で意図的にlong double版sqrtを呼ぶと、通常のdoubleのsqrtが誤差を生じたが要因不明。
• FPUのコントロールレジスタをいじって途中計算も64bitにさせる。

本番は色々試しているうち偶然3番が発動して解にたどり着いた。

void solve() {
	int i,j,k,l,r,x,y; string s;
	set<ll> R;
	
	FOR(i,100) {
		FOR(j,1000000) {
			ll a=j+(1LL<<29);
			a=a*a-50+i;
			ll d=(a-1)/4;
			
			//64bitなら以下のvolatile云々は不要でコメント通り1行で書ける。
			//32bit環境だと-msse2 -mfpmath=sseを付けるか、FPUのレジスタいじって80it→64bit精度にするか
			//以下のコードのように無理やりdoubleにキャストしてどうにかする
			
			//ll t=(ll)((-1 + sqrt(1+4*d)) / 2.0);
			volatile double c = 1 + 4*d;
			volatile ll b = c;
			ll t=(ll)((-1 + sqrt(b)) / 2.0);
			
			if(t*t+t>d || (t+1)*(t+1)+(t+1)<=d) {
				R.insert(d);
				if(R.size()==100000) {
					FORR(r,R) _P("%lld\n",r);
					return;
				}
			}
		}
	}
	return;
}

まとめ

勉強にはなったけど…。