找到几个数论算法的源代码

无心人 · 发表于 2008-11-30 21:44:28

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

http://www.mindspring.com/~pate/course/chap09.html APRCL.c 还少一个大数的库似乎在lip.h里定义的

/*
Author: Pate Williams (c) 1997
Lucas-Lehmer primality test.
See "Implementation of a New Primality Test",
Mathematics of Computation, Volume 48, Number 177,
January 1987, pages 103 - 121 by H. Cohen and
A. K. Lenstra.
*/
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"
#define BOUND 1000000l
#define t 25200l
int trial_division(long bound, verylong zn,
long *lm, long *lmSize,
long *lp, long *lpSize,
verylong *zfm,
verylong *zfp,
verylong *znp,
verylong *znm,
verylong *zrm,
verylong *zrp)
/* returns 0 if number is composite 1 if prime
2 if proved prime by trial division */
{
int factored = 0, flag = 0;
long b, p, r;
verylong za = 0, zs = 0, zn1 = 0;
*lmSize = *lpSize = 0;
zsadd(zn, - 1l, znm);
zsadd(zn, + 1l, znp);
zcopy(*znm, zrm);
zcopy(*znp, zrp);
zcopy(*znp, &zn1);
zpstart2();
p = zpnext();
zsqrt(zn, &zs, &za);
if (zscompare(zs, bound) > 0)
b = bound;
else {
flag = 1;
b = ztoint(zs);
}
while (!factored && p <= b) {
r = zsmod(*znp, p);
factored = r == 1;
if (!factored) {
if (r == 0) {
do {
zsdiv(*zrp, p, &za);
zcopy(za, zrp);
} while (zsmod(*zrp, p) == 0);
lp[*lpSize] = p;
*lpSize = *lpSize + 1;
}
if (r == 2) {
do {
zsdiv(*zrm, p, &za);
zcopy(za, zrm);
} while (zsmod(*zrm, p) == 0);
lm[*lmSize] = p;
*lmSize = *lmSize + 1;
}
}
p = zpnext();
}
zdiv(*znm, *zrm, zfm, &za);
zdiv(*znp, *zrp, zfp, &za);
zfree(&za);
zfree(&zs);
zfree(&zn1);
if (!factored && flag && p > b) factored = 2;
else factored = !factored;
return factored;
}
void ring_mul(long a, long r, long u,
verylong zn,
verylong *zx, verylong *zy,
verylong *zz)
{
verylong za = 0, zb = 0, zc = 0;
verylong zp0 = 0, zp1 = 0, zs0 = 0, zs1 = 0;
zmulmod(zx[0], zy[0], zn, &zp0);
zmulmod(zx[1], zy[1], zn, &zp1);
zaddmod(zx[0], zx[1], zn, &zs0);
zaddmod(zy[0], zy[1], zn, &zs1);
if (r == 1) {
zsmulmod(zp1, a, zn, &za);
zaddmod(zp0, za, zn, &zz[0]);
zmulmod(zs0, zs1, zn, &za);
zsubmod(za, zp0, zn, &zb);
zsubmod(zb, zp1, zn, &zz[1]);
}
else if (r == 3) {
zaddmod(zp0, zp1, zn, &zz[0]);
zsmulmod(zp1, u - 1, zn, &za);
zmulmod(zs0, zs1, zn, &zb);
zaddmod(zb, za, zn, &zc);
zsubmod(zc, zp0, zn, &zz[1]);
}
zfree(&za);
zfree(&zb);
zfree(&zc);
zfree(&zp0);
zfree(&zp1);
zfree(&zs0);
zfree(&zs1);
}
void ring_sqr(long u, verylong zn,
verylong *zx, verylong *zz)
{
verylong za = 0, zb = 0, zs = 0;
zsmulmod(zx[1], u, zn, &za);
zsmulmod(zx[0], 2l, zn, &zb);
zaddmod(za, zb, zn, &zs);
zmulmod(zx[0], zs, zn, &za);
zsadd(za, - 1l, &zb);
zmod(zb, zn, &zz[0]);
zmulmod(zx[1], zs, zn, &zz[1]);
zfree(&za);
zfree(&zb);
zfree(&zs);
}
void ring_pow(long a, long u,
verylong ze,
verylong zn,
verylong *zx,
verylong *zz)
/* calculates z = x ^ e in the ring
(Z / nZ)[T]/(T ^ 2 - u * T - 1) */
{
long i, r = zsmod(zn, 4l);
verylong za = 0, zb = 0;
verylong zA[2], zS[2], zT[2];
zA[0] = zA[1] = zS[0] = zS[1] = zT[0] = zT[1] = 0;
zcopy(ze, &za);
zzero(&zA[1]);
zone(&zA[0]);
zcopy(zx[0], &zS[0]);
zcopy(zx[1], &zS[1]);
while (zscompare(za, 0l) != 0) {
if (zodd(za)) {
ring_mul(a, r, u, zn, zA, zS, zz);
zcopy(zz[0], &zA[0]);
zcopy(zz[1], &zA[1]);
}
zrshift(za, 1l, &zb);
zcopy(zb, &za);
if (zscompare(za, 0l) != 0) {
ring_sqr(u, zn, zS, zT);
zcopy(zT[0], &zS[0]);
zcopy(zT[1], &zS[1]);
}
}
zfree(&za);
zfree(&zb);
for (i = 0; i < 2; i++) {
zfree(&zA[i]);
zfree(&zS[i]);
zfree(&zT[i]);
}
}
int test_for_nm(int *flag,
verylong zn,
verylong znm,
verylong zp,
verylong *zprod,
verylong **zbeta)
/* returns - 1 if test fails, 0 if n is composite,
1 otherwise */
{
int found = 0, value;
long i, l, p, q, x;
verylong za = 0, zb = 0, zo = 0, zx = 0;
zpstart2();
for (i = 0; !found && i < 50; i++) {
x = zpnext();
zdiv(znm, zp, &za, &zb);
zintoz(x, &zx);
zexpmod(zx, za, zn, &zb);
found = zscompare(zb, 1l) != 0;
}
if (found) {
zexpmod(zx, znm, zn, &za);
value = zscompare(za, 1l) == 0;
if (value) {
zcopy(*zprod, &zo);
zsadd(zb, - 1l, &za);
zmulmod(*zprod, za, zn, &zb);
zcopy(zb, zprod);
value = zscompare(*zprod, 0l) == 0;
if (value) {
zgcd(zo, zn, &za);
printf("factor = "); zwriteln(za);
value = 0;
}
else if (zprobprime(zp, 5l)) {
zpstart();
l = 1;
p = ztoint(zp);
do {
q = pow(p, l);
if (t % q == 0) {
flag[q] = 1;
for (i = 0; i < q - 1; i++) {
zsmul(znm, i, &za);
zsdiv(za, q, &zb);
zexpmod(zx, zb, zn, &zbeta[q][i]);
}
}
l++;
} while (q < t);
value = 1;
}
else value = 1;
}
}
else value = - 1;
if (!value) printf("composite by n - 1 test\n");
zfree(&za);
zfree(&zb);
zfree(&zo);
zfree(&zx);
return value;
}
int test_for_np(long a, long u,
verylong zn,
verylong znp,
verylong zp,
verylong *zprod)
/* returns
- 1 if test fails, 0 if n is composite,
+ 1 otherwise */
{
int factored = 0, found = 0, value;
long d, m;
verylong za = 0, zd = 0, ze = 0, zi = 0;
verylong zx[2], zz[2];
zx[0] = zx[1] = zz[0] = zz[1] = 0;
zdiv(znp, zp, &ze, &za);
for (m = 1; !factored && !found && m <= 50; m++) {
d = m * (m + u) - a;
zintoz(d, &zd);
zinvmod(zd, zn, &zi);
factored = zscompare(zi, 0l) == 0;
if (!factored) {
zintoz(m * m + a, &za);
zmulmod(za, zi, zn, &zx[0]);
zintoz((2 * m + u) * a, &za);
zmulmod(za, zi, zn, &zx[1]);
ring_pow(a, u, ze, zn, zx, zz);
found = zscompare(zz[0], 1l) != 0 ||
zscompare(zz[1], 0l) != 0;
}
}
if (!factored && found) {
ring_pow(a, u, znp, zn, zx, zz);
factored = zscompare(zz[0], 1l) == 1 &&
zscompare(zz[1], 0l) == 0;
if (!factored) {
if (zscompare(zx[0], 0l) != 0)
zmulmod(*zprod, zx[0], zn, &za);
else
zmulmod(*zprod, zx[1], zn, &za);
factored = zscompare(za, 0l) == 0;
if (factored) {
zgcd(*zprod, zn, &za);
printf("factor: ");
zwriteln(za);
}
else zcopy(za, zprod);
if (factored) value = 0; else value = 1;
}
}
else if (factored) value = 0;
else value = - 1;
zfree(&za);
zfree(&zd);
zfree(&ze);
zfree(&zi);
zfree(&zx[0]);
zfree(&zx[1]);
zfree(&zz[0]);
zfree(&zz[1]);
if (!value) printf("composite by n + 1 test\n");
return value;
}
void norm(long a, long u, verylong zn,
verylong *znorm, verylong *zx)
{
verylong za = 0, zb = 0, zc = 0;
zmulmod(zx[0], zx[0], zn, &za);
zmulmod(zx[0], zx[1], zn, &zb);
zsmulmod(zb, u, zn, &zc);
zaddmod(za, zc, zn, &zb);
zmulmod(zx[1], zx[1], zn, &za);
zsmulmod(za, - a, zn, &zc);
zaddmod(zb, zc, zn, znorm);
zfree(&za);
zfree(&zb);
zfree(&zc);
}
int Lucas_Lehmer(int *flag, long bound, long lmSize,
long *lm, long lpSize, long *lp,
verylong zn,
verylong zfm, verylong zfp,
verylong znm, verylong znp,
verylong zrm, verylong zrp,
verylong **zbeta)
/* returns
- 4 could not find suitable a,
- 3 could not find suitable u,
- 2 if n - 1 test fails,
- 1 n + 1 test fails, 0 if n is composite,
+ 1 if passed */
{
int f4_1, factored = 0, found = 0, test = 1;
long a, d, i, l, m, p, r, u;
verylong za = 0, zd = 0, ze = 0, zi = 0, zm = 0;
verylong zp = 0, zu = 0, zprod = 0, zx[2], zz[2];
zx[0] = zx[1] = zz[0] = zz[1] = 0;
zone(&zprod);
for (i = 0; test && i < lmSize; i++) {
p = lm[i];
zintoz(p, &zp);
test = test_for_nm(flag, zn, znm, zp, &zprod, zbeta);
}
if (test == 1) {
if (zcompare(zfm, zfp) > 0) zcopy(zfm, &za);
else zcopy(zfp, &za);
zmul(za, zfm, &zd);
zmul(zd, zfp, &ze);
zintoz(bound, &za);
zsexp(za, 3l, &zm);
zmul(ze, zm, &za);
f4_1 = zcompare(zn, za) < 0;
if (f4_1 && zscompare(zrm, 1l) != 0)
test = test_for_nm(flag, zn, znm, zrm, &zprod, zbeta);
}
if (test == 1) {
r = zsmod(zn, 4l);
if (r == 1) {
zrshift(znm, 1l, &ze);
zpstart2();
for (i = 1; !found && i <= 50; i++) {
a = zpnext();
zintoz(a, &za);
zexpmod(za, ze, zn, &zm);
found = zcompare(zm, znm) == 0;
}
if (!found) test = - 4;
else {
l = 1;
do {
p = pow(2, l);
zintoz(a, &za);
if (t % p == 0 && flag[p]) {
for (i = 0; i < p; i++) {
zsmul(znm, i, &zd);
zsdiv(zd, p, &ze);
zexpmod(za, ze, zn, &zbeta[p][i]);
}
}
l++;
} while (p < sqrt(t));
}
}
else {
a = 1;
for (i = 1; !found && i <= 50; i++) {
zintoz(i * i + 4l, &zu);
found = zjacobi(zu, zn) == - 1;
if (found) u = i;
}
if (!found) test = - 3;
else {
found = 0;
for (m = 1; !factored && !found && m <= 50; m++) {
d = m * (m + u) - a;
zintoz(d, &zd);
zinvmod(zd, zn, &zi);
factored = zscompare(zi, 0l) == 0;
if (!factored) {
found = 1;
zintoz(m * m + a, &za);
zmulmod(za, zi, zn, &zx[0]);
zintoz((2 * m + u) * a, &za);
zmulmod(za, zi, zn, &zx[1]);
}
}
if (!factored && !found) test = - 2;
else if (found) {
zone(&zx[0]);
zintoz(u, &zx[1]);
zrshift(znp, 1l, &ze);
ring_pow(a, u, ze, zn, zx, zz);
test = zcompare(zz[0], znm) == 0 &&
zscompare(zz[1], 0l) == 0;
}
else test = 0;
}
}
}
else if (test == - 1) test = - 2;
if (test == 1) {
for (i = 0; test && i < lpSize; i++) {
zintoz(lp[i], &zp);
test = test_for_np(a, u, zn, znp, zp, &zprod);
}
if (test == 1) {
if (f4_1) {
if (zscompare(zrp, 1l) == 0) test = 2;
else {
test = test_for_np(a, u, zn, znp, zrp, &zprod);
if (test) test = 2;
}
}
}
}
if (test == 1) {
zgcd(zprod, zn, &zd);
test = zscompare(zd, 1l) == 0;
}
zfree(&za);
zfree(&zd);
zfree(&ze);
zfree(&zi);
zfree(&zm);
zfree(&zp);
zfree(&zu);
zfree(&zprod);
zfree(&zx[0]);
zfree(&zx[1]);
zfree(&zz[0]);
zfree(&zz[1]);
return test;
}
int main(void)
{
char answer[256];
int flag[256], td;
long cols, rows;
long lm[256], lp[256], lmSize, lpSize;
long i, j, l, n = 0, p, q, r[32], tt;
verylong zn = 0;
verylong zfm = 0, zfp = 0, znm = 0, znp = 0;
verylong zrm = 0, zrp = 0;
verylong **zbeta;
zpstart2();
do {
p = zpnext();
if (t % p == 0) {
r[n++] = p;
l = 2;
do {
q = pow(p, l);
if (t % q == 0) r[n++] = q;
l++;
} while (q < sqrt(t));
}
} while (p < sqrt(t));
for (i = 0; i < n - 1; i++)
for (j = i + 1; j < n; j++)
if (r[i] > r[j])
tt = r[i], r[i] = r[j], r[j] = tt;
cols = rows = r[n - 1] + 1;
zbeta = calloc(rows, sizeof(verylong *));
if (!zbeta) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
exit(1);
}
for (i = 0; i < rows; i++) {
zbeta[i] = calloc(cols, sizeof(verylong));
if (!zbeta[i]) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
exit(1);
}
}
do {
printf("enter the number to be tested below:\n");
zread(&zn);
zwriteln(zn);
if (zprobprime(zn, 5l)) {
td = trial_division(BOUND, zn, lm, &lmSize,
lp, &lpSize, &zfm, &zfp, &znp,
&znm, &zrm, &zrp);
if (td == 1) {
for (i = 0; i < 256; i++) flag[i] = 0;
switch (Lucas_Lehmer(flag, BOUND, lmSize,
lm, lpSize, lp, zn,
zfm, zfp, znm, znp,
zrm, zrp, zbeta)) {
case - 4 :
printf("failed could not find suitable a\n");
break;
case - 3 :
printf("failed could not find suitable u\n");
break;
case - 2 :
printf("Lucas-Lehmer n - 1 test failed\n");
break;
case - 1 :
printf("Lucas-Lehmer n + 1 test failed\n");
break;
case + 0 :
break;
case + 1 :
printf("number passed Lucas-Lehmer test\n");
break;
case + 2 :
printf("number proven prime by Lucas-Lehmer test\n");
break;
}
}
else if (!td)
printf("number composite by trial division\n");
else
printf("number proven prime by trial division\n");
}
else printf("number composite by Miller-Rabin\n");
printf("another number (n or y) ? ");
scanf("%s", answer);
} while (tolower(answer[0]) == 'y');
zfree(&zn);
zfree(&zfm);
zfree(&zfp);
zfree(&znm);
zfree(&znp);
zfree(&zrm);
zfree(&zrp);
for (i = 0; i < rows; i++) {
for (j = 0; j < cols; j++)
zfree(&zbeta[i][j]);
free(zbeta[i]);
}
free(zbeta);
return 0;
}

复制代码

无心人 · 发表于 2008-11-30 21:47:47

*
Author: Pate Williams (c) 1997
Aldleman-Pomerance-Rumely-Cohen-Lenstra primality
test. See "Implementation of a New Primality Test",
Mathematics of Computation, Volume 48, Number 177,
January 1987, pages 103 - 121 by H. Cohen and
A. K. Lenstra. Table generation.
*/
#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"
#define t 25200l
struct power {long x, p;};
struct long_factor {long expon, prime;};
struct very_factor {int expon; verylong prime;};
void long_trial_division(long bound, verylong zn,
struct long_factor *f,
long *count)
{
int one = 0, pri = 0, value;
long e, p;
verylong za = 0, zb = 0;
*count = 0;
zcopy(zn, &za);
zpstart2();
do {
p = zpnext();
if (zsmod(za, p) == 0) {
e = 0;
do {
e++;
zsdiv(za, p, &zb);
zcopy(zb, &za);
} while (zsmod(za, p) == 0);
f[*count].expon = e;
f[*count].prime = p;
*count = *count + 1;
one = zscompare(za, 1l) == 0;
if (!one) pri = zprobprime(za, 5l);
}
} while (!one && !pri && p <= bound);
if (pri) {
f[*count].expon = 1;
f[*count].prime = ztoint(za);
*count = *count + 1;
value = zscompare(za, LONG_MAX) > 0;
}
if (value || p > bound) {
fprintf(stderr, "fatal error\ncan't factor number");
exit(1);
}
zfree(&za);
zfree(&zb);
}
void compute_e_of_t(verylong zt, long *Q,
verylong *ze, long *Q_count)
{
int found;
long e, i, lf_count, p, q;
struct long_factor lf[16];
verylong zq = 0, zr = 0;
*Q_count = 0;
long_trial_division(1000000l, zt, lf, &lf_count);
printf("t = "); zwriteln(zt);
printf("prime factorization of t\n");
for (i = 0; i < lf_count; i++) {
printf("%2ld", lf[i].prime);
if (lf[i].expon == 1) printf("\n");
else printf(" ^ %ld\n", lf[i].expon);
}
zintoz(2l, ze);
zpstart2();
do {
p = zpnext();
q = p - 1;
if (zsmod(zt, q) == 0) {
Q[*Q_count] = p;
*Q_count = *Q_count + 1;
e = 1;
found = 0;
for (i = 0; !found && i < lf_count; i++) {
if (p == lf[i].prime) {
found = 1;
e += lf[i].expon;
}
}
zintoz(p, &zq);
zsexp(zq, e, &zr);
zmul(*ze, zr, &zq);
zcopy(zq, ze);
}
} while (zscompare(zt, q) >= 0);
zfree(&zq);
zfree(&zr);
}
void very_trial_division(long bound, verylong zn,
struct very_factor *f,
long *count)
{
int e, one = 0, pri = 0;
long p;
verylong za = 0, zb = 0;
*count = 0;
zcopy(zn, &za);
zpstart2();
do {
p = zpnext();
if (zsmod(za, p) == 0) {
e = 0;
do {
e++;
zsdiv(za, p, &zb);
zcopy(zb, &za);
} while (zsmod(za, p) == 0);
f[*count].expon = e;
zintoz(p, &f[*count].prime);
*count = *count + 1;
one = zscompare(za, 1l) == 0;
if (!one) pri = zprobprime(za, 5l);
}
} while (!one && !pri && p <= bound);
if (pri) {
f[*count].expon = 1;
zcopy(za, &f[*count].prime);
*count = *count + 1;
}
if (p > bound) {
fprintf(stderr, "fatal error\ncan't factor number");
exit(1);
}
zfree(&za);
zfree(&zb);
}
void primitive_root(verylong zp, verylong *za)
{
long i, k;
struct very_factor f[16];
verylong zb = 0, ze = 0, zq = 0, zr = 0, zp1 = 0;
for (i = 0; i < 16; i++) f[i].prime = 0;
zsadd(zp, - 1l, &zp1);
very_trial_division(1000000l, zp1, f, &k);
zone(za);
L2:
i = 0;
zsadd(*za, 1l, &zb);
zcopy(zb, za);
L3:
zdiv(zp1, f[i].prime, &zq, &zr);
zexpmod(*za, zq, zp, &ze);
if (zscompare(ze, 1l) == 0) goto L2;
if (++i < k) goto L3;
zfree(&zb);
zfree(&ze);
zfree(&zq);
zfree(&zr);
zfree(&zp1);
for (i = 0; i < 16; i++) zfree(&f[i].prime);
}
int fcmp(const void *g1, const void *g2)
{
struct power *p1 = (struct power *) g1;
struct power *p2 = (struct power *) g2;
if (p1->p < p2->p) return - 1;
if (p1->p > p2->p) return + 1;
return 0;
}
void build_f_table(long q, verylong zg, verylong zq, FILE *file)
{
long f, q1 = q - 1, q2 = q - 2, x;
long *h = calloc(q, sizeof(long));
struct power *g, key, *ptr;
verylong zr = 0, zs = 0, zt = 0;
g = calloc(q, sizeof(struct power));
if (!g || !h) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
exit(1);
}
zone(&zr);
g[0].x = 0;
g[0].p = 1;
for (x = 1; x <= q2; x++) {
zsexpmod(zg, x, zq, &zs);
zsubmod(zr, zs, zq, &zt);
g[x].x = x;
g[x].p = ztoint(zs);
h[x] = ztoint(zt);
}
qsort(g, q1, sizeof(struct power), fcmp);
fwrite(&q, sizeof(q), 1, file);
for (x = 1; x <= q2; x++) {
key.x = x;
key.p = h[x];
ptr = bsearch(&key, g, q1, sizeof(struct power), fcmp);
if (!ptr) {
fprintf(stderr, "fatal error\ncan't find logarith\n");
exit(1);
}
f = ptr->x;
fwrite(&f, sizeof(f), 1, file);
}
free(g);
free(h);
zfree(&zr);
zfree(&zs);
zfree(&zt);
}
void generate_f_table(long q_count, long *q)
{
long i, i0, p;
FILE *file = fopen("f_table.dat", "w+b");
verylong zg = 0, zp = 0;
if (!file) {
fprintf(stderr, "fatal error\ncan't open file\n");
exit(1);
}
if (q[0] == 2) i0 = 1; else i0 = 0;
for (i = i0; i < q_count; i++) {
printf("\b\b%ld", i);
p = q[i];
zintoz(p, &zp);
primitive_root(zp, &zg);
build_f_table(p, zg, zp, file);
}
zfree(&zg);
zfree(&zp);
fclose(file);
}
void generate_j_table(long q_count, long q0)
{
long A, B, *a, *f, i, k, l, m, p, pk;
long pk1, pk3, q, q1, r, r0, x;
FILE *f_file = fopen("f_table.dat", "r+b");
FILE *j_file = fopen("j_table.dat", "w+b");
if (!f_file || !j_file) {
fprintf(stderr, "fatal error\ncan't open file(s)\n");
exit(1);
}
if (q0 == 2) r0 = 1; else r0 = 0;
for (r = r0; r < q_count; r++) {
printf("\b\b%ld", r);
fread(&q, sizeof(q), 1, f_file);
f = calloc(q, sizeof(long));
if (!f) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
exit(1);
}
for (i = 1; i <= q - 2; i++)
fread(&f[i], sizeof(f[i]), 1, f_file);
q1 = q - 1;
zpstart2();
do {
p = zpnext();
if (q1 % p == 0) {
k = 0;
do {
k++;
q1 /= p;
} while (q1 % p == 0);
}
pk = pow(p, k);
pk1 = pow(p, k - 1);
m = (p - 1) * pk1;
a = calloc(m, sizeof(long));
if (!a) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
exit(1);
}
if (pk != 2) {
for (i = 0; i < m; i++) a[i] = 0;
A = B = 1;
for (x = 1; x <= q - 2; x++) {
l = (A * x + B * f[x]) % pk;
if (l < m) a[l]++;
else
for (i = 1; i < p; i++)
a[((l - i) * pk1) % l]--;
}
fwrite(&m, sizeof(m), 1, j_file);
for (i = 0; i < m; i++)
fwrite(&a[i], sizeof(a[i]), 1, j_file);
}
else if (p == 2 && k >= 3) {
for (i = 0; i < m; i++) a[i] = 0;
A = 2, B = 1;
for (x = 1; x <= q - 2; x++) {
l = (A * x + B * f[x]) % pk;
if (l < m) a[l]++;
else
for (i = 1; i < p; i++)
a[((l - i) * pk1) % l]--;
}
fwrite(&m, sizeof(m), 1, j_file);
for (i = 0; i < m; i++)
fwrite(&a[i], sizeof(a[i]), 1, j_file);
for (i = 0; i < m; i++) a[i] = 0;
pk3 = pow(2, k - 3);
A = 3 * pk3, B = pk3;
for (x = 1; x <= q - 2; x++) {
l = (A * x + B * f[x]) % pk;
if (l < m) a[l]++;
else
for (i = 1; i < p; i++)
a[((l - i) * pk1) % l]--;
}
fwrite(&m, sizeof(m), 1, j_file);
for (i = 0; i < m; i++)
fwrite(&a[i], sizeof(a[i]), 1, j_file);
}
free(a);
} while (q1 != 1);
}
free(f);
fclose(f_file);
fclose(j_file);
}
int main(void)
{
long q_count, q[64];
verylong ze = 0, zt = 0;
zintoz(t, &zt);
compute_e_of_t(zt, q, &ze, &q_count);
printf("e(t) = "); zwriteln(ze);
printf("log10(e(t)) = %lf\n", zln(ze) / log(10.0));
printf("number of q primes = %ld\n", q_count);
generate_f_table(q_count, q);
printf("\ngenerated f_table\n");
generate_j_table(q_count, q[0]);
printf("\ngenerated j_table\n");
zfree(&ze);
zfree(&zt);
return 0;
}

复制代码

无心人 · 发表于 2008-11-30 21:49:49

/*
Author: Pate Williams (c) 1997
Algorithm 10.3.3 (Lenstra's ECM). See "A Course
in Computational Algebraic Number Theory" by
Henri Cohen page 488.
*/
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include "lip.h"
#define CURVES 50l
struct point {verylong zx, zy, zz;};
int partial_addition(verylong za,
verylong zn,
struct point P,
struct point Q,
struct point *R,
verylong *zd)
/* returns 0 if sum is found or 1 if divisor is found */
{
int value = 0;
verylong zb = 0, zc = 0, zl = 0, zs = 0, zt = 0;
verylong zx = 0, zy = 0, zy2 = 0;
if (zcompare(P.zx, Q.zx) == 0 &&
zscompare(P.zy, 0l) == 0 &&
zscompare(Q.zy, 0l) == 0) {
zzero(&R->zx);
zone(&R->zy);
return 0;
}
zsub(zn, Q.zy, &zb);
if (zcompare(P.zx, Q.zx) == 0 &&
zcompare(P.zy, zb) == 0) {
zzero(&R->zx);
zone(&R->zy);
zfree(&zb);
return 0;
}
if (zscompare(P.zx, 0l) == 0 &&
zscompare(P.zy, 1l) == 0 &&
zscompare(P.zz, 0l) == 0) {
/* O + Q = Q */
zcopy(Q.zx, &R->zx);
zcopy(Q.zy, &R->zy);
zcopy(Q.zz, &R->zz);
value = 0;
}
else if (zscompare(Q.zx, 0l) == 0 &&
zscompare(Q.zy, 1l) == 0 &&
zscompare(Q.zz, 0l) == 0) {
/* P + O = P */
zcopy(P.zx, &R->zx);
zcopy(P.zy, &R->zy);
zcopy(P.zz, &R->zz);
value = 0;
}
else {
/* P != O and Q != O */
zcopy(Q.zy, &zy2);
znegate(&zy2);
if (zcompare(P.zx, Q.zx) == 0 &&
zcompare(P.zy, zy2) == 0) {
zzero(&R->zx);
zone(&R->zy);
zzero(&R->zz);
}
else {
if (zcompare(P.zx, Q.zx) != 0) {
zsubmod(P.zx, Q.zx, zn, &zx);
zexteucl(zx, &zs, zn, &zt, zd);
if (zscompare(*zd, 1l) != 0) goto L1;
zsubmod(P.zy, Q.zy, zn, &zy);
zmulmod(zs, zy, zn, &zl);
}
else {
zaddmod(P.zy, Q.zy, zn, &zy);
zexteucl(zy, &zs, zn, &zt, zd);
if (zscompare(*zd, 1l) != 0) goto L1;
zmulmod(P.zx, P.zx, zn, &zb);
zsmulmod(zb, 3l, zn, &zc);
zaddmod(zc, za, zn, &zb);
zmulmod(zs, zb, zn, &zl);
}
zmulmod(zl, zl, zn, &zb);
zaddmod(P.zx, Q.zx, zn, &zc);
zsubmod(zb, zc, zn, &zx);
zcopy(zx, &R->zx);
zsubmod(zx, P.zx, zn, &zb);
zmulmod(zl, zb, zn, &zc);
zaddmod(zc, P.zy, zn, &zy);
znegate(&zy);
zcopy(zy, &R->zy);
zone(&R->zz);
goto L2;
L1:
value = 1;
L2:
}
}
zfree(&zb);
zfree(&zc);
zfree(&zl);
zfree(&zs);
zfree(&zt);
zfree(&zx);
zfree(&zy);
return value;
}
int multiply(long k,
verylong za,
verylong zn,
struct point P,
struct point *R,
verylong *zd)
{
int value = 0;
struct point A, S, T;
A.zx = A.zy = A.zz = S.zx = S.zy = S.zz = 0;
T.zx = T.zy = T.zz = 0;
zzero(&R->zx);
zone(&R->zy);
zzero(&R->zz);
zcopy(P.zx, &S.zx);
zcopy(P.zy, &S.zy);
zcopy(P.zz, &S.zz);
while (!value && k != 0) {
if (k & 1) {
value = partial_addition(za, zn, *R, S, &A, zd);
zcopy(A.zx, &R->zx);
zcopy(A.zy, &R->zy);
zcopy(A.zz, &R->zz);
}
k >>= 1;
if (!value && k != 0) {
value = partial_addition(za, zn, S, S, &T, zd);
zcopy(T.zx, &S.zx);
zcopy(T.zy, &S.zy);
zcopy(T.zz, &S.zz);
}
}
if (zscompare(R->zy, 0l) < 0) {
zadd(R->zy, zn, &A.zy);
zcopy(A.zy, &R->zy);
}
zfree(&A.zx);
zfree(&A.zy);
zfree(&A.zz);
zfree(&S.zx);
zfree(&S.zy);
zfree(&S.zz);
zfree(&T.zx);
zfree(&T.zy);
zfree(&T.zz);
return value;
}
int LenstrasECM(verylong *zN, verylong *zg)
{
int expon = 0, found;
long B = 2000000l, i, j, k = 0, l, *p, q, q1;
struct point x, y;
verylong za[CURVES], zb = 0, zd = 0;
for (i = 0; i < CURVES; i++) za[i] = 0;
x.zx = x.zy = x.zz = y.zx = y.zy = y.zz = 0;
zpstart2();
do { q = zpnext(); k++; } while (q < B);
p = calloc(k, sizeof(long));
if (!p) {
fprintf(stderr, "fatal error\ninsufficient memory\n");
exit(1);
}
zpstart2();
for (i = 0; i < k; i++) p[i] = zpnext();
for (i = 0; i < CURVES; i++)
zrandomb(*zN, &za[i]);
L2:
zone(&x.zx);
zone(&x.zy);
zzero(&x.zz);
i = - 1;
L3:
i++;
if (i == k) {
for (i = 0; i < CURVES; i++)
zrandomb(*zN, &za[i]);
goto L2;
}
q = p[i];
q1 = q;
l = B / q;
while (q1 <= l) q1 *= q;
found = 0;
for (j = 0; !found && j < CURVES; j++)
found = multiply(q1, za[j], *zN, x, &y, &zd);
if (!found) goto L3;
zcopy(y.zx, &x.zx);
zcopy(y.zy, &x.zy);
zcopy(y.zz, &x.zz);
zgcd(zd, *zN, zg);
if (zcompare(*zg, *zN) == 0) {
for (j = 0; j < CURVES; j++)
zrandomb(*zN, &za[j]);
goto L2;
}
do {
zdiv(*zN, *zg, &zb, &zd);
zcopy(zb, zN);
zmod(*zN, *zg, &zd);
expon++;
} while (zscompare(zd, 0l) == 0);
for (i = 0; i < CURVES; i++) zfree(&za[i]);
zfree(&zb);
zfree(&zd);
zfree(&x.zx);
zfree(&x.zy);
zfree(&x.zz);
zfree(&y.zx);
zfree(&y.zy);
zfree(&y.zz);
return expon;
}
int main(void)
{
int cant, expon, one, pri;
char answer[256];
verylong zN = 0, zd = 0, zg = 0;
do {
printf("enter the number to be factored below:\n");
zread(&zN);
if (zprobprime(zN, 5l))
printf("number is prime\n");
else {
printf("number is composite\n");
zintoz(6l, &zg);
zgcd(zN, zg, &zd);
if (zscompare(zd, 1l) != 0)
printf("number can't be factored using Lenstra's ECM\n");
else {
printf("factors:\n");
cant = 0, pri = 0;
do {
expon = LenstrasECM(&zN, &zg);
zwrite(zg);
if (expon == 1) printf("\n");
else printf(" ^ %d\n", expon);
one = zscompare(zN, 1l) == 0;
if (!one) pri = zprobprime(zN, 5l);
zintoz(6l, &zg);
zgcd(zN, zg, &zd);
if (zscompare(zd, 1l) != 0) {
printf("number can't be factored using Lenstra's ECM\n");
cant = 1;
}
} while (!cant && !one && !pri);
if (pri) zwriteln(zN);
if (cant) {
zwriteln(zN);
printf("last factor is composite\n");
}
}
}
printf("another number (n or y)? ");
scanf("%s", answer);
} while (tolower(answer[0]) == 'y');
zfree(&zN);
zfree(&zd);
zfree(&zg);
return 0;
}

复制代码

无心人 · 发表于 2008-11-30 21:54:49

需要freelip库 ftp://ftp.ox.ac.uk/pub/math/freelip/ 要休息了明天想办法下载下来这个库就补全了

xiugakei · 发表于 2008-12-1 11:56:44

斑竹果然是强人

无心人 · 发表于 2008-12-1 17:14:08

http://cr.yp.to/arith.html 几个算法，没仔细看

无心人 · 发表于 2008-12-1 17:20:29

http://topic.csdn.net/t/20011210/22/414003.html

freelip_1.1.tar.gz (197.02 KB, 下载次数: 18)

winxos · 发表于 2009-4-18 18:19:46

原帖由 无心人 于 2008-12-1 17:20 发表 http://topic.csdn.net/t/20011210/22/414003.html 602

我一直很诧异，为什么各位经常可以搞到这么多我看来这么离奇的地方的奇怪的东西？

gxqcn · 发表于 2009-4-18 18:56:58

说是无心实乃有心人。

只是呼吸 · 发表于 2010-2-13 12:20:28

我只能写VB的，对C还不会，就像上面的那两个几百行的程序，我把它复制下来，在C++6.0上却通不过编译器。这两个程序在编译时都指出是#include "lip.h"这一行出了问题，通不过。而我又不知道这是为什么。烦啊。

账号		自动登录	找回密码
密码			欢迎注册

[分享] 找到几个数论算法的源代码

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。