2024 iThome 鐵人賽

DAY 9

Python

為你自己讀 CPython 原始碼系列第 9 篇

Day 9 - 浮點數之小數點漂移記

16th鐵人賽 python python3 原始碼為你自己學

高見龍

2024-09-23 18:45:19

287 瀏覽

分享至

本文同步刊載於「為你自己學 Python - 浮點數之小數點漂移記」

浮點數之小數點漂移記

為你自己學 Python

什麼是「浮點數」？

首先，你可曾想過這個問題，以前我們在數學課都叫它「小數」，為什麼到了程式這邊卻叫做「浮點數（Floating Point Number）」？是在浮什麼？把它丟在水裡會浮起來嗎？

簡單來說，電腦無法直接處理我們日常看到的十進位小數點（像是 3.14 或 0.001），它需要用一種比較特別的方式來表示這些數字，浮點數的結構類似於科學記號表示法，例如 3.14 可以寫成 3.14 x 10^0，0.001 可寫成 1 x 10^-3。在電腦裡是用「有效數字（Mantissa）」和「指數（Exponent）」來表示。

想像一下，如果我們用固定小數點來表示 0.0000123 或是 123,000,000,000 之類的數字，當數字越大或越小，就會需要更多的 0，這可能會有點浪費空間。如果使用類似科學記號表示法的寫法就可以變成 1.23 x 10^-5 和 1.23 x 10^11，不需要一堆 0。

所以，浮點數的這個「浮」，是指小數點的位置可以自由變動，不用固定在哪一個位置上。這樣的設計讓電腦可以有效率地表示非常大或非常小的小數，不會因為數字太大或太小而造成溢位（Overflow）的問題，但缺點就是可能會有一些精準度上的問題。

浮點數的結構

在 Python 裡的浮點數也是個物件，它的結構如下：

// 檔案：Include/cpython/floatobject.h

typedef struct {
    PyObject_HEAD
    double ob_fval;
} PyFloatObject;

跟上個章節講到的整數 PyLongObject 比起來，這個 PyFloatObject 的結構更簡單一點，它只有一個 double 型態的成員變數 ob_fval 來存浮點數的值。有些程式語言，例如 C 語言本身就是，有分 32 位元的單精度浮點數（float）以及 64 位元的雙精度浮點數（double），但 CPython 的原始碼來看，它的浮點數就直接用借 C 語言的 double 來實作浮點數。

當執行 a = 3.14 時，CPython 會建立一個新的 PyFloatObject 物件，這會由 Objects/floatobject.c 裡所定義的 PyFloat_FromDouble() 函數完成，我把一些目前對我們來說比較不重要的條件編譯以及 { } 區塊拿掉，看起來會像這樣：

// 檔案：Objects/floatobject.c

PyObject *
PyFloat_FromDouble(double fval)
{
    PyFloatObject *op;

    op = PyObject_Malloc(sizeof(PyFloatObject));
    if (!op) {
        return PyErr_NoMemory();
    }

    _PyObject_Init((PyObject*)op, &PyFloat_Type);
    op->ob_fval = fval;
    return (PyObject *) op;
}

過程滿單純的，看起來是先呼叫 PyObject_Malloc() 函數要一塊記憶體，再 _PyObject_Init() 初始化，最後把傳入的 fval 指定給成員變數 ob_fval。
嗯...其實也沒這麼單純，這是因為我先省略了一些程式碼所以看起來比較簡單，待會我們看到浮點數的效能的部分再來補充。

正因為 CPython 的浮點數是用 C 語言裡的 double 實作的，所以有些 Python 的浮點數的設計以及運算結果都會跟 C 語言一樣。

關於浮點數

CPython 的浮點數等於 C 語言的 double，所以 CPython 的浮點數也跟 C 語言一樣都是按照 IEEE 754 標準來實作的。IEEE 754 對於雙精度浮點數是使用 64 個位元來表示一個數：

63  62         52                                      0
+---+----------+---------------------------------------+
| S |     E    |                 M                     |
+---+----------+---------------------------------------+

這裡我用幾個簡單的字母來表示：

S 是符號位元（Sign），佔 1 位元，0 代表正數，1 代表負數。
E 是指數位元（Exponent），佔 11 位元。
M 是尾數位元（Mantissa），佔 52 位元。

舉例來說，3.14 的二進位表示法算起來會是：

11.0010001111010111000010100011110101...

雖然這是一個除不盡的無限循環小數，但我們還是得用有限的資源或方式來表示它，所以會有不精準是很正常的。改用類似科學記號表示法的方式來表示，小數點往左移一位，變成：

1.10010001111010111000010100011110101... x 2^1

最後再補充一下 E 也就是指數位元的計算方式，這稍微複雜一點點。IEEE 754 是採用「指數偏移值」的方式計算，以 double 來說它的偏移值（Bias）是 2^10 - 1，也就是 1,023。現在我們知道 3.14 換成二進位會表示成 1.100100011.. x 2^1，這裡 2^1 的 1 就是指數位元 E，再加上偏移值 1,023 變成 1,024，換成二進位為 10000000000。使用偏移值來表示指數的原因是為了讓指數能夠用更簡單的方式來表達正的指數以及負的指數。

接著用 IEEE 754 的雙精度浮點數來表示的話，會是：

S 是 0，因為 3.14 是正數。
E 剛才算出來了，是 10000000000。
M 是小數點左邊的那一串 10010001111010111000010100011110101...

所以 3.14 轉換成 IEEE 754 的表示法會是：

0 10000000000 1001000111101011100001010001111010111000010100011110

從這裡也可以看的出來 M 部份根本沒辦法裝下所有的位數，所以大家常會說浮點數不精準，原因就是這樣。C 語言的 double 是照著 IEEE 754 標準來實作的，從原始碼就能看的出來 CPython 的浮點數其實就是 C 語言的 double，所以也有不精準的問題。

浮點數的運算

回來看一下 PyFloat_Type 結構的定義：

// 檔案：Objects/floatobject.c

PyTypeObject PyFloat_Type = {
    PyVarObject_HEAD_INIT(&PyType_Type, 0)
    "float",
    // ... 略 ...
    &float_as_number,                           /* tp_as_number */
    0,                                          /* tp_as_sequence */
    0,                                          /* tp_clear */
    0,                                          /* tp_as_mapping */
    // ... 略 ...
};

之前有介紹過 tp_as_xxxxx 這幾個成員變數的用途，如果要用來做四則運算的時候，看的就是 tp_as_number 的設定，順著 float_as_number 追過看看：

// 檔案：Objects/floatobject.c

static PyNumberMethods float_as_number = {
    float_add,          /* nb_add */
    float_sub,          /* nb_subtract */
    float_mul,          /* nb_multiply */
    float_rem,          /* nb_remainder */
    float_divmod,       /* nb_divmod */
    float_pow,          /* nb_power */
    (unaryfunc)float_neg, /* nb_negative */
    float_float,        /* nb_positive */
    (unaryfunc)float_abs, /* nb_absolute */
    (inquiry)float_bool, /* nb_bool */
    0,                  /* nb_invert */
    // ... 略 ...
};

光從名字就能猜到這些方法對應到就是四則運算之類的函數了，我們來追個加法 float_add 跟減法 float_add 看看：

// 檔案：Objects/floatobject.c

static PyObject *
float_add(PyObject *v, PyObject *w)
{
    double a,b;
    CONVERT_TO_DOUBLE(v, a);
    CONVERT_TO_DOUBLE(w, b);
    a = a + b;
    return PyFloat_FromDouble(a);
}

static PyObject *
float_sub(PyObject *v, PyObject *w)
{
    double a,b;
    CONVERT_TO_DOUBLE(v, a);
    CONVERT_TO_DOUBLE(w, b);
    a = a - b;
    return PyFloat_FromDouble(a);
}

滿簡單的，其實就是都先轉換成 double，算完再轉回 PyFloatObject 物件而已。

無限大！

在「為你自己學 Python」的數字與文章章節曾經提過「無限大」的概念，不管是正無限大還是負無限大，在 Python 其實都是一種浮點數。當浮點數大於它能表示的最大值的時候並不會出現溢位的問題，而是會變成無限大，反之如果小於能表示的最小值，會變成負的無限大。

不過現在如果知道 Python 的浮點數本質上就是 C 語言的 double 的話，Python 的無限大其實就是 C 語言裡的無限大。我們來寫個簡單的 C 語言程式證明看看：

#include <stdio.h>

int main() {
    // 挑一個接近 double 上限的數字
    double positive_float = 1e308;   // 正的最大值
    double negative_float = -1e308;  // 負的最大值

    // 故意再乘以 10 讓它超過範圍
    double result1 = positive_float * 10;
    double result2 = negative_float * 10;

    printf("結果為: %f\n", result1);
    printf("結果為: %f\n", result2);
}

結果分別會印出 inf 跟 -inf，這就是 C 語言的 double 在超過它能表示的範圍時的行為，所以在 Python 出現一樣的結果也不會太意外了。

不是數字！

NaN（Not a Number）是 IEEE 754 定義的一種特殊的浮點數，它代表一個不是數字的值。CPython 有定義了一個巨集用來判斷是否為 NaN 的巨集是：

#define Py_IS_NAN(X) isnan(X)

如果再往下追，會發現 isnan() 就是 C 語言本身的實作了，也就是這個 Py_IS_NAN 巨集本身就是直接呼叫 C 語言的 isnan 巨集進行判斷而已。

雖然說再追下去就已經不屬於 CPython 而是 C 語言的範圍了，但我還是很好奇想知道 C 語言的 isnan 巨集是怎麼判斷的：

#define isnan(x)                                                         \
    ( sizeof(x) == sizeof(float)  ? __inline_isnanf((float)(x))          \
    : sizeof(x) == sizeof(double) ? __inline_isnand((double)(x))         \
                                  : __inline_isnanl((long double)(x)))

看的出來會根據是 float 還是 double 來判斷要呼叫哪個函數，但如果再看這幾個函數就會發現這幾個函數的實作都差不多，我就挑 __inline_isnand() 看一下：

__header_always_inline int __inline_isnand(double __x) {
    return __x != __x;
}

就是回傳跟自身進行比較的結果，這是 IEEE 754 的定義，NaN 不等於任何值，包括它自己本身，所以只要自己不等於自己，就是 NaN。

浮點數的比較

跟一般的加減乘除相比，浮點數的比較就複雜得多了，連在程式碼裡都直接寫說浮點數的比較根本是惡夢！

// 檔案：Objects/floatobject.c

/* Comparison is pretty much a nightmare.  When comparing float to float,
 * we do it as straightforwardly (and long-windedly) as conceivable, so
 * that, e.g., Python x == y delivers the same result as the platform
 * C x == y when x and/or y is a NaN.
 * When mixing float with an integer type, there's no good *uniform* approach.
 * Converting the double to an integer obviously doesn't work, since we
 * may lose info from fractional bits.  Converting the integer to a double
 * also has two failure modes:  (1) an int may trigger overflow (too
 * large to fit in the dynamic range of a C double); (2) even a C long may have
 * more bits than fit in a C double (e.g., on a 64-bit box long may have
 * 63 bits of precision, but a C double probably has only 53), and then
 * we can falsely claim equality when low-order integer bits are lost by
 * coercion to double.  So this part is painful too.
 */

的確，浮點數本身的比較是一回事，可能還會遇到整數來跟浮點數比較的問題，甚至連不是數字的 NaN 或是無限大都會參一腳。我們來稍微看一下這到底是多麻煩，我們之前學過，像這種「比較」的東西，都是寫在 tp_richcompare 成員變數裡：

// 檔案：Objects/floatobject.c

static PyObject*
float_richcompare(PyObject *v, PyObject *w, int op)
{
    // ... 略 ...
}

這個函數大概有 170 行，裡面有一堆的 if 跟 else if，連判斷 NaN 的也有，難怪會說是惡夢了。

浮點數的效能

其實最一開始講到建立浮點數的 PyFloat_FromDouble() 函數，我刻意把某些條件編譯先拿掉了，它完整的程式碼是這樣：

// 檔案：Objects/floatobject.c

PyObject *
PyFloat_FromDouble(double fval)
{
    PyFloatObject *op;
#if PyFloat_MAXFREELIST > 0
    struct _Py_float_state *state = get_float_state();
    op = state->free_list;
    if (op != NULL) {
        state->free_list = (PyFloatObject *) Py_TYPE(op);
        state->numfree--;
        OBJECT_STAT_INC(from_freelist);
    }
    else
#endif

    // ... 略 ...
}

如果去找 PyFloat_MAXFREELIST 的定義會發現預設值是 100，這東西是為了提高效能而設計的。當浮點數不再使用時，它所佔用的記憶體空間不會立即被釋放，而是被放入一個空閒列表（Free List）中，來看看這個 Free List 的結構：

// 檔案：Include/internal/pycore_floatobject.h

struct _Py_float_state {
#if PyFloat_MAXFREELIST > 0
    /* Special free list
       free_list is a singly-linked list of available PyFloatObjects,
       linked via abuse of their ob_type members. */
    int numfree;
    PyFloatObject *free_list;
#endif
};

結構滿簡單的，有一個 numfree 來記錄目前有多少個空閒的空間，還有一個 free_list 來裝可以來裝浮點數物件。從註解也可看到這是一個單向鏈結串列的結構，透過「濫用」它們的 ob_type 成員把它們串在一起，這濫用不是我說的，是它註解裡這樣寫的。整個串列串起來大概的樣子會是：

+----------------+    +----------------+    +----------------+
| ob_fval = 3.14 | -> | ob_fval = 2.71 | -> | ob_fval = 1.41 | -> NULL
+----------------+    +----------------+    +----------------+

最後指向 NULL 代表這個串列的結束。

從 PyFloat_FromDouble() 函數的程式碼來看，Python 在建立浮點數的時候，會先翻一下 free_list 是不是存在，如果這個列表存在，接下來這行就是重點了：

state->free_list = (PyFloatObject *) Py_TYPE(op);

就以語法本身來看，看起來像是把一個新的 PyFloatObject 物件透過 Py_TYPE() 巨集轉換型態後再設定給 free_list，但事實上 Py_TYPE() 這個巨集是用來取得物件的 ob_type：

// 檔案：Include/object.h

static inline PyTypeObject* Py_TYPE(PyObject *ob) {
    return ob->ob_type;
}

這個函數會取得這顆物件身上的 ob_type，所以實際上這並不是在建立新的物件，以結果來說而是在移動 free_list 指針，使它指向鏈結串列的下一個元素，這就是剛剛講的「濫用」的部分。本質上是對這個列表進行 pop 的操作，會把這個鏈結串列最前面的記憶體位置給拿出來，讓指針指向下一個記憶體位置。這樣可以重複使用同樣大小的記憶體區塊，還能減少了記憶體碎片化的問題。

真有趣，我想都沒想過可以這樣做，不得不說，這設計實在十分巧妙！