← Back to branch summary

~ubuntu-branches/debian/sid/genius/sid

« back to all changes in this revision

Viewing changes to src/mpzextra.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel Holbach
  • Date: 2006-08-21 12:57:45 UTC
  • Revision ID: james.westby@ubuntu.com-20060821125745-sl9ks8v7fq324bdf
Tags: upstream-0.7.6.1
Import upstream version 0.7.6.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/mpzextra.c
 
1
/* GENIUS Calculator
 
2
 * Copyright (C) 1997-2003 George Lebl
 
3
 *
 
4
 * Author: George Lebl
 
5
 *
 
6
 * This program is free software; you can redistribute it and/or modify
 
7
 * it under the terms of the GNU General Public License as published by
 
8
 * the Free Software Foundation; either version 2 of the License, or
 
9
 * (at your option) any later version.
 
10
 *
 
11
 * This program 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 General Public License for more details.
 
15
 *
 
16
 * You should have received a copy of the GNU General Public License
 
17
 * along with this program; if not, write to the  Free Software
 
18
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
 
19
 * USA.
 
20
 */
 
21
 
 
22
#include "config.h"
 
23
 
 
24
#include <string.h>
 
25
#include <glib.h>
 
26
 
 
27
#include "calc.h" /* for evalnode_hook and i18n stuff */
 
28
 
 
29
#include "mpzextra.h"
 
30
 
 
31
extern gboolean interrupted;
 
32
 
 
33
/* The strong pseudoprime test code copied from GMP */
 
34
 
 
35
static int
 
36
strongpseudoprimesub (mpz_srcptr n, mpz_srcptr nm1, mpz_srcptr x, mpz_ptr y,
 
37
                      mpz_srcptr q, unsigned long int k)
 
38
{
 
39
  unsigned long int i;
 
40
 
 
41
  mpz_powm (y, x, q, n);
 
42
 
 
43
  if (mpz_cmp_ui (y, 1L) == 0 || mpz_cmp (y, nm1) == 0)
 
44
    return 1;
 
45
 
 
46
  for (i = 1; i < k; i++)
 
47
    {
 
48
      mpz_powm_ui (y, y, 2L, n);
 
49
      if (mpz_cmp (y, nm1) == 0)
 
50
        return 1;
 
51
      if (mpz_cmp_ui (y, 1L) == 0)
 
52
        return 0;
 
53
    }
 
54
  return 0;
 
55
}
 
56
 
 
57
 
 
58
gboolean
 
59
mympz_strong_pseudoprime_test (mpz_srcptr n, mpz_srcptr b)
 
60
{
 
61
  mpz_t nm1, y, q;
 
62
  unsigned long int k;
 
63
  int is_prime;
 
64
  unsigned long int size_n  = mpz_size (n) * GMP_LIMB_BITS;
 
65
 
 
66
  mpz_init2 (nm1, size_n + GMP_LIMB_BITS);
 
67
  mpz_sub_ui (nm1, n, 1L);
 
68
 
 
69
  mpz_init2 (y, 2 * size_n); /* mpz_powm_ui needs excessive memory!!! */
 
70
 
 
71
  mpz_init2 (q, size_n);
 
72
 
 
73
  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
 
74
  k = mpz_scan1 (nm1, 0L);
 
75
  mpz_tdiv_q_2exp (q, nm1, k);
 
76
 
 
77
  is_prime = strongpseudoprimesub (n, nm1, b, y, q, k);
 
78
 
 
79
  mpz_clear (nm1);
 
80
  mpz_clear (y);
 
81
  mpz_clear (q);
 
82
 
 
83
  return is_prime;
 
84
}
 
85
 
 
86
/* return true only if strong pseudoprime to all of 2,3,5 and 7 */
 
87
static gboolean
 
88
mympz_strong_pseudoprime_test_2_3_5_7 (mpz_srcptr n)
 
89
{
 
90
  mpz_t nm1, y, q, b;
 
91
  unsigned long int k;
 
92
  int is_prime;
 
93
  unsigned long int size_n  = mpz_size (n) * GMP_LIMB_BITS;
 
94
 
 
95
  mpz_init2 (nm1, size_n + GMP_LIMB_BITS);
 
96
  mpz_sub_ui (nm1, n, 1L);
 
97
 
 
98
  mpz_init2 (y, 2 * size_n); /* mpz_powm_ui needs excessive memory!!! */
 
99
 
 
100
  mpz_init2 (q, size_n);
 
101
 
 
102
  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
 
103
  k = mpz_scan1 (nm1, 0L);
 
104
  mpz_tdiv_q_2exp (q, nm1, k);
 
105
 
 
106
  is_prime = 0;
 
107
  mpz_init_set_ui (b, 2);
 
108
  if (strongpseudoprimesub (n, nm1, b, y, q, k)) {
 
109
          mpz_set_ui (b, 3);
 
110
          if (strongpseudoprimesub (n, nm1, b, y, q, k)) {
 
111
                  mpz_set_ui (b, 5);
 
112
                  if (strongpseudoprimesub (n, nm1, b, y, q, k)) {
 
113
                          mpz_set_ui (b, 7);
 
114
                          if (strongpseudoprimesub (n, nm1, b, y, q, k)) {
 
115
                                  is_prime = 1;
 
116
                          }
 
117
                  }
 
118
          }
 
119
  }
 
120
 
 
121
  mpz_clear (b);
 
122
  mpz_clear (nm1);
 
123
  mpz_clear (y);
 
124
  mpz_clear (q);
 
125
 
 
126
  return is_prime;
 
127
}
 
128
 
 
129
/* assuming generalized riemman hypothesis, if we test for
 
130
   every base from 2 to 2*lg(n)^2.
 
131
ref:
 
132
   Neil Koblitz, A Course in Number Theory and Cryptography, Springer, 1987 */
 
133
gboolean
 
134
mympz_miller_rabin_test_sure (mpz_srcptr n)
 
135
{
 
136
  mpz_t nm1, y, q, b, m;
 
137
  unsigned long int k;
 
138
  int is_prime;
 
139
  unsigned long int size_n  = mpz_size (n) * GMP_LIMB_BITS;
 
140
 
 
141
  mpz_init2 (nm1, size_n + GMP_LIMB_BITS);
 
142
  mpz_sub_ui (nm1, n, 1L);
 
143
 
 
144
  mpz_init2 (y, 2 * size_n); /* mpz_powm_ui needs excessive memory!!! */
 
145
 
 
146
  mpz_init2 (q, size_n);
 
147
 
 
148
  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
 
149
  k = mpz_scan1 (nm1, 0L);
 
150
  mpz_tdiv_q_2exp (q, nm1, k);
 
151
 
 
152
  /* compute 2*lg(n)^2 */
 
153
  mpz_init_set_ui (m, mpz_sizeinbase (n, 2));
 
154
  mpz_mul (m, m, m);
 
155
  mpz_mul_ui (m, m, 2);
 
156
 
 
157
  is_prime = 1;
 
158
  mpz_init_set_ui (b, 2);
 
159
  do {
 
160
          if ( ! strongpseudoprimesub (n, nm1, b, y, q, k)) {
 
161
                  is_prime = 0;
 
162
                  break;
 
163
          }
 
164
          mpz_add_ui (b, b, 1L);
 
165
          if (evalnode_hook != NULL) {
 
166
                  (*evalnode_hook)();
 
167
          }
 
168
          if (interrupted) {
 
169
                  is_prime = 0;
 
170
                  break;
 
171
          }
 
172
  } while (mpz_cmp (b, m) <= 0);
 
173
 
 
174
  mpz_clear (b);
 
175
  mpz_clear (m);
 
176
  mpz_clear (nm1);
 
177
  mpz_clear (y);
 
178
  mpz_clear (q);
 
179
 
 
180
  return is_prime;
 
181
}
 
182
 
 
183
/* we will really do 22+4 since we always test 2,3,5,7 as well,
 
184
   so the probability that a wrong result is returned is
 
185
   about 1/(4^26) which is about 8.9*10^-16.  Knuth says
 
186
   25 passes are reasonable (according to gmp source),
 
187
   and Knuth is never wrong. */
 
188
long int mympz_is_prime_miller_rabin_reps = 22;
 
189
 
 
190
gboolean
 
191
mympz_is_prime (mpz_srcptr n, int miller_rabin_reps)
 
192
{
 
193
        int ret;
 
194
        static mpz_t test;
 
195
        static gboolean inited_test = FALSE; 
 
196
        int sgn;
 
197
 
 
198
        sgn = mpz_sgn (n);
 
199
        if (sgn == 0)
 
200
                return 0;
 
201
        else if (sgn < 0) {
 
202
                mpz_t nn;
 
203
                mpz_init (nn);
 
204
                mpz_neg (nn, n);
 
205
                ret = mympz_is_prime (nn, miller_rabin_reps);
 
206
                mpz_clear (nn);
 
207
                return ret;
 
208
        }
 
209
 
 
210
        if (miller_rabin_reps < 0)
 
211
                miller_rabin_reps = mympz_is_prime_miller_rabin_reps;
 
212
 
 
213
        /* Use the probab prime for trial divisions and stuff
 
214
           and do one strong pseudoprime test for good meassure */
 
215
        ret = mpz_probab_prime_p (n, 1);
 
216
        /* if we are sure */
 
217
        if (ret == 2)
 
218
                return 1;
 
219
        else if (ret == 0)
 
220
                return 0;
 
221
 
 
222
        if ( ! mympz_strong_pseudoprime_test_2_3_5_7 (n))
 
223
                return 0;
 
224
 
 
225
        if ( ! inited_test) {
 
226
                /* set test to 25*10^9 */
 
227
                mpz_init_set_ui (test, 10);
 
228
                mpz_pow_ui (test, test, 9);
 
229
                mpz_mul_ui (test, test, 25);
 
230
        }
 
231
 
 
232
        /* if n < 25*10^9, we are now sure this
 
233
           is a prime since the only n less than that
 
234
           that is a composite and strong pseudoprime to 2,3,5,7
 
235
           is n = 3215031751
 
236
        ref:
 
237
           Neil Koblitz, A Course in Number Theory and Cryptography,
 
238
           Springer, 1987 */
 
239
        if (mpz_cmp (n, test) <= 0)
 
240
                return 1;
 
241
 
 
242
        if (evalnode_hook != NULL)
 
243
                (*evalnode_hook)();
 
244
        if (interrupted)
 
245
                return 0;
 
246
 
 
247
        return mpz_millerrabin (n, miller_rabin_reps-1);
 
248
}
 
249
 
 
250
static void
 
251
append_factor (GArray *fact, mpz_srcptr num)
 
252
{
 
253
        GelFactor f;
 
254
        int i;
 
255
        /* FIXME: implement faster search, this is sorted */
 
256
        /* We start at 1 since the 0 entry is always -1 or 1 */
 
257
        for (i = 1; i < fact->len; i++) {
 
258
                int cmp = mpz_cmp (g_array_index (fact, GelFactor, i).num, num);
 
259
                if (cmp > 0) {
 
260
                        /* must add */
 
261
                        mpz_init_set (f.num, num);
 
262
                        f.exp = 1;
 
263
                        g_array_insert_val (fact, i, f);
 
264
                        return;
 
265
                } else if (cmp == 0) {
 
266
                        /* already in the factorization */
 
267
                        g_array_index (fact, GelFactor, i).exp++;
 
268
                        return;
 
269
                }
 
270
        }
 
271
        /* not found, must add */
 
272
        mpz_init_set (f.num, num);
 
273
        f.exp = 1;
 
274
        g_array_append_val (fact, f);
 
275
}
 
276
 
 
277
static void
 
278
append_factor_uint (GArray *fact, unsigned long num)
 
279
{
 
280
        mpz_t znum;
 
281
        mpz_init_set_ui (znum, num);
 
282
        append_factor (fact, znum);
 
283
        mpz_clear (znum);
 
284
}
 
285
 
 
286
/* Pollard-Rho factorization code snarfed from the gmp examples */
 
287
 
 
288
static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
 
289
 
 
290
static void
 
291
factor_using_division (GArray *fact, mpz_t t, unsigned int limit)
 
292
{
 
293
  mpz_t q, r;
 
294
  unsigned long int f;
 
295
  int ai;
 
296
  unsigned *addv = add;
 
297
  unsigned int failures;
 
298
 
 
299
  mpz_init (q);
 
300
  mpz_init (r);
 
301
 
 
302
  f = mpz_scan1 (t, 0);
 
303
  mpz_div_2exp (t, t, f);
 
304
  while (f)
 
305
    {
 
306
      append_factor_uint (fact, 2);
 
307
      --f;
 
308
    }
 
309
 
 
310
  for (;;)
 
311
    {
 
312
      mpz_tdiv_qr_ui (q, r, t, 3);
 
313
      if (mpz_cmp_ui (r, 0) != 0)
 
314
        break;
 
315
      mpz_set (t, q);
 
316
      append_factor_uint (fact, 3);
 
317
    }
 
318
 
 
319
  for (;;)
 
320
    {
 
321
      mpz_tdiv_qr_ui (q, r, t, 5);
 
322
      if (mpz_cmp_ui (r, 0) != 0)
 
323
        break;
 
324
      mpz_set (t, q);
 
325
      append_factor_uint (fact, 5);
 
326
    }
 
327
 
 
328
  failures = 0;
 
329
  f = 7;
 
330
  ai = 0;
 
331
  while (mpz_cmp_ui (t, 1) != 0)
 
332
    {
 
333
      mpz_tdiv_qr_ui (q, r, t, f);
 
334
      if (mpz_cmp_ui (r, 0) != 0)
 
335
        {
 
336
          f += addv[ai];
 
337
          if (mpz_cmp_ui (q, f) < 0)
 
338
            break;
 
339
          ai = (ai + 1) & 7;
 
340
          failures++;
 
341
          if (failures > limit)
 
342
            break;
 
343
        }
 
344
      else
 
345
        {
 
346
          mpz_swap (t, q);
 
347
          append_factor_uint (fact, f);
 
348
          failures = 0;
 
349
        }
 
350
    }
 
351
 
 
352
  mpz_clear (q);
 
353
  mpz_clear (r);
 
354
}
 
355
 
 
356
static void
 
357
factor_using_pollard_rho (GArray *fact, mpz_t n, int a_int)
 
358
{
 
359
  mpz_t x, x1, y, P;
 
360
  mpz_t a;
 
361
  mpz_t g;
 
362
  mpz_t t1, t2;
 
363
  int k, l, c, i;
 
364
 
 
365
  mpz_init (g);
 
366
  mpz_init (t1);
 
367
  mpz_init (t2);
 
368
 
 
369
  mpz_init_set_si (a, a_int);
 
370
  mpz_init_set_si (y, 2);
 
371
  mpz_init_set_si (x, 2);
 
372
  mpz_init_set_si (x1, 2);
 
373
  k = 1;
 
374
  l = 1;
 
375
  mpz_init_set_ui (P, 1);
 
376
  c = 0;
 
377
 
 
378
  while (mpz_cmp_ui (n, 1) != 0)
 
379
    {
 
380
S2:
 
381
      if (evalnode_hook != NULL) {
 
382
              static int i = 0;
 
383
              if G_UNLIKELY ((i++ & RUN_HOOK_EVERY_MASK) == RUN_HOOK_EVERY_MASK) {
 
384
                      (*evalnode_hook)();
 
385
                      i = 0;
 
386
              }
 
387
      }
 
388
      if (interrupted) {
 
389
              mpz_set_ui (n, 1);
 
390
              continue;
 
391
      }
 
392
      mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 
393
      mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
 
394
      c++;
 
395
      if (c == 20)
 
396
        {
 
397
          c = 0;
 
398
          mpz_gcd (g, P, n);
 
399
          if (mpz_cmp_ui (g, 1) != 0)
 
400
            goto S4;
 
401
          mpz_set (y, x);
 
402
        }
 
403
      k--;
 
404
      if (k != 0)
 
405
        goto S2;
 
406
 
 
407
      mpz_gcd (g, P, n);
 
408
      if (mpz_cmp_ui (g, 1) != 0)
 
409
        goto S4;
 
410
 
 
411
      mpz_set (x1, x);
 
412
      k = l;
 
413
      l = 2 * l;
 
414
      for (i = 0; i < k; i++)
 
415
        {
 
416
          mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 
417
        }
 
418
      mpz_set (y, x);
 
419
      c = 0;
 
420
      goto S2;
 
421
S4:
 
422
      do
 
423
        {
 
424
          mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
 
425
          mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
 
426
        }
 
427
      while (mpz_cmp_ui (g, 1) == 0);
 
428
 
 
429
      if (!mympz_is_prime (g, -1))
 
430
        {
 
431
          do
 
432
            {
 
433
              mp_limb_t a_limb;
 
434
              mpn_random (&a_limb, (mp_size_t) 1);
 
435
              a_int = (int) a_limb;
 
436
            }
 
437
          while (a_int == -2 || a_int == 0);
 
438
 
 
439
          factor_using_pollard_rho (fact, g, a_int);
 
440
          break;
 
441
        }
 
442
      else
 
443
        {
 
444
          append_factor (fact, g);
 
445
        }
 
446
      mpz_div (n, n, g);
 
447
      mpz_mod (x, x, n);
 
448
      mpz_mod (x1, x1, n);
 
449
      mpz_mod (y, y, n);
 
450
      if (mympz_is_prime (n, -1))
 
451
        {
 
452
          append_factor (fact, n);
 
453
          break;
 
454
        }
 
455
    }
 
456
 
 
457
  mpz_clear (g);
 
458
  mpz_clear (P);
 
459
  mpz_clear (t2);
 
460
  mpz_clear (t1);
 
461
  mpz_clear (a);
 
462
  mpz_clear (x1);
 
463
  mpz_clear (x);
 
464
  mpz_clear (y);
 
465
}
 
466
 
 
467
static void
 
468
factor_number (GArray *fact, mpz_t t)
 
469
{
 
470
  unsigned int division_limit;
 
471
 
 
472
  /* Set the trial division limit according the size of t.  */
 
473
  division_limit = mpz_sizeinbase (t, 2);
 
474
  if (division_limit > 1000)
 
475
    division_limit = 1000 * 1000;
 
476
  else
 
477
    division_limit = division_limit * division_limit;
 
478
 
 
479
  factor_using_division (fact, t, division_limit);
 
480
 
 
481
  if (mpz_cmp_ui (t, 1) != 0)
 
482
    {
 
483
      if (mympz_is_prime (t, -1))
 
484
        append_factor (fact, t);
 
485
      else
 
486
        factor_using_pollard_rho (fact, t, 1);
 
487
    }
 
488
}
 
489
 
 
490
GArray *
 
491
mympz_pollard_rho_factorize (mpz_srcptr t)
 
492
{
 
493
        GArray *fact;
 
494
        mpz_t n;
 
495
        int sgn = mpz_sgn (t);
 
496
 
 
497
        fact = g_array_new (FALSE /* zero_terminated */,
 
498
                            FALSE /* clear */,
 
499
                            sizeof (GelFactor) /* element_size */);
 
500
 
 
501
        mpz_init_set (n, t);
 
502
 
 
503
        if (sgn == 0) {
 
504
                append_factor_uint (fact, 0);
 
505
        } if (sgn < 0) {
 
506
                GelFactor f;
 
507
                /* for negative numbers, add -1 to the factors
 
508
                   and factorize -t */
 
509
                mpz_neg (n, n);
 
510
 
 
511
                mpz_init_set_si (f.num, -1);
 
512
                f.exp = 1;
 
513
                g_array_append_val (fact, f);
 
514
 
 
515
                factor_number (fact, n);
 
516
        } else {
 
517
                append_factor_uint (fact, 1);
 
518
                /* just factor the number */
 
519
                factor_number (fact, n);
 
520
        }
 
521
 
 
522
        mpz_clear (n);
 
523
 
 
524
        if (interrupted) {
 
525
                mympz_factorization_free (fact);
 
526
                return NULL;
 
527
        }
 
528
 
 
529
        return fact;
 
530
 
 
531
}
 
532
 
 
533
void
 
534
mympz_factorization_free (GArray *fact)
 
535
{
 
536
        int i;
 
537
        for (i = 0; i < fact->len; i++) {
 
538
                mpz_clear (g_array_index (fact, GelFactor, i).num);
 
539
        }
 
540
        g_array_free (fact, TRUE /* free_segment */);
 
541
}