Ciao,
Il Mar, 20 Novembre 2018 10:34 pm, Adrien Prost-Boucle ha scritto:
> Mostly "use" existing mpz functions to reach a proof of concept,
> ... and I only finished the Lucas-Lehmer test for Mersenne numbers xD
> As you also suggested, a demo program would be a more appropriate.
To exploit the very special form of Mersenne's numbers, mpn are more
powerful I think.
Look at the attached source. It uses the same logic I implemented for the
function mpn_llriter in
https://gmplib.org/repo/gmp/file/17707/mpn/generic/strongfibo.c#l69
Ĝis,
m
--
http://bodrato.it/papers/
/* Tests primality of 2^exp-1, using the Lucas-Lehmer algorithm.
Copyright 2018 Marco Bodrato
This file is intendet for the GNU MP Library.
This is free software; you can redistribute it and/or modify it under
the terms of the GNU Affero General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at
your option) any later version.
This software is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero
General Public License for more details.
You should have received copies of the GNU Affero General Public
License with the software.
If not, see https://www.gnu.org/licenses/. */
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include "gmp.h"
#define MERSENNE_BITMAP 0x402008ee
/* (2 << (2/3) | 2 << (3/3) | 2 << (5/3) | 2 << (7/3) | */
/* 2 <<(13/3) | 2 <<(17/3) | 2 <<(19/3) | 2 <<(31/3) | */
/* 2 <<(61/3) | 2 <<(89/3) ) */
static int
lucas_lehmer (unsigned long int exp)
{
mpz_t e;
int ret;
mpz_init_set_ui (e, exp);
ret = mpz_probab_prime_p (e, 25) != 0;
mpz_clear (e);
if (ret != 0)
{
mp_size_t n;
mp_ptr value, square;
unsigned long int i;
mp_limb_t cy;
int c;
if (exp < 30 * 3)
return 2 & (MERSENNE_BITMAP >> exp / 3);
assert (exp > GMP_NUMB_BITS);
n = exp / GMP_NUMB_BITS + 1;
c = exp % GMP_NUMB_BITS;
assert (c > 0 && c < GMP_NUMB_BITS);
value = calloc (n, sizeof(mp_limb_t));
square = malloc (2 * n * sizeof(mp_limb_t));
#if GMP_NUMB_BITS > 30
*value = 1416317954; /* 37634^2 - 2 */
i = exp - 7;
#else
#if GMP_NUMB_BITS > 15
*value = 37634; /* 194^2 - 2 */
i = exp - 6;
#else
#if GMP_NUMB_BITS > 7
*value = 194; /* 14^2 - 2 */
i = exp - 5;
#else
#if GMP_NUMB_BITS > 3
*value = 14; /* 4^2 - 2 */
i = exp - 4;
#else
#if GMP_NUMB_BITS > 2
*value = 4;
i = exp - 3;
#else
#error Not implemented.
#endif /* 2 */
#endif /* 3 */
#endif /* 7 */
#endif /* 15 */
#endif /* 30 */
do
{
mp_limb_t cy;
mpn_sqr (square, value, n);
/* Next shift will not give a carry. */
cy = mpn_lshift (value, square + n, n, GMP_NUMB_BITS - c);
assert (cy == 0);
*value |= square[n - 1] >> c;
square[n - 1] &= GMP_NUMB_MASK >> GMP_NUMB_BITS - c;
/* Next add will not give a carry. */
cy = mpn_add_n (value, value, square, n);
assert (cy == 0);
cy = value[n - 1] >> c;
value[n - 1] &= GMP_NUMB_MASK >> GMP_NUMB_BITS - c;
/* Next add_1 will not give a carry. */
cy = mpn_add_1 (value, value, n, cy);
assert (cy == 0);
/* Subtract 2. */
cy = *value;
*value = (cy - 2) & GMP_NUMB_MASK;
/* Detect borrow, or fixed points. */
if (cy < 5) {
assert (n > 1);
if (mpn_zero_p (value + 1, n - 1))
{
/* Small value, will loop on -1 or 2. */
ret = 0;
break;
}
else if (cy < 2)
/* Take care of the borrow. */
mpn_sub_1 (value + 1, value + 1, n - 1, 1);
}
}
while (--i);
free (square);
if (ret != 0)
{
assert (exp & 1 == 1);
/* Check if value is +/- sqrt(2) i.e. +/- 2^((exp+1)/2) */
i = exp / 2 + 1;
cy = value [i / GMP_NUMB_BITS] ^= (mp_limb_t) 1 << i % GMP_NUMB_BITS;
/* The highest limb is not full, scan0 will not go out of bounds. */
assert (value[n - 1] != GMP_NUMB_MASK);
if ((cy == GMP_NUMB_MAX) && (mpn_scan0 (value, 0) == exp) ||
(cy == 0) && mpn_zero_p (value, n))
ret = 2;
else
ret = 0;
}
free (value);
}
return ret;
}
int
main (int argc, char **argv)
{
long exp;
exp = 0;
if (argc > 1)
{
char *end;
exp = strtol (argv[1], &end, 0);
if (*end || exp <= 0)
{
fprintf (stderr, "Invalid exponent: %s.\n", argv[1]);
return 1;
}
printf ("2^%lu-1 is ", exp);
if (lucas_lehmer (exp))
printf ("prime.\n");
else
printf ("composite.\n");
return 0;
}
fprintf (stderr, "WARNING: this program will loop forever.\n");
do {
if (lucas_lehmer (++exp))
printf ("%lu,\n", exp);
} while (1);
}_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
