Line data    Source code

       1             : #include "tommath_private.h"
       2             : #ifdef BN_S_MP_KARATSUBA_MUL_C
       3             : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
       4             : /* SPDX-License-Identifier: Unlicense */
       5             : 
       6             : /* c = |a| * |b| using Karatsuba Multiplication using
       7             :  * three half size multiplications
       8             :  *
       9             :  * Let B represent the radix [e.g. 2**MP_DIGIT_BIT] and
      10             :  * let n represent half of the number of digits in
      11             :  * the min(a,b)
      12             :  *
      13             :  * a = a1 * B**n + a0
      14             :  * b = b1 * B**n + b0
      15             :  *
      16             :  * Then, a * b =>
      17             :    a1b1 * B**2n + ((a1 + a0)(b1 + b0) - (a0b0 + a1b1)) * B + a0b0
      18             :  *
      19             :  * Note that a1b1 and a0b0 are used twice and only need to be
      20             :  * computed once.  So in total three half size (half # of
      21             :  * digit) multiplications are performed, a0b0, a1b1 and
      22             :  * (a1+b1)(a0+b0)
      23             :  *
      24             :  * Note that a multiplication of half the digits requires
      25             :  * 1/4th the number of single precision multiplications so in
      26             :  * total after one call 25% of the single precision multiplications
      27             :  * are saved.  Note also that the call to mp_mul can end up back
      28             :  * in this function if the a0, a1, b0, or b1 are above the threshold.
      29             :  * This is known as divide-and-conquer and leads to the famous
      30             :  * O(N**lg(3)) or O(N**1.584) work which is asymptopically lower than
      31             :  * the standard O(N**2) that the baseline/comba methods use.
      32             :  * Generally though the overhead of this method doesn't pay off
      33             :  * until a certain size (N ~ 80) is reached.
      34             :  */
      35         846 : mp_err s_mp_karatsuba_mul(const mp_int *a, const mp_int *b, mp_int *c)
      36             : {
      37           0 :    mp_int  x0, x1, y0, y1, t1, x0y0, x1y1;
      38           0 :    int     B;
      39         846 :    mp_err  err = MP_MEM; /* default the return code to an error */
      40             : 
      41             :    /* min # of digits */
      42         846 :    B = MP_MIN(a->used, b->used);
      43             : 
      44             :    /* now divide in two */
      45         846 :    B = B >> 1;
      46             : 
      47             :    /* init copy all the temps */
      48         846 :    if (mp_init_size(&x0, B) != MP_OKAY) {
      49           0 :       goto LBL_ERR;
      50             :    }
      51         846 :    if (mp_init_size(&x1, a->used - B) != MP_OKAY) {
      52           0 :       goto X0;
      53             :    }
      54         846 :    if (mp_init_size(&y0, B) != MP_OKAY) {
      55           0 :       goto X1;
      56             :    }
      57         846 :    if (mp_init_size(&y1, b->used - B) != MP_OKAY) {
      58           0 :       goto Y0;
      59             :    }
      60             : 
      61             :    /* init temps */
      62         846 :    if (mp_init_size(&t1, B * 2) != MP_OKAY) {
      63           0 :       goto Y1;
      64             :    }
      65         846 :    if (mp_init_size(&x0y0, B * 2) != MP_OKAY) {
      66           0 :       goto T1;
      67             :    }
      68         846 :    if (mp_init_size(&x1y1, B * 2) != MP_OKAY) {
      69           0 :       goto X0Y0;
      70             :    }
      71             : 
      72             :    /* now shift the digits */
      73         846 :    x0.used = y0.used = B;
      74         846 :    x1.used = a->used - B;
      75         846 :    y1.used = b->used - B;
      76             : 
      77             :    {
      78           0 :       int x;
      79           0 :       mp_digit *tmpa, *tmpb, *tmpx, *tmpy;
      80             : 
      81             :       /* we copy the digits directly instead of using higher level functions
      82             :        * since we also need to shift the digits
      83             :        */
      84         846 :       tmpa = a->dp;
      85         846 :       tmpb = b->dp;
      86             : 
      87         846 :       tmpx = x0.dp;
      88         846 :       tmpy = y0.dp;
      89       58374 :       for (x = 0; x < B; x++) {
      90       57528 :          *tmpx++ = *tmpa++;
      91       57528 :          *tmpy++ = *tmpb++;
      92             :       }
      93             : 
      94         846 :       tmpx = x1.dp;
      95       59220 :       for (x = B; x < a->used; x++) {
      96       58374 :          *tmpx++ = *tmpa++;
      97             :       }
      98             : 
      99         846 :       tmpy = y1.dp;
     100       59220 :       for (x = B; x < b->used; x++) {
     101       58374 :          *tmpy++ = *tmpb++;
     102             :       }
     103             :    }
     104             : 
     105             :    /* only need to clamp the lower words since by definition the
     106             :     * upper words x1/y1 must have a known number of digits
     107             :     */
     108         846 :    mp_clamp(&x0);
     109         846 :    mp_clamp(&y0);
     110             : 
     111             :    /* now calc the products x0y0 and x1y1 */
     112             :    /* after this x0 is no longer required, free temp [x0==t2]! */
     113         846 :    if (mp_mul(&x0, &y0, &x0y0) != MP_OKAY) {
     114           0 :       goto X1Y1;          /* x0y0 = x0*y0 */
     115             :    }
     116         846 :    if (mp_mul(&x1, &y1, &x1y1) != MP_OKAY) {
     117           0 :       goto X1Y1;          /* x1y1 = x1*y1 */
     118             :    }
     119             : 
     120             :    /* now calc x1+x0 and y1+y0 */
     121         846 :    if (s_mp_add(&x1, &x0, &t1) != MP_OKAY) {
     122           0 :       goto X1Y1;          /* t1 = x1 - x0 */
     123             :    }
     124         846 :    if (s_mp_add(&y1, &y0, &x0) != MP_OKAY) {
     125           0 :       goto X1Y1;          /* t2 = y1 - y0 */
     126             :    }
     127         846 :    if (mp_mul(&t1, &x0, &t1) != MP_OKAY) {
     128           0 :       goto X1Y1;          /* t1 = (x1 + x0) * (y1 + y0) */
     129             :    }
     130             : 
     131             :    /* add x0y0 */
     132         846 :    if (mp_add(&x0y0, &x1y1, &x0) != MP_OKAY) {
     133           0 :       goto X1Y1;          /* t2 = x0y0 + x1y1 */
     134             :    }
     135         846 :    if (s_mp_sub(&t1, &x0, &t1) != MP_OKAY) {
     136           0 :       goto X1Y1;          /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */
     137             :    }
     138             : 
     139             :    /* shift by B */
     140         846 :    if (mp_lshd(&t1, B) != MP_OKAY) {
     141           0 :       goto X1Y1;          /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
     142             :    }
     143         846 :    if (mp_lshd(&x1y1, B * 2) != MP_OKAY) {
     144           0 :       goto X1Y1;          /* x1y1 = x1y1 << 2*B */
     145             :    }
     146             : 
     147         846 :    if (mp_add(&x0y0, &t1, &t1) != MP_OKAY) {
     148           0 :       goto X1Y1;          /* t1 = x0y0 + t1 */
     149             :    }
     150         846 :    if (mp_add(&t1, &x1y1, c) != MP_OKAY) {
     151           0 :       goto X1Y1;          /* t1 = x0y0 + t1 + x1y1 */
     152             :    }
     153             : 
     154             :    /* Algorithm succeeded set the return code to MP_OKAY */
     155         846 :    err = MP_OKAY;
     156             : 
     157         846 : X1Y1:
     158         846 :    mp_clear(&x1y1);
     159         846 : X0Y0:
     160         846 :    mp_clear(&x0y0);
     161         846 : T1:
     162         846 :    mp_clear(&t1);
     163         846 : Y1:
     164         846 :    mp_clear(&y1);
     165         846 : Y0:
     166         846 :    mp_clear(&y0);
     167         846 : X1:
     168         846 :    mp_clear(&x1);
     169         846 : X0:
     170         846 :    mp_clear(&x0);
     171         846 : LBL_ERR:
     172         846 :    return err;
     173             : }
     174             : #endif