2003-06-16 Moritz Schulte <moritz@g10code.com>
[libgcrypt.git] / mpi / mpi-div.c
1 /* mpi-div.c  -  MPI functions
2  * Copyright (C) 1994, 1996, 1998, 2001, 2002, 2003 Free Software Foundation, Inc.
3  *
4  * This file is part of Libgcrypt.
5  *
6  * Libgcrypt is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as
8  * published by the Free Software Foundation; either version 2.1 of
9  * the License, or (at your option) any later version.
10  *
11  * Libgcrypt is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this program; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
19  *
20  * Note: This code is heavily based on the GNU MP Library.
21  *       Actually it's the same code with only minor changes in the
22  *       way the data is stored; this is to support the abstraction
23  *       of an optional secure memory allocation which may be used
24  *       to avoid revealing of sensitive data due to paging etc.
25  */
26
27 #include <config.h>
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include "mpi-internal.h"
31 #include "longlong.h"
32 #include "g10lib.h"
33
34
35 void
36 _gcry_mpi_fdiv_r( gcry_mpi_t rem, gcry_mpi_t dividend, gcry_mpi_t divisor )
37 {
38     int divisor_sign = divisor->sign;
39     gcry_mpi_t temp_divisor = NULL;
40
41     /* We need the original value of the divisor after the remainder has been
42      * preliminary calculated.  We have to copy it to temporary space if it's
43      * the same variable as REM.  */
44     if( rem == divisor ) {
45         temp_divisor = mpi_copy( divisor );
46         divisor = temp_divisor;
47     }
48
49     _gcry_mpi_tdiv_r( rem, dividend, divisor );
50
51     if( ((divisor_sign?1:0) ^ (dividend->sign?1:0)) && rem->nlimbs )
52         gcry_mpi_add( rem, rem, divisor);
53
54     if( temp_divisor )
55         mpi_free(temp_divisor);
56 }
57
58
59
60 /****************
61  * Division rounding the quotient towards -infinity.
62  * The remainder gets the same sign as the denominator.
63  * rem is optional
64  */
65
66 ulong
67 _gcry_mpi_fdiv_r_ui( gcry_mpi_t rem, gcry_mpi_t dividend, ulong divisor )
68 {
69     mpi_limb_t rlimb;
70
71     rlimb = _gcry_mpih_mod_1( dividend->d, dividend->nlimbs, divisor );
72     if( rlimb && dividend->sign )
73         rlimb = divisor - rlimb;
74
75     if( rem ) {
76         rem->d[0] = rlimb;
77         rem->nlimbs = rlimb? 1:0;
78     }
79     return rlimb;
80 }
81
82
83 void
84 _gcry_mpi_fdiv_q( gcry_mpi_t quot, gcry_mpi_t dividend, gcry_mpi_t divisor )
85 {
86     gcry_mpi_t tmp = mpi_alloc( mpi_get_nlimbs(quot) );
87     _gcry_mpi_fdiv_qr( quot, tmp, dividend, divisor);
88     mpi_free(tmp);
89 }
90
91 void
92 _gcry_mpi_fdiv_qr( gcry_mpi_t quot, gcry_mpi_t rem, gcry_mpi_t dividend, gcry_mpi_t divisor )
93 {
94     int divisor_sign = divisor->sign;
95     gcry_mpi_t temp_divisor = NULL;
96
97     if( quot == divisor || rem == divisor ) {
98         temp_divisor = mpi_copy( divisor );
99         divisor = temp_divisor;
100     }
101
102     _gcry_mpi_tdiv_qr( quot, rem, dividend, divisor );
103
104     if( (divisor_sign ^ dividend->sign) && rem->nlimbs ) {
105         gcry_mpi_sub_ui( quot, quot, 1 );
106         gcry_mpi_add( rem, rem, divisor);
107     }
108
109     if( temp_divisor )
110         mpi_free(temp_divisor);
111 }
112
113
114 /* If den == quot, den needs temporary storage.
115  * If den == rem, den needs temporary storage.
116  * If num == quot, num needs temporary storage.
117  * If den has temporary storage, it can be normalized while being copied,
118  *   i.e no extra storage should be allocated.
119  */
120
121 void
122 _gcry_mpi_tdiv_r( gcry_mpi_t rem, gcry_mpi_t num, gcry_mpi_t den)
123 {
124     _gcry_mpi_tdiv_qr(NULL, rem, num, den );
125 }
126
127 void
128 _gcry_mpi_tdiv_qr( gcry_mpi_t quot, gcry_mpi_t rem, gcry_mpi_t num, gcry_mpi_t den)
129 {
130     mpi_ptr_t np, dp;
131     mpi_ptr_t qp, rp;
132     mpi_size_t nsize = num->nlimbs;
133     mpi_size_t dsize = den->nlimbs;
134     mpi_size_t qsize, rsize;
135     mpi_size_t sign_remainder = num->sign;
136     mpi_size_t sign_quotient = num->sign ^ den->sign;
137     unsigned normalization_steps;
138     mpi_limb_t q_limb;
139     mpi_ptr_t marker[5];
140     int markidx=0;
141
142     /* Ensure space is enough for quotient and remainder.
143      * We need space for an extra limb in the remainder, because it's
144      * up-shifted (normalized) below.  */
145     rsize = nsize + 1;
146     mpi_resize( rem, rsize);
147
148     qsize = rsize - dsize;        /* qsize cannot be bigger than this.  */
149     if( qsize <= 0 ) {
150         if( num != rem ) {
151             rem->nlimbs = num->nlimbs;
152             rem->sign = num->sign;
153             MPN_COPY(rem->d, num->d, nsize);
154         }
155         if( quot ) {
156             /* This needs to follow the assignment to rem, in case the
157              * numerator and quotient are the same.  */
158             quot->nlimbs = 0;
159             quot->sign = 0;
160         }
161         return;
162     }
163
164     if( quot )
165         mpi_resize( quot, qsize);
166
167     /* Read pointers here, when reallocation is finished.  */
168     np = num->d;
169     dp = den->d;
170     rp = rem->d;
171
172     /* Optimize division by a single-limb divisor.  */
173     if( dsize == 1 ) {
174         mpi_limb_t rlimb;
175         if( quot ) {
176             qp = quot->d;
177             rlimb = _gcry_mpih_divmod_1( qp, np, nsize, dp[0] );
178             qsize -= qp[qsize - 1] == 0;
179             quot->nlimbs = qsize;
180             quot->sign = sign_quotient;
181         }
182         else
183             rlimb = _gcry_mpih_mod_1( np, nsize, dp[0] );
184         rp[0] = rlimb;
185         rsize = rlimb != 0?1:0;
186         rem->nlimbs = rsize;
187         rem->sign = sign_remainder;
188         return;
189     }
190
191
192     if( quot ) {
193         qp = quot->d;
194         /* Make sure QP and NP point to different objects.  Otherwise the
195          * numerator would be gradually overwritten by the quotient limbs.  */
196         if(qp == np) { /* Copy NP object to temporary space.  */
197             np = marker[markidx++] = mpi_alloc_limb_space(nsize,
198                                                           mpi_is_secure(quot));
199             MPN_COPY(np, qp, nsize);
200         }
201     }
202     else /* Put quotient at top of remainder. */
203         qp = rp + dsize;
204
205     count_leading_zeros( normalization_steps, dp[dsize - 1] );
206
207     /* Normalize the denominator, i.e. make its most significant bit set by
208      * shifting it NORMALIZATION_STEPS bits to the left.  Also shift the
209      * numerator the same number of steps (to keep the quotient the same!).
210      */
211     if( normalization_steps ) {
212         mpi_ptr_t tp;
213         mpi_limb_t nlimb;
214
215         /* Shift up the denominator setting the most significant bit of
216          * the most significant word.  Use temporary storage not to clobber
217          * the original contents of the denominator.  */
218         tp = marker[markidx++] = mpi_alloc_limb_space(dsize,mpi_is_secure(den));
219         _gcry_mpih_lshift( tp, dp, dsize, normalization_steps );
220         dp = tp;
221
222         /* Shift up the numerator, possibly introducing a new most
223          * significant word.  Move the shifted numerator in the remainder
224          * meanwhile.  */
225         nlimb = _gcry_mpih_lshift(rp, np, nsize, normalization_steps);
226         if( nlimb ) {
227             rp[nsize] = nlimb;
228             rsize = nsize + 1;
229         }
230         else
231             rsize = nsize;
232     }
233     else {
234         /* The denominator is already normalized, as required.  Copy it to
235          * temporary space if it overlaps with the quotient or remainder.  */
236         if( dp == rp || (quot && (dp == qp))) {
237             mpi_ptr_t tp;
238
239             tp = marker[markidx++] = mpi_alloc_limb_space(dsize, mpi_is_secure(den));
240             MPN_COPY( tp, dp, dsize );
241             dp = tp;
242         }
243
244         /* Move the numerator to the remainder.  */
245         if( rp != np )
246             MPN_COPY(rp, np, nsize);
247
248         rsize = nsize;
249     }
250
251     q_limb = _gcry_mpih_divrem( qp, 0, rp, rsize, dp, dsize );
252
253     if( quot ) {
254         qsize = rsize - dsize;
255         if(q_limb) {
256             qp[qsize] = q_limb;
257             qsize += 1;
258         }
259
260         quot->nlimbs = qsize;
261         quot->sign = sign_quotient;
262     }
263
264     rsize = dsize;
265     MPN_NORMALIZE (rp, rsize);
266
267     if( normalization_steps && rsize ) {
268         _gcry_mpih_rshift(rp, rp, rsize, normalization_steps);
269         rsize -= rp[rsize - 1] == 0?1:0;
270     }
271
272     rem->nlimbs = rsize;
273     rem->sign   = sign_remainder;
274     while( markidx )
275         mpi_free_limb_space(marker[--markidx]);
276 }
277
278 void
279 _gcry_mpi_tdiv_q_2exp( gcry_mpi_t w, gcry_mpi_t u, unsigned int count )
280 {
281     mpi_size_t usize, wsize;
282     mpi_size_t limb_cnt;
283
284     usize = u->nlimbs;
285     limb_cnt = count / BITS_PER_MPI_LIMB;
286     wsize = usize - limb_cnt;
287     if( limb_cnt >= usize )
288         w->nlimbs = 0;
289     else {
290         mpi_ptr_t wp;
291         mpi_ptr_t up;
292
293         RESIZE_IF_NEEDED( w, wsize );
294         wp = w->d;
295         up = u->d;
296
297         count %= BITS_PER_MPI_LIMB;
298         if( count ) {
299             _gcry_mpih_rshift( wp, up + limb_cnt, wsize, count );
300             wsize -= !wp[wsize - 1];
301         }
302         else {
303             MPN_COPY_INCR( wp, up + limb_cnt, wsize);
304         }
305
306         w->nlimbs = wsize;
307     }
308 }
309
310 /****************
311  * Check whether dividend is divisible by divisor
312  * (note: divisor must fit into a limb)
313  */
314 int
315 _gcry_mpi_divisible_ui(gcry_mpi_t dividend, ulong divisor )
316 {
317     return !_gcry_mpih_mod_1( dividend->d, dividend->nlimbs, divisor );
318 }
319
320
321 void
322 gcry_mpi_div (gcry_mpi_t quot, gcry_mpi_t rem, gcry_mpi_t dividend, gcry_mpi_t divisor, int round)
323 {
324   if (!round)
325     {
326       if (!rem)
327         {
328           gcry_mpi_t tmp = mpi_alloc (mpi_get_nlimbs(quot));
329           _gcry_mpi_tdiv_qr (quot, tmp, dividend, divisor);
330           mpi_free (tmp);
331         }
332       else
333         _gcry_mpi_tdiv_qr (quot, rem, dividend, divisor);
334     }
335   else if (round < 0)
336     {
337       if (!rem)
338         _gcry_mpi_fdiv_q (quot, dividend, divisor);
339       else if (!quot)
340         _gcry_mpi_fdiv_r (rem, dividend, divisor);
341       else
342         _gcry_mpi_fdiv_qr (quot, rem, dividend, divisor);
343     }
344   else
345     log_bug ("mpi rounding to ceiling not yet implemented\n");
346 }
347
348
349 void
350 gcry_mpi_mod (gcry_mpi_t rem, gcry_mpi_t dividend, gcry_mpi_t divisor)
351 {
352   _gcry_mpi_fdiv_r (rem, dividend, divisor);
353   rem->sign = 0;
354 }
355