cheng/mpir
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added a test for the sqrt code and made them work all the way up to

Browse Source
This commit is contained in:
wbhart 2009-12-08 01:54:07 +00:00
parent 4671fbc394
commit c3cf3b5716
3 changed files with 97 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
gmp-impl.h
View File

@ -784,6 +784,8 @@ void __gmp_default_free _PROTO ((void *, size_t));
int is_likely_prime_BPSW(mp_limb_t n);
mp_limb_t n_sqrt(mp_limb_t r);
void __gmpz_aorsmul_1 _PROTO ((REGPARM_3_1 (mpz_ptr w, mpz_srcptr u, mp_limb_t v, mp_size_t sub))) REGPARM_ATTR(1);
#define mpz_aorsmul_1(w,u,v,sub) __gmpz_aorsmul_1 (REGPARM_3_1 (w, u, v, sub))

58
mpz/likely_prime_p.c
View File

@ -1,7 +1,7 @@
/*
Copyright 2009 Jason Moxham
Copyright (C) 2008 Peter Shrimpton
Copyright (C) 2008, 2009 William Hart
Copyright (C) 2008 Peter Shrimpton
Copyright (C) 2008, 2009 William Hart
This file is part of the MPIR Library.
@ -52,16 +52,36 @@ mp_limb_t n_sqrt(mp_limb_t r)
mp_limb_t l;
} temp;
mp_limb_t bits32 = (r & GMP_LIMB_HIGHBIT);
mp_limb_t r2;
// algorithm can't handle 32 bits
if (bits32)
{
r2 = r;
r >>= 2;
}
temp.f = (float) r;
temp.l = (CNST_LIMB(0xbe6ec85e) - temp.l)>>1; // estimate of 1/sqrt(y)
x = temp.f;
z = (float) r*0.5;
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
res = is + ((is+1)*(is+1) <= r);
return res - (res*res > r);
res = is + ((is+1)*(is+1) <= r);
if (!bits32) return res - (res*res > r);
else
{
mp_limb_t sq;
res = res - (res*res > r);
res <<= 1;
sq = res*res;
res = res - ((sq > r2) || ((sq ^ r2) & GMP_LIMB_HIGHBIT));
sq = (res + 1)*(res + 1);
res = res + ((sq <= r2) && !((sq ^ r2) & GMP_LIMB_HIGHBIT));
return res;
}
#else
double x, z;
@ -70,18 +90,38 @@ mp_limb_t n_sqrt(mp_limb_t r)
mp_limb_t l;
} temp;
mp_limb_t bits64 = (r & GMP_LIMB_HIGHBIT);
mp_limb_t r2;
// algorithm can't handle 64 bits
if (bits64)
{
r2 = r;
r >>= 2;
}
temp.f = (double) r;
temp.l = (CNST_LIMB(0xbfcdd90a00000000) - temp.l)>>1; // estimate of 1/sqrt(y)
x = temp.f;
z = (double) r*0.5;
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
x = (1.5*x) - (x*x)*(x*z);
is = (mp_limb_t) (x*(double) r);
res = is + ((is+1)*(is+1) <= r);
return res - (res*res > r);
is = (mp_limb_t) (x*(double) r);
res = is + ((is+1)*(is+1) <= r);
if (!bits64) return res - (res*res > r);
else
{
mp_limb_t sq;
res = res - (res*res > r);
res <<= 1;
sq = res*res;
res = res - ((sq > r2) || ((sq ^ r2) & GMP_LIMB_HIGHBIT));
sq = (res + 1)*(res + 1);
res = res + ((sq <= r2) && !((sq ^ r2) & GMP_LIMB_HIGHBIT));
return res;
}
#endif
}

47
tests/mpz/t-likely_prime_p.c
View File

@ -26,6 +26,50 @@ MA 02110-1301, USA. */
#include "gmp-impl.h"
#include "tests.h"
void
check_sqrt (void)
{
gmp_randstate_t rands;
int i;
mpz_t x;
unsigned long bits;
mp_limb_t p, m, m2, s;
mpz_init (x);
gmp_randinit_default(rands);
for (i = 0; i < 2000000; i++)
{
do
{
mpz_set_ui(x, GMP_NUMB_BITS/2 - 1);
mpz_urandomm(x, rands, x);
bits = mpz_get_ui(x) + 2;
mpz_rrandomb(x, rands, bits);
p = mpz_getlimbn(x, 0);
m = p*p;
mpz_set_ui(x, MIN(p, 1000));
mpz_urandomm(x, rands, x);
m2 = m + mpz_get_ui(x) - MIN(p, 1000)/2;
} while ((m2 == m) || ((mp_limb_signed_t) (m2 ^ m) < (mp_limb_signed_t) 0) || (m2 == 0) || (m2 == 1));
s = n_sqrt(m2);
if (((m2 < m) && (s != p - 1)) || ((m2 > m) && (s != p)))
{
printf ("mpz_likely_prime_p\n");
printf ("n_sqrt is broken\n");
#if defined( _MSC_VER ) && defined( _WIN64 )
printf ("%llu is returned as n_sqrt(%llu)\n", s, m2);
#else
printf ("%lu is returned as n_sqrt(%lu)\n", s, m2);
#endif
abort();
}
}
mpz_clear (x);
gmp_randclear(rands);
}
void
check_rand (void)
{
@ -38,7 +82,7 @@ check_rand (void)
mpz_init (x);
gmp_randinit_default(rands);
for (i = 0; i < 200000; i++)
for (i = 0; i < 2000; i++)
{
do
{
@ -141,6 +185,7 @@ main (void)
{
tests_start ();
check_sqrt ();
check_rand ();
check_small ();