9f38a0a945687b205575731c366378a498369640
[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 #include <assert.h>
27
28 #include "mpi-internal.h"
29 #include "longlong.h"
30 #include "g10lib.h"
31
32
33 /* Context used with Barrett reduction.  */
34 struct barrett_ctx_s
35 {
36   gcry_mpi_t m;   /* The modulus - may not be modified. */
37   int m_copied;   /* If true, M needs to be released.  */
38   int k;
39   gcry_mpi_t y; 
40   gcry_mpi_t r1;  /* Helper MPI. */
41   gcry_mpi_t r2;  /* Helper MPI. */
42   gcry_mpi_t r3;  /* Helper MPI allocated on demand. */
43 };
44
45
46
47 void
48 _gcry_mpi_mod (gcry_mpi_t rem, gcry_mpi_t dividend, gcry_mpi_t divisor)
49 {
50   _gcry_mpi_fdiv_r (rem, dividend, divisor);
51   rem->sign = 0;
52 }
53
54
55 /* This function returns a new context for Barrett based operations on
56    the modulus M.  This context needs to be released using
57    _gcry_mpi_barrett_free.  If COPY is true M will be transferred to
58    the context and the user may change M.  If COPY is false, M may not
59    be changed until gcry_mpi_barrett_free has been called. */
60 mpi_barrett_t 
61 _gcry_mpi_barrett_init (gcry_mpi_t m, int copy)
62 {
63   mpi_barrett_t ctx;
64   gcry_mpi_t tmp;
65
66   mpi_normalize (m);
67   ctx = gcry_xcalloc (1, sizeof *ctx);
68
69   if (copy)
70     {
71       ctx->m = mpi_copy (m);
72       ctx->m_copied = 1;
73     }
74   else
75     ctx->m = m;
76
77   ctx->k = mpi_get_nlimbs (m);
78   tmp = mpi_alloc (ctx->k + 1);
79
80   /* Barrett precalculation: y = floor(b^(2k) / m). */
81   mpi_set_ui (tmp, 1);
82   mpi_lshift_limbs (tmp, 2 * ctx->k);
83   mpi_fdiv_q (tmp, tmp, m);
84
85   ctx->y  = tmp;
86   ctx->r1 = mpi_alloc ( 2 * ctx->k + 1 );
87   ctx->r2 = mpi_alloc ( 2 * ctx->k + 1 );
88
89   return ctx;
90 }
91
92 void
93 _gcry_mpi_barrett_free (mpi_barrett_t ctx)
94 {
95   if (ctx)
96     {
97       mpi_free (ctx->y);
98       mpi_free (ctx->r1);
99       mpi_free (ctx->r2);
100       if (ctx->r3)
101         mpi_free (ctx->r3);
102       if (ctx->m_copied)
103         mpi_free (ctx->m);
104       gcry_free (ctx);
105     }
106 }
107
108
109 /* R = X mod M
110
111    Using Barrett reduction.  Before using this function
112    _gcry_mpi_barrett_init must have been called to do the
113    precalculations.  CTX is the context created by this precalculation
114    and also conveys M.  If the Barret reduction could no be done a
115    starightforward reduction method is used.
116
117    We assume that these conditions are met:
118    Input:  x =(x_2k-1 ...x_0)_b
119            m =(m_k-1 ....m_0)_b   with m_k-1 != 0
120    Output: r = x mod m
121  */
122 void
123 _gcry_mpi_mod_barrett (gcry_mpi_t r, gcry_mpi_t x, mpi_barrett_t ctx)
124 {
125   gcry_mpi_t m = ctx->m;
126   int k = ctx->k;
127   gcry_mpi_t y = ctx->y;
128   gcry_mpi_t r1 = ctx->r1;
129   gcry_mpi_t r2 = ctx->r2;
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   /* 1. q1 = floor( x / b^k-1)
139    *    q2 = q1 * y
140    *    q3 = floor( q2 / b^k+1 )
141    * Actually, we don't need qx, we can work direct on r2 
142    */
143   mpi_set ( r2, x );
144   mpi_rshift_limbs ( r2, k-1 );
145   mpi_mul ( r2, r2, y );
146   mpi_rshift_limbs ( r2, k+1 );
147
148   /* 2. r1 = x mod b^k+1
149    *    r2 = q3 * m mod b^k+1
150    *    r  = r1 - r2
151    * 3. if r < 0 then  r = r + b^k+1
152    */
153   mpi_set ( r1, x );
154   if ( r1->nlimbs > k+1 ) /* Quick modulo operation.  */
155     r1->nlimbs = k+1;
156   mpi_mul ( r2, r2, m );
157   if ( r2->nlimbs > k+1 ) /* Quick modulo operation. */
158     r2->nlimbs = k+1;
159   mpi_sub ( r, r1, r2 );
160
161   if ( mpi_is_neg( r ) )
162     {
163       if (!ctx->r3)
164         {
165           ctx->r3 = mpi_alloc ( k + 2 );
166           mpi_set_ui (ctx->r3, 1);
167           mpi_lshift_limbs (ctx->r3, k + 1 );
168         }
169       mpi_add ( r, r, ctx->r3 );
170     }
171   
172   /* 4. while r >= m do r = r - m */
173   while ( mpi_cmp( r, m ) >= 0 )
174     mpi_sub ( r, r, m );
175
176 }
177
178
179 void
180 _gcry_mpi_mul_barrett (gcry_mpi_t w, gcry_mpi_t u, gcry_mpi_t v,
181                        mpi_barrett_t ctx)
182 {
183   gcry_mpi_mul (w, u, v);
184   mpi_mod_barrett (w, w, ctx);
185 }
186