Numerical Explorations

Win32 C Real Double Precision and Integer Matrix Add and Subtract by James Pate Williams, Jr.

void DblAdd(int l, double* a, double* b, double* c)
{
	_asm
	{
		mov eax, a; load a matrix address
		mov ebx, b; load a matrix address
		mov edi, c; load result address
		mov edx, l; load matrix length
		finit;
	ll:	fld qword ptr[eax];
		fld qword ptr[ebx];
		fadd;
		fstp qword ptr[edi];
		fwait;
		add edi, 8;
		add eax, 8;
		add ebx, 8;
		dec edx;
		jnz ll;
	}
}

void DblSub(int l, double* a, double* b, double* c)
{
	_asm
	{
		mov eax, a; load a matrix address
		mov ebx, b; load a matrix address
		mov edi, c; load result address
		mov edx, l; load matrix length
		finit;
	ll:	fld qword ptr[eax];
		fld qword ptr[ebx];
		fsub;
		fstp qword ptr[edi];
		fwait;
		add edi, 8;
		add eax, 8;
		add ebx, 8;
		dec edx;
		jnz ll;
	}
}

Win32 C Integer Matrix Operations: Add and Sub Using Assembly Language by James Pate Williams, Jr.

void Add(int l, int* a, int* b, int* c)
{
	_asm
	{
		mov eax, a; load a vector address
		mov ebx, b; load a vector address
		mov edi, c; load result address
		mov edx, l; load vector length
   ll : mov ecx, 0;
		add ecx, [eax];
		add ecx, [ebx];
		mov[edi], ecx;
		add edi, 4;
		add eax, 4;
		add ebx, 4;
		add ecx, 4;
		dec edx;
		jnz ll;
	}
}

void Sub(int l, int* a, int* b, int* c)
{
	_asm
	{
		mov eax, a; load a vector address
		mov ebx, b; load a vector address
		mov edi, c; load result address
		mov edx, l; load vector length
   ll : mov ecx, 0;
		add ecx, [eax];
		sub ecx, [ebx];
		mov[edi], ecx;
		add edi, 4;
		add eax, 4;
		add ebx, 4;
		add ecx, 4;
		dec edx;
		jnz ll;
	}
}

C Matrix Operations: Add, Subtract, Multiply, and Determinant Calculation by James Pate Williams, Jr.

Long Integer Addition Using Three Wide Character Stacks by James Pate Williams, Jr.

// LongIntSum.cpp : Defines the entry point for the application.
// Translated from the Pascal code found in the textbook
// "Applied Data Structures Using Pascal" by Guy J. Hale and
// Richard J. Easton pages 98 and 99.

#include "stdafx.h"
#include "CharArrayStack.h"
#include "LongIntSum.h"

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    AddDialog(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_LONGINTSUM, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_LONGINTSUM));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, NULL, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex;

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_LONGINTSUM));
    wcex.hCursor        = LoadCursor(NULL, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_LONGINTSUM);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE:  Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
			case ID_FILE_LONGINTADDITION:
				DialogBox(hInst, MAKEINTRESOURCE(IDD_DIALOG1), hWnd, AddDialog);
				break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

#define BUFFER_MAX 8192

WCHAR buffer1[BUFFER_MAX];
WCHAR buffer2[BUFFER_MAX];
WCHAR buffer3[BUFFER_MAX];
CHARARRAYSTACK stack1;
CHARARRAYSTACK stack2;
CHARARRAYSTACK stack3;

void Add(
	HWND hwnd,
	int id1,
	int id2,
	int id3)
{
	WCHAR ch1, ch2, ch3;
	int carry = 0, i, i1, i2, i3, count = 0;

	GetDlgItemText(hwnd, id1, buffer1, BUFFER_MAX);
	GetDlgItemText(hwnd, id2, buffer2, BUFFER_MAX);

	for (i = 0; i < (int)wcslen(buffer1); i++)
		PushCharArrayStack(&stack1, buffer1[i]);

	for (i = 0; i < (int)wcslen(buffer2); i++)
		PushCharArrayStack(&stack2, buffer2[i]);

	while (
		!EmptyCharArrayStack(&stack1) &&
		!EmptyCharArrayStack(&stack2))
	{
		PopCharArrayStack(&stack1, &ch1);
		PopCharArrayStack(&stack2, &ch2);
		i1 = (int)(ch1 - '0');
		i2 = (int)(ch2 - '0');
		i3 = i1 + i2 + carry;
		carry = i3 / 10;
		i3 %= 10;
		ch3 = (char)(i3 + '0');
		PushCharArrayStack(&stack3, ch3);
	}

	while (!EmptyCharArrayStack(&stack1))
	{
		PopCharArrayStack(&stack1, &ch1);
		i1 = (int)(ch1 - '0');
		i3 = i1 + carry;
		carry = i3 / 10;
		i3 %= 10;
		ch3 = (char)(i3 + '0');
		PushCharArrayStack(&stack3, ch3);
	}

	while (!EmptyCharArrayStack(&stack2))
	{
		PopCharArrayStack(&stack2, &ch2);
		i2 = (int)(ch2 - '0');
		i3 = i2 + carry;
		carry = i3 / 10;
		i3 %= 10;
		ch3 = (char)(i3 + '0');
		PushCharArrayStack(&stack3, ch3);
	}

	if (carry != 0)
	{
		ch3 = (char)(carry + '0');
		PushCharArrayStack(&stack3, ch3);
	}

	while (!EmptyCharArrayStack(&stack3))
	{
		PopCharArrayStack(&stack3, &ch3);
		buffer3[count++] = ch3;
	}

	buffer3[count] = '\0';
	SetDlgItemText(hwnd, id3, buffer3);
}

// Message handler for long int addition dialog box.
INT_PTR CALLBACK AddDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
	UNREFERENCED_PARAMETER(lParam);
	WORD word = LOWORD(wParam);

	switch (message)
	{
	case WM_INITDIALOG:
		InitCharArrayStack(&stack1);
		InitCharArrayStack(&stack2);
		InitCharArrayStack(&stack3);
		return (INT_PTR)TRUE;
	case WM_COMMAND:
		if (word == IDOK || word == IDCANCEL)
		{
			EndDialog(hDlg, word);
			return (INT_PTR)TRUE;
		}

		if (word == IDC_BUTTON1)
			Add(hDlg, IDC_EDIT1, IDC_EDIT2, IDC_EDIT3);

		break;
	}
	return (INT_PTR)FALSE;
}

#pragma once
#include "stdafx.h"

#define STACK_SIZE 8192

typedef struct myCharArrayStack {
	WCHAR data[STACK_SIZE];
	int top;
} CHARARRAYSTACK, *PCHARARRAYSTACK;

void InitCharArrayStack(
	PCHARARRAYSTACK stack);

BOOL EmptyCharArrayStack(
	PCHARARRAYSTACK stack);

BOOL FullCharArrayStack(
	PCHARARRAYSTACK stack);

BOOL PushCharArrayStack(
	PCHARARRAYSTACK stack,
	WCHAR data);

BOOL PopCharArrayStack(
	PCHARARRAYSTACK stack,
	WCHAR* data);

#include "stdafx.h"
#include "CharArrayStack.h"

void InitCharArrayStack(
	PCHARARRAYSTACK stack)
{
	int i;

	for (i = 0; i < stack->top; i++)
		stack->data[i] = '\0';

	stack->data[0] = 0;
	stack->top = -1;
}

BOOL EmptyCharArrayStack(
	PCHARARRAYSTACK stack)
{
	return stack->top == -1 ? TRUE : FALSE;
}

BOOL FullCharArrayStack(
	PCHARARRAYSTACK stack)
{
	return stack->top == STACK_SIZE - 1 ? TRUE : FALSE;
}

BOOL PushCharArrayStack(
	PCHARARRAYSTACK stack,
	WCHAR data)
{
	BOOL result = FALSE;

	if (!FullCharArrayStack(stack))
	{
		stack->top++;
		stack->data[stack->top] = 0;
		stack->data[stack->top] = data;
		result = TRUE;
	}

	return result;
}

BOOL PopCharArrayStack(
	PCHARARRAYSTACK stack,
	WCHAR* data)
{
	BOOL result = FALSE;

	if (!EmptyCharArrayStack(stack))
	{
		*data = stack->data[stack->top];
		stack->top--;
		result = TRUE;
	}

	return result;
}

History of My Checkers Game Software Development Using Machine Learning by James Pate Williams, Jr.

My history of developing artificially intelligent checkers applications began in the Winter Quarter of 1999 at Auburn University. I was taking a Machine Learning course taught by Professor Gerry V. Dozier. We used the textbook “Machine Learning” by Tom M. Mitchell. The first chapter of the now classic textbook is devoted to developing a framework for a checkers (draughts) game. Mitchell used an evaluation function with seven weights and a least mean squares training rule. Over the years I have expanded the number of weights to fourteen rules. I seem to recall my early efforts were in C and later Java. I began programming checkers and chess applications on a Palm Pilot in the Summer Semester of 2002 in the course Handheld Software Development taught by my PhD research advisor Professor Richard O. Chapman. This work was performed in Palm Pilot Operating System C programming language. I created a client server set of programs to operate over TCP/IP networks. I had a modem for my Palm Pilot and learned TCP/IP sockets programming on the Palm Pilot. I still have two volumes addressing Palm Pilot C programming. I created a C# checkers program beginning on October 19, 2017. I started creating a Win32 version of my C# application on October 27, 2022. Unfortunately, I cannot video capture a Win32 game using the Windows button followed by the letter G for Game Boy. The program has five dialogs and apparently Windows-G buttons do not record dialogs. You can find a video of the computer playing against the computer on my Facebook page. I show the main dialog of the Win32 application in this blog.

Latest Software Development Project Started on Saturday, ‎October ‎15, ‎2022, ‏‎12:07:19 AM by James Pate Williams, Jr.

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.

Generation of the Algebraic Integer Primes

Pollard Cubic Integer Factoring and Classical Shor Factoring Algorithms by James Pate Williams, Jr.

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.

Pollard’s Factoring with Cubic Integers 2^63-1

Six Large Integer Factoring Methods Using Arjen K. Lenstra’s Free Large Integer Package (Free LIP) by James Pate Wiliams, Jr.

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;
}

Lenstra’s Free LIP Lenstra’s ECM Factorization of the Eighth Fermat Number

Lenstra’s Free LIP Lenstra’s ECM Factorization of the Seventh Fermat Number

Three Factoring Algorithms Implemented Using the C Free Large Integer Package (Free LIP) by James Pate Williams, Jr.

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.

Trial Division 50-Digit Number Missing Two or More Prime Factors

/*
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);
}

Lenstra’s Elliptic Curve Method of Factoring Large Integers Using the Free Large Integer Package (FreeLIP) by James Pate Williams, Jr.

/*
	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;
}