Line data    Source code

       1             : #include "tommath_private.h"
       2             : #ifdef BN_S_MP_SQR_C
       3             : /* LibTomMath, multiple-precision integer library -- Tom St Denis */
       4             : /* SPDX-License-Identifier: Unlicense */
       5             : 
       6             : /* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
       7           0 : mp_err s_mp_sqr(const mp_int *a, mp_int *b)
       8             : {
       9           0 :    mp_int   t;
      10           0 :    int      ix, iy, pa;
      11           0 :    mp_err   err;
      12           0 :    mp_word  r;
      13           0 :    mp_digit u, tmpx, *tmpt;
      14             : 
      15           0 :    pa = a->used;
      16           0 :    if ((err = mp_init_size(&t, (2 * pa) + 1)) != MP_OKAY) {
      17           0 :       return err;
      18             :    }
      19             : 
      20             :    /* default used is maximum possible size */
      21           0 :    t.used = (2 * pa) + 1;
      22             : 
      23           0 :    for (ix = 0; ix < pa; ix++) {
      24             :       /* first calculate the digit at 2*ix */
      25             :       /* calculate double precision result */
      26           0 :       r = (mp_word)t.dp[2*ix] +
      27           0 :           ((mp_word)a->dp[ix] * (mp_word)a->dp[ix]);
      28             : 
      29             :       /* store lower part in result */
      30           0 :       t.dp[ix+ix] = (mp_digit)(r & (mp_word)MP_MASK);
      31             : 
      32             :       /* get the carry */
      33           0 :       u           = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
      34             : 
      35             :       /* left hand side of A[ix] * A[iy] */
      36           0 :       tmpx        = a->dp[ix];
      37             : 
      38             :       /* alias for where to store the results */
      39           0 :       tmpt        = t.dp + ((2 * ix) + 1);
      40             : 
      41           0 :       for (iy = ix + 1; iy < pa; iy++) {
      42             :          /* first calculate the product */
      43           0 :          r       = (mp_word)tmpx * (mp_word)a->dp[iy];
      44             : 
      45             :          /* now calculate the double precision result, note we use
      46             :           * addition instead of *2 since it's easier to optimize
      47             :           */
      48           0 :          r       = (mp_word)*tmpt + r + r + (mp_word)u;
      49             : 
      50             :          /* store lower part */
      51           0 :          *tmpt++ = (mp_digit)(r & (mp_word)MP_MASK);
      52             : 
      53             :          /* get carry */
      54           0 :          u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
      55             :       }
      56             :       /* propagate upwards */
      57           0 :       while (u != 0uL) {
      58           0 :          r       = (mp_word)*tmpt + (mp_word)u;
      59           0 :          *tmpt++ = (mp_digit)(r & (mp_word)MP_MASK);
      60           0 :          u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
      61             :       }
      62             :    }
      63             : 
      64           0 :    mp_clamp(&t);
      65           0 :    mp_exch(&t, b);
      66           0 :    mp_clear(&t);
      67           0 :    return MP_OKAY;
      68             : }
      69             : #endif