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        2869 : mp_err s_mp_karatsuba_sqr(const mp_int *a, mp_int *b)
      14             : {
      15           0 :    mp_int  x0, x1, t1, t2, x0x0, x1x1;
      16           0 :    int     B;
      17        2869 :    mp_err  err = MP_MEM;
      18             : 
      19             :    /* min # of digits */
      20        2869 :    B = a->used;
      21             : 
      22             :    /* now divide in two */
      23        2869 :    B = B >> 1;
      24             : 
      25             :    /* init copy all the temps */
      26        2869 :    if (mp_init_size(&x0, B) != MP_OKAY)
      27           0 :       goto LBL_ERR;
      28        2869 :    if (mp_init_size(&x1, a->used - B) != MP_OKAY)
      29           0 :       goto X0;
      30             : 
      31             :    /* init temps */
      32        2869 :    if (mp_init_size(&t1, a->used * 2) != MP_OKAY)
      33           0 :       goto X1;
      34        2869 :    if (mp_init_size(&t2, a->used * 2) != MP_OKAY)
      35           0 :       goto T1;
      36        2869 :    if (mp_init_size(&x0x0, B * 2) != MP_OKAY)
      37           0 :       goto T2;
      38        2869 :    if (mp_init_size(&x1x1, (a->used - B) * 2) != MP_OKAY)
      39           0 :       goto X0X0;
      40             : 
      41             :    {
      42           0 :       int x;
      43           0 :       mp_digit *dst, *src;
      44             : 
      45        2869 :       src = a->dp;
      46             : 
      47             :       /* now shift the digits */
      48        2869 :       dst = x0.dp;
      49      197961 :       for (x = 0; x < B; x++) {
      50      195092 :          *dst++ = *src++;
      51             :       }
      52             : 
      53        2869 :       dst = x1.dp;
      54      200830 :       for (x = B; x < a->used; x++) {
      55      197961 :          *dst++ = *src++;
      56             :       }
      57             :    }
      58             : 
      59        2869 :    x0.used = B;
      60        2869 :    x1.used = a->used - B;
      61             : 
      62        2869 :    mp_clamp(&x0);
      63             : 
      64             :    /* now calc the products x0*x0 and x1*x1 */
      65        2869 :    if (mp_sqr(&x0, &x0x0) != MP_OKAY)
      66           0 :       goto X1X1;           /* x0x0 = x0*x0 */
      67        2869 :    if (mp_sqr(&x1, &x1x1) != MP_OKAY)
      68           0 :       goto X1X1;           /* x1x1 = x1*x1 */
      69             : 
      70             :    /* now calc (x1+x0)**2 */
      71        2869 :    if (s_mp_add(&x1, &x0, &t1) != MP_OKAY)
      72           0 :       goto X1X1;           /* t1 = x1 - x0 */
      73        2869 :    if (mp_sqr(&t1, &t1) != MP_OKAY)
      74           0 :       goto X1X1;           /* t1 = (x1 - x0) * (x1 - x0) */
      75             : 
      76             :    /* add x0y0 */
      77        2869 :    if (s_mp_add(&x0x0, &x1x1, &t2) != MP_OKAY)
      78           0 :       goto X1X1;           /* t2 = x0x0 + x1x1 */
      79        2869 :    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        2869 :    if (mp_lshd(&t1, B) != MP_OKAY)
      84           0 :       goto X1X1;           /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */
      85        2869 :    if (mp_lshd(&x1x1, B * 2) != MP_OKAY)
      86           0 :       goto X1X1;           /* x1x1 = x1x1 << 2*B */
      87             : 
      88        2869 :    if (mp_add(&x0x0, &t1, &t1) != MP_OKAY)
      89           0 :       goto X1X1;           /* t1 = x0x0 + t1 */
      90        2869 :    if (mp_add(&t1, &x1x1, b) != MP_OKAY)
      91           0 :       goto X1X1;           /* t1 = x0x0 + t1 + x1x1 */
      92             : 
      93        2869 :    err = MP_OKAY;
      94             : 
      95        2869 : X1X1:
      96        2869 :    mp_clear(&x1x1);
      97        2869 : X0X0:
      98        2869 :    mp_clear(&x0x0);
      99        2869 : T2:
     100        2869 :    mp_clear(&t2);
     101        2869 : T1:
     102        2869 :    mp_clear(&t1);
     103        2869 : X1:
     104        2869 :    mp_clear(&x1);
     105        2869 : X0:
     106        2869 :    mp_clear(&x0);
     107        2869 : LBL_ERR:
     108        2869 :    return err;
     109             : }
     110             : #endif