This is an automated email from the git hooks/post-receive script. ppm-guest pushed a commit to annotated tag v0.06 in repository libmath-prime-util-perl.
commit 8f4796dfbec39a962b0ae093007b3ce097930e36 Author: Dana Jacobsen <d...@acm.org> Date: Fri Jun 8 23:05:11 2012 -0600 Add perfect square discriminator --- factor.c | 58 +++++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 49 insertions(+), 9 deletions(-) diff --git a/factor.c b/factor.c index dcfec63..656389a 100644 --- a/factor.c +++ b/factor.c @@ -178,6 +178,54 @@ int trial_factor(UV n, UV *factors, UV maxtrial) /* n^2 + a mod m */ #define sqraddmod(n, a, m) addmod(sqrmod(n,m), a,m) +/* Return 0 if n is not a perfect square. Set sqrtn to int(sqrt(n)) if so. + * + * A simplistic solution (but not unhelpful) is: + * + * return ( ((m&2)!= 0) || ((m&7)==5) || ((m&11) == 8) ) ? 0 : 1; + * + * The following Bloom filter cascade works very well indeed. Read all + * about it here: http://mersenneforum.org/showpost.php?p=110896 + */ +static int is_perfect_square(UV n, UV* sqrtn) +{ + UV m, lm; + m = n & 127; + if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a) return 0; + /* 82% of non-squares rejected here */ + +#if 0 + /* The big deal with this technique is that you do two total operations, + * one cheap (the & 127 above), one expensive (the modulo below) on n. + * The rest of the operations are 32-bit operations. This is a huge win + * if n is multiprecision. + * However, in this file we're doing native precision sqrt, so it just + * isn't expensive enough to justify this second filter set. + */ + lm = n % UVCONST(63*25*11*17*19*23*31); + m = lm % 63; + if ((m*0x3d491df7) & (m*0xc824a9f9) & 0x10f14008) return 0; + m = lm % 25; + if ((m*0x1929fc1b) & (m*0x4c9ea3b2) & 0x51001005) return 0; + m = 0xd10d829a*(lm%31); + if (m & (m+0x672a5354) & 0x21025115) return 0; + m = lm % 23; + if ((m*0x7bd28629) & (m*0xe7180889) & 0xf8300) return 0; + m = lm % 19; + if ((m*0x1b8bead3) & (m*0x4d75a124) & 0x4280082b) return 0; + m = lm % 17; + if ((m*0x6736f323) & (m*0x9b1d499) & 0xc0000300) return 0; + m = lm % 11; + if ((m*0xabf1a3a7) & (m*0x2612bf93) & 0x45854000) return 0; + /* 99.92% of non-squares are rejected now */ +#endif + m = sqrt(n); + if (n != (m*m)) + return 0; + + if (sqrtn != 0) *sqrtn = m; + return 1; +} /* Miller-Rabin probabilistic primality test * Returns 1 if probably prime relative to the bases, 0 if composite. @@ -360,15 +408,7 @@ int holf_factor(UV n, UV *factors, UV rounds) s = sqrt(n*i); /* TODO: overflow here */ if ( (s*s) != (n*i) ) s++; m = sqrmod(s, n); - /* Cheaper would be: - * if (m is probably not a perfect sequare) continue; - * f = sqrt(m); - * if (f*f == m) { yay } - */ - if ( ((m&2)!= 0) || ((m&7)==5) || ((m&11) == 8) ) - continue; - f = sqrt(m); /* expensive */ - if ( (f*f) == m ) { + if (is_perfect_square(m, &f)) { f = gcd_ui( (s>f) ? s-f : f-s, n); /* This should always succeed, but with overflow concerns.... */ if ((f == 1) || (f == n)) -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-perl/packages/libmath-prime-util-perl.git _______________________________________________ Pkg-perl-cvs-commits mailing list Pkg-perl-cvs-commits@lists.alioth.debian.org http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/pkg-perl-cvs-commits