Pythonの a//b と int(a/b) の違いが少しわかった。AtCoder ABC253-D

こんにちは、佐藤かえでです。

初めてQiitaの記事を書いています。
先日のAtCoderABC253-Dを解いていた際に、解法はあってるのになぜかACにならないということがありました。

だいたいこういうときは、floatとか小数点周りで誤差が出ているのでは？と思いいろいろ調べたのでここに記します。

# AtCoderABC253　D
まずそもそもの問題からです。
[D - FizzBuzz Sum Hard](https://atcoder.jp/contests/abc253/tasks/abc253_d)

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

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

詳しくは解説を見てください。
[解説](https://atcoder.jp/contests/abc253/editorial/4017)

## 間違いのコードと正解のコード
わたしが最初に提出した間違いのコードはこのようになります。
```python:間違い
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）のコードはこちらになります
```python:正解
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)
```
コード上での間違いの箇所は
```python: 間違いの箇所
##間違い
int( (n*1)*n/2 )

##正解
((N+1)*N)//2
```
です。

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

```python:間違いになる数字の調査
##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))
```
```shell:出力
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に対して以下のようなデータ構造与える方法です。

```math 
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)}
```

![IEEE754 double float](https://upload.wikimedia.org/wikipedia/commons/thumb/a/a9/IEEE_754_Double_Floating_Point_Format.svg/1920px-IEEE_754_Double_Floating_Point_Format.svg.png)


ここで注目したいのが、仮数部（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で表現できます。
```python:
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の挙動を観察すると

```python
A=9007199456067585
print(A/1)
print(A+1/1)
print(A+2/1)
print(A+3/1)
```
```shell:出力
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問題行けたのに・・・
コレのせいでレート下がりました・・・