Python中的长整型(Long)乘法C源码分析

最近学Python看的书是《Learning Python》第三版，在Chapter 2里有一个示例

print 'hello world'
print 2 ** 100

然后说了句“我将在这本书的后面解释print语句，以及为什么在Python中计算2的100次方不会溢出”。

I’ll explain the print statement, and why you can raise 2 to the power 100 in Python without overflowing, in later parts of this book.

Python中的长整型（long）和C语言的long有很大的区别（C语言的long对应Python里的plain integer），Python中的long可以实现无限精度（unlimited precision），很好奇这个在C代码中是怎么实现的，于是看了一下Python的C源码。

虽然求幂运算也有对应的算法，但是最终还是依赖于乘法来实现，所以在这里只研究一下Python的长整型乘法。长整型乘法在Python源码中的Objects目录的longobject.c中定义。

static PyObject *
long_mul(PyLongObject *v, PyLongObject *w)
{
    PyLongObject *a, *b, *z;

    if (!convert_binop((PyObject *)v, (PyObject *)w, &a, &b)) {
        Py_INCREF(Py_NotImplemented);
        return Py_NotImplemented;
    }

    z = k_mul(a, b);
    /* Negate if exactly one of the inputs is negative. */
    if (((a->ob_size ^ b->ob_size) < 0) && z)
        z->ob_size = -(z->ob_size);
    Py_DECREF(a);
    Py_DECREF(b);
    return (PyObject *)z;
}

其中调用了k_mul函数来进行运算，这个k_mul是什么呢？根据注释，这个函数是Karatsuba multiplication（一种比较高效的乘法算法）。

static PyLongObject *
k_mul(PyLongObject *a, PyLongObject *b)
{
    Py_ssize_t asize = ABS(Py_SIZE(a));
    Py_ssize_t bsize = ABS(Py_SIZE(b));
    /* 省略一些变量声明 */
    PyLongObject *t1, *t2, *t3;
    Py_ssize_t i;

    /* fiddle so that b is largest. */
    if (asize > bsize) {
        t1 = a;
        a = b;
        b = t1;
        i = asize;
        asize = bsize;
        bsize = i;
    }

    /* Use gradeschool math when either number is too small. */
    i = a == b ? KARATSUBA_SQUARE_CUTOFF : KARATSUBA_CUTOFF;
    if (asize <= i) {
        if (asize == 0)
            return _PyLong_New(0);
        else
            return x_mul(a, b);
    }
    /* 后面代码省略 */

可以看出，当其中有一个数比较小的时候（这里的小是相对的，整数的位数在70位以下，要知道2的100次方是31位），直接返回是x_mul函数的返回值。那么x_mul又是什么呢？根据注释，这个函数是Grade school multiplication（直接翻译难道是小学乘法？）。

static PyLongObject *
x_mul(PyLongObject *a, PyLongObject *b)
{
    PyLongObject *z;
    Py_ssize_t size_a = ABS(Py_SIZE(a));
    Py_ssize_t size_b = ABS(Py_SIZE(b));
    Py_ssize_t i;

    z = _PyLong_New(size_a + size_b);
    if (z == NULL)
        return NULL;

    memset(z->ob_digit, 0, Py_SIZE(z) * sizeof(digit));
    if (a == b) {
        /* 当a == b即计算a的平方时有更高效的算法 这里省略 */
    }
    else {      /* a is not the same as b -- gradeschool long mult */
        for (i = 0; i < size_a; ++i) {
            twodigits carry = 0;
            twodigits f = a->ob_digit[i];
            digit *pz = z->ob_digit + i;
            digit *pb = b->ob_digit;
            digit *pbend = b->ob_digit + size_b;

            SIGCHECK({
                    Py_DECREF(z);
                    return NULL;
                });

            while (pb < pbend) {
                carry += *pz + *pb++ * f;
                *pz++ = (digit)(carry & PyLong_MASK);
                carry >>= PyLong_SHIFT;
                assert(carry <= PyLong_MASK);
            }
            if (carry)
                *pz += (digit)(carry & PyLong_MASK);
            assert((carry >> PyLong_SHIFT) == 0);
        }
    }
    return long_normalize(z);
}

Karatsuba multiplication和Grade school multiplication的算法描述可以Google到，这里就不罗嗦了（其实是我看不懂）。

结论：C语言这种比较底层的高级语言不是一般人能够凌驾的，还是乖乖的用Python或者Java等做了大量封装的语言吧。