I am in the progress of translating (porting) my J. M. Pollard’s algorithm “Factoring with Cubic Integers” C# application to a Free LIP based vanilla C Windows 32-bit console application. The first phase of the method is to generate two factor bases namely a rational prime factor base and an algebraic integer prime factor base. I have included some preliminary results from this fast-moving computer programming task. I generated 2012 algebraic integer primes in about a minute and thirty seconds.
Back on Thursday, January 10, 2019, at 02:40 AM I started implementing J. M. Pollard’s factoring with cubic integers method and a classical variant of Shor’s quantum computer factoring algorithm. I have probably implemented at least three versions of Pollard’s magnificent work which led to the special number field sieve fast factoring method. I am working on translating my C# code of these two algorithms to Arjen K. Lenstra’s Free Large Integer Package in the C computer programming language. Appended to this electronic missive are factorizations of two small 19-decimal digit numbers 2^63-1 and 2^63+1. Pollard’s method requires that the exponent be divisible by three. Thus, we write the numbers as 2^ (21 * 3) plus or minus 1. You can factor numbers with different bases and small addends and subtrahends as long as the exponent is zero modulo 3. You will notice that Shor’s algorithm is very fast on these small numbers.
This C code uses singly linked lists as the key data structure. The factoring methods are trial division, Brent’s modification of the well-known Pollard rho method, Cohen-Pollard p – 1 algorithm first stage, Lehman’s method, Cohen’s algorithm for the Lenstra’s elliptic curve method, and the Free LIP elliptic curve method. I include a couple of sample factorizations by the Free LIP built-in Lenstra elliptic curve method. These include the seventh Fermat number 2^128+1 and the eighth Fermat number 2^256+1.
/*
Author: Pate Williams (c) 1997
Algorithm 8.4.1 (Lehman). See "A Course in
Computational Algebraic Number Theory" by
Henri Cohen page 425.
Author: Pate Williams (c) 1997 - 2022
Algorithm 8.5.2 (Pollard rho). See "A Course in
Computational Algebraic Number Theory" by Henri
Cohen page 429.
Author: Pate Williams (c) 1997 - 2022
Algorithm 8.8.2 (p - 1 First Stage). See "A Course
in Computational Algebraic Number Theory" by Henri
Cohen page 439.
Author: Pate Williams (c) 1997 - 2022
Algorithm 10.3.3 (Lenstra's ECM). See "A Course
in Computational Algebraic Number Theory" by
Henri Cohen page 488.
*/
#include "LIP_data_structures.h"
#include <crtdbg.h>
#ifndef CLK_TCK
#define CLK_TCK CLOCKS_PER_SEC
#endif
typedef struct myThreadData
{
int zacount;
verylong zn;
PFACTORNODE current;
PFACTORNODE first;
PFACTORNODE(*tdFunction)(
verylong*,
PFACTORNODE,
int*);
} THREADDATA, *PTHREADDATA;
DWORD WINAPI TdThread(LPVOID lpParam)
{
clock_t time0 = clock();
PTHREADDATA td = (PTHREADDATA)lpParam;
PFACTORNODE head = td->first;
td->current = td->tdFunction(
&td->zn, td->first, &td->zacount);
td->first = td->current;
PFACTORNODE last = find_last(td->first);
quick_sort(td->first, last);
if (td->first != NULL)
{
while (td->first != NULL)
{
if (zprobprime(td->first->factor.zfactor, 20))
{
zwrite(td->first->factor.zfactor);
if (td->first->factor.exponent > 1)
printf(" ^ %d ", td->first->factor.exponent);
printf("\n");
}
td->first = td->first->next;
}
}
PFACTORNODE ptr = head;
while (ptr != NULL)
{
ptr = delete_factor_node(head);
head = ptr;
}
double time = (clock() - time0) / (double)CLK_TCK;
printf("total time required: %f seconds\n", time);
return TRUE;
}
int main()
{
int i = 0, zacount = 0;
verylong zN = 0, zf = 0, zn = 0, zo = 0;
verylong zaddend = 0, zbase = 0, ztemp = 0;
verylong zq = 0, zr = 0, zs = 0;
PFACTORNODE current = NULL, first = NULL;
while (1)
{
char optionStr[256], numStr[256] = { 0 }, option, *bPtr, *ePtr;
int base, expon, addend, iExpon = 0, positive;
printf("Menu\n");
printf("1 Trial Division\n");
printf("2 Brent-Pollard rho Method\n");
printf("3 Cohen-Pollard p - 1 Method\n");
printf("4 Lehman's Method\n");
printf("5 Cohen Lenstra's Elliptic Curve Method\n");
printf("6 FreeLIP Lenstra's Elliptic Curve Method\n");
printf("7 Exit\n");
printf("Enter choice (1 - 7): ");
scanf_s("%s", optionStr, 256);
option = optionStr[0];
zone(&zo);
if (option == '7')
break;
if (option < '1' || option > '6')
{
printf("Unknown choice, please reenter choice\n");
continue;
}
printf("enter the number below (0 to quit):\n");
scanf_s("%s", numStr, 256);
bPtr = strchr(numStr, '^');
if (bPtr != NULL)
{
*bPtr = '\0';
base = atoi(numStr);
ePtr = strchr(bPtr + 1, '+');
if (ePtr != NULL)
positive = 1;
else {
ePtr = strchr(bPtr + 1, '-');
if (ePtr != NULL)
positive = 0;
}
if (ePtr != NULL)
*ePtr = '\0';
expon = atoi(bPtr + 1);
if (ePtr != NULL) {
addend = atoi(ePtr + 1);
if (positive == 0)
addend = -addend;
}
else
addend = 0;
zintoz(base, &zbase);
zintoz(addend, &zaddend);
zsexp(zbase, expon, &ztemp);
zadd(ztemp, zaddend, &zN);
numStr[0] = '\0';
}
else
zstrtoz(numStr, &zN);
zcopy(zN, &zn);
if (zmcomposite(zn, 20) == 0)
{
printf("Number is prime\n");
continue;
}
first = insert_first_factor_node(zo, 0);
zacount = 1;
DWORD dwMilliseconds;
PTHREADDATA pThreadData = (PTHREADDATA)
malloc(sizeof(THREADDATA));
if (pThreadData == NULL)
exit(0);
memset(pThreadData, 0, sizeof(THREADDATA));
pThreadData->current = NULL;
pThreadData->first = first;
pThreadData->zacount = 1;
zcopy(zn, &pThreadData->zn);
printf("Wait milliseconds = ");
scanf_s("%ld", &dwMilliseconds);
if (option == '1')
{
pThreadData->tdFunction = trial_division;
DWORD threadId;
HANDLE hMutex = CreateMutex(
NULL, FALSE, NULL);
if (hMutex == NULL)
exit(0);
HANDLE thread = CreateThread(
NULL,
0,
(LPTHREAD_START_ROUTINE)TdThread,
pThreadData,
0,
&threadId);
WaitForSingleObject(thread, dwMilliseconds);
CloseHandle(thread);
CloseHandle(hMutex);
}
else if (option == '2')
{
pThreadData->tdFunction = Brent_Pollard_Method;
DWORD threadId;
HANDLE hMutex = CreateMutex(
NULL, FALSE, NULL);
if (hMutex == NULL)
exit(0);
HANDLE thread = CreateThread(
NULL,
0,
(LPTHREAD_START_ROUTINE)TdThread,
pThreadData,
0,
&threadId);
if (thread == NULL)
exit(0);
WaitForSingleObject(thread, dwMilliseconds);
CloseHandle(thread);
CloseHandle(hMutex);
}
else if (option == '3')
{
pThreadData->tdFunction = Cohen_Pollard_Method;
DWORD threadId;
HANDLE hMutex = CreateMutex(
NULL, FALSE, NULL);
if (hMutex == NULL)
exit(0);
HANDLE thread = CreateThread(
NULL,
0,
(LPTHREAD_START_ROUTINE)TdThread,
pThreadData,
0,
&threadId);
WaitForSingleObject(thread, dwMilliseconds);
CloseHandle(thread);
CloseHandle(hMutex);
}
else if (option == '4')
{
pThreadData->tdFunction = Lehman;
DWORD threadId;
HANDLE hMutex = CreateMutex(
NULL, FALSE, NULL);
if (hMutex == NULL)
exit(0);
HANDLE thread = CreateThread(
NULL,
0,
(LPTHREAD_START_ROUTINE)TdThread,
pThreadData,
0,
&threadId);
WaitForSingleObject(thread, dwMilliseconds);
CloseHandle(thread);
CloseHandle(hMutex);
}
else if (option == '5')
{
first = CohenEllipticCurveMethod(
numStr, base, expon, addend);
}
else if (option == '6')
{
pThreadData->tdFunction = LenstraEllipticCurveMethod;
DWORD threadId;
HANDLE hMutex = CreateMutex(
NULL, FALSE, NULL);
if (hMutex == NULL)
exit(0);
HANDLE thread = CreateThread(
NULL,
0,
(LPTHREAD_START_ROUTINE)TdThread,
pThreadData,
0,
&threadId);
WaitForSingleObject(thread, dwMilliseconds);
CloseHandle(thread);
CloseHandle(hMutex);
}
zfree(&pThreadData->zn);
free(pThreadData);
}
zfree(&zaddend);
zfree(&zbase);
zfree(&ztemp);
zfree(&zN);
zfree(&zf);
zfree(&zn);
zfree(&zo);
zfree(&zq);
zfree(&zr);
zfree(&zs);
return 1;
}
#pragma once
#include "..\LIPFactoring\lip.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#define my_B 1000000l
#define my_NUMBER_PRIMES 78498l
#define my_BOUND 100000000l
/* https://www.tutorialspoint.com/data_structures_algorithms/linked_list_program_in_c.htm */
/* https://www.geeksforgeeks.org/quicksort-on-singly-linked-list */
typedef struct my_factor
{
long exponent; // secondary storage key
verylong zfactor; // primary storage key
} FACTOR, *PFACTOR;
typedef struct my_factor_node
{
FACTOR factor;
struct my_factor_node* next;
} FACTORNODE, *PFACTORNODE;
/* basic singly linked list functions
** the c means complete comparison including
** the exponents of the factor or factor_node
** s denotes a shallow comparison */
int factor_z_compare_c(FACTOR lt, FACTOR rt);
int factor_z_compare_s(FACTOR lt, FACTOR rt);
int factor_node_z_compare_c(FACTORNODE lt, FACTORNODE rt);
int factor_node_z_compare_s(FACTORNODE lt, FACTORNODE rt);
PFACTORNODE insert_first_factor_node(verylong zf, int exponent);
PFACTORNODE insert_factor_node(PFACTORNODE current, FACTOR factor);
PFACTORNODE delete_factor_node(PFACTORNODE first);
PFACTORNODE linear_search_factor_c(PFACTORNODE first, FACTOR factor);
PFACTORNODE linear_search_factor_s(PFACTORNODE first, FACTOR factor);
void reverse(PFACTORNODE* first_ref);
// old school factoring algorithms
PFACTORNODE trial_division(verylong *zN, PFACTORNODE fn, int *count);
int trial_division_1(long B, verylong *zN, verylong *zf);
PFACTORNODE Brent_Pollard_Method(verylong *zN, PFACTORNODE fn, int *count);
PFACTORNODE Cohen_Pollard_Method(verylong *zN, PFACTORNODE fn, int *count);
PFACTORNODE Lehman(verylong *zN, PFACTORNODE fn, int *count);
PFACTORNODE find_last(PFACTORNODE first);
void quick_sort(PFACTORNODE first, PFACTORNODE last);
/* Lenstra's Elliptic Curve Method define and point structure */
#define CURVES 1024l
typedef struct ec_point
{
verylong zx, zy, zz;
} ECPOINT, *PECPOINT;
PFACTORNODE CohenEllipticCurveMethod(
char *numStr, int base, int expon, int addend);
PFACTORNODE LenstraEllipticCurveMethod(
verylong *zN, PFACTORNODE first, int *count);
#include "LIP_data_structures.h"
/* complete factor comparison */
int factor_z_compare_c(FACTOR lt, FACTOR rt)
{
verylong zlt, zrt;
zcopy(lt.zfactor, &zlt);
zcopy(rt.zfactor, &zrt);
int zc = zcompare(zlt, zrt);
if (zc > 0)
return +1;
if (zc < 0)
return -1;
if (lt.exponent < rt.exponent)
return -1;
if (lt.exponent > rt.exponent)
return -1;
return 0;
}
/* shallow factor comparison */
int factor_z_compare_s(FACTOR lt, FACTOR rt)
{
verylong zlt = 0, zrt = 0;
zcopy(lt.zfactor, &zlt);
zcopy(rt.zfactor, &zrt);
int zc = zcompare(zlt, zrt);
if (zc > 0)
return +1;
if (zc < 0)
return -1;
return 0;
}
/* complete factor_node comparison */
int factor_node_z_compare_c(FACTORNODE lt, FACTORNODE rt)
{
verylong zlt, zrt;
zcopy(lt.factor.zfactor, &zlt);
zcopy(rt.factor.zfactor, &zrt);
int zc = zcompare(zlt, zrt);
if (zc > 0)
return +1;
if (zc < 0)
return -1;
if (lt.factor.exponent < rt.factor.exponent)
return -1;
if (lt.factor.exponent > rt.factor.exponent)
return -1;
return 0;
}
/* shallow factor node comparison */
int factor_node_z_compare_s(FACTORNODE lt, FACTORNODE rt)
{
verylong zlt, zrt;
zcopy(lt.factor.zfactor, &zlt);
zcopy(rt.factor.zfactor, &zrt);
int zc = zcompare(zlt, zrt);
if (zc > 0)
return +1;
if (zc < 0)
return -1;
return 0;
}
/* returns the first factor node */
PFACTORNODE insert_first_factor_node(verylong zf, int exponent)
{
PFACTORNODE first = (PFACTORNODE) malloc(sizeof(FACTORNODE));
memset(first, 0, sizeof(FACTORNODE));
zcopy(zf, &first->factor.zfactor);
first->factor.exponent = exponent;
first->next = NULL;
return first;
}
/* insert a factor node in order */
PFACTORNODE insert_factor_node(PFACTORNODE first, FACTOR factor)
{
if (first == NULL)
return NULL;
PFACTORNODE next = (PFACTORNODE)malloc(sizeof(FACTORNODE));
memset(next, 0, sizeof(FACTORNODE));
next->factor.exponent = factor.exponent;
zcopy(factor.zfactor, &next->factor.zfactor);
next->next = first;
return next;
}
PFACTORNODE delete_factor_node(PFACTORNODE first)
{
if (first == NULL)
return NULL;
PFACTORNODE toDelete = first;
first = first->next;
zfree(&toDelete->factor.zfactor);
free(toDelete);
return first;
}
PFACTORNODE linear_search_factor_c(PFACTORNODE first, FACTOR factor)
{
int found = 0;
PFACTORNODE ptr = first;
if (ptr == NULL)
return NULL;
do
{
found = factor_z_compare_c(ptr->factor, factor);
if (found == 0)
return ptr;
else
ptr = ptr->next;
} while (ptr != NULL);
return NULL;
}
void reverse(PFACTORNODE* first_ref)
{
PFACTORNODE prev = NULL;
PFACTORNODE current = *first_ref;
PFACTORNODE next;
while (current != NULL) {
next = current->next;
current->next = prev;
prev = current;
current = next;
}
*first_ref = prev;
}
PFACTORNODE linear_search_factor_s(PFACTORNODE first, FACTOR factor)
{
int found = 0;
PFACTORNODE ptr = first;
if (ptr == NULL)
return NULL;
do
{
found = factor_z_compare_s(ptr->factor, factor);
if (found == 0)
return ptr;
else
ptr = ptr->next;
} while (ptr != NULL);
return NULL;
}
PFACTORNODE partition(PFACTORNODE first, PFACTORNODE last)
{
if (first == last || first == NULL || last == NULL)
return first;
PFACTORNODE pivot = first;
PFACTORNODE front = first;
FACTOR temp;
while (front != NULL && front != last)
{
if (zcompare(
front->factor.zfactor, last->factor.zfactor) < 0) {
pivot = first;
temp = first->factor;
first->factor = front->factor;
front->factor = temp;
first = first->next;
}
front = front->next;
}
temp = first->factor;
first->factor = last->factor;
last->factor = temp;
return pivot;
}
PFACTORNODE find_last(PFACTORNODE first)
{
while (first->next != NULL)
{
first = first->next;
}
return first;
}
void quick_sort(PFACTORNODE first, PFACTORNODE last)
{
if (first == last) {
return;
}
PFACTORNODE pivot = partition(first, last);
if (pivot != NULL && pivot->next != NULL) {
quick_sort(pivot->next, last);
}
if (pivot != NULL && first != pivot) {
quick_sort(first, pivot);
}
}
PFACTORNODE do_Brent_Pollard(verylong *zN, PFACTORNODE first, int *count)
{
int e, one, pr;
long c, i, k, l;
verylong zP = 0, za = 0, zb = 0, zg = 0, zn = 0;
verylong zx = 0, zx1 = 0, zy = 0;
FACTOR factor;
PFACTORNODE current = NULL;
zcopy(*zN, &zn);
do {
c = 0, k = l = 1;
zone(&zP);
zintoz(2l, &zy);
zintoz(2l, &zx);
zintoz(2l, &zx1);
L2:
zsq(zx, &za);
zsadd(za, 1l, &zb);
zmod(zb, zn, &zx);
zsub(zx1, zx, &za);
zmulmod(za, zP, zn, &zb);
zcopy(zb, &zP);
if (++c == 20) {
zgcd(zP, zn, &zg);
if (zscompare(zg, 1l) > 0)
goto L4;
zcopy(zx, &zy);
c = 0;
}
if (--k != 0)
goto L2;
zgcd(zP, zn, &zg);
if (zscompare(zg, 1l) > 0)
goto L4;
zcopy(zx, &zx1);
k = l, l *= 2;
for (i = 0; i < k; i++)
{
zsq(zx, &za);
zsadd(za, 1l, &zb);
zmod(zb, zn, &zx);
}
zcopy(zx, &zy);
c = 0;
goto L2;
L4:
do {
zsq(zy, &za);
zsadd(za, 1l, &zb);
zmod(zb, zn, &zy);
zsub(zx1, zy, &za);
zgcd(za, zn, &zg);
} while (zscompare(zg, 1l) == 0);
if (zcompare(zg, zn) == 0)
{
fprintf(stderr, "fatal error\nBrent's method failed\n");
exit(1);
}
if (!zprobprime(zg, 20l))
{
zcopy(zg, &za);
if ((current = trial_division(&zg, first, count)) == NULL) {
fprintf(stderr, "fatal error\ncould not trial divide\n");
exit(1);
}
first = current;
zcopy(za, &zg);
zdiv(zn, zg, &za, &zb);
zcopy(za, &zn);
}
else
{
e = 0;
do {
zdiv(zn, zg, &za, &zb);
zcopy(za, &zn);
zmod(zn, zg, &za);
e++;
} while (zscompare(za, 0l) == 0);
if (first != NULL)
{
factor.exponent = e;
factor.zfactor = 0;
zcopy(zg, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
}
}
one = zscompare(zn, 1l) == 0;
if (!one)
{
pr = zprobprime(zn, 20l);
if (pr)
{
factor.exponent = 1;
factor.zfactor = 0;
zcopy(zn, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
break;
}
}
} while (!one);
zfree(&zP);
zfree(&za);
zfree(&zb);
zfree(&zg);
zfree(&zn);
zfree(&zx);
zfree(&zx1);
zfree(&zy);
return first;
}
PFACTORNODE Brent_Pollard_Method(
verylong *zN,
PFACTORNODE first,
int *count)
{
verylong zq = 0, zr = 0;
FACTOR factor;
PFACTORNODE current = NULL, ptr;
*count = 0;
while (1)
{
current = do_Brent_Pollard(zN, first, count);
first = current;
ptr = first;
while (ptr != NULL)
{
zdiv(*zN, ptr->factor.zfactor, &zq, &zr);
while (zscompare(zr, 0l))
{
zdiv(zq, ptr->factor.zfactor, &zq, &zr);
}
if (zscompare(zq, 1l))
goto out_loop;
if (zprobprime(zq, 20l))
goto out_loop;
ptr = ptr->next;
}
if (zscompare(*zN, 1l) == 0)
return first;
if (zprobprime(*zN, 20l) > 0)
break;
}
out_loop:
zfree(&zq);
zfree(&zr);
factor.exponent = 1;
factor.zfactor = 0;
zcopy(*zN, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
return first;
}
PFACTORNODE first_stage(long k, verylong *zN, long x0, long *p,
verylong *zx, PFACTORNODE first, int *count)
{
long c = 0, i = -1, j = i, l, q, q1;
static verylong zg = 0, zq1 = 0, zt = 0, zx1 = 0, zy = 0;
FACTOR factor;
PFACTORNODE current = NULL;
*count = 0;
zzero(&zg);
zsadd(zg, x0, zx);
zcopy(*zx, &zy);
L2:
i++;
if (i > k) {
zsadd(*zx, -1l, &zx1);
zgcd(zx1, *zN, &zg);
if (zscompare(zg, 1l) == 0) return NULL;
else {
i = j;
zcopy(zy, zx);
goto L5;
}
}
else {
q = p[i];
q1 = q;
l = my_B / q;
}
while (q1 <= l) q1 *= q;
zzero(&zt);
zsadd(zt, q1, &zq1);
zcopy(*zx, &zt);
zexpmod(zt, zq1, *zN, zx);
if (++c < 20) goto L2;
zsadd(*zx, -1l, &zx1);
zgcd(zx1, *zN, &zg);
if (zscompare(zg, 1l) == 0) {
c = 0;
j = i;
zcopy(*zx, &zy);
goto L2;
}
else {
i = j;
zcopy(zy, zx);
}
L5:
i++;
q = p[i];
q1 = q;
L6:
zzero(&zt);
zsadd(zt, q, &zq1);
zcopy(*zx, &zt);
zexpmod(zt, zq1, *zN, zx);
zsadd(*zx, -1l, &zx1);
zgcd(zx1, *zN, &zg);
if (zscompare(zg, 1l) == 0) {
q1 *= q;
if (q1 <= my_B) goto L6; else goto L5;
}
else
{
if (zcompare(zg, *zN) < 0) {
factor.exponent = 1;
factor.zfactor = 0;
zcopy(zg, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
zcopy(*zN, &zq1);
zdiv(zq1, zg, zN, &zx1);
return first;
}
if (zcompare(zg, *zN) == 0)
{
first = NULL;
}
}
zfree(&zg);
zfree(&zq1);
zfree(&zt);
zfree(&zx1);
zfree(&zy);
return first;
}
PFACTORNODE Cohen_Pollard_Method(
verylong *zN,
PFACTORNODE first,
int *count)
{
long cnt, i, *p = (long *)malloc(my_NUMBER_PRIMES * sizeof(long));
verylong zx = 0;
PFACTORNODE current = NULL;
zpstart2();
for (i = 0; i < my_NUMBER_PRIMES; i++)
p[i] = zpnext();
i = 0;
int j = p[0];
do {
while ((current = first_stage(my_NUMBER_PRIMES, zN,
j, p, &zx, first, &cnt)) != NULL && i < my_NUMBER_PRIMES)
{
i++;
j = p[i];
first = current;
count += cnt;
}
j = p[++i];
} while (
!zscompare(*zN, 1l) &&
!zprobprime(*zN, 20) &&
i < my_NUMBER_PRIMES);
FACTOR factor;
factor.exponent = 1;
factor.zfactor = *zN;
current = insert_factor_node(first, factor);
first = current;
free(p);
return first;
}
int square_test(verylong zn, verylong *zq)
{
int square = 1;
long q11[11], q63[63], q64[64], q65[65];
long k, r, t;
verylong zd = 0;
for (k = 0; k < 11; k++) q11[k] = 0;
for (k = 0; k < 6; k++) q11[(k * k) % 11] = 1;
for (k = 0; k < 63; k++) q63[k] = 0;
for (k = 0; k < 32; k++) q63[(k * k) % 63] = 1;
for (k = 0; k < 64; k++) q64[k] = 0;
for (k = 0; k < 32; k++) q64[(k * k) % 64] = 1;
for (k = 0; k < 65; k++) q65[k] = 0;
for (k = 0; k < 33; k++) q65[(k * k) % 65] = 1;
t = zsmod(zn, 64l);
r = zsmod(zn, 45045);
if (q64[t] == 0) square = 0;
if (square && q63[r % 63] == 0) square = 0;
if (square && q65[r % 65] == 0) square = 0;
if (square && q11[r % 11] == 0) square = 0;
if (square) {
zsqrt(zn, zq, &zd);
if (zscompare(zd, 0l) != 0) square = 0;
}
zfree(&zd);
return square;
}
int trial_division_1(long B, verylong *zN, verylong *zf)
{
int e = 0;
long p;
verylong za = 0, zb = 0;
zcopy(*zN, &za);
zpstart2();
do {
p = zpnext();
if (zsmod(za, p) == 0) {
do {
zsdiv(za, p, &zb);
zcopy(zb, &za);
e++;
} while (zsmod(za, p) == 0);
zintoz(p, zf);
zcopy(za, zN);
}
} while (!e && p <= B);
zfree(&za);
zfree(&zb);
return e;
}
int do_Lehman(verylong *zN, verylong *zf)
{
long B = pow(zdoub(*zN), 1.0 / 3.0), a, k = 0, m, r;
int e = trial_division_1(B, zN, zf);
verylong za = 0, zb = 0, zc = 0, zh = 0, zl = 0;
verylong zr = 0, zs = 0;
if (!e) {
L2:
k++;
if (k > B) goto L4;
if (!(k & 1)) {
zone(&zr);
m = 2;
}
else {
zsadd(*zN, k, &zr);
m = 4;
}
zsmul(*zN, k, &za);
zlshift(za, 2l, &zl);
zintoz(B, &zh);
zsq(zh, &zb);
zadd(zl, zb, &zh);
zsqrt(zl, &za, &zb);
zsq(za, &zs);
while (zcompare(zs, zl) < 0) {
zsadd(za, 1l, &zb);
zcopy(zb, &za);
zsq(za, &zs);
}
r = zsmod(zr, m);
while (!e && zcompare(zs, zh) <= 0) {
a = zsmod(za, m);
if (a == r) {
zsub(zs, zl, &zc);
if (zscompare(zc, 0l) == 0 || square_test(zc, &zb)) {
zadd(za, zb, &zc);
zgcd(zc, *zN, zf);
do {
zdiv(*zN, *zf, &za, &zb);
zcopy(za, zN);
zmod(za, *zf, &zb);
e++;
} while (zscompare(zb, 0l) == 0);
}
}
if (!e) {
zsadd(za, 1l, &zb);
zcopy(zb, &za);
zsq(za, &zs);
}
}
if (!e) goto L2;
}
L4:
zfree(&za);
zfree(&zb);
zfree(&zc);
zfree(&zh);
zfree(&zl);
zfree(&zr);
zfree(&zs);
return e;
}
PFACTORNODE Lehman(
verylong *zN,
PFACTORNODE first,
int *count)
{
int e = 0;
verylong zf = 0;
FACTOR factor;
PFACTORNODE current = NULL;
while (1)
{
e = do_Lehman(zN, &zf);
if (e)
{
factor.exponent = e;
factor.zfactor = 0;
zcopy(zf, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
}
else if (zprobprime(*zN, 20l))
{
factor.exponent = 1;
factor.zfactor = 0;
zcopy(*zN, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
break;
}
else
break;
}
zfree(&zf);
return first;
}
PFACTORNODE trial_division(
verylong *zN,
PFACTORNODE first,
int *count)
{
int e, one = 0, pr = 0;
long B, p;
verylong za = 0, zb = 0, zl = 0, zo = 0, zp = 0;
FACTOR factor;
PFACTORNODE current = NULL;
zsqrt(*zN, &za, &zb);
zone(&zo);
zcopy(*zN, &zl);
if (zscompare(za, my_BOUND) > 0)
B = my_BOUND;
else
B = ztoint(za);
zcopy(*zN, &za);
zpstart2();
do {
p = zpnext();
if (zsmod(za, p) == 0)
{
e = 0;
do {
zsdiv(za, p, &zb);
zcopy(zb, &za);
e++;
} while (zsmod(za, p) == 0);
zintoz(p, &zp);
factor.exponent = e;
factor.zfactor = 0;
zcopy(zp, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
zcopy(za, zN);
one = zscompare(*zN, 1l) == 0;
pr = zprobprime(*zN, 20l);
if (!one && pr)
{
factor.exponent = 1;
zcopy(*zN, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
break;
}
else
{
factor.exponent = 1;
factor.zfactor = 0;
zcopy(*zN, &factor.zfactor);
current = insert_factor_node(first, factor);
first = current;
*count = *count + 1;
}
}
} while (!one && p <= B);
PFACTORNODE last = (PFACTORNODE)malloc(sizeof(FACTORNODE));
last->factor.exponent = 1;
last->factor.zfactor = 0;
zcopy(zl, &last->factor.zfactor);
current = insert_factor_node(first, last->factor);
first = current;
*count = *count + 1;
zfree(&za);
zfree(&zb);
zfree(&zl);
zfree(&zo);
zfree(&zp);
return first;
}
int partial_addition(verylong za,
verylong zn,
ECPOINT P,
ECPOINT Q,
PECPOINT R,
verylong *zd)
/* returns 0 if sum is found or 1 if divisor is found */
{
int result = -1, value = 0;
verylong zb = 0, zc = 0, zl = 0, zs = 0;
verylong zt = 0, 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);
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);
zfree(&zy2);
return value;
}
int multiply(long k,
verylong za,
verylong zn,
ECPOINT P,
PECPOINT R,
verylong *zd)
{
int value = 0;
ECPOINT A = { 0, 0, 0 };
ECPOINT S = { 0, 0, 0 };
ECPOINT T = { 0, 0, 0 };
zone(&A.zx);
zone(&A.zy);
zone(&A.zz);
zone(&S.zx);
zone(&S.zy);
zone(&S.zz);
zone(&T.zx);
zone(&T.zy);
zone(&T.zz);
R = (PECPOINT)malloc(sizeof(ECPOINT));
if (R == NULL)
exit(0);
memset(R, 0, sizeof(ECPOINT));
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;
}
/* the following definition limits the L3 loop */
int LenstrasECM(verylong *zN, verylong *zg)
{
int expon = 0, found;
long B = 100000000l, i, j, k = 0, l, *p, q, q1;
long giveUp = z2log(*zN);
ECPOINT 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) {
expon = -1;
goto L4;
}
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 == giveUp) {
expon = 0;
goto L4;
}
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;
}
if (!zprobprime(*zg, 5l))
goto L4;
expon = 0;
do {
zdiv(*zN, *zg, &zb, &zd);
zcopy(zb, zN);
zmod(*zN, *zg, &zd);
expon++;
} while (zscompare(zd, 0l) == 0);
L4:
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;
}
PFACTORNODE CohenEllipticCurveMethod(
char *numStr, int base, int iExpon, int addend)
{
PFACTORNODE current, first;
FACTOR factor;
int result;
int cant, expon, one, pri;
verylong zn = 0, zN = 0, zd = 0, zg = 0;
if (strlen(numStr) != 0)
zstrtoz(numStr, &zN);
else
{
static verylong zb = 0, zp = 0;
zintoz(base, &zb);
zsexp(zb, iExpon, &zp);
zsadd(zp, addend, &zN);
}
zcopy(zN, &zn);
first = insert_first_factor_node(zN, 1);
cant = pri = 0;
do {
expon = LenstrasECM(&zN, &zg);
if (expon == -1)
{
result = -1;
goto L1;
}
if (zprobprime(zg, 20l))
{
factor.exponent = expon;
factor.zfactor = zg;
current = insert_factor_node(first, factor);
first = current;
}
one = zscompare(zN, 1l) == 0;
if (zprobprime(zN, 20l))
{
factor.exponent = 1;
factor.zfactor = zN;
current = insert_factor_node(first, factor);
first = current;
break;
}
} while (!one);
L1:
zfree(&zn);
zfree(&zN);
zfree(&zd);
zfree(&zg);
return first;
}
PFACTORNODE LenstraEllipticCurveMethod(
verylong *zN, PFACTORNODE first, int *count)
{
long curvebnd = 1024, e = 0, exp = 0, phase1bnd = 64;
long grow = 16, info = 2, nbtests = 20, s = 1;
verylong zf = 0, zn = 0, zq = 0, zr = 0;
FACTOR factor;
PFACTORNODE current = NULL;
zcopy(*zN, &zn);
while (exp != -1)
{
exp = zfecm(zn, &zf, s, &curvebnd, &phase1bnd,
grow, nbtests, info, (FILE *)NULL);
if (exp == 1 || exp == -1)
{
if (zscompare(zf, 1) != 0)
{
*count = *count + 1;
e = 0;
zdiv(zn, zf, &zq, &zr);
while (zscompare(zr, 0l) == 0)
{
e++;
zcopy(zq, &zn);
zdiv(zn, zf, &zq, &zr);
}
factor.exponent = e;
factor.zfactor = zf;
current = insert_factor_node(first, factor);
first = current;
}
else
{
if (zprobprime(zn, 20))
{
factor.exponent = 1;
factor.zfactor = zn;
current = insert_factor_node(first, factor);
first = current;
}
}
}
else
break;
}
zfree(&zf);
zfree(&zn);
zfree(&zq);
zfree(&zr);
return first;
}
This set of experiments is dedicated to Alfred J. Menezes, Arjen K. Lenstra, H. W. Lenstra, Jr., J. M. Pollard, Henri Cohen, Hans Riesel, H. T. Lau, Charles Petzold, Richard O. Chapman, Gerry V. Dozier, Homer Carlisle, Fay A. Riddle, Brooks Shelhorse, among many others whose names have faded from my memory.
The factoring algorithms are trial division, Brent’s modification of the famous Pollard rho Method, and Lenstra’s Elliptic Curve Method. I used repetitions of the digit string “1234567890” of varying length. I chose these strings because of the richness of prime factors.
These algorithms are not guaranteed to completely factor these relatively large numbers. I gave up on using the elliptic curve method for numbers of greater than 40-digits. I have to give up on Brent’s modification of the Pollard rho method for numbers of more than 70-digits.
/*
Author: Pate Williams (c) 1997 - 2022
Algorithm 8.5.2 (Pollard rho). See "A Course in
Computational Algebraic Number Theory" by Henri
Cohen page 429.
Author: Pate Williams (c) 1997 - 2022
Algorithm 10.3.3 (Lenstra's ECM). See "A Course
in Computational Algebraic Number Theory" by
Henri Cohen page 488.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "lip.h"
typedef struct fact {
verylong zf;
int expon;
} FACTOR, PFACTOR;
int BOUND = 100000000;
void add_factor(int e, verylong zp,
FACTOR f[], int *count)
{
int done = 0, found = 0, i, j = 0;
for (i = 0; !done && i < *count; i++)
{
found = zcompare(zp, f[i].zf);
if (found == 0)
{
done = 1;
j = i;
}
}
if (found == 0)
{
for (i = *count - 2; i >= j; i--)
{
f[i + 1].expon = f[i].expon;
zcopy(f[i].zf, &f[i + 1].zf);
}
*count = *count + 1;
f[j].expon = e;
zcopy(zp, &f[j].zf);
}
else if (found == 0)
f[j].expon++;
else
{
f[*count].expon = e;
zcopy(zp, &f[*count].zf);
*count = *count + 1;
}
}
int trial_division(verylong *zN, FACTOR f[], int *count)
{
int e, one = 0, pr = 0;
long B, p;
verylong za = 0, zb = 0, zp = 0;
zsqrt(*zN, &za, &zb);
if (zscompare(za, BOUND) > 0)
B = BOUND;
else
B = ztoint(za);
zcopy(*zN, &za);
zpstart2();
do {
p = zpnext();
if (zsmod(za, p) == 0)
{
e = 0;
do {
zsdiv(za, p, &zb);
zcopy(zb, &za);
e++;
} while (zsmod(za, p) == 0);
zintoz(p, &zp);
add_factor(e, zp, f, count);
zcopy(za, zN);
one = zscompare(*zN, 1l) == 0;
if (!one)
pr = zprobprime(*zN, 5l);
}
} while (!one && !pr && p <= B);
if (pr)
add_factor(1, *zN, f, count);
zfree(&za);
zfree(&zb);
zfree(&zp);
return p <= B;
}
void Brent(verylong *zN, FACTOR f[], int *count)
{
int e, one, pr;
long c, i, k, l;
verylong zP = 0, za = 0, zb = 0, zg = 0, zn = 0;
verylong zx = 0, zx1 = 0, zy = 0;
zcopy(*zN, &zn);
do {
c = 0, k = l = 1;
zone(&zP);
zintoz(2l, &zy);
zintoz(2l, &zx);
zintoz(2l, &zx1);
L2:
zsq(zx, &za);
zsadd(za, 1l, &zb);
zmod(zb, zn, &zx);
zsub(zx1, zx, &za);
zmulmod(za, zP, zn, &zb);
zcopy(zb, &zP);
if (++c == 20) {
zgcd(zP, zn, &zg);
if (zscompare(zg, 1l) > 0) goto L4;
zcopy(zx, &zy);
c = 0;
}
if (--k != 0)
goto L2;
zgcd(zP, zn, &zg);
if (zscompare(zg, 1l) > 0)
goto L4;
zcopy(zx, &zx1);
k = l, l *= 2;
for (i = 0; i < k; i++)
{
zsq(zx, &za);
zsadd(za, 1l, &zb);
zmod(zb, zn, &zx);
}
zcopy(zx, &zy);
c = 0;
goto L2;
L4:
do {
zsq(zy, &za);
zsadd(za, 1l, &zb);
zmod(zb, zn, &zy);
zsub(zx1, zy, &za);
zgcd(za, zn, &zg);
} while (zscompare(zg, 1l) == 0);
if (zcompare(zg, zn) == 0) {
fprintf(stderr, "fatal error\nBrent's method failed\n");
exit(1);
}
if (!zprobprime(zg, 20l))
{
zcopy(zg, &za);
if (!trial_division(&zg, f, count)) {
fprintf(stderr, "fatal error\ncould not trial divide\n");
exit(1);
}
zcopy(za, &zg);
zdiv(zn, zg, &za, &zb);
zcopy(za, &zn);
}
else
{
e = 0;
do {
zdiv(zn, zg, &za, &zb);
zcopy(za, &zn);
zmod(zn, zg, &za);
e++;
} while (zscompare(za, 0l) == 0);
add_factor(e, zg, f, count);
}
one = zscompare(zn, 1l) == 0;
if (!one) pr = zprobprime(zn, 5l);
} while (!one && !pr);
if (!one)
add_factor(1, zn, f, count);
zfree(&zP);
zfree(&za);
zfree(&zb);
zfree(&zg);
zfree(&zn);
zfree(&zx);
zfree(&zx1);
zfree(&zy);
}
int partition(FACTOR a[], int n, int lo, int hi)
{
int pivotIndex = lo + (hi - lo) / 2;
FACTOR x = a[pivotIndex];
FACTOR t = x;
a[pivotIndex] = a[hi];
a[hi] = t;
int storeIndex = lo;
for (int i = lo; i < hi; i++)
{
if (zcompare(a[i].zf, x.zf) < 0)
{
t = a[i];
a[i] = a[storeIndex];
a[storeIndex++] = t;
}
}
t = a[storeIndex];
a[storeIndex] = a[hi];
a[hi] = t;
return storeIndex;
}
void do_quick_sort(FACTOR a[], int n, int p, int r)
{
if (p < r)
{
int q = partition(a, n, p, r);
do_quick_sort(a, n, p, q - 1);
do_quick_sort(a, n, q + 1, r);
}
}
void quick_sort(FACTOR a[], int n)
{
do_quick_sort(a, n, 0, n - 1);
}
#ifndef CLK_TCK
#define CLK_TCK CLOCKS_PER_SEC
#endif
#define CURVES 1024l
struct point { verylong zx, zy, zz; };
struct factor { int expon; verylong prime; };
static char cexp[32][32], cfac[32][512];
static int cnt, primeCount;
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);
zfree(&zy2);
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;
}
/* the following definition limits the L3 loop */
int LenstrasECM(verylong *zN, verylong *zg)
{
int expon = 0, found = 0;
long B = 2000000l, i, j, k = 0, l, *p, q, q1;
long giveUp = z2log(*zN);
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) {
expon = -1;
goto L4;
}
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 == giveUp) {
expon = 0;
goto L4;
}
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;
}
if (!zprobprime(*zg, 5l))
goto L4;
expon = 0;
do {
zdiv(*zN, *zg, &zb, &zd);
zcopy(zb, zN);
zmod(*zN, *zg, &zd);
expon++;
} while (zscompare(zd, 0l) == 0);
L4:
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 ECMVLConvert(verylong zx, char *str, int length)
{
verylong zy = 0;
int a, result = 1, i = 0;
while (zscompare(zx, 0) > 0)
{
zsdiv(zx, 10, &zy);
a = zsmod(zx, 10);
zcopy(zy, &zx);
if (i < length)
str[i++] = a + '0';
else
{
result = 0;
break;
}
}
str[i] = '\0';
_strrev(str);
zfree(&zy);
return result;
}
int ECMLast(char exp[32][32], char fac[32][512], struct factor *f, int *cnt)
{
int i;
for (i = 0; i < *cnt; i++)
{
if (ECMVLConvert(f[i].prime, fac[i], 512) == 1)
sprintf(exp[i], "%d", f[i].expon);
else
return 0;
}
return 1;
}
int DoECM(char exp[32][32], char fac[32][512], char *numStr,
int *cnt, int base, int iExpon, int addend)
{
struct factor f[32];
int i, j, result;
int cant, expon, one, pri;
verylong zn = 0, zN = 0, zd = 0, zg = 0;
if (strlen(numStr) != 0)
zstrtoz(numStr, &zN);
else
{
static verylong zb = 0, zp = 0;
zintoz(base, &zb);
zsexp(zb, iExpon, &zp);
zsadd(zp, addend, &zN);
}
zcopy(zN, &zn);
for (i = 0; i < 32; i++)
{
f[i].expon = 0;
f[i].prime = 0;
}
*cnt = 0;
f[*cnt].expon = 1;
zcopy(zN, &f[*cnt].prime);
*cnt = *cnt + 1;
if (z2log(zn) > 512)
{
result = 3;
goto L1;
}
if (zprobprime(zN, 5l))
{
result = 2;
goto L1;
}
cant = pri = 0;
do {
expon = LenstrasECM(&zN, &zg);
if (expon == -1)
{
result = -1;
goto L1;
}
if (zprobprime(zg, 5l))
{
f[*cnt].expon = expon;
zcopy(zg, &f[*cnt].prime);
*cnt = *cnt + 1;
}
one = zscompare(zN, 1l) == 0;
if (zprobprime(zN, 5l))
{
f[*cnt].expon = 1;
zcopy(zN, &f[*cnt].prime);
*cnt = *cnt + 1;
break;
}
} while (!one);
/* selection sort the factors */
for (i = 1; i < *cnt - 1; i++)
{
for (j = i + 1; j < *cnt; j++)
{
if (zcompare(f[i].prime, f[j].prime) > 0)
{
int temp;
verylong zt = 0;
zcopy(f[i].prime, &zt);
zcopy(f[j].prime, &f[i].prime);
zcopy(zt, &f[j].prime);
temp = f[i].expon;
f[i].expon = f[j].expon;
f[j].expon = temp;
}
}
}
result = ECMLast(exp, fac, f, cnt);
L1:
if (result == 2)
{
int temp = ECMLast(exp, fac, f, cnt);
if (temp == 0)
result = 3;
}
zfree(&zn);
zfree(&zN);
zfree(&zd);
zfree(&zg);
return result;
}
int PrimeCheck(char fac[32][512], int cnt)
{
int i, result = 0;
verylong zp = 0;
for (i = 0; i < cnt; i++)
{
zstrtoz(fac[i], &zp);
if (zprobprime(zp, 20l))
result++;
}
return result;
}
int ECMMethod(char *numStr, int base, int expon, int addend)
{
int i, result;
cnt = 0;
for (i = 0; i < 32; i++)
{
cexp[i][0] = '\0';
cfac[i][0] = '\0';
}
result = DoECM(cexp, cfac, numStr, &cnt, base, expon, addend);
primeCount = PrimeCheck(cfac, cnt);
return result;
}
int main()
{
int i = 0, zacount = 0;
verylong zN = 0, zf = 0, zn = 0, zo = 0;
verylong zq = 0, zr = 0, zs = 0;
FACTOR za[32] = { 0 };
zintoz(1, &zo);
while (1) {
clock_t time0;
char numStr[256] = { 0 }, option[256] = { 0 }, *bPtr, *ePtr;
double time;
int c, base, expon, addend, iExpon = 0, positive, result;
printf("Menu\n");
printf("1 Trial Division\n");
printf("2 Pollard rho\n");
printf("3 Lenstra's Elliptic Curve Method\n");
printf("4 Exit\n");
printf("Enter choice (1 - 4): ");
c = getchar();
i = 0;
while (c != '\n' && i < 256)
{
option[i++] = c;
c = getchar();
}
option[i] = '\0';
if (option[0] == '4')
break;
if (option[0] < '1' || option[0] > '3')
{
printf("Unknown choice, please reenter choice\n");
continue;
}
printf("enter the number below (0 to quit):\n");
i = 0;
c = getchar();
while (c != '\n' && i < 256)
{
numStr[i++] = c;
c = getchar();
}
numStr[i] = '\0';
bPtr = strchr(numStr, '^');
if (bPtr != NULL)
{
*bPtr = '\0';
base = atoi(numStr);
ePtr = strchr(bPtr + 1, '+');
if (ePtr != NULL)
positive = 1;
else
{
ePtr = strchr(bPtr + 1, '-');
if (ePtr != NULL)
positive = 0;
}
if (ePtr != NULL)
*ePtr = '\0';
expon = atoi(bPtr + 1);
if (ePtr != NULL) {
addend = atoi(ePtr + 1);
if (positive == 0)
addend = -addend;
}
else
addend = 0;
numStr[0] = '\0';
}
else
{
addend = atoi(numStr);
if (addend == 0)
break;
}
if (strlen(numStr) != 0)
zstrtoz(numStr, &zN);
else
{
static verylong zb = 0, zp = 0;
zintoz(base, &zb);
zsexp(zb, iExpon, &zp);
zsadd(zp, addend, &zN);
}
zcopy(zN, &zn);
if (zmcomposite(zn, 20) == 0)
{
printf("Number is prime\n");
continue;
}
zacount = 0;
for (i = 0; i < 32; i++)
{
if (za[i].zf != 0)
zfree(&za[i].zf);
za[i].expon = 0;
}
time0 = clock();
if (option[0] == '1')
{
printf("Trial division factors:\n");
trial_division(&zN, za, &zacount);
}
else if (option[0] == '2')
{
printf("Pollard rho factors:\n");
Brent(&zN, za, &zacount);
}
else if (option[0] == '3')
{
printf("Lenstra's ECM factors:\n");
result = ECMMethod(numStr, base, expon, addend);
if (result == 0)
printf("can't completely factor this number\n");
else if (result == 1)
{
printf("%s\nnumber is composite factors:\n", cfac[0]);
for (i = 1; i < cnt; i++)
if (strcmp(cexp[i], "1") == 0)
printf("\t%s\n", cfac[i]);
else
printf("\t%s ^ %s\n", cfac[i], cexp[i]);
}
else if (result == 2)
printf("%s\nnumber is prime:\n", cfac[0]);
else if (result == 3)
printf("number or prime factor is too large\n");
}
if (option[0] < '3' && zacount > 0)
{
quick_sort(za, zacount);
for (i = 0; i < zacount; i++)
{
zwrite(za[i].zf);
if (za[i].expon > 1)
printf(" ^ %d ", za[i].expon);
printf("\n");
}
for (i = 0; i < zacount; i++)
if (za[i].zf != 0)
zfree(&za[i].zf);
zacount = 0;
}
time = (clock() - time0) / (double)CLK_TCK;
printf("total time required: %f seconds\n", time);
}
zfree(&zN);
zfree(&zn);
zfree(&zf);
zfree(&zo);
zfree(&zq);
zfree(&zr);
}
/*
Author: Pate Williams (c) 1997 - 2009
Algorithm 10.3.3 (Lenstra's ECM). See "A Course
in Computational Algebraic Number Theory" by
Henri Cohen page 488.
*/
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "lip.h"
#ifndef CLK_TCK
#define CLK_TCK CLOCKS_PER_SEC
#endif
/*
if you creating a Windows only program
the following function can be edited
out or completely removed since strrev
is defined in Windows Visual Studio C
*/
char * strrev(char *s)
{
int i, length = strlen(s), j = length - 1;
char *t = malloc((length + 1) * sizeof(char));
if (t == NULL)
return NULL;
for (i = 0; i < length; i++, j--)
t[j] = s[i];
t[length] = '\0';
strcpy(s, t);
free(t);
return s;
}
#define CURVES 1024l
struct point {verylong zx, zy, zz;};
struct factor {int expon; verylong prime;};
static char cexp[32][32], cfac[32][512];
static int cnt, primeCount;
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);
zfree(&zy2);
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;
}
/* the following definition limits the L3 loop */
int LenstrasECM(verylong *zN, verylong *zg)
{
int expon = 0, found;
long B = 2000000l, i, j, k = 0, l, *p, q, q1;
long giveUp = z2log(*zN);
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) {
expon = - 1;
goto L4;
}
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 == giveUp) {
expon = 0;
goto L4;
}
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;
}
if (!zprobprime(*zg, 5l))
goto L4;
expon = 0;
do {
zdiv(*zN, *zg, &zb, &zd);
zcopy(zb, zN);
zmod(*zN, *zg, &zd);
expon++;
} while (zscompare(zd, 0l) == 0);
L4:
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 ECMVLConvert(verylong zx, char *str, int length)
{
verylong zy = 0;
int a, result = 1, i = 0;
while (zscompare(zx, 0) > 0)
{
zsdiv(zx, 10, &zy);
a = zsmod(zx, 10);
zcopy(zy, &zx);
if (i < length)
str[i++] = a + '0';
else
{
result = 0;
break;
}
}
str[i] = '\0';
strrev(str);
zfree(&zy);
return result;
}
int ECMLast(char exp[32][32], char fac[32][512], struct factor *f, int *cnt)
{
int i;
for (i = 0; i < *cnt; i++)
{
if (ECMVLConvert(f[i].prime, fac[i], 512) == 1)
sprintf(exp[i], "%d", f[i].expon);
else
return 0;
}
return 1;
}
int DoECM(char exp[32][32], char fac[32][512], char *numStr,
int *cnt, int base, int iExpon, int addend)
{
struct factor f[32];
int i, j, result;
int cant, expon, one, pri;
verylong zn = 0, zN = 0, zd = 0, zg = 0;
if (strlen(numStr) != 0)
zstrtoz(numStr, &zN);
else
{
static verylong zb = 0, zp = 0;
zintoz(base, &zb);
zsexp(zb, iExpon, &zp);
zsadd(zp, addend, &zN);
}
zcopy(zN, &zn);
for (i = 0; i < 32; i++)
{
f[i].expon = 0;
f[i].prime = 0;
}
*cnt = 0;
f[*cnt].expon = 1;
zcopy(zN, &f[*cnt].prime);
*cnt = *cnt + 1;
if (z2log(zn) > 512)
{
result = 3;
goto L1;
}
if (zprobprime(zN, 5l))
{
result = 2;
goto L1;
}
cant = pri = 0;
do {
expon = LenstrasECM(&zN, &zg);
if (expon == - 1)
{
result = - 1;
goto L1;
}
if (zprobprime(zg, 5l))
{
f[*cnt].expon = expon;
zcopy(zg, &f[*cnt].prime);
*cnt = *cnt + 1;
}
one = zscompare(zN, 1l) == 0;
if (zprobprime(zN, 5l))
{
f[*cnt].expon = 1;
zcopy(zN, &f[*cnt].prime);
*cnt = *cnt + 1;
break;
}
} while (!one);
/* selection sort the factors */
for (i = 1; i < *cnt - 1; i++)
{
for (j = i + 1; j < *cnt; j++)
{
if (zcompare(f[i].prime, f[j].prime) > 0)
{
int temp;
verylong zt = 0;
zcopy(f[i].prime, &zt);
zcopy(f[j].prime, &f[i].prime);
zcopy(zt, &f[j].prime);
temp = f[i].expon;
f[i].expon = f[j].expon;
f[j].expon = temp;
}
}
}
result = ECMLast(exp, fac, f, cnt);
L1:
if (result == 2)
{
int temp = ECMLast(exp, fac, f, cnt);
if (temp == 0)
result = 3;
}
zfree(&zn);
zfree(&zN);
zfree(&zd);
zfree(&zg);
return result;
}
int PrimeCheck(char fac[32][512], int cnt)
{
int i, result = 0;
verylong zp = 0;
for (i = 0; i < cnt; i++)
{
zstrtoz(fac[i], &zp);
if (zprobprime(zp, 5l))
result++;
}
return result;
}
int ECMMethod(char *numStr, int base, int expon, int addend)
{
int i, result;
cnt = 0;
for (i = 0; i < 32; i++)
{
cexp[i][0] = '\0';
cfac[i][0] = '\0';
}
result = DoECM(cexp, cfac, numStr, &cnt, base, expon, addend);
primeCount = PrimeCheck(cfac, cnt);
return result;
}
int main(int argc, char *argv[])
{
while (1) {
clock_t time0;
char numStr[1024] = {0}, *bPtr, *ePtr;
double time;
int c, base, expon, addend, i = 0, positive, result;
printf("enter the number below (0 to quit):\n");
c = getchar();
while (c != '\n' && i < 1023)
{
numStr[i++] = c;
c = getchar();
}
numStr[i] = '\0';
bPtr = strchr(numStr, '^');
if (bPtr != NULL)
{
*bPtr = '\0';
base = atoi(numStr);
ePtr = strchr(bPtr + 1, '+');
if (ePtr != NULL)
positive = 1;
else {
ePtr = strchr(bPtr + 1, '-');
if (ePtr != NULL)
positive = 0;
}
if (ePtr != NULL)
*ePtr = '\0';
expon = atoi(bPtr + 1);
if (ePtr != NULL) {
addend = atoi(ePtr + 1);
if (positive == 0)
addend = - addend;
}
else
addend = 0;
numStr[0] = '\0';
}
else
{
addend = atoi(numStr);
if (addend == 0)
break;
}
printf("\n");
time0 = clock();
result = ECMMethod(numStr, base, expon, addend);
if (result == 0)
printf("can't completely factor this number\n");
else if (result == 1)
{
printf("%s\nnumber is composite factors:\n", cfac[0]);
for (i = 1; i < cnt; i++)
if (strcmp(cexp[i], "1") == 0)
printf("\t%s\n", cfac[i]);
else
printf("\t%s ^ %s\n", cfac[i], cexp[i]);
}
else if (result == 2)
printf("%s\nnumber is prime:\n", cfac[0]);
else if (result == 3)
printf("number or prime factor is too large\n");
time = (clock() - time0) / (double) CLK_TCK;
printf("total time required: %f seconds\n", time);
printf("\n");
}
return 0;
}
/*
Author: James Pate Williams, Jr.
Solution of 8-puzzle using A* search
*/
import java.util.*;
import java.awt.*;
import java.awt.event.*;
import javax.swing.*;
class Puzzle8Panel extends JPanel {
char[] s;
int deltaX, deltaY;
int x0, x1, y0, y1;
public Puzzle8Panel(int iWidth, int iHeight, char[] solution) {
x0 = iWidth / 8;
x1 = 7 * x0;
y0 = iHeight / 8;
y1 = 7 * y0;
deltaX = (x1 - x0) / 3;
deltaY = (y1 - y0) / 3;
s = new char[9];
s = solution;
}
public void sevenSegmentDisplay(Graphics g, char digit,
int x, int y, int xUnit, int yUnit) {
int xInc = xUnit / 10;
int yInc = yUnit / 10;
int xPoint1 = x + xInc;
int yPoint1 = y + yInc;
int xPoint2 = x + xUnit - xInc;
int yPoint2 = yPoint1;
int xPoint3 = xPoint1;
int yPoint3 = y + yUnit / 2;
int xPoint4 = xPoint2;
int yPoint4 = yPoint3;
int xPoint5 = xPoint3;
int yPoint5 = y + yUnit - yInc;
int xPoint6 = xPoint4;
int yPoint6 = yPoint5;
g.setColor(Color.white);
switch (digit) {
case '0':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint2, yPoint2, xPoint6, yPoint6);
g.drawLine(xPoint6, yPoint6, xPoint5, yPoint5);
g.drawLine(xPoint5, yPoint5, xPoint1, yPoint1);
break;
case '1':
g.drawLine(xPoint1, yPoint1, xPoint5, yPoint5);
break;
case '2':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint2, yPoint2, xPoint4, yPoint4);
g.drawLine(xPoint4, yPoint4, xPoint3, yPoint3);
g.drawLine(xPoint3, yPoint3, xPoint5, yPoint5);
g.drawLine(xPoint5, yPoint5, xPoint6, yPoint6);
break;
case '3':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint2, yPoint2, xPoint4, yPoint4);
g.drawLine(xPoint4, yPoint4, xPoint3, yPoint3);
g.drawLine(xPoint4, yPoint4, xPoint6, yPoint6);
g.drawLine(xPoint6, yPoint6, xPoint5, yPoint5);
break;
case '4':
g.drawLine(xPoint1, yPoint1, xPoint3, yPoint3);
g.drawLine(xPoint3, yPoint3, xPoint4, yPoint4);
g.drawLine(xPoint4, yPoint4, xPoint2, yPoint2);
g.drawLine(xPoint4, yPoint4, xPoint6, yPoint6);
break;
case '5':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint1, yPoint1, xPoint3, yPoint3);
g.drawLine(xPoint3, yPoint3, xPoint4, yPoint4);
g.drawLine(xPoint4, yPoint4, xPoint6, yPoint6);
g.drawLine(xPoint6, yPoint6, xPoint5, yPoint5);
break;
case '6':
g.drawLine(xPoint1, yPoint1, xPoint3, yPoint3);
g.drawLine(xPoint3, yPoint3, xPoint4, yPoint4);
g.drawLine(xPoint4, yPoint4, xPoint6, yPoint6);
g.drawLine(xPoint6, yPoint6, xPoint5, yPoint5);
g.drawLine(xPoint5, yPoint5, xPoint3, yPoint3);
break;
case '7':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint2, yPoint2, xPoint6, yPoint6);
break;
case '8':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint2, yPoint2, xPoint6, yPoint6);
g.drawLine(xPoint6, yPoint6, xPoint5, yPoint5);
g.drawLine(xPoint5, yPoint5, xPoint1, yPoint1);
g.drawLine(xPoint3, yPoint3, xPoint4, yPoint4);
break;
case '9':
g.drawLine(xPoint1, yPoint1, xPoint2, yPoint2);
g.drawLine(xPoint2, yPoint2, xPoint6, yPoint6);
g.drawLine(xPoint1, yPoint1, xPoint3, yPoint3);
g.drawLine(xPoint3, yPoint3, xPoint4, yPoint4);
break;
}
}
public void paintComponent(Graphics g) {
int i, j, x, y;
int xUnit = deltaY / 9;
int yUnit = deltaY / 15;
super.paintComponent(g);
y = y0;
for (i = 0; i < 3; i++) {
x = x0;
for (j = 0; j < 3; j++) {
g.setColor(Color.white);
g.drawRect(x, y, deltaX, deltaY);
g.setColor(Color.black);
g.fillRect(x, y, deltaX, deltaY);
sevenSegmentDisplay(g, s[3 * i + j], x, y, deltaX, deltaY);
x += deltaX;
}
y += deltaY;
}
}
public void setSolution(char[] solution) {
s = new char[9];
s = solution;
}
}
class Puzzle8Frame extends JFrame implements Runnable {
boolean next;
int iHeight, iWidth;
JButton jButton1 = new JButton();
JPanel jPanel = new JPanel();
BorderLayout borderLayout = new BorderLayout();
Puzzle8Panel puzzle8Panel;
// step 3 - percentage size the window
void setDesktopSize(JFrame frame, int wPerc, int hPerc) {
Dimension screen = Toolkit.getDefaultToolkit().getScreenSize();
iWidth = screen.width * wPerc / 100;
iHeight = screen.height * hPerc / 100;
frame.setSize(iWidth, iHeight);
}
// step 4 - center the window
void centerOnScreen(JFrame frame) {
Dimension screen = Toolkit.getDefaultToolkit().getScreenSize();
Dimension window = frame.getSize();
int iCenterX = screen.width / 2;
int iCenterY = screen.height / 2;
frame.setLocation(iCenterX - window.width / 2, iCenterY - window.height / 2);
}
public Puzzle8Frame(char[] solution) {
String title = "Puzzle8 by James Pate Williams, Jr. (c) 2001";
next = false;
jButton1.setToolTipText("Perform one iteration of algorithm");
jButton1.setText("Next");
jButton1.setVerticalAlignment(SwingConstants.BOTTOM);
jButton1.addActionListener(new java.awt.event.ActionListener() {
public void actionPerformed(ActionEvent e) {
jButton1_actionPerformed(e);
}
});
this.getContentPane().setLayout(borderLayout);
setTitle(title);
addWindowListener(new WindowAdapter() {
public void windowClosing(WindowEvent event) {
System.exit(0);
}
});
setDesktopSize(this, 100, 100);
centerOnScreen(this);
Container contentPane = getContentPane();
contentPane.add(jPanel, BorderLayout.SOUTH);
jPanel.add(jButton1, BorderLayout.CENTER);
puzzle8Panel = new Puzzle8Panel(iWidth, iHeight, solution);
contentPane.add(puzzle8Panel, BorderLayout.CENTER);
this.show();
(new Thread(this)).run();
}
public boolean getNext() {
return next;
}
public void setNext(boolean n) {
next = n;
}
void jButton1_actionPerformed(ActionEvent e) {
next = true;
}
public void run() {
Thread.yield();
}
public void draw(char[] solution) {
puzzle8Panel.setSolution(solution);
puzzle8Panel.paintComponent(getGraphics());
}
}
class Puzzle {
int g, nodesExpanded;
int[][] board;
Date date;
Random random;
public static final int MaxMoves = 100;
public Puzzle() {
boolean found;
int digit, i, j, k;
int[] placed = new int[9];
date = new Date();
random = new Random(date.getTime());
for (i = 0; i < 9; i++)
placed[i] = 0;
board = new int[3][3];
g = nodesExpanded = 0;
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
board[i][j] = 0;
for (i = 0; i < 9; i++) {
found = false;
do {
digit = random.nextInt(9);
found = placed[digit] == 0;
if (found)
placed[digit] = 1;
} while (!found);
do {
j = random.nextInt(3);
k = random.nextInt(3);
found = board[j][k] == 0;
if (found)
board[j][k] = digit;
} while (!found);
}
}
int getNodesExpanded() {
return nodesExpanded;
}
int expand(int[][] square, int[][][] tempSquare) {
int b = - 1, col = - 1, i, j, k, row = - 1;
for (i = 0; i < 4; i++)
for (j = 0; j < 3; j++)
for (k = 0; k < 3; k++)
tempSquare[i][j][k] = square[j][k];
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
if (square[i][j] == 0) {
row = i;
col = j;
break;
}
}
}
if (row == 0 && col == 0) {
tempSquare[0][0][0] = tempSquare[0][0][1];
tempSquare[0][0][1] = 0;
tempSquare[1][0][0] = tempSquare[1][1][0];
tempSquare[1][1][0] = 0;
b = 2;
}
else if (row == 0 && col == 1) {
tempSquare[0][0][1] = tempSquare[0][0][0];
tempSquare[0][0][0] = 0;
tempSquare[1][0][1] = tempSquare[1][1][1];
tempSquare[1][1][1] = 0;
tempSquare[2][0][1] = tempSquare[2][0][2];
tempSquare[2][0][2] = 0;
b = 3;
}
else if (row == 0 && col == 2) {
tempSquare[0][0][2] = tempSquare[0][0][1];
tempSquare[0][0][1] = 0;
tempSquare[1][0][2] = tempSquare[1][1][2];
tempSquare[1][1][2] = 0;
b = 2;
}
else if (row == 1 && col == 0) {
tempSquare[0][1][0] = tempSquare[0][0][0];
tempSquare[0][0][0] = 0;
tempSquare[1][1][0] = tempSquare[1][1][1];
tempSquare[1][1][1] = 0;
tempSquare[2][1][0] = tempSquare[2][2][0];
tempSquare[2][2][0] = 0;
b = 3;
}
else if (row == 1 && col == 1) {
tempSquare[0][1][1] = tempSquare[0][1][0];
tempSquare[0][1][0] = 0;
tempSquare[1][1][1] = tempSquare[1][0][1];
tempSquare[1][0][1] = 0;
tempSquare[2][1][1] = tempSquare[2][1][2];
tempSquare[2][1][2] = 0;
tempSquare[3][1][1] = tempSquare[3][2][1];
tempSquare[3][2][1] = 0;
b = 4;
}
else if (row == 1 && col == 2) {
tempSquare[0][1][2] = tempSquare[0][0][2];
tempSquare[0][0][2] = 0;
tempSquare[1][1][2] = tempSquare[1][1][1];
tempSquare[1][1][1] = 0;
tempSquare[2][1][2] = tempSquare[2][2][2];
tempSquare[2][2][2] = 0;
b = 3;
}
else if (row == 2 && col == 0) {
tempSquare[0][2][0] = tempSquare[0][1][0];
tempSquare[0][1][0] = 0;
tempSquare[1][2][0] = tempSquare[1][2][1];
tempSquare[1][2][1] = 0;
b = 2;
}
else if (row == 2 && col == 1) {
tempSquare[0][2][1] = tempSquare[0][2][0];
tempSquare[0][2][0] = 0;
tempSquare[1][2][1] = tempSquare[1][1][1];
tempSquare[1][1][1] = 0;
tempSquare[2][2][1] = tempSquare[2][2][2];
tempSquare[2][2][2] = 0;
b = 3;
}
else if (row == 2 && col == 2) {
tempSquare[0][2][2] = tempSquare[0][2][1];
tempSquare[0][2][1] = 0;
tempSquare[1][2][2] = tempSquare[1][1][2];
tempSquare[1][1][2] = 0;
b = 2;
}
return b;
}
int heuristic(int[][] square) {
return ManhattenDistance(square);
}
boolean move() {
int b, count, i, j, k, min;
int[] f = new int[4], index = new int[4];
int[][][] tempSquare = new int[4][3][3];
if (board[0][0] == 1 && board[0][1] == 2 && board[0][2] == 3 &&
board[1][0] == 8 && board[1][1] == 0 && board[1][2] == 4 &&
board[2][0] == 7 && board[2][1] == 6 && board[2][2] == 5)
return true;
b = expand(board, tempSquare);
for (i = 0; i < b; i++) {
f[i] = g + heuristic(tempSquare[i]);
for (j = 0; j < 3; j++)
for (k = 0; k < 3; k++)
board[j][k] = tempSquare[i][j][k];
}
// find the node of minimum f
min = f[0];
for (i = 1; i < b; i++)
if (f[i] < min)
min = f[i];
for (count = i = 0; i < b; i++)
if (f[i] == min)
index[count++] = i;
i = index[random.nextInt(count)];
// increment the cost of the path thus far
g++;
nodesExpanded += b;
for (j = 0; j < 3; j++)
for (k = 0; k < 3; k++)
board[j][k] = tempSquare[i][j][k];
return false;
}
int outOfPlace(int[][] square) {
int i, j, oop = 0;
int[][] goal = new int[3][3];
goal[0][0] = 1;
goal[0][1] = 2;
goal[0][2] = 3;
goal[1][0] = 8;
goal[1][1] = 0;
goal[1][2] = 4;
goal[2][0] = 7;
goal[2][1] = 6;
goal[2][2] = 5;
for (i = 0; i < 3; i++)
for (j = 0; j < 3; j++)
if (square[i][j] != goal[i][j])
oop++;
return oop;
}
int ManhattenDistance(int[][] square) {
// city block or Manhatten distance heuristic
int md = 0;
if (square[0][0] == 1)
md += 0;
else if (square[0][0] == 2)
md += 1;
else if (square[0][0] == 3)
md += 2;
else if (square[0][0] == 4)
md += 3;
else if (square[0][0] == 5)
md += 4;
else if (square[0][0] == 6)
md += 3;
else if (square[0][0] == 7)
md += 2;
else if (square[0][0] == 8)
md += 1;
if (square[0][1] == 1)
md += 1;
else if (square[0][1] == 2)
md += 0;
else if (square[0][1] == 3)
md += 1;
else if (square[0][1] == 4)
md += 2;
else if (square[0][1] == 5)
md += 3;
else if (square[0][1] == 6)
md += 2;
else if (square[0][1] == 7)
md += 3;
else if (square[0][1] == 8)
md += 2;
if (square[0][2] == 1)
md += 2;
else if (square[0][2] == 2)
md += 1;
else if (square[0][2] == 3)
md += 0;
else if (square[0][2] == 4)
md += 1;
else if (square[0][2] == 5)
md += 2;
else if (square[0][2] == 6)
md += 3;
else if (square[0][2] == 7)
md += 4;
else if (square[0][2] == 8)
md += 3;
if (square[1][0] == 1)
md += 1;
else if (square[1][0] == 2)
md += 2;
else if (square[1][0] == 3)
md += 3;
else if (square[1][0] == 4)
md += 2;
else if (square[1][0] == 5)
md += 3;
else if (square[1][0] == 6)
md += 2;
else if (square[1][0] == 7)
md += 1;
else if (square[1][0] == 8)
md += 0;
if (square[1][1] == 1)
md += 2;
else if (square[1][1] == 2)
md += 1;
else if (square[1][1] == 3)
md += 2;
else if (square[1][1] == 4)
md += 1;
else if (square[1][1] == 5)
md += 2;
else if (square[1][1] == 6)
md += 1;
else if (square[1][1] == 7)
md += 2;
else if (square[1][1] == 8)
md += 1;
if (square[1][2] == 1)
md += 3;
else if (square[1][2] == 2)
md += 2;
else if (square[1][2] == 3)
md += 1;
else if (square[1][2] == 4)
md += 0;
else if (square[1][2] == 5)
md += 1;
else if (square[1][2] == 6)
md += 2;
else if (square[1][2] == 7)
md += 3;
else if (square[1][2] == 8)
md += 2;
if (square[2][0] == 1)
md += 2;
else if (square[2][0] == 2)
md += 3;
else if (square[2][0] == 3)
md += 4;
else if (square[2][0] == 4)
md += 3;
else if (square[2][0] == 5)
md += 2;
else if (square[2][0] == 6)
md += 1;
else if (square[2][0] == 7)
md += 0;
else if (square[2][0] == 8)
md += 1;
if (square[2][1] == 1)
md += 3;
else if (square[2][1] == 2)
md += 2;
else if (square[2][1] == 3)
md += 3;
else if (square[2][1] == 4)
md += 2;
else if (square[2][1] == 5)
md += 1;
else if (square[2][1] == 6)
md += 0;
else if (square[2][1] == 7)
md += 1;
else if (square[2][1] == 8)
md += 2;
if (square[2][2] == 1)
md += 4;
else if (square[2][2] == 2)
md += 3;
else if (square[2][2] == 3)
md += 2;
else if (square[2][2] == 4)
md += 1;
else if (square[2][2] == 5)
md += 0;
else if (square[2][2] == 6)
md += 1;
else if (square[2][2] == 7)
md += 2;
else if (square[2][2] == 8)
md += 3;
return md;
}
int solve(char[][] solution) {
boolean found;
int i, j, k, m = 0;
do {
for (i = k = 0; i < 3; i++)
for (j = 0; j < 3; j++)
solution[m][k++] = (char) (board[i][j] + '0');
found = move();
m++;
} while (!found && m < MaxMoves);
for (i = k = 0; i < 3; i++)
for (j = 0; j < 3; j++)
solution[m][k++] = (char) (board[i][j] + '0');
return m;
}
}
class Puzzle8 implements Runnable {
char[][] solution = null;
int moves;
public Puzzle8 () {
solution = new char[Puzzle.MaxMoves + 1][9];
do {
Puzzle puzzle = new Puzzle();
moves = puzzle.solve(solution);
} while (moves == Puzzle.MaxMoves);
}
public void run() {
Puzzle8Frame puzzle8Frame = new Puzzle8Frame(solution[0]);
System.out.println("moves = " + moves);
for (int i = 0; i < moves; i++) {
while (!puzzle8Frame.getNext())
Thread.yield();
puzzle8Frame.setNext(false);
puzzle8Frame.draw(solution[i]);
}
puzzle8Frame.draw(solution[moves]);
}
static void main(String[] arg) {
(new Puzzle8()).run();
}
}
/*
Author: James Pate Williams, Jr.
Stochastic L-System to draw a plant. See "Simulating Plant Growth" by
Marco Grubert, Crossroads, the ACM Student Magazine Issue 8.2, Winter
2001, pages 20-27.
*/
import java.awt.*;
import java.awt.event.*;
import java.lang.*;
import java.util.*;
import javax.swing.*;
class Cursor {
private double factor, rTheta;
private int height, theta, width, x, x0, y, y0;
public Cursor() {
}
public Cursor(int w, int h) {
width = w;
height = h;
factor = Math.PI / 180.0;
theta = 270;
rTheta = factor * theta;
x0 = width / 2;
y0 = height;
x = 0;
y = 0;
}
public Cursor(int iw, int ih, int iTheta, int ix, int iy) {
width = iw;
height = ih;
theta = iTheta;
factor = Math.PI / 180.0;
rTheta = factor * theta;
x0 = width / 2;
y0 = height;
x = ix;
y = iy;
}
public void F(Graphics g, int n) {
double nX = n * Math.cos(rTheta) + x, nY = n * Math.sin(rTheta) + y;
int X = x, Y = y;
x = (int) nX;
y = (int) nY;
g.drawLine(X + x0, Y + y0, x + x0, y + y0);
}
public void f(Graphics g, int n) {
double nX = n * Math.cos(rTheta) + x, nY = n * Math.sin(rTheta) + y;
x = (int) nX;
y = (int) nY;
g.translate(x + x0, y + y0);
}
public void plus(int alpha) {
theta = (theta - alpha) % 360;
if (theta < 0)
theta += 360;
rTheta = factor * theta;
}
public void minus(int alpha) {
theta = (theta + alpha) % 360;
if (theta < 0)
theta += 360;
rTheta = factor * theta;
}
public void exclamation() {
theta = (theta + 180) % 360;
if (theta < 0)
theta += 360;
rTheta = factor * theta;
}
public int getWidth() {
return width;
}
public int getHeight() {
return height;
}
public int getTheta() {
return theta;
}
public int getX() {
return x;
}
public int getY() {
return y;
}
public int getX0() {
return x0;
}
public int getY0() {
return y0;
}
}
class PlantPanel extends JPanel {
int depth;
Cursor cursor = null;
Random random = null;
public PlantPanel(int iDepth, int iWidth, int iHeight, int seed) {
depth = iDepth;
cursor = new Cursor(iWidth, iHeight);
random = new Random(seed);
}
public void lSystem1(Graphics g, int d, int depth, int n, int width, int height,
int theta, int x, int y) {
int alpha = 20;
Cursor cursor = new Cursor(width, height, theta, x, y);
Cursor stackCursor = null;
if (d == depth)
return;
cursor.F(g, n);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
stackCursor = new Cursor(cursor.getWidth(), cursor.getHeight(),
cursor.getTheta(), cursor.getX(), cursor.getY());
stackCursor.plus(alpha);
stackCursor.F(g, n);
theta = stackCursor.getTheta();
x = stackCursor.getX();
y = stackCursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
cursor.F(g, n);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
stackCursor = new Cursor(cursor.getWidth(), cursor.getHeight(),
cursor.getTheta(), cursor.getX(), cursor.getY());
stackCursor.minus(alpha);
stackCursor.F(g, n);
theta = stackCursor.getTheta();
x = stackCursor.getX();
y = stackCursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
cursor.F(g, n);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
}
public void lSystem2(Graphics g, int d, int depth, int n, int width, int height,
int theta, int x, int y) {
int alpha = 20;
Cursor cursor = new Cursor(width, height, theta, x, y);
Cursor stackCursor = null;
if (d == depth)
return;
cursor.F(g, n);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
stackCursor = new Cursor(cursor.getWidth(), cursor.getHeight(),
cursor.getTheta(), cursor.getX(), cursor.getY());
stackCursor.plus(alpha);
stackCursor.F(g, n);
theta = stackCursor.getTheta();
x = stackCursor.getX();
y = stackCursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
cursor.F(g, n);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
}
public void lSystem3(Graphics g, int d, int depth, int n, int width, int height,
int theta, int x, int y) {
int alpha = 20;
Cursor cursor = new Cursor(width, height, theta, x, y);
Cursor stackCursor = null;
if (d == depth)
return;
cursor.F(g, n);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
stackCursor = new Cursor(cursor.getWidth(), cursor.getHeight(),
cursor.getTheta(), cursor.getX(), cursor.getY());
stackCursor.minus(alpha);
stackCursor.F(g, n);
theta = stackCursor.getTheta();
x = stackCursor.getX();
y = stackCursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
}
public void lSystem0(Graphics g, int d, int depth, int n, int width, int height,
int theta, int x, int y) {
double r = random.nextDouble();
if (d == depth)
return;
if (r <= 0.33)
lSystem1(g, d, depth, n, width, height, theta, x, y);
else if (r <= 0.66)
lSystem2(g, d, depth, n, width, height, theta, x, y);
else
lSystem3(g, d, depth, n, width, height, theta, x, y);
theta = cursor.getTheta();
x = cursor.getX();
y = cursor.getY();
lSystem0(g, d + 1, depth, n / 2, width, height, theta, x, y);
}
public void lSystem(Graphics g, int d, int depth, int n, int width, int height,
int theta, int x, int y) {
lSystem0(g, d, depth, n, width, height, theta, x, y);
}
public void paintComponent(Graphics g) {
int n = 120;
int width = cursor.getWidth(), height = cursor.getHeight();
int theta = cursor.getTheta(), x0 = cursor.getX0(), y0 = cursor.getY0();
lSystem(g, 0, depth, n, width, height, theta, 0, 0);
}
}
class PlantFrame extends JFrame {
int iHeight, iWidth;
PlantPanel plantPanel;
// step 3 - percentage size the window
void setDesktopSize(JFrame frame, int wPerc, int hPerc) {
Dimension screen = Toolkit.getDefaultToolkit().getScreenSize();
iWidth = screen.width * wPerc / 100;
iHeight = screen.height * hPerc / 100;
frame.setSize(iWidth, iHeight);
}
// step 4 - center the window
void centerOnScreen(JFrame frame) {
Dimension screen = Toolkit.getDefaultToolkit().getScreenSize();
Dimension window = frame.getSize();
int iCenterX = screen.width / 2;
int iCenterY = screen.height / 2;
frame.setLocation(iCenterX - window.width / 2, iCenterY - window.height / 2);
}
public PlantFrame(int depth, int seed) {
String title = "Plant by James Pate Williams, Jr. (c) 2001";
setTitle(title);
addWindowListener(new WindowAdapter() {
public void windowClosing(WindowEvent event) {
System.exit(0);
}
});
setDesktopSize(this, 100, 100);
centerOnScreen(this);
Container contentPane = getContentPane();
plantPanel = new PlantPanel(depth, iWidth, iHeight, seed);
contentPane.add(plantPanel, BorderLayout.CENTER);
this.show();
}
}
class Plant {
public static void main(String[] args) {
if (args.length != 2) {
System.out.println("usage: java Plant depth seed");
System.out.println("where 5 <= depth <= 8 and seed >= 1");
System.exit(0);
}
int depth = (new Integer(args[0])). intValue();
if (depth < 5 || depth > 8) {
System.out.println("depth not in range 5 <= depth <= 8!");
System.exit(0);
}
int seed = (new Integer(args[1])).intValue();
if (seed < 1) {
System.out.println("seed must be >= 1");
System.exit(0);
}
PlantFrame plantFrame = new PlantFrame(depth, seed);
}
}
/*
Author: James Pate Williams, Jr. (c) 2000
Min-conflicts hill-climbing applied to N-Queens problem.
Wallace and Freuder's algorithm with an unbiased
preprocessor.
*/
#include <stdlib.h>
#include <time.h>
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;
const bool Debug = true;
const int MaxRepairs = 1000000;
int conflicts(int a, int i, int n, vector<int> &Q) {
int b, c = 0, j;
for (j = 0; j < n; j++) {
b = Q[j];
if (i != j && a != - 1 && b != - 1) {
if (a == b) c++;
if (i - j == a - b || i - j == b - a) c++;
}
}
return c;
}
int constraintsViolated(vector<int> &Q) {
int a, b, c = 0, i, j, n;
n = Q.size();
for (i = 0; i < n; i++) {
a = Q[i];
for (j = 0; j < n; j++) {
b = Q[j];
if (i != j && a != - 1 && b != - 1) {
if (a == b) c++;
if (i - j == a - b || i - j == b - a) c++;
}
}
}
return c;
}
int mchcAlgorithm(int n, vector<int> &compoundLabel) {
bool commit, firstLoop = true;
double p = 0.1, r;
int c, i, index, j, minConflict, v;
int repairs = 0;
vector<int> conflict, domain, position;
// preprocessing
for (i = 0; i < n; i++) {
conflict.erase(conflict.begin(), conflict.end());
for (j = 0; j < n; j++)
domain.push_back(j);
for (commit = false, j = 0; !commit && j < n; j++) {
v = domain[rand() % domain.size()];
vector<int>::iterator vIterator = find(domain.begin(), domain.end(), v);
domain.erase(vIterator);
compoundLabel[i] = v;
c = conflicts(v, i, i, compoundLabel);
conflict.push_back(c);
if (c == 0) commit = true;
}
if (!commit) {
minConflict = conflict[0];
for (j = 1; j < n; j++)
if (conflict[j] < minConflict) minConflict = conflict[j];
position.erase(position.begin(), position.end());
for (j = 0; j < n; j++)
if (minConflict == conflict[j]) position.push_back(j);
compoundLabel[i] = position[rand() % position.size()];
}
}
for (;;) {
conflict.erase(conflict.begin(), conflict.end());
for (i = 0, j = 0; i < n; i++) {
c = conflicts(compoundLabel[i], i, n, compoundLabel);
conflict.push_back(c);
j += c;
}
if (firstLoop) {
cout << "conflicts after preprocessor: " << j << endl;
firstLoop = false;
}
if (j == 0) break;
position.erase(position.begin(), position.end());
for (i = 0; i < n; i++)
if (conflict[i] >= 1) position.push_back(i);
index = position[rand() % position.size()];
r = (double) rand() / RAND_MAX;
if (r <= p)
compoundLabel[index] = rand() % n;
else {
for (i = 0; i < n; i++) {
compoundLabel[index] = i;
conflict[i] = conflicts(i, index, n, compoundLabel);
repairs++;
if (repairs == MaxRepairs) return repairs;
}
minConflict = conflict[0];
for (i = 1; i < n; i++)
if (conflict[i] < minConflict) minConflict = conflict[i];
position.erase(position.begin(), position.end());
for (i = 0; i < n; i++)
if (minConflict == conflict[i]) position.push_back(i);
compoundLabel[index] = position[rand() % position.size()];
}
}
return repairs;
}
void printSolution(vector<int> &solution) {
char hyphen[256];
int column, i, i4, n = solution.size(), row;
if (n <= 10) {
for (i = 0; i < n; i++) {
i4 = i * 4;
hyphen[i4 + 0] = '-';
hyphen[i4 + 1] = '-';
hyphen[i4 + 2] = '-';
hyphen[i4 + 3] = '-';
}
i4 = i * 4;
hyphen[i4 + 0] = '-';
hyphen[i4 + 1] = '\n';
hyphen[i4 + 2] = '\0';
for (row = 0; row < n; row++) {
column = solution[row];
cout << hyphen;
for (i = 0; i < column; i++)
cout << "| ";
cout << "| Q ";
for (i = column + 1; i < n; i++)
cout << "| ";
cout << '|' << endl;
}
cout << hyphen;
}
else
for (row = 0; row < n; row++)
cout << row << ' ' << solution[row] << endl;
}
int main(int argc, char *argv[]) {
if (argc != 2) {
cout << "usage: " << argv[0] << " numberOfQueens" << endl;
exit(1);
}
int n = atoi(argv[1]), repairs;
double seconds;
clock_t time0, time1;
vector<int> solution(n, - 1);
srand(time(NULL));
time0 = clock();
repairs = mchcAlgorithm(n, solution);
time1 = clock();
time1 -= time0;
seconds = (double) time1 / CLOCKS_PER_SEC;
if (Debug) printSolution(solution);
cout << "Repairs: " << repairs << endl;
cout << "seconds: " << seconds << endl;
return 0;
}
/*
* NQueensBruteForce.cpp
* NQueensBruteForce
*
* Brute force determination of the total number of
* solutions of the n-queens constraint satisfaction
* problem (CSP). Order n! complexity algorithm.
*
* Created by James Pate Williams, Jr. on 12/10/07.
* Copyright 2007. All rights reserved.
*
*/
#include "stdafx.h"
#include <stdlib.h>
#include <time.h>
#include <algorithm>
#include <functional>
#include <iostream>
#include <vector>
int ConstraintViolations(std::vector<int> &Q, int &cvCount) {
int a, b, cv = 0, i, j, n;
n = Q.size();
for (i = 0; i < n; i++) {
a = Q[i];
for (j = 0; j < n; j++) {
b = Q[j];
if (i != j && a != -1 && b != -1) {
if (a == b)
cv++;
if (i - j == a - b || i - j == b - a)
cv++;
}
}
}
cvCount++;
return cv;
}
void PrintSolution(std::vector<int> &solution) {
char hyphen[256];
int column, i, i4, n = solution.size(), row;
if (n <= 10) {
for (i = 0; i < n; i++) {
i4 = i * 4;
hyphen[i4 + 0] = '-';
hyphen[i4 + 1] = '-';
hyphen[i4 + 2] = '-';
hyphen[i4 + 3] = '-';
}
i4 = i * 4;
hyphen[i4 + 0] = '-';
hyphen[i4 + 1] = '\n';
hyphen[i4 + 2] = '\0';
for (row = 0; row < n; row++) {
column = solution[row];
std::cout << hyphen;
for (i = 0; i < column; i++)
std::cout << "| ";
std::cout << "| Q ";
for (i = column + 1; i < n; i++)
std::cout << "| ";
std::cout << '|' << std::endl;
}
std::cout << hyphen;
}
else
for (row = 0; row < n; row++)
std::cout << row << ' ' << solution[row] << std::endl;
}
int main(int argc, char * const argv[])
{
if (argc != 3) {
std::cout << "usage: " << argv[0] << " numberOfQueens print";
std::cout << std::endl;
exit(-1);
}
int numberOfQueens = atoi(argv[1]);
if (numberOfQueens <= 0) {
std::cout << "numberOfQueens must be positive" << std::endl;
exit(-1);
}
int count = 0, i, cvCount = 0, print = atoi(argv[2]), space = 0;
clock_t clock0;
std::vector<int> queenVector;
clock0 = clock();
for (i = 0; i < numberOfQueens; i++)
queenVector.push_back(i);
do {
int cv = ConstraintViolations(queenVector, cvCount);
if (cv == 0) {
count++;
if (print == 1) {
PrintSolution(queenVector);
std::cout << std::endl;
}
}
space++;
} while (std::next_permutation(
queenVector.begin(),
queenVector.end(),
std::less<int>()));
double seconds = (double)(clock() - clock0) / CLOCKS_PER_SEC;
std::cout << "search space size = " << space << std::endl;
std::cout << "number of con checks = " << cvCount << std::endl;
std::cout << "number of solutions = " << count << std::endl;
std::cout << "total time (seconds) = " << seconds << std::endl;
return 0;
}
#pragma once
// Insertion sort, heap sort, and quick sort algorithms
// From "Algortihms" by Thomas H. Cormen, et. al.
// Selection sort, Singleton's sort, and Tree Sort 3, etc.
// From "Sorting and Sort Systems" by Harold Lorin
class Sorting {
public:
// Translated from Lorin's PL/I procedure
static void BasicExchange(long long tosort[], int number);
// Translated from Lorin's PL/I procedure
static void StandardExchange(long long tosort[], int number);
static void InsertionSort(long long a[], int n);
static void SelectionSort(long long a[], int n);
// CACM Algorithm 175
static void SimpleShifting(long long tosort[], int number);
// CACM Algorithm 201 1963
static void ShellSort(long long tosort[], int number);
static void HeapSort(long long a[], int n);
static void QuickSort(long long a[], int n);
// CACM Algorithm 347 Ricahrd C. Singleton 1968
// Translated from Lorin's PL/I procedure
static void SingletonsSort(long long a[], int n);
// CACM Algorithm 245 Robert W. Floyd 1964
// Translated from Lorin's PL/I procedure
static void TreeSort3(long long a[], int n);
// see "Algortihms" by Thomas H. Cormen, et. al.
// A is the array to be sorted, B is the sorted
// array and k is maximum number in A
static void CountingSort(
long long A[],
long long B[],
long long C[],
int n, long long k);
private:
// five helper methods for Heap Sort
static int Parent(int i);
static int Left(int i);
static int Right(int i);
static void BuildHeap(long long a[], int n, int& heapSize);
static void Heapify(long long a[], int i, int n, int& heapSize);
/* helper functions for SETHEAP
static void SWOPHEAP(long long A[], int n, long long& in, long long& out);
static void INHEAP(long long A[], int& n, long long& in);
static void OUTHEAP(long long A[], int& n, long long& out);
static void SETHEAP(long long A[], int n);*/
// helper methods for Quick Sort
static int Partition(long long a[], int n, int lo, int hi);
static void DoQuickSort(long long a[], int n, int p, int r);
// helper function for Singleton's sort
static void DoSingletonsSort(long long a[], int ii, int jj);
// helper function for TreeSort3
static void SiftUp(long long a[], int i, int n);
};
#include "Sorting.h"
#include <math.h>
void Sorting::BasicExchange(long long tosort[], int number)
{
int adjust = 2 * (number / 2);
int elimit, olimit;
if (adjust < number)
{
elimit = adjust;
olimit = adjust - 1;
}
else
{
elimit = number - 2;
olimit = number - 1;
}
odd:
int passw = 1;
int limit = olimit;
int oddeve = 1;
pass:
int excount = 0;
for (int i = oddeve; i <= limit; i += 2)
{
if (tosort[i] > tosort[i + 1])
{
long long temp = tosort[i];
tosort[i] = tosort[i + 1];
tosort[i + 1] = temp;
excount = 1;
}
}
if (excount == 0)
return;
if (passw == 1)
{
passw = 0;
oddeve = 0;
limit = elimit;
goto pass;
}
goto odd;
}
void Sorting::StandardExchange(long long tosort[], int number)
{
pass:
int excount = 0;
for (int i = 0; i < number - 1; i++)
{
if (tosort[i] > tosort[i + 1])
{
long long temp = tosort[i];
tosort[i] = tosort[i + 1];
tosort[i + 1] = temp;
excount = 1;
}
}
if (excount == 0)
return;
goto pass;
}
void Sorting::InsertionSort(long long a[], int n)
{
for (int j = 1; j < n; j++)
{
long long key = a[j];
int i = j - 1;
while (i >= 0 && a[i] > key)
{
a[i + 1] = a[i];
i--;
}
a[i + 1] = key;
}
}
void Sorting::SelectionSort(long long a[], int n)
{
for (int i = 0; i < n - 1; i++)
{
for (int j = i + 1; j < n; j++)
{
if (a[i] > a[j])
{
long long t = a[i];
a[i] = a[j];
a[j] = t;
}
}
}
}
void Sorting::SimpleShifting(long long tosort[], int number)
{
for (int i = 0; i < number - 1; i++)
{
for (int j = i; j >= 0; j--)
{
if (tosort[j] > tosort[j + 1])
{
long long temp = tosort[j];
tosort[j] = tosort[j + 1];
tosort[j + 1] = temp;
}
}
}
}
void Sorting::ShellSort(long long tosort[], int number)
{
int lognumber = log2(number);
int distance = (int)pow(2, lognumber) - 1;
while (distance > 0)
{
int limit = number - distance;
for (int j = 0; j < limit; j++)
{
for (int i = j; i >= 0; i--)
{
if (tosort[i + distance] < tosort[i])
{
long long temp = tosort[i];
tosort[i] = tosort[i + distance];
tosort[i + distance] = temp;
}
}
}
distance /= 2;
}
}
int Sorting::Parent(int i)
{
return i / 2;
}
int Sorting::Left(int i)
{
return 2 * i;
}
int Sorting::Right(int i)
{
return 2 * i + 1;
}
void Sorting::Heapify(long long a[], int i, int n, int& heapSize)
{
int l = Left(i);
int r = Right(i);
int largest = -1;
if (l < heapSize && a[l] > a[i])
largest = l;
else
largest = i;
if (r < heapSize && a[r] > a[largest])
largest = r;
if (largest != i)
{
long long t = a[i];
a[i] = a[largest];
a[largest] = t;
Heapify(a, largest, n, heapSize);
}
}
void Sorting::BuildHeap(long long a[], int n, int& heapSize)
{
for (int i = n / 2; i >= 0; i--)
Heapify(a, i, n, heapSize);
}
void Sorting::HeapSort(long long a[], int n)
{
int heapSize = n;
BuildHeap(a, n, heapSize);
for (int i = n - 1; i >= 0; i--)
{
long long t = a[0];
a[0] = a[i];
a[i] = t;
heapSize--;
Heapify(a, 0, n, heapSize);
}
}
void SWOPHEAP(long long A[], int n,
long long& in, long long& out)
{
if (in <= A[0])
out = in;
else
{
int i = 0;
A[n] = in;
out = A[0];
scan:
int j = i + i;
if (j <= n - 1)
{
long long temp = A[j];
long long temp1 = A[j + 1];
if (temp1 < temp)
{
temp = temp1;
j++;
}
if (temp < in)
{
A[i] = temp;
i = j;
goto scan;
}
}
A[i] = in;
}
}
void INHEAP(long long A[], int& n, long long& in)
{
int i = n;
n++;
scan:
if (i > 0)
{
int j = i / 2;
if (in < A[j])
{
A[i] = A[j];
i = j;
goto scan;
}
}
A[i] = in;
}
void OUTHEAP(long long A[], int& n, long long& out)
{
SWOPHEAP(A, n - 1, A[n - 1], out);
n = n - 2;
}
void SETHEAP(long long A[], int n)
{
int j = 0;
L:
INHEAP(A, j, A[j + 1]);
if (j < n - 2)
goto L;
}
void HEAP(long long A[], int n)
{
SETHEAP(A, n);
for (int i = n - 1; i >= 1; i--)
OUTHEAP(A, i, A[i]);
}
int Sorting::Partition(long long a[], int n, int lo, int hi)
{
int pivotIndex = lo + (hi - lo) / 2;
long long x = a[pivotIndex];
long long t = x;
a[pivotIndex] = a[hi];
a[hi] = t;
int storeIndex = lo;
for (int i = lo; i < hi; i++)
{
if (a[i] < x)
{
t = a[i];
a[i] = a[storeIndex];
a[storeIndex++] = t;
}
}
t = a[storeIndex];
a[storeIndex] = a[hi];
a[hi] = t;
return storeIndex;
}
void Sorting::DoQuickSort(long long a[], int n, int p, int r)
{
if (p < r)
{
int q = Partition(a, n, p, r);
DoQuickSort(a, n, p, q - 1);
DoQuickSort(a, n, q + 1, r);
}
}
void Sorting::QuickSort(long long a[], int n)
{
DoQuickSort(a, n, 0, n - 1);
}
void QuickerSort(long long tosort[], int number)
{
const int limit = 20;
int highlim = 0, lowlim, origin, partind, pivot;
long long exch, temp = 0;
int partablow[limit] = { 0 }, parttabhigh[limit] = { 0 };
origin = partind = 0;
testsize:
if (number - 1 - origin > 0)
{
pivot = (number - 1 + origin) / 2;
temp = tosort[pivot];
tosort[pivot] = tosort[origin];
highlim = number - 1;
for (lowlim = origin + 1; lowlim <= highlim; lowlim++)
{
if (tosort[lowlim] > temp)
{
for (highlim = highlim; highlim >= lowlim; highlim--)
{
if (tosort[highlim] < temp)
{
exch = tosort[lowlim];
tosort[lowlim] = tosort[highlim];
tosort[highlim] = exch;
highlim = lowlim - 1;
goto limsmet;
}
}
highlim = lowlim - 1;
goto limsmet;
}
}
}
limsmet:
tosort[origin] = tosort[highlim];
tosort[highlim] = temp;
if (2 * highlim > origin + number - 1)
{
partablow[partind] = origin;
parttabhigh[partind] = highlim - 1;
origin = highlim + 1;
}
else
{
partablow[partind] = highlim + 1;
parttabhigh[partind] = number - 1;
number = highlim - 1;
}
partind++;
goto testsize;
if (origin == number - 1)
goto setpart;
if (tosort[origin] > tosort[number - 1])
{
exch = tosort[origin];
tosort[origin] = tosort[number - 1];
tosort[number - 1] = exch;
}
setpart:
partind--;
if (partind > -1)
{
origin = partablow[partind];
number = parttabhigh[partind];
goto testsize;
}
}
void Sorting::DoSingletonsSort(long long a[], int ii, int jj)
{
int m = 0;
int i = ii;
int j = jj;
int iu[50] = { 0 };
int il[50] = { 0 };
int ij = 0, l = 1, k = 0;
long long t = 0, tt = 0;
Label1:
if (i >= j)
goto Label8;
goto Label2;
Label2:
k = i;
ij = (j + i) / 2;
t = a[ij];
if (a[i] <= t)
goto Label3;
a[ij] = a[i];
a[i] = t;
t = a[ij];
goto Label3;
Label3:
l = j;
if (a[j] >= t)
goto Label5;
a[ij] = a[j];
a[j] = t;
t = a[ij];
if (a[i] <= t)
goto Label5;
a[ij] = a[i];
a[i] = t;
t = a[ij];
goto Label5;
Label4:
a[l] = a[k];
a[k] = tt;
goto Label5;
Label5:
l--;
if (a[l] > t)
goto Label5;
tt = a[l];
goto Label6;
Label6:
k++;
if (a[k] < t)
goto Label6;
if (k <= l)
goto Label4;
if (l - i <= j - k)
goto Label7;
il[m] = i;
iu[m] = l;
i = k;
m++;
goto Label9;
Label7:
il[m] = k;
iu[m] = j;
j = l;
m++;
goto Label9;
Label8:
m--;
if (m == -1)
return;
i = il[m];
j = iu[m];
goto Label9;
Label9:
if (j - i > 10)
goto Label2;
if (i == ii)
goto Label1;
i--;
goto Label10;
Label10:
i++;
if (i == j)
goto Label8;
t = a[i + 1];
if (a[i] <= t)
goto Label10;
k = i;
goto Label11;
Label11:
a[k + 1] = a[k];
k--;
if (t < a[k])
goto Label11;
a[k + 1] = t;
goto Label10;
}
void Sorting::SingletonsSort(long long a[], int n)
{
DoSingletonsSort(a, 0, n - 1);
}
void Sorting::SiftUp(long long a[], int i, int n)
{
long long copy = a[i];
int j;
while (true)
{
j = i + i;
if (j <= n)
{
if (j < n)
{
if (a[j + 1] > a[j])
j++;
}
if (a[j] > copy)
{
a[i] = a[j];
i = j;
continue;
}
}
break;
}
a[i] = copy;
}
void Sorting::TreeSort3(long long a[], int n)
{
for (int i = n / 2; i >= 1; i--)
SiftUp(a, i, n - 1);
for (int i = n - 1; i >= 1; i--)
{
SiftUp(a, 0, i);
long long t = a[0];
a[0] = a[i];
a[i] = t;
}
}
void Sorting::CountingSort(
long long A[],
long long B[],
long long C[],
int n, long long k)
{
for (long long i = 0; i <= k; i++)
C[i] = 0;
for (int j = 0; j < n; j++)
{
B[j] = 0;
C[A[j]]++;
}
for (long long i = 1; i <= k; i++)
C[i] += C[i - 1];
for (int i = n - 1; i >= 0; i--)
{
long long j = A[i];
C[j]--;
B[C[j]] = A[i];
}
}
#pragma once
#include <limits.h>
#include <vector>
using namespace std;
class SampleStatistics {
public:
inline static long long Mean(vector<long long> x) {
long long sum = 0;
for (int i = 0; i < (int)x.size(); i++)
sum += (long long)x[i];
return sum / (int)x.size();
};
inline static long long Max(vector<long long> x)
{
long long max = LLONG_MIN;
for (int i = 0; i < (int)x.size(); i++)
if (x[i] > max)
max = x[i];
return max;
};
inline static long long Min(vector<long long> x)
{
long long min = LLONG_MAX;
for (int i = 0; i < (int)x.size(); i++)
if (x[i] < min)
min = x[i];
return min;
};
inline static long long Median(vector<long long> x)
{
unsigned int n0 = (unsigned int)x.size();
unsigned int n1 = n0 / 2 - 1;
unsigned int n2 = n0 / 2;
long long median = 0, lt = x[n1], rt = x[n2];
if ((n0 & 1) == 1)
median = x[n0 / 2];
else
median = (lt + rt) / 2;
return median;
};
inline static long long Variance(vector<long long> x)
{
long long mean = Mean(x), sumSq = 0;
int n = (int)x.size();
int n1 = n - 1;
for (int i = 0; i < n; i++)
{
long long delta = x[i] - mean;
sumSq += delta * delta;
}
return sumSq / n1;
};
};
#include "SampleStatistics.h"
#include "Sorting.h"
#include <stdlib.h>
#include <algorithm>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std::chrono;
using namespace std;
// https://www.geeksforgeeks.org/measure-execution-time-function-cpp/
int GetSubOption(char title[])
{
char line[128];
cout << title << endl;
cout << "1 Basic Exchange" << endl;
cout << "2 Standard Exchange" << endl;
cout << "3 Insertion Sort" << endl;
cout << "4 Selection Sort" << endl;
cout << "5 Simple Sifting (Shuttle Sort)" << endl;
cout << "6 Shell Sort" << endl;
cout << "7 Heap Sort" << endl;
cout << "8 Quick Sort" << endl;
cout << "9 Singleton's Sort" << endl;
cout << "10 Tree Sort 3" << endl;
cout << "11 Counting Sort" << endl;
cout << "12 Exit Submenu" << endl;
cout << "13 Exit Program" << endl;
cout << "Suboption Number: ";
cin.getline(line, 128);
int so = 0;
if (strlen(line) > 1)
{
so = 10 * (line[0] - '0');
so += line[1] - '0';
}
else
so = line[0] - '0';
return so;
}
int DoSort(char option) {
char line[128];
char title[128] = { 0 };
char name[11][64] = { { 0 } };
int subOption;
strcpy_s(name[0], 64, "Basic Exchange");
strcpy_s(name[1], 64, "Standard Exchange");
strcpy_s(name[2], 64, "Insertion Sort");
strcat_s(name[3], 64, "Selection Sort");
strcpy_s(name[4], 64, "Simple Shifting");
strcpy_s(name[5], 64, "Shell Sort");
strcat_s(name[6], 64, "Heap Sort");
strcat_s(name[7], 64, "Quick Sort");
strcat_s(name[8], 64, "Singleton's Sort");
strcat_s(name[9], 64, "Tree Sort 3");
strcpy_s(name[10], 64, "Counting Sort");
vector<long long> runtime;
if (option == '1')
{
strcpy_s(title, 128, "Single Sorting Tests Submenu");
}
else {
strcpy_s(title, 128, "Statistical Sorting Tests Submenu");
}
subOption = GetSubOption(title);
if (subOption < 1 || subOption > 13)
{
cout << "Illegal Suboption Exiting Application!" << endl;
exit(-1);
}
if (subOption == 12)
return subOption;
if (subOption == 13)
exit(1);
int index = subOption - 1, maximum = 0, n = 0, seed = 1;
cout << name[index] << endl;
do {
cout << "PRNG Seed = ";
cin.getline(line, 128);
seed = atoi(line);
} while (seed < 1);
do {
cout << "Number of Samples = ";
cin.getline(line, 128);
n = atoi(line);
} while (n < 2);
do {
cout << "Maximum Sample Value = ";
cin.getline(line, 128);
maximum = atoi(line);
} while (maximum < 10);
srand(seed);
int number = 0;
if (option == '2')
{
do {
cout << "Number of Experiments = ";
cin.getline(line, 128);
number = atoi(line);
} while (number < 1);
}
else
{
number = 1;
}
long long k = 0;
long long* B = new long long[n];
long long* sample = new long long[n];
long long* shadow = new long long[n];
for (int i = 0; i < number; i++)
{
long long* C = NULL;
for (int j = 0; j < n; j++)
{
sample[j] = rand() % maximum;
if (option == '1') {
shadow[j] = sample[j];
}
}
if (subOption == 11)
{
for (int j = 0; j < n; j++)
if (sample[j] > k)
k = sample[j];
C = new long long[(unsigned int)(k + 1)];
}
auto start = high_resolution_clock::now();
switch (subOption) {
case 1:
Sorting::BasicExchange(sample, n);
break;
case 2:
Sorting::StandardExchange(sample, n);
break;
case 3:
Sorting::InsertionSort(sample, n);
break;
case 4:
Sorting::SelectionSort(sample, n);
break;
case 5:
Sorting::SimpleShifting(sample, n);
break;
case 6:
Sorting::ShellSort(sample, n);
break;
case 7:
Sorting::HeapSort(sample, n);
break;
case 8:
Sorting::QuickSort(sample, n);
break;
case 9:
Sorting::SingletonsSort(sample, n);
break;
case 10:
Sorting::TreeSort3(sample, n);
break;
case 11:
Sorting::CountingSort(sample, B, C, n, k);
break;
default:
cout << "Unknown Suboption" << endl;
break;
}
auto stop = high_resolution_clock::now();
auto duration = duration_cast<microseconds>(stop - start);
long long rt = (long long)duration.count();
runtime.push_back(rt);
if (subOption == 11)
{
if (C != NULL)
delete[] C;
}
}
if (option == '1')
{
for (int i = 0; i < n; i++)
{
cout << setw(2) << shadow[i] << ' ';
if (subOption != 11)
cout << setw(2) << sample[i] << endl;
else
cout << setw(2) << B[i] << endl;
}
}
else if (option == '2') {
sort(runtime.begin(), runtime.end());
cout << "Runtimes in Microseconds" << endl;
if (runtime.size() <= 15)
{
for (int i = 0; i < (int)runtime.size(); i++)
cout << setw(4) << runtime[i] << endl;
}
cout << "Minimum Runtime = "
<< setw(4) << fixed << setprecision(0)
<< SampleStatistics::Min(runtime) << endl;
cout << "Maximum Runtime = "
<< setw(4) << fixed << setprecision(0)
<< SampleStatistics::Max(runtime) << endl;
cout << "Mean Runtime = "
<< setw(4) << fixed << setprecision(0)
<< SampleStatistics::Mean(runtime) << endl;
cout << "Median Runtime = "
<< setw(4) << fixed << setprecision(0) <<
SampleStatistics::Median(runtime) << endl;
}
delete[] B;
delete[] sample;
delete[] shadow;
return subOption;
}
int main() {
while (true) {
char line[128];
cout << "Menu" << endl;
cout << "1 Single Sorting Tests" << endl;
cout << "2 Statistical Tests" << endl;
cout << "3 Exit" << endl;;
cout << "Option number: ";
cin.getline(line, 128);
if (line[0] == '3')
break;
if (line[0] == '1' || line[0] == '2') {
while (true)
{
int doSort = DoSort(line[0]);
if (doSort == 12)
break;
else if (doSort == 13)
exit(1);
}
}
}
return 0;
}
Numerical Explorations 