Lenstra's Solving pell equations with index calculus example code

Lenstra.c

//Full guide : https://leetarxiv.substack.com/p/solving-pell-equations-with-index
#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <stdbool.h>
#include <math.h>
#include <gmp.h>
#define STB_DS_IMPLEMENTATION
#include "stb_ds.h"
#include "ff_asm_primes.h"
#define PRIMES_COUNT 5000
//clear && gcc Lenstra.c -lm -lgmp -o m.o && ./m.o
typedef struct smooth_number_struct *SmoothNumber;
struct smooth_number_struct
{
	mpz_t x;
	mpz_t y;
	mpz_t D;
	mpz_t norm;
	int normSign;
	int *primes;
	int *powers;
};

double Find_Log_MPZ_Double(mpz_t x, double base){if(mpz_cmp_ui(x, 0) == 0){return 0;}signed long int ex;const double di = mpz_get_d_2exp(&ex, x);return ( (log(di) + log(base) * (double) ex) /log(base));}
double PomeranceSmoothnessBound(mpz_t N)
{
	//Find log2 of N
	double logBase2 = Find_Log_MPZ_Double(N, 2);	
	
	//Find natural log
	double ln_N = logBase2 * log(2.0);
	double ln_ln_N = log(ln_N);
	
	//Find bound
	double B = exp(sqrt(ln_N * ln_ln_N));
	return B;
}

double OdlyzkoBound(mpz_t N)
{
	double logBase2 = Find_Log_MPZ_Double(N, 2);	
	double c = 1 / (4 * log(2));
	c = sqrt(c);
	double B = c * sqrt(logBase2 * log(logBase2));
	B = pow(2, B);
	return B;
}

SmoothNumber CreateSmoothNumber()
{
	SmoothNumber smoothNumber = malloc(sizeof(struct smooth_number_struct));
	mpz_init(smoothNumber->x);
	mpz_init(smoothNumber->y);
	mpz_init(smoothNumber->D);
	mpz_init(smoothNumber->norm);
	smoothNumber->normSign = 0;
	smoothNumber->primes = NULL;
	smoothNumber->powers = NULL;
	return smoothNumber;
}

void DestroySmoothNumber(SmoothNumber smoothNumber)
{
	if(smoothNumber)
	{
		mpz_clear(smoothNumber->x);
		mpz_clear(smoothNumber->y);
		mpz_clear(smoothNumber->D);
		mpz_clear(smoothNumber->norm);
		arrfree(smoothNumber->primes);
		arrfree(smoothNumber->powers);
		free(smoothNumber);
	}
}

int SmoothNumber_LegendreSymbol(mpz_t a, mpz_t p)
{
	if(mpz_cmp_ui(p,2) == 0)
	{
	        if(mpz_even_p(a)){return 0;}	
		mpz_t mod8; mpz_init(mod8);mpz_mod_ui(mod8, a, 8);
		int res;
		switch(mpz_get_ui(mod8))
		{
			case 1: case 7: res =  1; break;
			case 3: case 5: res = -1; break;
			default: res = 0; // Shouldn't happen since a is odd
		}
		mpz_clear(mod8);
		return res;
	}
	else
	{
		return mpz_legendre(a, p);
	}
}

bool SmoothNumber_Factorize(mpz_t integer, mpz_t D, int bSmooth, int **primes, int **powers)
{
	//Factorize while ensuring legendre check
	mpz_t remainder, primeHolder;mpz_inits(remainder,primeHolder,NULL);
	mpz_set(remainder, integer);
	bool isSmooth = true;

	for(int i = 0; i < PRIMES_COUNT; i++)
	{
		int primeNumber = first5239[i];
		//Stop if primeNumber is smooth
		if(primeNumber > bSmooth){break;}
		int count = 0;
		while(mpz_divisible_ui_p(remainder, primeNumber))
		{
			mpz_divexact_ui(remainder, remainder, primeNumber);
			count++;
		}
		if(count > 0)
		{
			//Check legendre symbol == 1
			mpz_set_ui(primeHolder,primeNumber);
			if(SmoothNumber_LegendreSymbol(D, primeHolder) == 1)
			{
				arrput(*primes, primeNumber);
				arrput(*powers, count);
			}
			else
			{
				isSmooth = false;
				break;
			}
			//gmp_printf("%lu^%d ", primeNumber, count);
		}

		// Early exit if fully factored
		if(mpz_cmp_ui(remainder, 1) == 0){break;}
	}
	//Check if remaining factor is 1 or smooth
	if(isSmooth && mpz_cmp_ui(remainder, 1) != 0)
	{
		if(mpz_cmp_ui(remainder, bSmooth) <= 0)
		{
			//gmp_printf("%Zd^1 ", remainder);
			arrput(*primes,  mpz_get_ui(remainder));
			arrput(*powers, 1);
		}
		else
		{
			//gmp_printf("[non-smooth factor: %Zd] ", remainder);
			isSmooth = false;
		}
	}
	mpz_clears(remainder,primeHolder,NULL);
	return isSmooth;
}

void SmoothNumber_FindNorm(SmoothNumber smoothNumber)
{
	mpz_t temp;
	mpz_inits(temp, NULL);
	//Set norm to x^2
	mpz_mul(smoothNumber->norm, smoothNumber->x, smoothNumber->x);
	//Set temp to y^2
	mpz_mul(temp, smoothNumber->y, smoothNumber->y);
	//multiply temp by D
	mpz_mul(temp, smoothNumber->D, temp);
	
	//Set norm to x^2 - Dy^2
	mpz_sub(smoothNumber->norm, smoothNumber->norm, temp);
	
	//Find sign of norm
	smoothNumber->normSign = mpz_sgn(smoothNumber->norm);
	//Set norm to it's absolute value
	mpz_abs(smoothNumber->norm, smoothNumber->norm);
	
	mpz_clears(temp, NULL);
}

void PrintSmoothNumber(SmoothNumber smoothNumber)
{
	gmp_printf("(%Zd + %Zd√%Zd): Norm: %d * %Zd : ",smoothNumber->x,smoothNumber->y, smoothNumber->D,smoothNumber->normSign, smoothNumber->norm);
	for(size_t i = 0; i < arrlen(smoothNumber->primes); i++)
	{
		printf("(%2d ^%2d)",smoothNumber->primes[i],smoothNumber->powers[i]);
	}
	printf("\n");
}

void TestSmoothNumber()
{
	SmoothNumber smoothNumber = CreateSmoothNumber();
	mpz_set_ui(smoothNumber->x, 4351);
	mpz_set_ui(smoothNumber->y, 2);
	mpz_set_ui(smoothNumber->D, 4729494);
	SmoothNumber_FindNorm(smoothNumber);
	
	int smoothnessFactor = 50;
	bool isSmooth = SmoothNumber_Factorize(smoothNumber->norm, smoothNumber->D, smoothnessFactor, &(smoothNumber->primes), &(smoothNumber->powers));
	if(isSmooth)
	{
		PrintSmoothNumber(smoothNumber);
	}

	DestroySmoothNumber(smoothNumber);
}

void TestLegendreSymbol()
{
	mpz_t D, p;
	mpz_inits(D,p, NULL);
	//Set D
	mpz_set_ui(D, 4729494);
	for(int i = 0; i < 15; i++)
	{
		mpz_set_ui(p, first5239[i]);
		if(SmoothNumber_LegendreSymbol(D, p) == 1)
		{
			printf("%d, ", first5239[i]);
		}
	}
	printf("\n");
	mpz_clears(D,p,NULL);
}


SmoothNumber* SmoothNumber_FindSmoothNumbers(mpz_t D, int bSmooth, int maxFoundCount, int maxY, int maxX)
{
	SmoothNumber *acceptedNumbers = NULL;
	mpz_t sqrtD, ySquareD, xSearch,xSearchBound,gcd;
	mpz_inits(sqrtD, ySquareD, xSearch, xSearchBound, gcd, NULL);
	
	//Set squareRoot of D and variables
	int foundCount = 0;
	mpz_sqrt(sqrtD, D);
	for(int y = 1; y < maxY; y++)
	{
		//Set xSearchBound to ySquareD
		mpz_mul_ui(xSearchBound, sqrtD, y);
		//Set ySquareD
		mpz_mul_ui(ySquareD, D, y);
		mpz_mul_ui(ySquareD, ySquareD, y);
		for(int x = -maxX; x < maxX; x++)
		{
			//Search x^2 values close to ySquareD
			mpz_set(xSearch, xSearchBound);
			if(x < 0)
			{
				mpz_sub_ui(xSearch, xSearch, abs(x));
			}
			else
			{
				mpz_add_ui(xSearch, xSearch, x);	
			}

			//Ensure x is coprime to y
			mpz_gcd_ui(gcd, xSearch, y);
			if(mpz_cmp_ui(gcd, 1) == 0)
			{
				//gmp_printf("\t%Zd ^2 - %d^2 * %Zd = ",xSearch, y, D);
				mpz_mul(xSearch, xSearch, xSearch);
				mpz_sub(xSearch, xSearch, ySquareD);
				
				//Find absolute value
				int sign = mpz_sgn(xSearch);
				//Set norm to it's absolute value
				mpz_abs(xSearch, xSearch);
				//gmp_printf("%Zd\n",xSearch);
		
				//Test for smoothness
				int *primes = NULL;
				int *powers = NULL;
				bool isSmooth = SmoothNumber_Factorize(xSearch, D, bSmooth, &primes, &powers);
				if(isSmooth)
				{
					SmoothNumber smoothNumber = CreateSmoothNumber();
					mpz_set(smoothNumber->x, xSearch);
					mpz_set_ui(smoothNumber->y, y);
					mpz_set(smoothNumber->D, D);
					SmoothNumber_FindNorm(smoothNumber);
					smoothNumber->primes = primes;
					smoothNumber->powers = powers;
					arrput(acceptedNumbers, smoothNumber);
				}
				else
				{
					arrfree(primes);
					arrfree(powers);
				}
			}
		}
	}
	
	mpz_clears(sqrtD, ySquareD, xSearch,xSearchBound,gcd,NULL);
	return acceptedNumbers;
}

void PrintMatrixForm(SmoothNumber *smoothNumbers)
{
	
}
void TestGenerateSmoothNumbers()
{
	mpz_t D;mpz_inits(D, NULL);
	mpz_set_ui(D, 4729494);
		
	int smoothness = 50;
	int maxFoundCount = 20;
	int maxY = 100;
	int maxX = 200;
	SmoothNumber *smoothNumbers = SmoothNumber_FindSmoothNumbers(D, smoothness,maxFoundCount,maxY,maxX);
	//for(size_t i = 0; i < arrlen(smoothNumbers); i++){PrintSmoothNumber(smoothNumbers[i]);printf("\n");	}
	
	PrintMatrixForm(smoothNumbers);
	
	for(size_t i = 0; i < arrlen(smoothNumbers); i++){DestroySmoothNumber(smoothNumbers[i]);}
	arrfree(smoothNumbers);
	mpz_clears(D,NULL);
}



void TestSmoothnessBound()
{
	mpz_t D;mpz_inits(D, NULL);
	mpz_set_ui(D, 4729494);
		
	double smoothnessBound = PomeranceSmoothnessBound(D);
	double smoothnessBoun = OdlyzkoBound(D);
	printf("Rough B:%.3f\n Better B:%.3f \n",smoothnessBound,smoothnessBoun);
	mpz_clears(D,NULL);
}

int main()
{
	//TestLegendreSymbol();
	//TestSmoothNumber();
	//TestGenerateSmoothNumbers();
	TestSmoothnessBound();
	return 0;
}