Find Complex Primitive Root of Eisenstein Number

ComplexPrimitiveRoot.c

#include <stdio.h>
#include <stdint.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <stdbool.h>
#include <time.h>
#include <math.h>
#include <gmp.h>
#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpzi.h>
#include <flint/fmpq.h>
#include <flint/fmpz_factor.h>
#define STB_DS_IMPLEMENTATION
#include "stb_ds.h"
#define INFINITE_LOOP_CHECK_SUM 1000
#define MPF_DEFAULT_PRECISION 512
// clear && gcc Eisenstein2.c -o m.o -lm -lgmp -lmpfr -lflint && ./m.o 

typedef struct discrete_log_system_struct *System;
typedef struct eisenstein_integer_struct *fmpzE;
struct eisenstein_integer_struct
{
	fmpz_t a;
	fmpz_t b;
	fmpz_t norm;
	int exponent;
};
struct discrete_log_system_struct
{
	fmpz_t integerPrime;
	fmpz_t integerPrimeMinusOne;
	fmpz_t integerPrimitiveRoot;
	fmpz_t rootOfUnity;
	fmpz_t twoInverse;
	fmpz_t heegnerNumber;
	fmpz_t heegnerNumberRoot;
	fmpz_factor_t normFactors;
	fmpzE *complexFactors;
	fmpzE complexPrime;
	fmpzE complexPrimitiveRoot;
		
};


fmpzE fmpzE_init()
{
	fmpzE e = malloc(sizeof(struct eisenstein_integer_struct));
	e->exponent = 0;
	fmpz_init(e->a);
	fmpz_init(e->b);
	fmpz_init(e->norm);
	return e;
}

void fmpzE_clear(fmpzE e)
{
	fmpz_clear(e->a);
	fmpz_clear(e->b);
	fmpz_clear(e->norm);
	free(e);
}

void fmpzE_add(fmpzE r, fmpzE a, fmpzE b)
{
	fmpz_add(r->a, a->a, b->a);
	fmpz_add(r->b, a->b, b->b);
}

void fmpzE_sub(fmpzE r, fmpzE a, fmpzE b)
{
	fmpz_sub(r->a, a->a, b->a);
	fmpz_sub(r->b, a->b, b->b);
}

void fmpzE_mul(fmpzE r, fmpzE a, fmpzE b)
{
	fmpz_t temp0;fmpz_init(temp0);
	fmpz_mul(r->a, a->a, b->a);
	fmpz_mul(temp0, a->b, b->b);
	fmpz_sub(r->a,r->a,temp0);
	
	
	fmpz_set(r->b, temp0);
	fmpz_neg(r->b, r->b);	
	fmpz_mul(temp0, a->a, b->b);
	fmpz_add(r->b,r->b,temp0);
	fmpz_mul(temp0, a->b, b->a);
	fmpz_add(r->b,r->b,temp0);
	fmpz_clear(temp0);
}

void fmpzE_conjugate(fmpzE result, fmpzE e)
{
	fmpz_sub(result->a, e->a, e->b);
	fmpz_neg(result->b, e->b);
}

void fmpzE_norm(fmpz_t norm, fmpz_t a, fmpz_t b)
{
	fmpz_t temp0;
	fmpz_init(temp0);
	fmpz_mul(temp0, a,a);//a^2
	fmpz_mul(norm, b,b);//b^2
	fmpz_add(norm, norm, temp0);
	
	fmpz_mul(temp0, a, b);//ab
	fmpz_sub(norm, norm, temp0);
	fmpz_clear(temp0);
}

void fmpz_set_mpf_round(fmpz_t result, mpf_t x)
{
	mpf_t f_tmp,frac;
	mpf_init(f_tmp);mpf_init(frac);
	
	mpf_floor(f_tmp, x);
	mpf_sub(frac, x, f_tmp);
	if(mpf_cmp_d(frac, 0.5) < 0) 
	{
		fmpz_set_mpf(result, f_tmp); 
	}
	else
	{
		mpf_ceil(f_tmp, x); 
		fmpz_set_mpf(result, f_tmp);
	}
	mpf_clear(f_tmp);mpf_clear(frac);
	
}

void fmpzE_set_fmpzE(fmpzE result, fmpzE e)
{
	fmpz_set(result->a, e->a);
	fmpz_set(result->b, e->b);
	fmpz_set(result->norm, e->norm);
}
void fmpzE_divide(fmpzE r, fmpzE a, fmpzE b)
{
	fmpzE bConjugate = fmpzE_init();
	fmpzE num = fmpzE_init();
	fmpzE_conjugate(bConjugate, b);
	
	fmpzE_mul(num, a, bConjugate);
	fmpzE_norm(b->norm, b->a, b->b);
	
	mpf_t f_NumA, f_NumB, f_norm;
	mpf_init(f_NumA);mpf_init(f_NumB);mpf_init(f_norm);

	fmpz_get_mpf(f_NumA, num->a);
	fmpz_get_mpf(f_NumB, num->b);
	fmpz_get_mpf(f_norm, b->norm);

	mpf_div(f_NumA, f_NumA, f_norm);
	mpf_div(f_NumB, f_NumB, f_norm);

	
	fmpz_set_mpf_round(r->a, f_NumA);
	fmpz_set_mpf_round(r->b, f_NumB);
	
	mpf_clear(f_NumA);mpf_clear(f_NumB);mpf_clear(f_norm);
	fmpzE_clear(num);
	fmpzE_clear(bConjugate);
}

void fmpzE_mod(fmpzE r, fmpzE a, fmpzE b)
{
	fmpzE q = fmpzE_init();
	fmpzE temp = fmpzE_init();
	
	fmpzE_divide(q, a, b);  
	fmpzE_mul(temp, q, b);   
	fmpzE_sub(r, a, temp);  

	fmpzE_clear(q);
	fmpzE_clear(temp);
}

void fmpzE_gcd(fmpzE result, fmpzE alpha, fmpzE beta)
{
	fmpzE a = fmpzE_init();
	fmpzE b = fmpzE_init();
	fmpzE r = fmpzE_init();
	fmpzE zero = fmpzE_init();
	fmpz_set(a->a, alpha->a);fmpz_set(a->b, alpha->b);
	fmpz_set(b->a, beta->a);fmpz_set(b->b, beta->b);
	fmpz_zero(zero->a);fmpz_zero(zero->b);
	int infiniteLoopCheck = 0;
	while(!(fmpz_equal(b->a, zero->a) && fmpz_equal(b->b, zero->b)))
	{
		//r = a % b
		fmpzE_mod(r, a, b);

		//a = b
		fmpz_set(a->a, b->a);fmpz_set(a->b, b->b);

		//b = r
		fmpz_set(b->a, r->a);fmpz_set(b->b, r->b);
    		
    		if(infiniteLoopCheck >= INFINITE_LOOP_CHECK_SUM)
		{
			printf("Infinite loop in fmpzE_gcd\n");
			assert(infiniteLoopCheck < INFINITE_LOOP_CHECK_SUM);
		}
		infiniteLoopCheck += 1;
    	
    	}

	//result = a
	fmpz_set(result->a, a->a);
	fmpz_set(result->b, a->b);

	fmpzE_clear(a);
	fmpzE_clear(b);
	fmpzE_clear(r);
	fmpzE_clear(zero);
}

void fmpzE_powm_fmpzE(fmpzE result, fmpzE alpha, fmpz_t exponent, fmpzE modulo)
{
	//Case 0^0 = 0
	if(fmpz_is_zero(alpha->a) && fmpz_is_zero(alpha->b) && fmpz_is_zero(exponent))
	{
		fmpz_set_ui(result->a, 0);
		fmpz_set_ui(result->b, 0);
	}
	else
	{
		fmpzE temp0 = fmpzE_init();
		//Initialize result to 1
		fmpz_set_ui(result->a, 1);fmpz_set_ui(result->b, 0);
		int size = fmpz_sizeinbase(exponent, 2);
		for(int i = size - 1; i >= 0; i--)
		{
			fmpzE_mul(temp0, result, result);	
			fmpzE_set_fmpzE(result, temp0);
			
			fmpzE_mod(temp0, result, modulo);	
			fmpzE_set_fmpzE(result, temp0);
			if(fmpz_tstbit(exponent, i))
			{
				fmpzE_mul(temp0, result, alpha);	
				fmpzE_set_fmpzE(result, temp0);
			
				fmpzE_mod(temp0, result, modulo);	
				fmpzE_set_fmpzE(result, temp0);
			}
		}
		fmpzE_clear(temp0);
	}
}

void fmpzE_FindRandomPrimitiveRoot(System system)
{
	fmpz_t temp0;
	fmpz_init(temp0);
	fmpzE result = fmpzE_init();
	flint_rand_t state;
	flint_rand_init(state);
	bool isPrimitiveRoot = false;
	int infiniteLoopCheck = 0;
	while(isPrimitiveRoot == false)
	{
		//Generate random candidate
		fmpz_randm(system->complexPrimitiveRoot->a, state, system->integerPrime); 
		//fmpz_randm(system->complexPrimitiveRoot->b, state, system->integerPrime); 
		
		//Test random candidate
		isPrimitiveRoot = true;
		for(slong i = 0; i < system->normFactors->num; i++)
		{
			fmpz_divexact(temp0, system->integerPrimeMinusOne, &system->normFactors->p[i]);
			fmpzE_powm_fmpzE(result, system->complexPrimitiveRoot, temp0,  system->complexPrime);
			fmpzE_norm(result->norm, result->a, result->b);
			if(fmpz_cmp_ui(result->norm, 1) == 0)
			{
				isPrimitiveRoot = false;
				break;
			}
		}
		if(infiniteLoopCheck >= INFINITE_LOOP_CHECK_SUM)
		{
			printf("Infinite loop in fmpzE_gcd\n");
			assert(infiniteLoopCheck < INFINITE_LOOP_CHECK_SUM);
		}
		infiniteLoopCheck += 1;

	}

	fmpzE_clear(result);
	fmpz_clear(temp0);
	flint_rand_clear(state);
}

void fmpzE_GCDSplit(fmpz_t prime, fmpz_t heegnerNumber, fmpzE complexFactor)
{
	fmpz_t temp2;
	fmpz_init(temp2);
	int foundSquareRoot = fmpz_sqrtmod(temp2, heegnerNumber, prime);
	assert(foundSquareRoot == 1);
	fmpzE pE = fmpzE_init();
	fmpzE gcdTemp = fmpzE_init();
	fmpzE divTemp = fmpzE_init();
	fmpz_set(pE->a, prime);fmpz_set_ui(pE->b, 0);
	fmpz_add_ui(gcdTemp->a, temp2, 1);fmpz_set_ui(gcdTemp->b, 2);

	//Set complex factor to gcd result
	fmpz_set(pE->a, prime);fmpz_set_ui(pE->b, 0);
	fmpzE_gcd(complexFactor, gcdTemp, pE);
	
	fmpzE_clear(pE);fmpzE_clear(gcdTemp);fmpzE_clear(divTemp);fmpz_clear(temp2);
}
 
void fmpzE_print(fmpzE e)
{
	printf("(");fmpz_print(e->a);printf(" + ");fmpz_print(e->b);printf(")\n");
}

void fmpzE_printPrimeFactors(fmpzE *primeFactors)
{
	for(size_t i = 0; i < arrlen(primeFactors); i++)
	{
		printf("(");fmpz_print(primeFactors[i]->a);printf(" + ");fmpz_print(primeFactors[i]->b);printf(") ^ %d, ",primeFactors[i]->exponent);	
	}
}

void fmpzE_freePrimeFactors(fmpzE *primeFactors)
{
	for(size_t i = 0; i < arrlen(primeFactors); i++)
	{
		fmpzE_clear(primeFactors[i]);
	}
	arrfree(primeFactors);
}



bool TestPrimitiveRoot(fmpz_t primitiveRoot, fmpz_t primeMinusOne, fmpz_t prime, fmpz_factor_t factors)
{
	fmpz_t temp0;
	fmpz_init(temp0);
	bool isPrimitiveRoot = true;
	for(slong i = 0; i < factors->num; i++)
	{
		fmpz_divexact(temp0, primeMinusOne, &factors->p[i]);
		fmpz_powm(temp0, primitiveRoot, temp0, prime);
		if(fmpz_cmp_ui(temp0, 1) == 0)
		{
			isPrimitiveRoot = false;
			break;
		}
	}
	fmpz_clear(temp0);
	return isPrimitiveRoot;
}


bool FindTV_Cornacchia(fmpz_t T, fmpz_t V, fmpz_t heegner, fmpz_t root, fmpz_t prime)
{
	//Solve T^2 + |D| * V^2 = p
	fmpz_t a0, a1, remainder, rootP,D,rhs;
	fmpz_init(a0);fmpz_init(a1);fmpz_init(remainder);fmpz_init(D);fmpz_init(rootP);fmpz_init(rhs);
	
	//Find square root of prime
	fmpz_root(rootP, prime, 2);
	fmpz_set(a0, prime);
	fmpz_set(a1, root);
	//Ensure we are working with a nageative heegner number
	assert(fmpz_cmp_ui(heegner, 0) < 0);
	fmpz_mul_si(D, heegner, -1);
	
	while(fmpz_cmp(a1, rootP) > 0)
	{
		//Find remainder = a0 mod a1
		//Ensure a1 != 0
		if(fmpz_cmp_ui(a1, 0) == 0){break;}
		fmpz_mod(remainder, a0, a1);
		fmpz_set(a0, a1);
		fmpz_set(a1, remainder);
	}
	bool a1not0 = (fmpz_cmp_ui(a1, 0) != 0);
	bool isRhsSquare = false;
	bool isRhsDivisibleByD = false;
	//no solution if a1 = 0
	if(a1not0 == true)
	{
		assert(fmpz_cmp_ui(a1, 0) != 0);
		fmpz_set(T, a1);
		//Compute rhs = (p - T^2) / |D|
		fmpz_mul(rhs, T, T);
		fmpz_sub(rhs, prime, rhs);
		
		//rhs must be divisible by |D|
		isRhsDivisibleByD = fmpz_divisible(rhs, D);
		
		if(isRhsDivisibleByD)	
		{
			fmpz_fdiv_q(rhs, rhs, D);
			//rhs must be a perfect square
			isRhsSquare = fmpz_is_square(rhs);
			if(isRhsSquare)
			{
				fmpz_sqrt(V, rhs);
			}
		}
	}

	
	fmpz_clear(a0);fmpz_clear(a1);fmpz_clear(remainder);fmpz_clear(D);fmpz_clear(rootP);fmpz_clear(rhs);
	return (isRhsSquare == true) && (isRhsDivisibleByD == true) && (a1not0 == true);
}

void FindRandomIntegerPrimitiveRoot(System system, fmpz_t primitiveRoot,fmpz_factor_t factors, fmpz_t prime)
{
	fmpz_t primeMinusOne,temp0;
	flint_rand_t state;
	flint_rand_init(state);
	fmpz_init(primeMinusOne);
	fmpz_init(temp0);
	fmpz_sub_ui(primeMinusOne, prime, 1);
	fmpz_factor(factors, primeMinusOne);
	
	bool isPrimitiveRoot = false;
	fmpz_set_si(temp0, -3);
	fmpz_set_ui(system->integerPrimitiveRoot, 2);
	while(isPrimitiveRoot == false)
	{
		fmpz_nextprime(system->integerPrimitiveRoot, system->integerPrimitiveRoot, 0);
		isPrimitiveRoot = TestPrimitiveRoot(system->integerPrimitiveRoot, primeMinusOne, prime, factors);
		if(isPrimitiveRoot == true)
		{
			int legendre = fmpz_jacobi(temp0, system->integerPrimitiveRoot);	
			if(legendre != 1)
			{
				isPrimitiveRoot = false;
			}	
			else
			{
				fmpzE_GCDSplit(system->integerPrimitiveRoot, system->heegnerNumber, system->complexPrimitiveRoot);
				fmpzE_norm(system->complexPrimitiveRoot->norm, system->complexPrimitiveRoot->a, system->complexPrimitiveRoot->b);
				if(fmpz_cmp(system->complexPrimitiveRoot->norm, system->integerPrimitiveRoot) != 0)
				{
					isPrimitiveRoot = false;
				}
			}
		}
	}

	fmpz_clear(temp0);
	fmpz_clear(primeMinusOne);
	flint_rand_clear(state);
}



System CreateSystem(fmpz_t prime, fmpz_factor_t factors)
{
	size_t mpfPrecision = fmpz_sizeinbase(prime, 2) * 10;
	mpf_set_default_prec(mpfPrecision);
	assert(factors->num > 0);
	
	fmpz_t temp0, temp1, temp2;
	fmpz_init(temp0);fmpz_init(temp1);fmpz_init(temp2);
	System system = malloc(sizeof(struct discrete_log_system_struct));
	/*Initialize fmpz_t*/			
	fmpz_init(system->integerPrime);
	fmpz_init(system->integerPrimeMinusOne);
	fmpz_init(system->integerPrimitiveRoot);
	fmpz_init(system->heegnerNumber);
	fmpz_init(system->heegnerNumberRoot);
	fmpz_init(system->rootOfUnity);
	fmpz_init(system->twoInverse);
	fmpz_factor_init(system->normFactors);
	/*Initialize fmpzE*/			
	system->complexPrime = fmpzE_init();
	system->complexPrimitiveRoot = fmpzE_init();
	system->complexFactors = NULL;
	
	/*Ensure prime splits*/			
	fmpz_set_si(system->heegnerNumber, -3);
	int legendre = fmpz_jacobi(system->heegnerNumber, prime);
	if(legendre != 1)
	{
		printf("Error: prime does not split root -3\n");
		assert(legendre == 1);			
	}
	/*Set p, p-1 and 2 inverse*/			
	fmpz_set(system->integerPrime, prime);
	fmpz_sub_ui(system->integerPrimeMinusOne, prime,1);
	fmpz_set_ui(system->twoInverse, 2);
	fmpz_invmod(system->twoInverse, system->twoInverse, system->integerPrime);
	
	/*Find a non-trivial root of unity*/	
	int foundSquareRoot = fmpz_sqrtmod(system->heegnerNumberRoot, system->heegnerNumber, system->integerPrime);
	assert(foundSquareRoot == 1);
	
	//Try either heegnerNumberRoots in roots of unity to find T and V
	bool TV_solutionExists = false;
	int solutionSearch = 0;
	while(TV_solutionExists == false && solutionSearch < 2)
	{
		//Set root of unity
		fmpz_add_si(system->rootOfUnity, system->heegnerNumberRoot, -1);
		fmpz_mul(system->rootOfUnity, system->rootOfUnity, system->twoInverse);
		fmpz_mod(system->rootOfUnity, system->rootOfUnity, system->integerPrime);
		/*Find T and V, we use norm as a temp var*/
		fmpz_mul_ui(system->complexPrime->norm,system->integerPrime, 4);
		//Solve for U^2 + 3V^2 = 4p
		TV_solutionExists = FindTV_Cornacchia(system->complexPrime->a, system->complexPrime->b, system->heegnerNumber, system->heegnerNumberRoot, system->complexPrime->norm);
		
		if(TV_solutionExists == false)
		{
			//Try other heegnerNumberRoot
			fmpz_sub(system->heegnerNumberRoot,system->integerPrime, system->heegnerNumberRoot);	
		}
		solutionSearch += 1;
	}
	
	assert(TV_solutionExists == true);
	//Solve for T and norm
	fmpz_add(system->complexPrime->a, system->complexPrime->a, system->complexPrime->b);
	fmpz_divexact_ui(system->complexPrime->a, system->complexPrime->a, 2);
	fmpzE_norm(system->complexPrime->norm, system->complexPrime->a, system->complexPrime->b);
	assert(fmpz_cmp(system->complexPrime->norm, system->integerPrime) == 0);
	
	/*Set factors of p-1*/
	fmpz_set_ui(temp0, 1);
	for(int i = 0; i < factors->num; i++)
	{
		fmpz_powm(temp1, &factors->p[i], &factors->exp[i], system->integerPrime);
		fmpz_mul(temp0, temp0, temp1);
		assert(factors->exp[i] == 1);
		int mod = fmpz_mod_ui(temp1, &factors->p[i], 3);
		_fmpz_factor_append(system->normFactors, &factors->p[i], legendre);	
		
		if(mod == 1)
		{
			//Find eisenstein form
			fmpzE alpha = fmpzE_init();
			fmpzE_GCDSplit(&factors->p[i], system->heegnerNumber, alpha);
			fmpzE_norm(alpha->norm,alpha->a,alpha->b);
			//printf("Mu: ");fmpzE_print(alpha); printf("N: ");fmpz_print(alpha->norm);printf(" | ");fmpz_print(&factors->p[i]);printf("\n");
			arrput(system->complexFactors, alpha);
		}
	}
	if(fmpz_cmp(temp0, system->integerPrimeMinusOne) != 0)
	{
		printf("Manual factors don't match prime - 1: ");
		fmpz_print(temp0);
		printf("\n");
		assert(fmpz_cmp(temp0, system->integerPrimeMinusOne) == 0);
	}
	
	/*Set primitive roots*/
	FindRandomIntegerPrimitiveRoot(system, system->integerPrimitiveRoot, system->normFactors, system->integerPrime);
	
	fmpz_clear(temp0);fmpz_clear(temp1);fmpz_clear(temp2);
	return system;
}

void PrintSystem(System system)
{
	printf("Integer Prime: ");fmpz_print(system->integerPrime);printf("\n");
	printf("Integer Prime-1: ");fmpz_print(system->integerPrimeMinusOne);printf("\n");
	printf("Integer Primitive Root: ");fmpz_print(system->integerPrimitiveRoot);printf("\n");
	printf("Norm factors:\n");
	for(slong i = 0; i < system->normFactors->num; i++)
	{
		printf("(");
		fmpz_print(&system->normFactors->p[i]); printf(" ^ ");
		fmpz_print(&system->normFactors->exp[i]); printf(")");
	}
	printf("\n");
	printf("Complex factors:\n");
	for(size_t i = 0; i < arrlen(system->complexFactors); i++)
	{
		printf("(");fmpz_print(system->complexFactors[i]->norm);printf(") ");
		fmpzE_print(system->complexFactors[i]);
	}
	printf("\n");
	printf("Heegner Number: ");fmpz_print(system->heegnerNumber);printf("\n");
	printf("RootOfUnity (S): ");fmpz_print(system->rootOfUnity);printf("\n");
	printf("T: ");fmpz_print(system->complexPrime->a);printf("\n");
	printf("V: ");fmpz_print(system->complexPrime->b);printf("\n");
	printf("Complex prime: (");fmpz_print(system->complexPrime->a);printf(" + ");fmpz_print(system->complexPrime->b);printf("√-3) Norm = ");fmpz_print(system->complexPrime->norm);printf("\n");
	printf("Complex PrimitiveRoot: (");fmpz_print(system->complexPrimitiveRoot->a);printf(" + ");fmpz_print(system->complexPrimitiveRoot->b);printf("√-3) Norm = ");fmpz_print(system->complexPrimitiveRoot->norm);printf("\n");
}

void DestroySystem(System system)
{
	for(size_t i = 0; i < arrlen(system->complexFactors); i++)
	{
		fmpzE_clear(system->complexFactors[i]);
	}
	arrfree(system->complexFactors);
	fmpz_factor_clear(system->normFactors);
	fmpz_clear(system->heegnerNumberRoot);
	fmpz_clear(system->twoInverse);
	fmpz_clear(system->rootOfUnity);
	fmpz_clear(system->integerPrimeMinusOne);
	fmpz_clear(system->integerPrime);
	fmpz_clear(system->integerPrimitiveRoot);
	fmpz_clear(system->heegnerNumber);
	fmpzE_clear(system->complexPrime);
	fmpzE_clear(system->complexPrimitiveRoot);
	free(system);
}

void TestSmallSystem()
{
	fmpz_factor_t factors;
	fmpz_t prime, primeMinusOne;
	fmpz_init(prime);fmpz_init(primeMinusOne);
	fmpz_factor_init(factors);
	fmpz_set_ui(prime, 20959);
	fmpz_sub_ui(primeMinusOne, prime, 1);
	fmpz_factor(factors, primeMinusOne);
	

	System system = CreateSystem(prime, factors);
	PrintSystem(system);

	

	DestroySystem(system);
	fmpz_factor_clear(factors);
	fmpz_clear(prime);fmpz_clear(primeMinusOne);	
}

void TestBitcoinSystem()
{
	char *primeNumberHexadecimal = "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F";	
	char *factor4Base10 = "205115282021455665897114700593932402728804164701536103180137503955397371";	
	fmpz_factor_t factors;
	fmpz_t prime, primeMinusOne,temp0;
	fmpz_init(prime);fmpz_init(primeMinusOne);fmpz_init(temp0);
	fmpz_factor_init(factors);
	fmpz_set_str(prime, primeNumberHexadecimal, 16);
	
	fmpz_set_ui(temp0,2);_fmpz_factor_append(factors, temp0, 1);	
	fmpz_set_ui(temp0,3);_fmpz_factor_append(factors, temp0, 1);	
	fmpz_set_ui(temp0,7);_fmpz_factor_append(factors, temp0, 1);		
	fmpz_set_ui(temp0,13441);_fmpz_factor_append(factors, temp0, 1);	
	fmpz_set_str(temp0,factor4Base10, 10);_fmpz_factor_append(factors, temp0, 1);	

	
	System system = CreateSystem(prime, factors);
	PrintSystem(system);
	fmpz_factor_clear(factors);
	fmpz_clear(prime);fmpz_clear(primeMinusOne);fmpz_clear(temp0);
	DestroySystem(system);
}


int main()
{
	TestSmallSystem();
	TestBitcoinSystem();
	flint_cleanup();
	return 0;
}