/* Find highest digit where a and b differ: */ i = size_a; while (--i >= 0 && a->ob_digit[i] == b->ob_digit[i]) ; if (i < 0) return _PyLong_New(0); if (a->ob_digit[i] < b->ob_digit[i]) { sign = -1; { PyLongObject *temp = a; a = b; b = temp; } } size_a = size_b = i+1; } z = _PyLong_New(size_a); if (z == NULL) return NULL; for (i = 0; i < size_b; ++i) { /* The following assumes unsigned arithmetic works module 2**N for some N>SHIFT. */ borrow = a->ob_digit[i] - b->ob_digit[i] - borrow; z->ob_digit[i] = borrow & MASK; borrow >>= SHIFT; borrow &= 1; /* Keep only one sign bit */ } for (; i < size_a; ++i) { borrow = a->ob_digit[i] - borrow; z->ob_digit[i] = borrow & MASK; borrow >>= SHIFT; borrow &= 1; /* Keep only one sign bit */ } assert(borrow == 0); if (sign < 0) z->ob_size = -(z->ob_size); return long_normalize(z);} static PyObject *long_add(PyLongObject *v, PyLongObject *w){ PyLongObject *a, *b, *z; CONVERT_BINOP((PyObject *)v, (PyObject *)w, &a, &b); if (a->ob_size < 0) { if (b->ob_size < 0) { z = x_add(a, b); if (z != NULL && z->ob_size != 0) z->ob_size = -(z->ob_size); } else z = x_sub(b, a); } else { if (b->ob_size < 0) z = x_sub(a, b); else z = x_add(a, b); } Py_DECREF(a); Py_DECREF(b); return (PyObject *)z;} static PyObject *long_sub(PyLongObject *v, PyLongObject *w){ PyLongObject *a, *b, *z; CONVERT_BINOP((PyObject *)v, (PyObject *)w, &a, &b); if (a->ob_size < 0) { if (b->ob_size < 0) z = x_sub(a, b); else z = x_add(a, b); if (z != NULL && z->ob_size != 0) z->ob_size = -(z->ob_size); } else { if (b->ob_size < 0) z = x_add(a, b); else z = x_sub(a, b); } Py_DECREF(a); Py_DECREF(b); return (PyObject *)z;} /* Grade school multiplication, ignoring the signs. * Returns the absolute value of the product, or NULL if error. */static PyLongObject *x_mul(PyLongObject *a, PyLongObject *b){ PyLongObject *z; Py_ssize_t size_a = ABS(a->ob_size); Py_ssize_t size_b = ABS(b->ob_size); Py_ssize_t i; z = _PyLong_New(size_a + size_b); if (z == NULL) return NULL; memset(z->ob_digit, 0, z->ob_size * sizeof(digit)); if (a == b) { /* Efficient squaring per HAC, Algorithm 14.16: * http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf * Gives slightly less than a 2x speedup when a == b, * via exploiting that each entry in the multiplication * pyramid appears twice (except for the size_a squares). */ for (i = 0; i < size_a; ++i) { twodigits carry; twodigits f = a->ob_digit[i]; digit *pz = z->ob_digit + (i << 1); digit *pa = a->ob_digit + i + 1; digit *paend = a->ob_digit + size_a; SIGCHECK({ Py_DECREF(z); return NULL; }) carry = *pz + f * f; *pz++ = (digit)(carry & MASK); carry >>= SHIFT; assert(carry <= MASK); /* Now f is added in twice in each column of the * pyramid it appears. Same as adding f<<1 once. */ f <<= 1; while (pa < paend) { carry += *pz + *pa++ * f; *pz++ = (digit)(carry & MASK); carry >>= SHIFT; assert(carry <= (MASK << 1)); } if (carry) { carry += *pz; *pz++ = (digit)(carry & MASK); carry >>= SHIFT; } if (carry) *pz += (digit)(carry & MASK); assert((carry >> SHIFT) == 0); } } 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 & MASK); carry >>= SHIFT; assert(carry <= MASK); } if (carry) *pz += (digit)(carry & MASK); assert((carry >> SHIFT) == 0); } } return long_normalize(z);} /* A helper for Karatsuba multiplication (k_mul). Takes a long "n" and an integer "size" representing the place to split, and sets low and high such that abs(n) == (high << size) + low, viewing the shift as being by digits. The sign bit is ignored, and the return values are >= 0. Returns 0 on success, -1 on failure.*/static intkmul_split(PyLongObject *n, Py_ssize_t size, PyLongObject **high, PyLongObject **low){ PyLongObject *hi, *lo; Py_ssize_t size_lo, size_hi; const Py_ssize_t size_n = ABS(n->ob_size); size_lo = MIN(size_n, size); size_hi = size_n - size_lo; if ((hi = _PyLong_New(size_hi)) == NULL) return -1; if ((lo = _PyLong_New(size_lo)) == NULL) { Py_DECREF(hi); return -1; } memcpy(lo->ob_digit, n->ob_digit, size_lo * sizeof(digit)); memcpy(hi->ob_digit, n->ob_digit + size_lo, size_hi * sizeof(digit)); *high = long_normalize(hi); *low = long_normalize(lo); return 0;} static PyLongObject *k_lopsided_mul(PyLongObject *a, PyLongObject *b); /* Karatsuba multiplication. Ignores the input signs, and returns the * absolute value of the product (or NULL if error). * See Knuth Vol. 2 Chapter 4.3.3 (Pp. 294-295). */static PyLongObject *k_mul(PyLongObject *a, PyLongObject *b){ Py_ssize_t asize = ABS(a->ob_size); Py_ssize_t bsize = ABS(b->ob_size); PyLongObject *ah = NULL; PyLongObject *al = NULL; PyLongObject *bh = NULL; PyLongObject *bl = NULL; PyLongObject *ret = NULL; PyLongObject *t1, *t2, *t3; Py_ssize_t shift; /* the number of digits we split off */ Py_ssize_t i; /* (ah*X+al)(bh*X+bl) = ah*bh*X*X + (ah*bl + al*bh)*X + al*bl * Let k = (ah+al)*(bh+bl) = ah*bl + al*bh + ah*bh + al*bl * Then the original product is * ah*bh*X*X + (k - ah*bh - al*bl)*X + al*bl * By picking X to be a power of 2, "*X" is just shifting, and it's * been reduced to 3 multiplies on numbers half the size. */ /* We want to split based on the larger number; 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); } /* If a is small compared to b, splitting on b gives a degenerate * case with ah==0, and Karatsuba may be (even much) less efficient * than "grade school" then. However, we can still win, by viewing * b as a string of "big digits", each of width a->ob_size. That * leads to a sequence of balanced calls to k_mul. */ if (2 * asize <= bsize) return k_lopsided_mul(a, b); /* Split a & b into hi & lo pieces. */ shift = bsize >> 1; if (kmul_split(a, shift, &ah, &al) < 0) goto fail; assert(ah->ob_size > 0); /* the split isn't degenerate */ if (a == b) { bh = ah; bl = al; Py_INCREF(bh); Py_INCREF(bl); } else if (kmul_split(b, shift, &bh, &bl) < 0) goto fail; /* The plan: * 1. Allocate result space (asize + bsize digits: that's always * enough). * 2. Compute ah*bh, and copy into result at 2*shift. * 3. Compute al*bl, and copy into result at 0. Note that this * can't overlap with #2. * 4. Subtract al*bl from the result, starting at shift. This may * underflow (borrow out of the high digit), but we don't care: * we're effectively doing unsigned arithmetic mod * BASE**(sizea + sizeb), and so long as the *final* result fits, * borrows and carries out of the high digit can be ignored. * 5. Subtract ah*bh from the result, starting at shift. * 6. Compute (ah+al)*(bh+bl), and add it into the result starting * at shift. */ /* 1. Allocate result space. */ ret = _PyLong_New(asize + bsize); if (ret == NULL) goto fail;#ifdef Py_DEBUG /* Fill with trash, to catch reference to uninitialized digits. */ memset(ret->ob_digit, 0xDF, ret->ob_size * sizeof(digit));#endif /* 2. t1 <- ah*bh, and copy into high digits of result. */ if ((t1 = k_mul(ah, bh)) == NULL) goto fail; assert(t1->ob_size >= 0); assert(2*shift + t1->ob_size <= ret->ob_size); memcpy(ret->ob_digit + 2*shift, t1->ob_digit, t1->ob_size * sizeof(digit)); /* Zero-out the digits higher than the ah*bh copy. */ i = ret->ob_size - 2*shift - t1->ob_size; if (i) memset(ret->ob_digit + 2*shift + t1->ob_size, 0, i * sizeof(digit)); /* 3. t2 <- al*bl, and copy into the low digits. */ if ((t2 = k_mul(al, bl)) == NULL) { Py_DECREF(t1); goto fail; } assert(t2->ob_size >= 0); assert(t2->ob_size <= 2*shift); /* no overlap with high digits */ memcpy(ret->ob_digit, t2->ob_digit, t2->ob_size * sizeof(digit)); /* Zero out remaining digits. */ i = 2*shift - t2->ob_size; /* number of uninitialized digits */ if (i) memset(ret->ob_digit + t2->ob_size, 0, i * sizeof(digit)); /* 4 & 5. Subtract ah*bh (t1) and al*bl (t2). We do al*bl first * because it's fresher in cache. */ i = ret->ob_size - shift; /* # digits after shift */ (void)v_isub(ret->ob_digit + shift, i, t2->ob_digit, t2->ob_size); Py_DECREF(t2); (void)v_isub(ret->ob_digit + shift, i, t1->ob_digit, t1->ob_size); Py_DECREF(t1); /* 6. t3 <- (ah+al)(bh+bl), and add into result. */ if ((t1 = x_add(ah, al)) == NULL) goto fail; Py_DECREF(ah); Py_DECREF(al); ah = al = NULL; if (a == b) { t2 = t1; Py_INCREF(t2); } else if ((t2 = x_add(bh, bl)) == NULL) { Py_DECREF(t1); goto fail; } Py_DECREF(bh); Py_DECREF(bl); bh = bl = NULL; t3 = k_mul(t1, t2); Py_DECREF(t1); Py_DECREF(t2); if (t3 == NULL) goto fail; assert(t3->ob_size >= 0); /* Add t3. It's not obvious why we can't run out of room here. * See the (*) comment after this function. */ (void)v_iadd(ret->ob_digit + shift, i, t3->ob_digit, t3->ob_size); Py_DECREF(t3); return long_normalize(ret); fail: Py_XDECREF(ret); Py_XDECREF(ah); Py_XDECREF(al); Py_XDECREF(bh); Py_XDECREF(bl); return NULL;} /* (*) Why adding t3 can't "run out of room" above. Let f(x) mean the floor of x and c(x) mean the ceiling of x. Some factsto start with: 1. For any integer i, i = c(i/2) + f(i/2). In particular, bsize = c(bsize/2) + f(bsize/2).2. shift = f(bsize/2)3. asize <= bsize4. Since we call k_lopsided_mul if asize*2 <= bsize, asize*2 > bsize in this routine, so asize > bsize/2 >= f(bsize/2) in this routine. We allocated asize + bsize result digits, and add t3 into them at an offsetof shift. This leaves asize+bsize-shift allocated digit positions for t3to fit into, = (by #1 and #2) asize + f(bsize/2) + c(bsize/2) - f(bsize/2) =asize + c(bsize/2) available digit positions. bh has c(bsize/2) digits, and bl at most f(size/2) digits. So bh+hl hasat most c(bsize/2) digits + 1 bit. If asize == bsize, ah has c(bsize/2) digits, else ah has at most f(bsize/2)digits, and al has at most f(bsize/2) digits in any case. So ah+al has atmost (asize == bsize ? c(bsize/2) : f(bsize/2)) digits + 1 bit. The product (ah+al)*(bh+bl) therefore has at most c(bsize/2) + (asize == bsize ? c(bsize/2) : f(bsize/2)) digits + 2 bits and we have asize + c(bsize/2) available digit positions. We need to showthis is always enough. An instance of c(bsize/2) cancels out in both, sothe question reduces to whether asize digits is enough to hold(asize == bsize ? c(bsize/2) : f(bsize/2)) digits + 2 bits. If asize < bsize,then we're asking whether asize digits >= f(bsize/2) digits + 2 bits. By #4,asize is at least f(bsize/2)+1 digits, so this in turn reduces to whether 1digit is enough to hold 2 bits. This is so since SHIFT=15 >= 2. Ifasize == bsize, then we're asking whether bsize digits is enough to holdc(bsize/2) digits + 2 bits, or equivalently (by #1) whether f(bsize/2) digitsis enough to hold 2 bits. This is so if bsize >= 2, which holds becausebsize >= KARATSUBA_CUTOFF >= 2. Note that since there's always enough room for (ah+al)*(bh+bl), and that'sclearly >= each of ah*bh and al*bl, there's always enough room to subtractah*bh and al*bl too.*/ /* b has at least twice the digits of a, and a is big enough that Karatsuba * would pay off *if* the inputs had balanced sizes. View b as a sequence * of slices, each with a->ob_size digits, and multiply the slices by a, * one at a time. This gives k_mul balanced inputs to work with, and is * also cache-friendly (we compute one double-width slice of the result * at a time, then move on, never bactracking except for the helpful * single-width slice overlap between successive partial sums). */static PyLongObject *k_lopsided_mul(PyLongObject *a, PyLongObject *b){ const Py_ssize_t asize = ABS(a->ob_size); Py_ssize_t bsize = ABS(b->ob_size); Py_ssize_t nbdone; /* # of b digits already multiplied */ PyLongObject *ret; PyLongObject *bslice = NULL; assert(asize > KARATSUBA_CUTOFF); assert(2 * asize <= bsize); /* Allocate result space, and zero it out. */ ret = _PyLong_New(asize + bsize); if (ret == NULL) return NULL; memset(ret->ob_digit, 0, ret->ob_size * sizeof(digit)); /* Successive slices of b are copied into bslice. */ bslice = _PyLong_New(asize); if (bslice == NULL) goto fail; nbdone = 0; while (bsize > 0) { PyLongObject *product; const Py_ssize_t nbtouse = MIN(bsize, asize); /* Multiply the next slice of b by a. */ memcpy(bslice->ob_digit, b->ob_digit + nbdone, nbtouse * sizeof(digit)); bslice->ob_size = nbtouse; product = k_mul(a, bslice); if (product == NULL) goto fail; /* Add into result. */ (void)v_iadd(ret->ob_digit + nbdone, ret->ob_size - nbdone, product->ob_digit, product->ob_size); Py_DECREF(product); bsize -= nbtouse; nbdone += nbtouse; } Py_DECREF(bslice); return long_normalize(ret); fail: Py_DECREF(ret); Py_XDECREF(bslice); return NULL;} 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. */Lines 2001–2500 of 3458. Previous lines Load more
/* Find highest digit where a and b differ: */ i = size_a; while (--i >= 0 && a->ob_digit[i] == b->ob_digit[i]) ; if (i < 0) return _PyLong_New(0); if (a->ob_digit[i] < b->ob_digit[i]) { sign = -1; { PyLongObject *temp = a; a = b; b = temp; } } size_a = size_b = i+1; } z = _PyLong_New(size_a); if (z == NULL) return NULL; for (i = 0; i < size_b; ++i) { /* The following assumes unsigned arithmetic works module 2**N for some N>SHIFT. */ borrow = a->ob_digit[i] - b->ob_digit[i] - borrow; z->ob_digit[i] = borrow & MASK; borrow >>= SHIFT; borrow &= 1; /* Keep only one sign bit */ } for (; i < size_a; ++i) { borrow = a->ob_digit[i] - borrow; z->ob_digit[i] = borrow & MASK; borrow >>= SHIFT; borrow &= 1; /* Keep only one sign bit */ } assert(borrow == 0); if (sign < 0) z->ob_size = -(z->ob_size); return long_normalize(z); } static PyObject * long_add(PyLongObject *v, PyLongObject *w) { PyLongObject *a, *b, *z; CONVERT_BINOP((PyObject *)v, (PyObject *)w, &a, &b); if (a->ob_size < 0) { if (b->ob_size < 0) { z = x_add(a, b); if (z != NULL && z->ob_size != 0) z->ob_size = -(z->ob_size); } else z = x_sub(b, a); } else { if (b->ob_size < 0) z = x_sub(a, b); else z = x_add(a, b); } Py_DECREF(a); Py_DECREF(b); return (PyObject *)z; } static PyObject * long_sub(PyLongObject *v, PyLongObject *w) { PyLongObject *a, *b, *z; CONVERT_BINOP((PyObject *)v, (PyObject *)w, &a, &b); if (a->ob_size < 0) { if (b->ob_size < 0) z = x_sub(a, b); else z = x_add(a, b); if (z != NULL && z->ob_size != 0) z->ob_size = -(z->ob_size); } else { if (b->ob_size < 0) z = x_add(a, b); else z = x_sub(a, b); } Py_DECREF(a); Py_DECREF(b); return (PyObject *)z; } /* Grade school multiplication, ignoring the signs. * Returns the absolute value of the product, or NULL if error. */ static PyLongObject * x_mul(PyLongObject *a, PyLongObject *b) { PyLongObject *z; Py_ssize_t size_a = ABS(a->ob_size); Py_ssize_t size_b = ABS(b->ob_size); Py_ssize_t i; z = _PyLong_New(size_a + size_b); if (z == NULL) return NULL; memset(z->ob_digit, 0, z->ob_size * sizeof(digit)); if (a == b) { /* Efficient squaring per HAC, Algorithm 14.16: * http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf * Gives slightly less than a 2x speedup when a == b, * via exploiting that each entry in the multiplication * pyramid appears twice (except for the size_a squares). */ for (i = 0; i < size_a; ++i) { twodigits carry; twodigits f = a->ob_digit[i]; digit *pz = z->ob_digit + (i << 1); digit *pa = a->ob_digit + i + 1; digit *paend = a->ob_digit + size_a; SIGCHECK({ Py_DECREF(z); return NULL; }) carry = *pz + f * f; *pz++ = (digit)(carry & MASK); carry >>= SHIFT; assert(carry <= MASK); /* Now f is added in twice in each column of the * pyramid it appears. Same as adding f<<1 once. */ f <<= 1; while (pa < paend) { carry += *pz + *pa++ * f; *pz++ = (digit)(carry & MASK); carry >>= SHIFT; assert(carry <= (MASK << 1)); } if (carry) { carry += *pz; *pz++ = (digit)(carry & MASK); carry >>= SHIFT; } if (carry) *pz += (digit)(carry & MASK); assert((carry >> SHIFT) == 0); } } 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 & MASK); carry >>= SHIFT; assert(carry <= MASK); } if (carry) *pz += (digit)(carry & MASK); assert((carry >> SHIFT) == 0); } } return long_normalize(z); } /* A helper for Karatsuba multiplication (k_mul). Takes a long "n" and an integer "size" representing the place to split, and sets low and high such that abs(n) == (high << size) + low, viewing the shift as being by digits. The sign bit is ignored, and the return values are >= 0. Returns 0 on success, -1 on failure. */ static int kmul_split(PyLongObject *n, Py_ssize_t size, PyLongObject **high, PyLongObject **low) { PyLongObject *hi, *lo; Py_ssize_t size_lo, size_hi; const Py_ssize_t size_n = ABS(n->ob_size); size_lo = MIN(size_n, size); size_hi = size_n - size_lo; if ((hi = _PyLong_New(size_hi)) == NULL) return -1; if ((lo = _PyLong_New(size_lo)) == NULL) { Py_DECREF(hi); return -1; } memcpy(lo->ob_digit, n->ob_digit, size_lo * sizeof(digit)); memcpy(hi->ob_digit, n->ob_digit + size_lo, size_hi * sizeof(digit)); *high = long_normalize(hi); *low = long_normalize(lo); return 0; } static PyLongObject *k_lopsided_mul(PyLongObject *a, PyLongObject *b); /* Karatsuba multiplication. Ignores the input signs, and returns the * absolute value of the product (or NULL if error). * See Knuth Vol. 2 Chapter 4.3.3 (Pp. 294-295). */ static PyLongObject * k_mul(PyLongObject *a, PyLongObject *b) { Py_ssize_t asize = ABS(a->ob_size); Py_ssize_t bsize = ABS(b->ob_size); PyLongObject *ah = NULL; PyLongObject *al = NULL; PyLongObject *bh = NULL; PyLongObject *bl = NULL; PyLongObject *ret = NULL; PyLongObject *t1, *t2, *t3; Py_ssize_t shift; /* the number of digits we split off */ Py_ssize_t i; /* (ah*X+al)(bh*X+bl) = ah*bh*X*X + (ah*bl + al*bh)*X + al*bl * Let k = (ah+al)*(bh+bl) = ah*bl + al*bh + ah*bh + al*bl * Then the original product is * ah*bh*X*X + (k - ah*bh - al*bl)*X + al*bl * By picking X to be a power of 2, "*X" is just shifting, and it's * been reduced to 3 multiplies on numbers half the size. */ /* We want to split based on the larger number; 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); } /* If a is small compared to b, splitting on b gives a degenerate * case with ah==0, and Karatsuba may be (even much) less efficient * than "grade school" then. However, we can still win, by viewing * b as a string of "big digits", each of width a->ob_size. That * leads to a sequence of balanced calls to k_mul. */ if (2 * asize <= bsize) return k_lopsided_mul(a, b); /* Split a & b into hi & lo pieces. */ shift = bsize >> 1; if (kmul_split(a, shift, &ah, &al) < 0) goto fail; assert(ah->ob_size > 0); /* the split isn't degenerate */ if (a == b) { bh = ah; bl = al; Py_INCREF(bh); Py_INCREF(bl); } else if (kmul_split(b, shift, &bh, &bl) < 0) goto fail; /* The plan: * 1. Allocate result space (asize + bsize digits: that's always * enough). * 2. Compute ah*bh, and copy into result at 2*shift. * 3. Compute al*bl, and copy into result at 0. Note that this * can't overlap with #2. * 4. Subtract al*bl from the result, starting at shift. This may * underflow (borrow out of the high digit), but we don't care: * we're effectively doing unsigned arithmetic mod * BASE**(sizea + sizeb), and so long as the *final* result fits, * borrows and carries out of the high digit can be ignored. * 5. Subtract ah*bh from the result, starting at shift. * 6. Compute (ah+al)*(bh+bl), and add it into the result starting * at shift. */ /* 1. Allocate result space. */ ret = _PyLong_New(asize + bsize); if (ret == NULL) goto fail; #ifdef Py_DEBUG /* Fill with trash, to catch reference to uninitialized digits. */ memset(ret->ob_digit, 0xDF, ret->ob_size * sizeof(digit)); #endif /* 2. t1 <- ah*bh, and copy into high digits of result. */ if ((t1 = k_mul(ah, bh)) == NULL) goto fail; assert(t1->ob_size >= 0); assert(2*shift + t1->ob_size <= ret->ob_size); memcpy(ret->ob_digit + 2*shift, t1->ob_digit, t1->ob_size * sizeof(digit)); /* Zero-out the digits higher than the ah*bh copy. */ i = ret->ob_size - 2*shift - t1->ob_size; if (i) memset(ret->ob_digit + 2*shift + t1->ob_size, 0, i * sizeof(digit)); /* 3. t2 <- al*bl, and copy into the low digits. */ if ((t2 = k_mul(al, bl)) == NULL) { Py_DECREF(t1); goto fail; } assert(t2->ob_size >= 0); assert(t2->ob_size <= 2*shift); /* no overlap with high digits */ memcpy(ret->ob_digit, t2->ob_digit, t2->ob_size * sizeof(digit)); /* Zero out remaining digits. */ i = 2*shift - t2->ob_size; /* number of uninitialized digits */ if (i) memset(ret->ob_digit + t2->ob_size, 0, i * sizeof(digit)); /* 4 & 5. Subtract ah*bh (t1) and al*bl (t2). We do al*bl first * because it's fresher in cache. */ i = ret->ob_size - shift; /* # digits after shift */ (void)v_isub(ret->ob_digit + shift, i, t2->ob_digit, t2->ob_size); Py_DECREF(t2); (void)v_isub(ret->ob_digit + shift, i, t1->ob_digit, t1->ob_size); Py_DECREF(t1); /* 6. t3 <- (ah+al)(bh+bl), and add into result. */ if ((t1 = x_add(ah, al)) == NULL) goto fail; Py_DECREF(ah); Py_DECREF(al); ah = al = NULL; if (a == b) { t2 = t1; Py_INCREF(t2); } else if ((t2 = x_add(bh, bl)) == NULL) { Py_DECREF(t1); goto fail; } Py_DECREF(bh); Py_DECREF(bl); bh = bl = NULL; t3 = k_mul(t1, t2); Py_DECREF(t1); Py_DECREF(t2); if (t3 == NULL) goto fail; assert(t3->ob_size >= 0); /* Add t3. It's not obvious why we can't run out of room here. * See the (*) comment after this function. */ (void)v_iadd(ret->ob_digit + shift, i, t3->ob_digit, t3->ob_size); Py_DECREF(t3); return long_normalize(ret); fail: Py_XDECREF(ret); Py_XDECREF(ah); Py_XDECREF(al); Py_XDECREF(bh); Py_XDECREF(bl); return NULL; } /* (*) Why adding t3 can't "run out of room" above. Let f(x) mean the floor of x and c(x) mean the ceiling of x. Some facts to start with: 1. For any integer i, i = c(i/2) + f(i/2). In particular, bsize = c(bsize/2) + f(bsize/2). 2. shift = f(bsize/2) 3. asize <= bsize 4. Since we call k_lopsided_mul if asize*2 <= bsize, asize*2 > bsize in this routine, so asize > bsize/2 >= f(bsize/2) in this routine. We allocated asize + bsize result digits, and add t3 into them at an offset of shift. This leaves asize+bsize-shift allocated digit positions for t3 to fit into, = (by #1 and #2) asize + f(bsize/2) + c(bsize/2) - f(bsize/2) = asize + c(bsize/2) available digit positions. bh has c(bsize/2) digits, and bl at most f(size/2) digits. So bh+hl has at most c(bsize/2) digits + 1 bit. If asize == bsize, ah has c(bsize/2) digits, else ah has at most f(bsize/2) digits, and al has at most f(bsize/2) digits in any case. So ah+al has at most (asize == bsize ? c(bsize/2) : f(bsize/2)) digits + 1 bit. The product (ah+al)*(bh+bl) therefore has at most c(bsize/2) + (asize == bsize ? c(bsize/2) : f(bsize/2)) digits + 2 bits and we have asize + c(bsize/2) available digit positions. We need to show this is always enough. An instance of c(bsize/2) cancels out in both, so the question reduces to whether asize digits is enough to hold (asize == bsize ? c(bsize/2) : f(bsize/2)) digits + 2 bits. If asize < bsize, then we're asking whether asize digits >= f(bsize/2) digits + 2 bits. By #4, asize is at least f(bsize/2)+1 digits, so this in turn reduces to whether 1 digit is enough to hold 2 bits. This is so since SHIFT=15 >= 2. If asize == bsize, then we're asking whether bsize digits is enough to hold c(bsize/2) digits + 2 bits, or equivalently (by #1) whether f(bsize/2) digits is enough to hold 2 bits. This is so if bsize >= 2, which holds because bsize >= KARATSUBA_CUTOFF >= 2. Note that since there's always enough room for (ah+al)*(bh+bl), and that's clearly >= each of ah*bh and al*bl, there's always enough room to subtract ah*bh and al*bl too. */ /* b has at least twice the digits of a, and a is big enough that Karatsuba * would pay off *if* the inputs had balanced sizes. View b as a sequence * of slices, each with a->ob_size digits, and multiply the slices by a, * one at a time. This gives k_mul balanced inputs to work with, and is * also cache-friendly (we compute one double-width slice of the result * at a time, then move on, never bactracking except for the helpful * single-width slice overlap between successive partial sums). */ static PyLongObject * k_lopsided_mul(PyLongObject *a, PyLongObject *b) { const Py_ssize_t asize = ABS(a->ob_size); Py_ssize_t bsize = ABS(b->ob_size); Py_ssize_t nbdone; /* # of b digits already multiplied */ PyLongObject *ret; PyLongObject *bslice = NULL; assert(asize > KARATSUBA_CUTOFF); assert(2 * asize <= bsize); /* Allocate result space, and zero it out. */ ret = _PyLong_New(asize + bsize); if (ret == NULL) return NULL; memset(ret->ob_digit, 0, ret->ob_size * sizeof(digit)); /* Successive slices of b are copied into bslice. */ bslice = _PyLong_New(asize); if (bslice == NULL) goto fail; nbdone = 0; while (bsize > 0) { PyLongObject *product; const Py_ssize_t nbtouse = MIN(bsize, asize); /* Multiply the next slice of b by a. */ memcpy(bslice->ob_digit, b->ob_digit + nbdone, nbtouse * sizeof(digit)); bslice->ob_size = nbtouse; product = k_mul(a, bslice); if (product == NULL) goto fail; /* Add into result. */ (void)v_iadd(ret->ob_digit + nbdone, ret->ob_size - nbdone, product->ob_digit, product->ob_size); Py_DECREF(product); bsize -= nbtouse; nbdone += nbtouse; } Py_DECREF(bslice); return long_normalize(ret); fail: Py_DECREF(ret); Py_XDECREF(bslice); return NULL; } 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. */