@@ -960,6 +960,175 @@ mpz_sqr_basic_limbs(mp_limb *result, const mp_limb *x, size_t n)
960960 }
961961}
962962963+/*
964+ * Karatsuba Multiplication
965+ *
966+ * Splits inputs into 2 parts: A = A1*B^half + A0, B = B1*B^half + B0
967+ * Computes: z0 = A0*B0, z2 = A1*B1, z1 = (A0+A1)*(B0+B1) - z0 - z2
968+ * Result: z2*B^(2*half) + z1*B^half + z0
969+ *
970+ * Complexity: O(n^1.585) - trades 4 multiplications for 3 plus additions
971+ */
972+973+#define KARATSUBA_THRESHOLD 32
974+975+static inline mrb_bool
976+should_use_karatsuba(size_t n)
977+{
978+return n >= KARATSUBA_THRESHOLD;
979+}
980+981+/* Calculate scratch space needed for Karatsuba */
982+static size_t
983+karatsuba_scratch_size(size_t n)
984+{
985+if (n < KARATSUBA_THRESHOLD) {
986+return 0;
987+ }
988+989+size_t half = (n + 1) / 2;
990+991+/*
992+ * Per level storage:
993+ * - 2 evaluation results: (A0+A1), (B0+B1), each up to (half+1) limbs
994+ * - 3 products: z0, z1, z2, each up to 2*(half+1) limbs
995+ */
996+size_t eval_len = half + 1;
997+size_t prod_len = 2 * eval_len;
998+999+size_t eval_size = 2 * eval_len; /* 2 evaluation temps */
1000+size_t prod_size = 3 * prod_len; /* 3 products */
1001+size_t current_level = eval_size + prod_size;
1002+1003+/* Recursive scratch - sequential calls reuse same buffer */
1004+size_t sub_scratch = karatsuba_scratch_size(eval_len);
1005+1006+return current_level + sub_scratch + 8; /* +8 safety margin */
1007+}
1008+1009+/* Forward declaration for recursive calls */
1010+static void
1011+mpz_mul_karatsuba_limbs(mp_limb *result,
1012+const mp_limb *x, size_t x_len,
1013+const mp_limb *y, size_t y_len,
1014+mp_limb *scratch);
1015+1016+/*
1017+ * Karatsuba multiplication on raw limb arrays.
1018+ *
1019+ * result must have space for x_len + y_len limbs.
1020+ * scratch must have karatsuba_scratch_size(max(x_len, y_len)) limbs.
1021+ */
1022+static void
1023+mpz_mul_karatsuba_limbs(mp_limb *result,
1024+const mp_limb *x, size_t x_len,
1025+const mp_limb *y, size_t y_len,
1026+mp_limb *scratch)
1027+{
1028+size_t min_len = (x_len < y_len) ? x_len : y_len;
1029+size_t max_len = (x_len > y_len) ? x_len : y_len;
1030+1031+/* Base case - use schoolbook */
1032+if (!should_use_karatsuba(min_len)) {
1033+mpz_mul_basic_limbs(result, x, x_len, y, y_len);
1034+return;
1035+ }
1036+1037+/*
1038+ * Split: x = x1*B^half + x0, y = y1*B^half + y0
1039+ * where B = base^half
1040+ */
1041+size_t half = (max_len + 1) / 2;
1042+1043+/* Determine actual lengths of each part */
1044+size_t x0_len = (x_len > half) ? half : x_len;
1045+size_t x1_len = (x_len > half) ? x_len - half : 0;
1046+size_t y0_len = (y_len > half) ? half : y_len;
1047+size_t y1_len = (y_len > half) ? y_len - half : 0;
1048+1049+const mp_limb *x0 = x;
1050+const mp_limb *x1 = x + half;
1051+const mp_limb *y0 = y;
1052+const mp_limb *y1 = y + half;
1053+1054+/* Allocate scratch space */
1055+size_t eval_len = half + 1;
1056+size_t prod_len = 2 * eval_len;
1057+1058+size_t offset = 0;
1059+mp_limb *sum_x = scratch + offset; offset += eval_len; /* x0 + x1 */
1060+mp_limb *sum_y = scratch + offset; offset += eval_len; /* y0 + y1 */
1061+mp_limb *z0 = scratch + offset; offset += prod_len; /* x0 * y0 */
1062+mp_limb *z2 = scratch + offset; offset += prod_len; /* x1 * y1 */
1063+mp_limb *z1 = scratch + offset; offset += prod_len; /* (x0+x1)*(y0+y1) */
1064+mp_limb *recursive_scratch = scratch + offset;
1065+1066+/* Compute sum_x = x0 + x1 */
1067+mpn_zero(sum_x, eval_len);
1068+mpn_copyi(sum_x, x0, x0_len);
1069+if (x1_len > 0) {
1070+mpn_add(sum_x, sum_x, eval_len, x1, x1_len);
1071+ }
1072+size_t sum_x_len = eval_len;
1073+while (sum_x_len > 1 && sum_x[sum_x_len - 1] == 0) sum_x_len--;
1074+1075+/* Compute sum_y = y0 + y1 */
1076+mpn_zero(sum_y, eval_len);
1077+mpn_copyi(sum_y, y0, y0_len);
1078+if (y1_len > 0) {
1079+mpn_add(sum_y, sum_y, eval_len, y1, y1_len);
1080+ }
1081+size_t sum_y_len = eval_len;
1082+while (sum_y_len > 1 && sum_y[sum_y_len - 1] == 0) sum_y_len--;
1083+1084+/* z0 = x0 * y0 */
1085+mpn_zero(z0, prod_len);
1086+if (x0_len > 0 && y0_len > 0) {
1087+mpz_mul_karatsuba_limbs(z0, x0, x0_len, y0, y0_len, recursive_scratch);
1088+ }
1089+1090+/* z2 = x1 * y1 */
1091+mpn_zero(z2, prod_len);
1092+if (x1_len > 0 && y1_len > 0) {
1093+mpz_mul_karatsuba_limbs(z2, x1, x1_len, y1, y1_len, recursive_scratch);
1094+ }
1095+1096+/* z1 = (x0 + x1) * (y0 + y1) */
1097+mpn_zero(z1, prod_len);
1098+mpz_mul_karatsuba_limbs(z1, sum_x, sum_x_len, sum_y, sum_y_len, recursive_scratch);
1099+1100+/* z1 = z1 - z0 - z2 */
1101+mpn_sub(z1, z1, prod_len, z0, prod_len);
1102+mpn_sub(z1, z1, prod_len, z2, prod_len);
1103+1104+/*
1105+ * Combine: result = z2*B^(2*half) + z1*B^half + z0
1106+ */
1107+size_t result_len = x_len + y_len;
1108+mpn_zero(result, result_len);
1109+1110+/* Add z0 at position 0 */
1111+size_t z0_actual_len = prod_len;
1112+while (z0_actual_len > 0 && z0[z0_actual_len - 1] == 0) z0_actual_len--;
1113+if (z0_actual_len > 0) {
1114+mpn_copyi(result, z0, z0_actual_len);
1115+ }
1116+1117+/* Add z1 at position half */
1118+size_t z1_actual_len = prod_len;
1119+while (z1_actual_len > 0 && z1[z1_actual_len - 1] == 0) z1_actual_len--;
1120+if (z1_actual_len > 0) {
1121+limb_add_at(result, result_len, z1, z1_actual_len, half);
1122+ }
1123+1124+/* Add z2 at position 2*half */
1125+size_t z2_actual_len = prod_len;
1126+while (z2_actual_len > 0 && z2[z2_actual_len - 1] == 0) z2_actual_len--;
1127+if (z2_actual_len > 0 && 2 * half < result_len) {
1128+limb_add_at(result, result_len, z2, z2_actual_len, 2 * half);
1129+ }
1130+}
1131+9631132/*
9641133 * Toom-3 (Toom-Cook 3-way) Multiplication
9651134 *
@@ -980,11 +1149,15 @@ should_use_toom3(size_t n)
9801149return n >= TOOM3_THRESHOLD;
9811150}
9821151983-/* Calculate scratch space needed for Toom-3 */
1152+/* Calculate scratch space needed for Toom-3 (including Karatsuba at base) */
9841153static size_t
9851154toom3_scratch_size(size_t n)
9861155{
9871156if (!should_use_toom3(n)) {
1157+/* For Karatsuba range, return Karatsuba scratch size */
1158+if (should_use_karatsuba(n)) {
1159+return karatsuba_scratch_size(n);
1160+ }
9881161return 0;
9891162 }
9901163@@ -1158,14 +1331,19 @@ mpz_mul_toom3(mpz_ctx_t *ctx, mp_limb *result,
11581331mp_limb *scratch)
11591332{
11601333/*
1161- * Base case - use schoolbook.
1334+ * Base case - use Karatsuba or schoolbook.
11621335 * Toom-3 requires both operands to be large enough to avoid
11631336 * buffer overflow when writing at offset 4*third.
11641337 */
11651338size_t min_len = (x_len < y_len) ? x_len : y_len;
11661339size_t n = (x_len > y_len) ? x_len : y_len;
11671340if (!should_use_toom3(min_len)) {
1168-mpz_mul_basic_limbs(result, x, x_len, y, y_len);
1341+if (should_use_karatsuba(min_len)) {
1342+mpz_mul_karatsuba_limbs(result, x, x_len, y, y_len, scratch);
1343+ }
1344+else {
1345+mpz_mul_basic_limbs(result, x, x_len, y, y_len);
1346+ }
11691347return;
11701348 }
11711349@@ -2232,14 +2410,45 @@ mpz_mul(mpz_ctx_t *ctx, mpz_t *ww, mpz_t *u, mpz_t *v)
22322410 }
2233241122342412/*
2235- * Toom-3 requires both operands to be reasonably similar in size.
2236- * Use schoolbook for cases where the smaller operand is below threshold.
2413+ * Use schoolbook for small operands below Karatsuba threshold.
22372414 */
2238-if (!should_use_toom3(min_sz)) {
2415+if (!should_use_karatsuba(min_sz)) {
22392416mpz_mul_basic(ctx, ww, u, v);
22402417return;
22412418 }
224224192420+/*
2421+ * Karatsuba for medium-sized operands (KARATSUBA_THRESHOLD <= min_sz < TOOM3_THRESHOLD).
2422+ */
2423+if (!should_use_toom3(min_sz)) {
2424+size_t result_size = u->sz + v->sz;
2425+mpz_realloc(ctx, ww, result_size);
2426+2427+size_t scratch_size = karatsuba_scratch_size(max_sz);
2428+scratch_size += (scratch_size >> 3) + 16; /* safety margin */
2429+size_t pool_state = pool_save(ctx);
2430+mp_limb *scratch = NULL;
2431+2432+if (MPZ_HAS_POOL(ctx)) {
2433+scratch = pool_alloc(MPZ_POOL(ctx), scratch_size);
2434+ }
2435+2436+if (scratch) {
2437+mpz_mul_karatsuba_limbs(ww->p, u->p, u->sz, v->p, v->sz, scratch);
2438+pool_restore(ctx, pool_state);
2439+ }
2440+else {
2441+scratch = (mp_limb*)mrb_malloc(MPZ_MRB(ctx), scratch_size * sizeof(mp_limb));
2442+mpz_mul_karatsuba_limbs(ww->p, u->p, u->sz, v->p, v->sz, scratch);
2443+mrb_free(MPZ_MRB(ctx), scratch);
2444+ }
2445+2446+ww->sz = result_size;
2447+ww->sn = u->sn * v->sn;
2448+trim(ww);
2449+return;
2450+ }
2451+22432452/*
22442453 * Toom-3 writes at offset 4*third with products up to 2*(third+3)+16 limbs.
22452454 * Maximum write position: 4*n/3 + 2*n/3 + 22 = 2*n + 22, where n = max size.