diff options
Diffstat (limited to 'src/benchmarks/cfrac/cfrac.c')
| -rw-r--r-- | src/benchmarks/cfrac/cfrac.c | 268 |
1 files changed, 268 insertions, 0 deletions
diff --git a/src/benchmarks/cfrac/cfrac.c b/src/benchmarks/cfrac/cfrac.c new file mode 100644 index 0000000..9185d18 --- /dev/null +++ b/src/benchmarks/cfrac/cfrac.c @@ -0,0 +1,268 @@ +#include <string.h> +#include <stdio.h> +#include <math.h> /* for findk */ + +#if defined(_WIN32) +#include <windows.h> +#endif + +#include "pdefs.h" + +#ifdef __STDC__ +#include <stdlib.h> +#endif +#include "precision.h" +#include "pfactor.h" + +#ifdef __STDC__ +extern unsigned *pfactorbase(precision n, unsigned k, + unsigned *m, unsigned aborts); +extern double pomeranceLpow(double n, double alpha); +#else +extern unsigned *pfactorbase(); +extern double pomeranceLpow(); +#endif + +int verbose = 0; +int debug = 0; + +extern unsigned cfracNabort; +extern unsigned cfracTsolns; +extern unsigned cfracPsolns; +extern unsigned cfracT2solns; +extern unsigned cfracFsolns; + + +extern unsigned short primes[]; +extern unsigned primesize; + +/* + * Return the value of "f(p,d)" from Knuth's exercise 28 + */ +float pfKnuthEx28(p, d) + unsigned p; + precision d; +{ + register float res; + precision k = pUndef; + + (void) pparm(d); + if (p == 2) { + if (peven(d)) { + pset(&k, phalf(d)); + if (peven(k)) { + res = 2.0/3.0 + pfKnuthEx28(2,k)/2.0; /* eliminate powers of 2 */ + } else { /* until only one 2 left in d. */ + res = 1.0/3.0; /* independent of (the now odd) k. Wow! */ + } + } else { /* d now odd */ + pset(&k, phalf(d)); + if (podd(k)) { + res = 1.0/3.0; /* f(2,4k+3): d%8 == 3 or 7 */ + } else { + if (podd(phalf(k))) { + res = 2.0/3.0; /* f(2,8k+5): d%8 == 5 */ + } else { + res = 4.0/3.0; /* f(2,8k+1): d%8 == 1 */ + } + } + } + } else { /* PART 3: p odd, d could still be even (OK) */ + pset(&k, utop(p)); + if peq(ppowmod(d, phalf(psub(k, pone)), k), pone) { + res = (float) (p+p) / (((float) p)*p-1.0); /* beware int overflow! */ + } else { + res = 0.0; + } + } + + pdestroy(k); + pdestroy(d); + if (debug > 1) { + fprintf(stdout, "f(%u,", p); + fprintf(stdout, "d) = %9.7f\n", res); + } + return res; +} + +float cfrac_logf(unsigned p, precision n, unsigned k) +{ + register float res; + + (void) pparm(n); + +#if 0 /* old code for non-float machines; not worth the cost */ + pset(&r, utop(k)); + log2sqrtk = plogb(pipow(r, q >> 1), ptwo); + fplog2p = (f(p,pmul(r,n),q) * plogb(pipow(utop(p),q),ptwo)+(q>>1))/q; +#endif + + res = pfKnuthEx28(p, pmul(itop(k),n)) * log((double) p); + /* res -= log((double) k) * 0.5; */ + + pdestroy(n); + return res; +} + +/* + * Find the best value of k for the given n and m. + * + * Input/Output: + * n - the number to factor + * m - pointer to size of factorbase (0 = select "best" size) + * aborts - the number of early aborts + */ +unsigned findk(n, m, aborts, maxk) + precision n; + register unsigned *m; + unsigned aborts, maxk; +{ + unsigned k, bestk = 0, count, bestcount = 0, maxpm; + float sum, max = -1.0E+15; /* should be small enough */ + unsigned *p; + register unsigned i; + register unsigned short *primePtr; + + (void) pparm(n); + + for (k = 1; k < maxk; k++) { /* maxk should best be m+m? */ + if (debug) { + fputs("kN = ", stdout); + fputp(stdout, pmul(utop(k), n)); putc('\n', stdout); + } + count = *m; + p = pfactorbase(n, k, &count, aborts); + if (p == (unsigned *) 0) { + fprintf(stderr, "couldn't compute factor base in findk\n"); + exit(1); + } + + maxpm = p[count-1]; + + sum = 0.0; + primePtr = primes; + while (*primePtr <= maxpm) { + sum += cfrac_logf((unsigned) *primePtr++, n, k); + } + sum -= log((double) k) * 0.5; + if (verbose > 2) fprintf(stdout, "%u: %5.2f", k, sum); + if (debug) fprintf(stdout, " log(k)/2=%5.2f", log((double) k) * 0.5); + if (verbose > 2) { + fputs("\n", stdout); + fflush(stdout); + } + if (sum > max) { + max = sum; + bestk = k; + bestcount = count; + } +#ifndef IGNOREFREE + free(p); +#endif + } + + *m = bestcount; + pdestroy(n); + return bestk; +} + +extern char *optarg; +extern int optind; + +char *progName; + +extern int getopt(); + +int main(argc, argv) + int argc; + char *argv[]; +{ + unsigned m = 0, k = 0; + unsigned maxCount = 1<<30, count, maxk = 0; + int ch; + precision n = pUndef, f = pUndef; + unsigned aborts = 3; + unsigned *p; + double d; + + progName = *argv; + + while ((ch = getopt(argc, argv, "a:k:i:dv")) != EOF) switch (ch) { + case 'a': + aborts = atoi(optarg); + break; + case 'k': + maxk = atoi(optarg); + break; + case 'i': + maxCount = atoi(optarg); + break; + case 'd': + debug++; + break; + case 'v': + verbose++; + break; + default: +usage: fprintf(stderr, + "usage: %s [-dv] [-a aborts ] [-k maxk ] [-i maxCount ] n [[ m ] k ]\n", + progName); + return 1; + } + argc -= optind; + argv += optind; + + if (argc == 0) { + argc = 1; + static char* argvx[2] = { "17545186520507317056371138836327483792736", NULL }; + argv = argvx; + } + + if (argc < 1 || argc > 3) goto usage; + + pset(&n, atop(*argv++)); --argc; + if (argc) { m = atoi(*argv++); --argc; } + if (argc) { k = atoi(*argv++); --argc; } + + if (k == 0) { + if (maxk == 0) { + maxk = m / 2 + 5; + if (verbose) fprintf(stdout, "maxk = %u\n", maxk); + } + k = findk(n, &m, aborts, maxk); + if (verbose) { + fprintf(stdout, "k = %u\n", k); + } + } + + count = maxCount; + + pcfracInit(m, k, aborts); + + pset(&f, pcfrac(n, &count)); + count = maxCount - count; + if (verbose) { + putc('\n', stdout); + fprintf(stdout, "Iterations : %u\n", count); + fprintf(stdout, "Early Aborts : %u\n", cfracNabort); + fprintf(stdout, "Total Partials : %u\n", cfracTsolns); + fprintf(stdout, "Used Partials : %u\n", cfracT2solns); + fprintf(stdout, "Full Solutions : %u\n", cfracPsolns); + fprintf(stdout, "Factor Attempts: %u\n", cfracFsolns); + } + + if (f != pUndef) { + fputp(stdout, n); + fputs(" = ", stdout); + fputp(stdout, f); + fputs(" * ", stdout); + pdivmod(n, f, &n, pNull); + fputp(stdout, n); + putc('\n', stdout); + } + + pdestroy(f); + pdestroy(n); + + return 0; +} |