AtCoderABC253　D

まずそもそもの問題からです。
D - FizzBuzz Sum Hard

問題文
1 以上 N 以下の整数であって、A の倍数でも B の倍数でもないものの総和を求めてください。

解法としては、
1~Nまでの総和から、同区間のAの倍数の総和、同区間のBの倍数の総和を引いて、
重複するAとBの最小公倍数の倍数の総和を足せばokです。

詳しくは解説を見てください。
解説

間違いのコードと正解のコード

わたしが最初に提出した間違いのコードはこのようになります。

間違い

import fractions
 
def wa(N):
  Ans = int((N+1)*N/2)
  return Ans
 
N,A,B=map(int,input().split())
S=0
 
ax=int(N/A)
bx=int(N/B)
G=int( A * B /fractions.gcd(A, B))
gx=int( N/(A * B /fractions.gcd(A, B)))
 
S=wa(N) -wa(ax)*A -wa(bx)*B + wa(gx)*G
 
print(S)

正解（AC）のコードはこちらになります

正解

import fractions
 
def wa(N):
  Ans = ((N+1)*N)//2
  return Ans
 
N,A,B=map(int,input().split())
S=0
 
ax=N//A
bx=N//B
G= A * B //fractions.gcd(A, B)
gx=N//G
 
S=wa(N) -wa(ax)*A -wa(bx)*B + wa(gx)*G
 
print(S)

コード上での間違いの箇所は

間違いの箇所

##間違い
int( (n*1)*n/2 )

##正解
((N+1)*N)//2

です。

Pythonはint(N)でNの少数を切り捨てできるので、//と同じ結果になるはずです。
そこで、一体どんなNをwa(N)に入力すると、違いが発生するのかをAtCoderのコードテストで検証しました。

間違いになる数字の調査

##Python3だと遅くて終わらないのでPyPy3で実行
for i in range(10**9):
  if(  (i+1)*i//2 != int((i+1)*i/2) ):
    print(i,((i+1)*i)//2 , int((i+1)*i/2))

出力

134217729 9007199456067585 9007199456067584
134217730 9007199590285315 9007199590285316
134217733 9007199992938511 9007199992938512
134217734 9007200127156245 9007200127156244
134217737 9007200529809453 9007200529809452
...

答えが9007199456067584より大きい数字で起きることがあることがわかりました。
ここでおそらく型の問題だなあと思いました。

Floatのこと　IEEE754 double

PythonのFloatはC++でいうところのDoubleです。
これはIEEEで決まっていて、64bitに対して以下のようなデータ構造与える方法です。

Float = (-1)^{Sign}\times1.fraction\times10^{exponent-11 1111 1111}\\
-----------------
\\
ex)\\
sign=0\\
fraction = 1000 ... 0000\\
exponent = 100 0000 0011\\
Float = (-1)^{0}\times1.1000...0000\times10^{100 0000 0011 - 11 1111 1111}\\=24_{(10)}

ここで注目したいのが、仮数部（fraction）が52bitということです。
試しに$2^{53}$を計算してみると　(仮数部は1.nを表しているため)

$2^{53}=9,007,199,254,740,992$

デジャブ感がありますね。
なぜこの数字を超えると、バグが出るのかそれを考えてみます。

Floatの精度

結論から言うと
$2^{53}-1$より大きい数字は1刻みではなく、2刻み、4刻みでしか表現できません。
これがFloatの精度です。

原理としては単純です。
仮数部の最大値である1.111...1111に対して、$2^{52}$をかけると、指数部の小数点が右へ52回移動するため、ちょうど整数になります。
この数字が$2^{53}-1=9,007,199,254,740,991$となります。

今回の問題の場合、$n(n+1)$は偶数のため、$2^{54}$までは正しくfloatで表現できます。

for i in range(10**9):
  if( (i+1)*i > 2**54 ):
    break
  elif(  (i+1)*i//2 != int((i+1)*i/2) ):
    print(i,((i+1)*i)//2 , int((i+1)*i/2))

試しにコレを回すと、何も出力されません。2**55にすると出てきます（それはそう）

そして、問題となった9007199456067585について　/1をすることでfloatの挙動を観察すると

A=9007199456067585
print(A/1)
print(A+1/1)
print(A+2/1)
print(A+3/1)

出力

9007199456067584.0
9007199456067584.0
9007199456067586.0
9007199456067588.0

となります。
なんでこのような出力になるのかは厳密にはわかりませんがPythonは小さい方に丸められるらしいので、
...86/1の際は丸めで84へ、...87/1の際も丸めで86へっていうことでしょうか。
...86/1の丸めが納得いきませんが多分そうなんでしょう。

まとめ

なんでこんな丸め方になるのかはわからないのですが、とりあえずint(a/b)とa//bの違いについて
基本はa//bを使うっていうことがわかりました。
a//bは素直にintで計算してくれます。

ううう・・・D問題行けたのに・・・
コレのせいでレート下がりました・・・