1/*	$NetBSD: bn_s_mp_sqr.c,v 1.2 2017/01/28 21:31:47 christos Exp $	*/
2
3#include <tommath.h>
4#ifdef BN_S_MP_SQR_C
5/* LibTomMath, multiple-precision integer library -- Tom St Denis
6 *
7 * LibTomMath is a library that provides multiple-precision
8 * integer arithmetic as well as number theoretic functionality.
9 *
10 * The library was designed directly after the MPI library by
11 * Michael Fromberger but has been written from scratch with
12 * additional optimizations in place.
13 *
14 * The library is free for all purposes without any express
15 * guarantee it works.
16 *
17 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
18 */
19
20/* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
21int s_mp_sqr (mp_int * a, mp_int * b)
22{
23  mp_int  t;
24  int     res, ix, iy, pa;
25  mp_word r;
26  mp_digit u, tmpx, *tmpt;
27
28  pa = a->used;
29  if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
30    return res;
31  }
32
33  /* default used is maximum possible size */
34  t.used = 2*pa + 1;
35
36  for (ix = 0; ix < pa; ix++) {
37    /* first calculate the digit at 2*ix */
38    /* calculate double precision result */
39    r = ((mp_word) t.dp[2*ix]) +
40        ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
41
42    /* store lower part in result */
43    t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
44
45    /* get the carry */
46    u           = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
47
48    /* left hand side of A[ix] * A[iy] */
49    tmpx        = a->dp[ix];
50
51    /* alias for where to store the results */
52    tmpt        = t.dp + (2*ix + 1);
53
54    for (iy = ix + 1; iy < pa; iy++) {
55      /* first calculate the product */
56      r       = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
57
58      /* now calculate the double precision result, note we use
59       * addition instead of *2 since it's easier to optimize
60       */
61      r       = ((mp_word) *tmpt) + r + r + ((mp_word) u);
62
63      /* store lower part */
64      *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
65
66      /* get carry */
67      u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
68    }
69    /* propagate upwards */
70    while (u != ((mp_digit) 0)) {
71      r       = ((mp_word) *tmpt) + ((mp_word) u);
72      *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
73      u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
74    }
75  }
76
77  mp_clamp (&t);
78  mp_exch (&t, b);
79  mp_clear (&t);
80  return MP_OKAY;
81}
82#endif
83
84/* Source: /cvs/libtom/libtommath/bn_s_mp_sqr.c,v  */
85/* Revision: 1.4  */
86/* Date: 2006/12/28 01:25:13  */
87