Line data    Source code

       1             : #include "tommath_private.h"
       2             : #ifdef BN_S_MP_TOOM_SQR_C
       3             : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
       4             : /* SPDX-License-Identifier: Unlicense */
       5             : 
       6             : /* squaring using Toom-Cook 3-way algorithm */
       7             : 
       8             : /*
       9             :    This file contains code from J. Arndt's book  "Matters Computational"
      10             :    and the accompanying FXT-library with permission of the author.
      11             : */
      12             : 
      13             : /* squaring using Toom-Cook 3-way algorithm */
      14             : /*
      15             :    Setup and interpolation from algorithm SQR_3 in
      16             : 
      17             :      Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."
      18             :      18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.
      19             : 
      20             : */
      21           0 : mp_err s_mp_toom_sqr(const mp_int *a, mp_int *b)
      22             : {
      23           0 :    mp_int S0, a0, a1, a2;
      24           0 :    mp_digit *tmpa, *tmpc;
      25           0 :    int B, count;
      26           0 :    mp_err err;
      27             : 
      28             : 
      29             :    /* init temps */
      30           0 :    if ((err = mp_init(&S0)) != MP_OKAY) {
      31           0 :       return err;
      32             :    }
      33             : 
      34             :    /* B */
      35           0 :    B = a->used / 3;
      36             : 
      37             :    /** a = a2 * x^2 + a1 * x + a0; */
      38           0 :    if ((err = mp_init_size(&a0, B)) != MP_OKAY)                   goto LBL_ERRa0;
      39             : 
      40           0 :    a0.used = B;
      41           0 :    if ((err = mp_init_size(&a1, B)) != MP_OKAY)                   goto LBL_ERRa1;
      42           0 :    a1.used = B;
      43           0 :    if ((err = mp_init_size(&a2, B + (a->used - (3 * B)))) != MP_OKAY) goto LBL_ERRa2;
      44             : 
      45           0 :    tmpa = a->dp;
      46           0 :    tmpc = a0.dp;
      47           0 :    for (count = 0; count < B; count++) {
      48           0 :       *tmpc++ = *tmpa++;
      49             :    }
      50           0 :    tmpc = a1.dp;
      51           0 :    for (; count < (2 * B); count++) {
      52           0 :       *tmpc++ = *tmpa++;
      53             :    }
      54           0 :    tmpc = a2.dp;
      55           0 :    for (; count < a->used; count++) {
      56           0 :       *tmpc++ = *tmpa++;
      57           0 :       a2.used++;
      58             :    }
      59           0 :    mp_clamp(&a0);
      60           0 :    mp_clamp(&a1);
      61           0 :    mp_clamp(&a2);
      62             : 
      63             :    /** S0 = a0^2;  */
      64           0 :    if ((err = mp_sqr(&a0, &S0)) != MP_OKAY)                       goto LBL_ERR;
      65             : 
      66             :    /** \\S1 = (a2 + a1 + a0)^2 */
      67             :    /** \\S2 = (a2 - a1 + a0)^2  */
      68             :    /** \\S1 = a0 + a2; */
      69             :    /** a0 = a0 + a2; */
      70           0 :    if ((err = mp_add(&a0, &a2, &a0)) != MP_OKAY)                  goto LBL_ERR;
      71             :    /** \\S2 = S1 - a1; */
      72             :    /** b = a0 - a1; */
      73           0 :    if ((err = mp_sub(&a0, &a1, b)) != MP_OKAY)                    goto LBL_ERR;
      74             :    /** \\S1 = S1 + a1; */
      75             :    /** a0 = a0 + a1; */
      76           0 :    if ((err = mp_add(&a0, &a1, &a0)) != MP_OKAY)                  goto LBL_ERR;
      77             :    /** \\S1 = S1^2;  */
      78             :    /** a0 = a0^2; */
      79           0 :    if ((err = mp_sqr(&a0, &a0)) != MP_OKAY)                       goto LBL_ERR;
      80             :    /** \\S2 = S2^2;  */
      81             :    /** b = b^2; */
      82           0 :    if ((err = mp_sqr(b, b)) != MP_OKAY)                           goto LBL_ERR;
      83             : 
      84             :    /** \\ S3 = 2 * a1 * a2  */
      85             :    /** \\S3 = a1 * a2;  */
      86             :    /** a1 = a1 * a2; */
      87           0 :    if ((err = mp_mul(&a1, &a2, &a1)) != MP_OKAY)                  goto LBL_ERR;
      88             :    /** \\S3 = S3 << 1;  */
      89             :    /** a1 = a1 << 1; */
      90           0 :    if ((err = mp_mul_2(&a1, &a1)) != MP_OKAY)                     goto LBL_ERR;
      91             : 
      92             :    /** \\S4 = a2^2;  */
      93             :    /** a2 = a2^2; */
      94           0 :    if ((err = mp_sqr(&a2, &a2)) != MP_OKAY)                       goto LBL_ERR;
      95             : 
      96             :    /** \\ tmp = (S1 + S2)/2  */
      97             :    /** \\tmp = S1 + S2; */
      98             :    /** b = a0 + b; */
      99           0 :    if ((err = mp_add(&a0, b, b)) != MP_OKAY)                      goto LBL_ERR;
     100             :    /** \\tmp = tmp >> 1; */
     101             :    /** b = b >> 1; */
     102           0 :    if ((err = mp_div_2(b, b)) != MP_OKAY)                         goto LBL_ERR;
     103             : 
     104             :    /** \\ S1 = S1 - tmp - S3  */
     105             :    /** \\S1 = S1 - tmp; */
     106             :    /** a0 = a0 - b; */
     107           0 :    if ((err = mp_sub(&a0, b, &a0)) != MP_OKAY)                    goto LBL_ERR;
     108             :    /** \\S1 = S1 - S3;  */
     109             :    /** a0 = a0 - a1; */
     110           0 :    if ((err = mp_sub(&a0, &a1, &a0)) != MP_OKAY)                  goto LBL_ERR;
     111             : 
     112             :    /** \\S2 = tmp - S4 -S0  */
     113             :    /** \\S2 = tmp - S4;  */
     114             :    /** b = b - a2; */
     115           0 :    if ((err = mp_sub(b, &a2, b)) != MP_OKAY)                      goto LBL_ERR;
     116             :    /** \\S2 = S2 - S0;  */
     117             :    /** b = b - S0; */
     118           0 :    if ((err = mp_sub(b, &S0, b)) != MP_OKAY)                      goto LBL_ERR;
     119             : 
     120             : 
     121             :    /** \\P = S4*x^4 + S3*x^3 + S2*x^2 + S1*x + S0; */
     122             :    /** P = a2*x^4 + a1*x^3 + b*x^2 + a0*x + S0; */
     123             : 
     124           0 :    if ((err = mp_lshd(&a2, 4 * B)) != MP_OKAY)                    goto LBL_ERR;
     125           0 :    if ((err = mp_lshd(&a1, 3 * B)) != MP_OKAY)                    goto LBL_ERR;
     126           0 :    if ((err = mp_lshd(b, 2 * B)) != MP_OKAY)                      goto LBL_ERR;
     127           0 :    if ((err = mp_lshd(&a0, 1 * B)) != MP_OKAY)                    goto LBL_ERR;
     128           0 :    if ((err = mp_add(&a2, &a1, &a2)) != MP_OKAY)                  goto LBL_ERR;
     129           0 :    if ((err = mp_add(&a2, b, b)) != MP_OKAY)                      goto LBL_ERR;
     130           0 :    if ((err = mp_add(b, &a0, b)) != MP_OKAY)                      goto LBL_ERR;
     131           0 :    if ((err = mp_add(b, &S0, b)) != MP_OKAY)                      goto LBL_ERR;
     132             :    /** a^2 - P  */
     133             : 
     134             : 
     135           0 : LBL_ERR:
     136           0 :    mp_clear(&a2);
     137           0 : LBL_ERRa2:
     138           0 :    mp_clear(&a1);
     139           0 : LBL_ERRa1:
     140           0 :    mp_clear(&a0);
     141           0 : LBL_ERRa0:
     142           0 :    mp_clear(&S0);
     143             : 
     144           0 :    return err;
     145             : }
     146             : 
     147             : #endif