Pandas DataFrame (Series) における sum の違い

お前は何を言っているんだというタイトルだと思いますが、 DataFrame (Series) に対する `sum` に対するこの結果を見てみましょう。

```python
import math
import pandas as pd

df = pd.DataFrame({"a": [1e-7] * int(1e7)})

df["a"].sum()
1.0000000000000173

sum(df["a"])
0.99999999975017

math.fsum(df["a"])
1.0
```

いかがでしたでしょうか。みたいな薄っぺらい商材ブログの文体とは異なり、ストイックな追求をここではしていく。

## 前提知識

そもそも単純だと思われるsum関数において、いくつかの加算アルゴリズムが存在しているという事実に注意が必要だ（と偉そうに書いていてこの記事を書くために調べて初めてきちんと理解した）。具体的には [Kahan summation](https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%8F%E3%83%B3%E3%81%AE%E5%8A%A0%E7%AE%97%E3%82%A2%E3%83%AB%E3%82%B4%E3%83%AA%E3%82%BA%E3%83%A0), [pairwise summation](https://en.wikipedia.org/wiki/Pairwise_summation) が今回の題材である。

## いきなり結論

- pandas Series.sum() は内部で numpy.sum() を呼び、それは pairwise summation を使う。
  - 一方で GroupBy sum() や Rolling sum() method は、kahan summation を使うため、上記と結果が一致しないことがある
- build-in sum関数(CPython)は単純ループ加算（呼び名不明）を使う
- math.fsum は kahan summation を使う

## 調査環境

```python
>>> import sys
>>> sys.version
'3.9.7 (default, Dec  5 2021, 20:50:26) \n[GCC 9.3.0]'
>>> import numpy as np
>>> import pandas as pd
>>> np.__version__
'1.21.4'
>>> pd.__version__
'1.3.5'
```

## 個々の呼び出し元のソースコードの確認
 
## pandas Series の sum method

（注意: 長文です）

`breakpoint()` の step-in で雑にコードをたどると 
pandas/core/generic.py(10708)sum()
pandas/core/generic.py(10437)sum()
pandas/core/generic.py(104390)_min_count_stat_function()
pandas/core/series.py(4359)_reduce()
pandas/core/nanops.py(83)_f()
pandas/core/nanops.py(395)new_func()
pandas/core/nanops.py(546)nansum() -> (591) ndarray.sum()
numpy/core/_methods.py(46)_sum()

```python
# numpy/core/_methods.py(48)
return umr_sum(a, axis, dtype, out, keepdims, initial, where)
```
ここで debug 上のstepで `_methods.py(48)_sum()->1.0000000000000173` と値が返ってきてしまった。

さて、 `umr_sum()` とは事実上 `numpy.core.umath.add.reduce` を呼び出していることが同ファイルから伺える。

```python
# numpy/core/_methods.py (10, 20)
from numpy.core import umath as um
...
umr_sum = um.add.reduce
```

`numpy/core/umath.py` を見ると `__all__` で `"add"` が定義されているので、 `from . import _multiarray_umath` を遡ろうとするが、 `numpy/core/_multiarray_umath` というファイルは存在せず `numpy/core/_multiarray_umath.cpython-39-x86_64-linux-gnu.so` のバイナリファイルがあった（調査環境は Ubuntu 20.04 on WSL2）。

一方で `numpy.ufunc` の method だということもわかる。

```sh
Pdb) p umr_sum
<built-in method reduce of numpy.ufunc object at 0x7fe743c56a90>
```


さて、行き詰ってしまったので `numpy.core.umath.add` を探してみると `generate_umath.py(293)` 付近になにやら怪しいものを発見

```sh
$ find -type f -exec grep -nH 'numpy.core.umath.add' {} \;
./numpy/core/code_generators/generate_umath.py:293:          docstrings.get('numpy.core.umath.add'),
```


C言語のコード生成で作っている

```python

# This dictionary describes all the ufunc implementations, generating
# all the function names and their corresponding ufunc signatures.  TD is
# an object which expands a list of character codes into an array of
# TypeDescriptions.
defdict = {
'add':
    Ufunc(2, 1, Zero,
          docstrings.get('numpy.core.umath.add'),
          'PyUFunc_AdditionTypeResolver',
          TD(notimes_or_obj, simd=[('avx2', ints)], dispatch=[('loops_arithm_fp', 'fdFD')]),
          [TypeDescription('M', FullTypeDescr, 'Mm', 'M'),
           TypeDescription('m', FullTypeDescr, 'mm', 'm'),
           TypeDescription('M', FullTypeDescr, 'mM', 'M'),
          ],
          TD(O, f='PyNumber_Add'),
          ),
# 略
# L1206
if __name__ == "__main__":
    filename = __file__
    code = make_code(defdict, filename)
    with open('__umath_generated.c', 'w') as fid:
        fid.write(code)
```

`make_code()` でc言語のコードを作成する DSL 的パートであろう。注目するのは `PyNumber_Add` と `dispatch=[('loops_arithm_fp', 'fdFD')]` の `'loops_arithm_fp` だ。

`loops_arithm_fp` をファイル名で調べると [numpy/core/src/umath/loops_utils.h.src](https://github.com/numpy/numpy/blob/v1.21.4/numpy/core/src/umath/loops_utils.h.src) というファイルが見つかる。そ、そしてついに [L78](https://github.com/numpy/numpy/blob/v1.21.4/numpy/core/src/umath/loops_utils.h.src#L78) に `@TYPE@_pairwise_sum(char *a, npy_intp n, npy_intp stride)` という目的っぽいC言語の関数を見つける！！！（ここまで2.5時間）

---

<details>
<summary>以下ボツとなった調査プロセス（供養）。探索的調査方法が誰かの参考になるため残しておく。</summary>
<div>
さて、ここで行き詰まるのも勿体ないので `numpy/core/include/numpy/*.h` あたりを読むと何やら `PyArray_Sum()` なるものを `numpy/core/include/numpy/__multiarray_api.h` にて発見できた。とはいえローカルファイル上ではソースが存在しないため、ここから github のソースコードを探すとにする

[numpy/core/src/multiarray/calculation.c#L521](https://github.com/numpy/numpy/blob/v1.21.4/numpy/core/src/multiarray/calculation.c#L521) にて `PyArray_Sum()` が定義され、その中で `PyArray_GenericReduceFunction()` を呼び出している。それは [numpy/core/src/multiarray/number.c#L228](https://github.com/numpy/numpy/blob/v1.21.4/numpy/core/src/multiarray/number.c#L228) で定義され、その中で

<pre>
 meth = PyObject_GetAttrString(op, "reduce");
    if (meth && PyCallable_Check(meth)) {
        ret = PyObject_Call(meth, args, kwds);
    }
</pre>

とある。ここで `op = numpy.add` と推測されることから、結局は `numpy.add.reduce` を呼び出しているのであろうか。コードリーディング力が低すぎて循環してしまった。。。（分かる方コメントを頂けると助かります）。

ということで stackoverflow に助けを借りにいくと、[今回の調べたいことそのもの](https://stackoverflow.com/questions/30998305/weird-numpy-sum-behavior-when-adding-zeros) の回答を見つけた。ずばり [numpy/core/src/umath/loops_utils.h.src#L78](https://github.com/numpy/numpy/blob/v1.21.4/numpy/core/src/umath/loops_utils.h.src#L78) にて `@TYPE@_pairwise_sum(...)` で定義されている。
</div>
</details>

---


この `@TYPE@_pairwise_sum(...)` を読むと、次の特徴がわかる

1. 配列要素数が8未満の場合、単純な足し算ループ
1. 配列要素数が `PW_BLOCKSIZE` (=128) 以下の場合、　`r[0]` から `r[7]` までの8つの配列を用意し、r[0] には i = 0, 8, 16, ..., r[1] には i = 1, 9, 17, ..., [...], r[7] には i =  7, 15, 23, ... が詰め込まれ、最後に `res = ((r[0] + r[1]) + (r[2] + r[3])) + ((r[4] + r[5]) + (r[6] + r[7]))` の順番で足し合わされる。
1. 配列要素が 128 超の場合、　前半と後半の約半分に（前半部分が8の倍数になるように）分け、再帰的に足し合わせる。


### Groupby.Sum

pandas/core/groupby/groupby.py(1838)sum()
pandas/core/groupby/groupby.py(1352)_agg_general()
pandas/core/groupby/generic.py(1056)_cython_agg_general()
pandas/core/internals/managers.py(530)get_numeric_data()
pandas/core/internals/managers.py(539)_combine()
pandas/core/groupby/generic.py(1066)_cython_agg_general()
pandas/core/internals/managers.py(1213)grouped_reduce()
pandas/core/internals/blocks.py(374)apply()
pandas/core/groupby/generic.py(1066)array_func()
pandas/core/groupby/ops.py(991)_cython_operation()
pandas/core/groupby/ops.py(631)cython_operation()
pandas/core/groupby/ops.py(488)_cython_op_ndim_compat()
pandas/core/groupby/ops.py(525)_call_cython_op()
(580) の func の call で結果が返っている

```
(Pdb) p func
<cyfunction group_add at 0x7f2f9a657ac0>
```

```sh
$ find -type f -exec grep -nH 'group_add' {} \;
# 略
/pandas/_libs/groupby.pyx:482:def group_add
```

[/pandas/_libs/groupby.pyx を確認](https://github.com/pandas-dev/pandas/blob/v1.3.5/pandas/_libs/groupby.pyx#L250) すると

> For floats, use Kahan summation to reduce floating-point

とあるので、 kahan summation が使われていることがわかる。実装はこの部分を読めば良い（解説略）

### Rolling.sum()

pandas/core/window/rolling.py(1670)sum()
pandas/core/window/rolling.py(1221)sum()
pandas/core/window/rolling.py(482)_apply()
pandas/core/window/rolling.py(411)_apply_blockwise()
pandas/core/window/rolling.py(392)_apply_series()
pandas/core/window/rolling.py(514)homogeneous_func()
pandas/core/window/rolling.py(520)calc()
pandas/core/window/rolling.py(531)calc()

```
(Pdb) p func
<built-in function roll_sum>
```

```sh
$ find -type f -exec grep -nH 'roll_sum' {} \;
# 略
pandas/_libs/window/aggregations.pyx:118:def roll_sum
```
[pandas/_libs/window/aggregations.pyx を確認](https://github.com/pandas-dev/pandas/blob/v1.3.5/pandas/_libs/window/aggregations.pyx#L86) すると `""" remove a value from the sum calc using Kahan summation """` と記述があることが確認できた。



## math.fsum

[CPython 上の math.fsum](https://github.com/python/cpython/blob/3363e1cb05d0d19ed172ea63606d8cb6268747fc/Modules/mathmodule.c#L1446) を見てみると（少し拡張された） kahan summation であることがわかる。

```c
// L1474-1485
for (i = j = 0; j < n; j++) {       /* for y in partials */
    y = p[j];
    if (fabs(x) < fabs(y)) {
        t = x; x = y; y = t;
    }
    hi = x + y;
    yr = hi - x;
    lo = y - yr;
    if (lo != 0.0)
        p[i++] = lo;
    x = hi;
}
```

## sum 関数

[bltinmodule.c](https://github.com/python/cpython/blob/135cabd328504e1648d17242b42b675cdbd0193b/Python/bltinmodule.c#L2407) の `builtin_sum_impl(...)` の中の [L2564](https://github.com/python/cpython/blob/135cabd328504e1648d17242b42b675cdbd0193b/Python/bltinmodule.c#L2564-L2595) あたりの for 文。単純なループでの順番での足し上げ。

```C
    for(;;) {
        item = PyIter_Next(iter);
        if (item == NULL) {
            /* error, or end-of-sequence */
            if (PyErr_Occurred()) {
                Py_DECREF(result);
                result = NULL;
            }
            break;
        }
        temp = PyNumber_Add(result, item);
        Py_DECREF(result);
        Py_DECREF(item);
        result = temp;
        if (result == NULL)
            break;
    }
    Py_DECREF(iter);
    return result;
```

`PyNumber_Add` は、[内部](https://github.com/python/cpython/blob/db671b010cb4662011c6e6567a1f4fcd509102b8/Objects/abstract.c#L1071)で `BINARY_OP1(v, w, NB_SLOT(nb_add), "+")` を call する。以降の解説は [Python処理系入門 〜1 + 1 で学ぶ処理系解読の基礎〜](https://yigarashi.hatenablog.com/entry/cpython#:~:text=3.3.%20PyNumber_Add%20%E3%82%92%E8%BF%BD%E3%81%86)を参考にすると良い。

## 再現

### pairwise summation by numpy 

python, C言語でも結果が再現できない。原因がわかる方、コメントもらえると大変うれしいです。

<details>
<summary>python code (再現できず)</summary><div>

```python
def my_pairwise(lst):
    PW_BLOCKSIZE = 128
    n = len(lst)

    if n < 8:
        out = 0.0
        for i in range(n):
            out += lst[i]
        return out
    elif n <= PW_BLOCKSIZE:
        r = [0.0] * 8
        r[0] = lst[0]
        r[1] = lst[1]
        r[2] = lst[2]
        r[3] = lst[3]
        r[4] = lst[4]
        r[5] = lst[5]
        r[6] = lst[6]
        r[7] = lst[7]
        for i in range(8, n - (n % 8), 8):
            r[0] += lst[0 + i]
            r[1] += lst[1 + i]
            r[2] += lst[2 + i]
            r[3] += lst[3 + i]
            r[4] += lst[4 + i]
            r[5] += lst[5 + i]
            r[6] += lst[6 + i]
            r[7] += lst[7 + i]
        res = ((r[0] + r[1]) + (r[2] + r[3])) +\
              ((r[4] + r[5]) + (r[6] + r[7]))
        
        # rest of 8-blocks
        for i in  range(n - (n % 8), n, 1):
            res += lst[i]
        
        return res
    else:
        n2 = n // 2
        n2 -= n2 % 8
        return my_pairwise(lst[:n2]) + my_pairwise(lst[n2:])

```

```python
out = my_pairwise([1e-7] * int(1e7))
"%0.20f" % out
# '0.99999999999999988898'
# 再現できない
```

</div></details>

<details>
<summary>c言語 (こちらも再現できず)</summary><div>

https://paiza.io/projects/9ELMbLsizKFydjQTTMgaug?language=c

```C
#include <stdio.h>

#define PW_BLOCKSIZE 128
#define IN_SIZE 10000000

double my_pairwise_sum(double *a, int n, int stride);


int main(void){
    double in_array[IN_SIZE];
    for (int i = 0; i < IN_SIZE; i++) {
        in_array[i] = 1.0 / IN_SIZE;
    }
    double out = my_pairwise_sum(in_array, IN_SIZE, 0);
    printf("%.20f\n", out);
}

double my_pairwise_sum(double *a, int n, int stride)
{
    if (n < 8) {
        int i;
        double res = 0.;
        for (i = 0; i < n; i++) {
            res += *(a + i * stride);
        }
        return res;
    }
    else if (n <= PW_BLOCKSIZE) {
        int i;
        double r[8], res;

        /*
         * sum a block with 8 accumulators
         * 8 times unroll reduces blocksize to 16 and allows vectorization with
         * avx without changing summation ordering
         */
        r[0] = *(a + 0 * stride);
        r[1] = *(a + 1 * stride);
        r[2] = *(a + 2 * stride);
        r[3] = *(a + 3 * stride);
        r[4] = *(a + 4 * stride);
        r[5] = *(a + 5 * stride);
        r[6] = *(a + 6 * stride);
        r[7] = *(a + 7 * stride);

        for (i = 8; i < n - (n % 8); i += 8) {
            /* small blocksizes seems to mess with hardware prefetch */
            r[0] += *(a + (i + 0) * stride);
            r[1] += *(a + (i + 1) * stride);
            r[2] += *(a + (i + 2) * stride);
            r[3] += *(a + (i + 3) * stride);
            r[4] += *(a + (i + 4) * stride);
            r[5] += *(a + (i + 5) * stride);
            r[6] += *(a + (i + 6) * stride);
            r[7] += *(a + (i + 7) * stride);
        }

        /* accumulate now to avoid stack spills for single peel loop */
        res = ((r[0] + r[1]) + (r[2] + r[3])) +
              ((r[4] + r[5]) + (r[6] + r[7]));

        /* do non multiple of 8 rest */
        for (; i < n; i++) {
            res += *(a + i * stride);
        }
        return res;
    }
    else {
        /* divide by two but avoid non-multiples of unroll factor */
        int n2 = n / 2;

        n2 -= n2 % 8;
        return my_pairwise_sum(a, n2, stride) +
               my_pairwise_sum(a + n2 * stride, n - n2, stride);
    }
}
// 0.99999999999999988898 
```

</div></details>

### fsum

<details>
<summary>python code</summary><div>

```python
def my_fsum(lst):
    res = 0.0 
    c = 0.0
 
    for l in lst:
        y = l - c
        t = res + y
        c = (t - res) - y
        res = t
 
    return res        
```

```python
1.0
# 一致
```

</div></details>


### basic sum

<details>
<summary>python code</summary><div>

```python
def my_sum(lst):
    out = 0.0
    for l in lst:
        out += l
    return out
```

```python
l =  [1/1e7] * int(1e7)
my_sum(l)
# 0.99999999975017
```

</div></details>


## その他おまけ

```python
# https://github.com/python/cpython/blob/c602c1be439e295fed9ebab47e895ef1d9df28be/Lib/statistics.py#L153
>>> import statistics
>>> statistics._sum([1e50, 1, -1e50] * 1000)
(<class 'float'>, Fraction(1000, 1), 3000)
```

```python
import numpy as np
from math import fsum
from statistics import _sum

for i in range(1, 21):
    lst = [0.1] * i
    print(f"{sum(lst)=:.18f}, {np.sum(lst)=:.18f}, {fsum(lst)=:.18f}, {_sum(lst)[1]=}")
```

```python
sum(lst)=0.100000000000000006, np.sum(lst)=0.100000000000000006, fsum(lst)=0.100000000000000006, _sum(lst)[1]=Fraction(3602879701896397, 36028797018963968)
sum(lst)=0.200000000000000011, np.sum(lst)=0.200000000000000011, fsum(lst)=0.200000000000000011, _sum(lst)[1]=Fraction(3602879701896397, 18014398509481984)
sum(lst)=0.300000000000000044, np.sum(lst)=0.300000000000000044, fsum(lst)=0.300000000000000044, _sum(lst)[1]=Fraction(10808639105689191, 36028797018963968)
sum(lst)=0.400000000000000022, np.sum(lst)=0.400000000000000022, fsum(lst)=0.400000000000000022, _sum(lst)[1]=Fraction(3602879701896397, 9007199254740992)
sum(lst)=0.500000000000000000, np.sum(lst)=0.500000000000000000, fsum(lst)=0.500000000000000000, _sum(lst)[1]=Fraction(18014398509481985, 36028797018963968)
sum(lst)=0.599999999999999978, np.sum(lst)=0.599999999999999978, fsum(lst)=0.600000000000000089, _sum(lst)[1]=Fraction(10808639105689191, 18014398509481984)
sum(lst)=0.699999999999999956, np.sum(lst)=0.699999999999999956, fsum(lst)=0.700000000000000067, _sum(lst)[1]=Fraction(25220157913274779, 36028797018963968)
sum(lst)=0.799999999999999933, np.sum(lst)=0.800000000000000044, fsum(lst)=0.800000000000000044, _sum(lst)[1]=Fraction(3602879701896397, 4503599627370496)
sum(lst)=0.899999999999999911, np.sum(lst)=0.900000000000000022, fsum(lst)=0.900000000000000022, _sum(lst)[1]=Fraction(32425917317067573, 36028797018963968)
sum(lst)=0.999999999999999889, np.sum(lst)=1.000000000000000000, fsum(lst)=1.000000000000000000, _sum(lst)[1]=Fraction(18014398509481985, 18014398509481984)
sum(lst)=1.099999999999999867, np.sum(lst)=1.100000000000000089, fsum(lst)=1.100000000000000089, _sum(lst)[1]=Fraction(39631676720860367, 36028797018963968)
sum(lst)=1.199999999999999956, np.sum(lst)=1.200000000000000178, fsum(lst)=1.200000000000000178, _sum(lst)[1]=Fraction(10808639105689191, 9007199254740992)
sum(lst)=1.300000000000000044, np.sum(lst)=1.300000000000000266, fsum(lst)=1.300000000000000044, _sum(lst)[1]=Fraction(46837436124653161, 36028797018963968)
sum(lst)=1.400000000000000133, np.sum(lst)=1.400000000000000355, fsum(lst)=1.400000000000000133, _sum(lst)[1]=Fraction(25220157913274779, 18014398509481984)
sum(lst)=1.500000000000000222, np.sum(lst)=1.500000000000000444, fsum(lst)=1.500000000000000000, _sum(lst)[1]=Fraction(54043195528445955, 36028797018963968)
sum(lst)=1.600000000000000311, np.sum(lst)=1.600000000000000089, fsum(lst)=1.600000000000000089, _sum(lst)[1]=Fraction(3602879701896397, 2251799813685248)
sum(lst)=1.700000000000000400, np.sum(lst)=1.700000000000000178, fsum(lst)=1.700000000000000178, _sum(lst)[1]=Fraction(61248954932238749, 36028797018963968)
sum(lst)=1.800000000000000488, np.sum(lst)=1.800000000000000266, fsum(lst)=1.800000000000000044, _sum(lst)[1]=Fraction(32425917317067573, 18014398509481984)
sum(lst)=1.900000000000000577, np.sum(lst)=1.900000000000000355, fsum(lst)=1.900000000000000133, _sum(lst)[1]=Fraction(68454714336031543, 36028797018963968)
sum(lst)=2.000000000000000444, np.sum(lst)=2.000000000000000444, fsum(lst)=2.000000000000000000, _sum(lst)[1]=Fraction(18014398509481985, 9007199254740992)
```

## 公式ドキュメントでの言及

[15. 浮動小数点演算、その問題と制限](https://docs.python.org/ja/3/tutorial/floatingpoint.html)
を読むと良い。

```python
x = 0.1
x.as_integer_ratio()
# (3602879701896397, 36028797018963968)

x.hex()
# '0x1.999999999999ap-4'
```



## 最後に

後半が疲れからかグダグダになった感は否めない。


## 参考

https://qiita.com/jo7ueb/items/25ede2bd48c4a2b322e2
https://qiita.com/yucatio/items/c03def631bc5eaaa9887
https://github.com/numpy/numpy/pull/10635/
https://yigarashi.hatenablog.com/entry/cpython
https://note.nkmk.me/python-float-hex/