3d6248b2275fb2d67ad6113bc3c76d46d4ddc25c
[libgcrypt.git] / mpi / mpi-mod.c
1 /* mpi-mod.c -  Modular reduction
2    Copyright (C) 1998, 1999, 2001, 2002, 2003,
3                  2007  Free Software Foundation, Inc.
4
5    This file is part of Libgcrypt.
6
7    Libgcrypt is free software; you can redistribute it and/or modify
8    it under the terms of the GNU Lesser General Public License as
9    published by the Free Software Foundation; either version 2.1 of
10    the License, or (at your option) any later version.
11
12    Libgcrypt is distributed in the hope that it will be useful,
13    but WITHOUT ANY WARRANTY; without even the implied warranty of
14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15    GNU Lesser General Public License for more details.
16
17    You should have received a copy of the GNU Lesser General Public
18    License along with this program; if not, write to the Free Software
19    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
20    USA.  */
21
22
23 #include <config.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26
27 #include "mpi-internal.h"
28 #include "longlong.h"
29 #include "g10lib.h"
30
31
32 /* Context used with Barrett reduction.  */
33 struct barrett_ctx_s
34 {
35   gcry_mpi_t m;   /* The modulus - may not be modified. */
36   int m_copied;   /* If true, M needs to be released.  */
37   int k;
38   gcry_mpi_t y;
39   gcry_mpi_t r1;  /* Helper MPI. */
40   gcry_mpi_t r2;  /* Helper MPI. */
41   gcry_mpi_t r3;  /* Helper MPI allocated on demand. */
42 };
43
44
45
46 void
47 _gcry_mpi_mod (gcry_mpi_t rem, gcry_mpi_t dividend, gcry_mpi_t divisor)
48 {
49   _gcry_mpi_fdiv_r (rem, dividend, divisor);
50   rem->sign = 0;
51 }
52
53
54 /* This function returns a new context for Barrett based operations on
55    the modulus M.  This context needs to be released using
56    _gcry_mpi_barrett_free.  If COPY is true M will be transferred to
57    the context and the user may change M.  If COPY is false, M may not
58    be changed until gcry_mpi_barrett_free has been called. */
59 mpi_barrett_t
60 _gcry_mpi_barrett_init (gcry_mpi_t m, int copy)
61 {
62   mpi_barrett_t ctx;
63   gcry_mpi_t tmp;
64
65   mpi_normalize (m);
66   ctx = gcry_xcalloc (1, sizeof *ctx);
67
68   if (copy)
69     {
70       ctx->m = mpi_copy (m);
71       ctx->m_copied = 1;
72     }
73   else
74     ctx->m = m;
75
76   ctx->k = mpi_get_nlimbs (m);
77   tmp = mpi_alloc (ctx->k + 1);
78
79   /* Barrett precalculation: y = floor(b^(2k) / m). */
80   mpi_set_ui (tmp, 1);
81   mpi_lshift_limbs (tmp, 2 * ctx->k);
82   mpi_fdiv_q (tmp, tmp, m);
83
84   ctx->y  = tmp;
85   ctx->r1 = mpi_alloc ( 2 * ctx->k + 1 );
86   ctx->r2 = mpi_alloc ( 2 * ctx->k + 1 );
87
88   return ctx;
89 }
90
91 void
92 _gcry_mpi_barrett_free (mpi_barrett_t ctx)
93 {
94   if (ctx)
95     {
96       mpi_free (ctx->y);
97       mpi_free (ctx->r1);
98       mpi_free (ctx->r2);
99       if (ctx->r3)
100         mpi_free (ctx->r3);
101       if (ctx->m_copied)
102         mpi_free (ctx->m);
103       gcry_free (ctx);
104     }
105 }
106
107
108 /* R = X mod M
109
110    Using Barrett reduction.  Before using this function
111    _gcry_mpi_barrett_init must have been called to do the
112    precalculations.  CTX is the context created by this precalculation
113    and also conveys M.  If the Barret reduction could no be done a
114    straightforward reduction method is used.
115
116    We assume that these conditions are met:
117    Input:  x =(x_2k-1 ...x_0)_b
118            m =(m_k-1 ....m_0)_b   with m_k-1 != 0
119    Output: r = x mod m
120  */
121 void
122 _gcry_mpi_mod_barrett (gcry_mpi_t r, gcry_mpi_t x, mpi_barrett_t ctx)
123 {
124   gcry_mpi_t m = ctx->m;
125   int k = ctx->k;
126   gcry_mpi_t y = ctx->y;
127   gcry_mpi_t r1 = ctx->r1;
128   gcry_mpi_t r2 = ctx->r2;
129   int sign;
130
131   mpi_normalize (x);
132   if (mpi_get_nlimbs (x) > 2*k )
133     {
134       mpi_mod (r, x, m);
135       return;
136     }
137
138   sign = x->sign;
139   x->sign = 0;
140
141   /* 1. q1 = floor( x / b^k-1)
142    *    q2 = q1 * y
143    *    q3 = floor( q2 / b^k+1 )
144    * Actually, we don't need qx, we can work direct on r2
145    */
146   mpi_set ( r2, x );
147   mpi_rshift_limbs ( r2, k-1 );
148   mpi_mul ( r2, r2, y );
149   mpi_rshift_limbs ( r2, k+1 );
150
151   /* 2. r1 = x mod b^k+1
152    *    r2 = q3 * m mod b^k+1
153    *    r  = r1 - r2
154    * 3. if r < 0 then  r = r + b^k+1
155    */
156   mpi_set ( r1, x );
157   if ( r1->nlimbs > k+1 ) /* Quick modulo operation.  */
158     r1->nlimbs = k+1;
159   mpi_mul ( r2, r2, m );
160   if ( r2->nlimbs > k+1 ) /* Quick modulo operation. */
161     r2->nlimbs = k+1;
162   mpi_sub ( r, r1, r2 );
163
164   if ( mpi_has_sign ( r ) )
165     {
166       if (!ctx->r3)
167         {
168           ctx->r3 = mpi_alloc ( k + 2 );
169           mpi_set_ui (ctx->r3, 1);
170           mpi_lshift_limbs (ctx->r3, k + 1 );
171         }
172       mpi_add ( r, r, ctx->r3 );
173     }
174
175   /* 4. while r >= m do r = r - m */
176   while ( mpi_cmp( r, m ) >= 0 )
177     mpi_sub ( r, r, m );
178
179   x->sign = sign;
180 }
181
182
183 void
184 _gcry_mpi_mul_barrett (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v,
185                        mpi_barrett_t ctx)
186 {
187   gcry_mpi_mul (w, u, v);
188   mpi_mod_barrett (w, w, ctx);
189 }