Line data    Source code

       1             : #include "tommath_private.h"
       2             : #ifdef BN_S_MP_KARATSUBA_SQR_C
       3             : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
       4             : /* SPDX-License-Identifier: Unlicense */
       5             : 
       6             : /* Karatsuba squaring, computes b = a*a using three
       7             :  * half size squarings
       8             :  *
       9             :  * See comments of karatsuba_mul for details.  It
      10             :  * is essentially the same algorithm but merely
      11             :  * tuned to perform recursive squarings.
      12             :  */
      13           0 : mp_err s_mp_karatsuba_sqr(const mp_int *a, mp_int *b)
      14             : {
      15             :    mp_int  x0, x1, t1, t2, x0x0, x1x1;
      16             :    int     B;
      17           0 :    mp_err  err = MP_MEM;
      18             : 
      19             :    /* min # of digits */
      20           0 :    B = a->used;
      21             : 
      22             :    /* now divide in two */
      23           0 :    B = B >> 1;
      24             : 
      25             :    /* init copy all the temps */
      26           0 :    if (mp_init_size(&x0, B) != MP_OKAY)
      27           0 :       goto LBL_ERR;
      28           0 :    if (mp_init_size(&x1, a->used - B) != MP_OKAY)
      29           0 :       goto X0;
      30             : 
      31             :    /* init temps */
      32           0 :    if (mp_init_size(&t1, a->used * 2) != MP_OKAY)
      33           0 :       goto X1;
      34           0 :    if (mp_init_size(&t2, a->used * 2) != MP_OKAY)
      35           0 :       goto T1;
      36           0 :    if (mp_init_size(&x0x0, B * 2) != MP_OKAY)
      37           0 :       goto T2;
      38           0 :    if (mp_init_size(&x1x1, (a->used - B) * 2) != MP_OKAY)
      39           0 :       goto X0X0;
      40             : 
      41             :    {
      42             :       int x;
      43             :       mp_digit *dst, *src;
      44             : 
      45           0 :       src = a->dp;
      46             : 
      47             :       /* now shift the digits */
      48           0 :       dst = x0.dp;
      49           0 :       for (x = 0; x < B; x++) {
      50           0 :          *dst++ = *src++;
      51             :       }
      52             : 
      53           0 :       dst = x1.dp;
      54           0 :       for (x = B; x < a->used; x++) {
      55           0 :          *dst++ = *src++;
      56             :       }
      57             :    }
      58             : 
      59           0 :    x0.used = B;
      60           0 :    x1.used = a->used - B;
      61             : 
      62           0 :    mp_clamp(&x0);
      63             : 
      64             :    /* now calc the products x0*x0 and x1*x1 */
      65           0 :    if (mp_sqr(&x0, &x0x0) != MP_OKAY)
      66           0 :       goto X1X1;           /* x0x0 = x0*x0 */
      67           0 :    if (mp_sqr(&x1, &x1x1) != MP_OKAY)
      68           0 :       goto X1X1;           /* x1x1 = x1*x1 */
      69             : 
      70             :    /* now calc (x1+x0)**2 */
      71           0 :    if (s_mp_add(&x1, &x0, &t1) != MP_OKAY)
      72           0 :       goto X1X1;           /* t1 = x1 - x0 */
      73           0 :    if (mp_sqr(&t1, &t1) != MP_OKAY)
      74           0 :       goto X1X1;           /* t1 = (x1 - x0) * (x1 - x0) */
      75             : 
      76             :    /* add x0y0 */
      77           0 :    if (s_mp_add(&x0x0, &x1x1, &t2) != MP_OKAY)
      78           0 :       goto X1X1;           /* t2 = x0x0 + x1x1 */
      79           0 :    if (s_mp_sub(&t1, &t2, &t1) != MP_OKAY)
      80           0 :       goto X1X1;           /* t1 = (x1+x0)**2 - (x0x0 + x1x1) */
      81             : 
      82             :    /* shift by B */
      83           0 :    if (mp_lshd(&t1, B) != MP_OKAY)
      84           0 :       goto X1X1;           /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */
      85           0 :    if (mp_lshd(&x1x1, B * 2) != MP_OKAY)
      86           0 :       goto X1X1;           /* x1x1 = x1x1 << 2*B */
      87             : 
      88           0 :    if (mp_add(&x0x0, &t1, &t1) != MP_OKAY)
      89           0 :       goto X1X1;           /* t1 = x0x0 + t1 */
      90           0 :    if (mp_add(&t1, &x1x1, b) != MP_OKAY)
      91           0 :       goto X1X1;           /* t1 = x0x0 + t1 + x1x1 */
      92             : 
      93           0 :    err = MP_OKAY;
      94             : 
      95           0 : X1X1:
      96           0 :    mp_clear(&x1x1);
      97           0 : X0X0:
      98           0 :    mp_clear(&x0x0);
      99           0 : T2:
     100           0 :    mp_clear(&t2);
     101           0 : T1:
     102           0 :    mp_clear(&t1);
     103           0 : X1:
     104           0 :    mp_clear(&x1);
     105           0 : X0:
     106           0 :    mp_clear(&x0);
     107           0 : LBL_ERR:
     108           0 :    return err;
     109             : }
     110             : #endif