Skip to content

Instantly share code, notes, and snippets.

@MurageKibicho

MurageKibicho/Eisenstein.c

Created December 6, 2025 12:12
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save MurageKibicho/3e4c5394d36caaa5e7a08f19dff86d8a to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save MurageKibicho/3e4c5394d36caaa5e7a08f19dff86d8a to your computer and use it in GitHub Desktop.
Download ZIP
Arithmetic with Eisenstein numbers
Raw
Eisenstein.c
//https://github.com/RasmusFL/EisensteinIntegers/blob/main/EisensteinIntegers/eisenstein_integers.cpp
#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 100000
int HeegnerNumbers[] = {-1, -2, -3, -7, -11, -19, -43, -67, -163};
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;
};
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_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_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);
}
struct discrete_log_system_struct
{
fmpz_t prime;
fmpz_t primeMinusOne;
fmpz_t primitiveRoot;
fmpz_t heegnerNumberField;
fmpz_t heegnerNumberRoot;
fmpz_t T;
fmpz_t V;
int B0;
int B1;
};
double PomeranceSmoothnessBound(fmpz_t N)
{
//Find log2 of N
double logBase2 = fmpz_sizeinbase(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;
}
void PrintFactors(fmpz_factor_t factors)
{
for(slong i = 0; i < factors->num; i++)
{
int size = fmpz_sizeinbase(&factors->p[i], 2);
//printf("(|%d| ", size);
fmpz_print(&factors->p[i]); printf(" ^ ");
fmpz_print(&factors->exp[i]); printf(", ");
}
printf("\n");
}
void 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);
}
//no solution if a1 = 0
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|
assert(fmpz_divisible(rhs, D));
fmpz_fdiv_q(rhs, rhs, D);
//rhs must be a perfect square
assert(fmpz_is_square(rhs));
fmpz_sqrt(V, rhs);
fmpz_clear(a0);fmpz_clear(a1);fmpz_clear(remainder);fmpz_clear(D);fmpz_clear(rootP);fmpz_clear(rhs);
}
int FindLargestFactorSize(fmpz_factor_t factors, int *largePrimeCount, int bound)
{
int largest = 0;
*largePrimeCount = 0;
for(slong i = 0; i < factors->num; i++)
{
int size = fmpz_sizeinbase(&factors->p[i], 2);
if(size > largest)
{
largest = size;
}
if(size > bound)
{
*largePrimeCount += 1;
}
//fmpz_print(&factors->p[i]); printf(" ^ ");
//fmpz_print(&factors->exp[i]); printf(")");
}
return largest;
}
void GaussReduce(fmpz_t x1, fmpz_t y1, fmpz_t x2, fmpz_t y2)
{
//Shortest solution is stored in x1,y1
fmpz_t mu, dot12, dot11, temp0,tmp2,tmp3;
fmpz_init(mu);fmpz_init(tmp2);fmpz_init(tmp3); fmpz_init(dot12); fmpz_init(dot11); fmpz_init(temp0);
int infiniteLoopCheck = 0;
while(1)
{
//mu = round((v1·v2)/(v1·v1))
fmpz_mul(dot12, x1, x2);fmpz_mul(temp0, y1, y2);fmpz_add(dot12, dot12, temp0);
fmpz_mul(dot11, x1, x1);fmpz_mul(temp0, y1, y1);fmpz_add(dot11, dot11, temp0);
fmpz_tdiv_q(mu, dot12, dot11);
// v2 = v2 - mu * v1
fmpz_mul(temp0, mu, x1); fmpz_sub(x2, x2, temp0);
fmpz_mul(temp0, mu, y1); fmpz_sub(y2, y2, temp0);
// ||v2|| < ||v1||, swap
fmpz_mul(temp0, x1, x1);
fmpz_mul(tmp2, y1, y1); fmpz_add(temp0, temp0, tmp2);
fmpz_mul(tmp2, x2, x2);
fmpz_mul(tmp3, y2, y2); fmpz_add(tmp2, tmp2, tmp3);
if(fmpz_cmp(tmp2, temp0) < 0){fmpz_swap(x1, x2);fmpz_swap(y1, y2);}
else{break;}
infiniteLoopCheck += 1;
if(infiniteLoopCheck > INFINITE_LOOP_CHECK_SUM)
{
printf("Gauss reduce may be in infinite loop\n");
assert(infiniteLoopCheck < INFINITE_LOOP_CHECK_SUM);
}
}
fmpz_clear(mu);fmpz_clear(tmp2);fmpz_clear(tmp3); fmpz_clear(dot12); fmpz_clear(dot11); fmpz_clear(temp0);
}
System CreateSystem(fmpz_t primitiveRoot, fmpz_t primeNumber)
{
System system = malloc(sizeof(struct discrete_log_system_struct));
fmpz_init(system->prime);
fmpz_init(system->heegnerNumberField);
fmpz_init(system->heegnerNumberRoot);
fmpz_init(system->primitiveRoot);
fmpz_init(system->primeMinusOne);
fmpz_init(system->T);
fmpz_init(system->V);
fmpz_set(system->prime, primeNumber);
fmpz_set(system->primitiveRoot, primitiveRoot);
fmpz_sub_ui(system->primeMinusOne, system->prime, 1);
//Find number field
bool foundField = false;
for(size_t i = 0; i < sizeof(HeegnerNumbers) / sizeof(int); i++)
{
fmpz_set_si(system->heegnerNumberField, HeegnerNumbers[i]);
int legendreSymbol = fmpz_jacobi(system->heegnerNumberField, system->prime);
if(legendreSymbol == 1)
{
foundField = true;
break;
}
}
assert(foundField == true);
int foundSquareRoot = fmpz_sqrtmod(system->heegnerNumberRoot, system->heegnerNumberField, system->prime);
assert(foundSquareRoot == 1);
FindTV_Cornacchia(system->T, system->V, system->heegnerNumberField, system->heegnerNumberRoot, system->prime);
system->B0 = (int) ceil(log2(PomeranceSmoothnessBound(system->prime)));
system->B1 = 1.5 * system->B0;
return system;
}
void DestroySystem(System system)
{
fmpz_clear(system->T);
fmpz_clear(system->V);
fmpz_clear(system->heegnerNumberRoot);
fmpz_clear(system->heegnerNumberField);
fmpz_clear(system->primeMinusOne);
fmpz_clear(system->prime);
fmpz_clear(system->primitiveRoot);
free(system);
}
void PrintSystem(System system)
{
printf("Prime: ");fmpz_print(system->prime);printf("\n");
printf("Primitive Root: ");fmpz_print(system->primitiveRoot);printf("\n");
printf("Heegner Number(r): ");fmpz_print(system->heegnerNumberField);printf("\n");
printf("Heegner Root(S): ");fmpz_print(system->heegnerNumberRoot);printf("\n");
printf("T: ");fmpz_print(system->T);printf("\n");
printf("V: ");fmpz_print(system->V);printf("\n");
printf("Bound0: %3d bits\n", system->B0);
printf("Bound1: %3d bits\n", system->B1);
}
fmpzE *PrimeFactorization(fmpzE alpha)
{
fmpzE *primeFactors = NULL;
fmpzE pE = fmpzE_init();
fmpzE gcdTemp = fmpzE_init();
fmpzE divTemp = fmpzE_init();
fmpzE alphaCopy = fmpzE_init();
fmpz_t mod, neg3,neg3Root;
fmpz_factor_t normFactors;fmpz_factor_init(normFactors);
fmpz_init(mod);fmpz_init(neg3);fmpz_init(neg3Root);
fmpz_set_si(neg3, -3);
fmpzE_norm(alpha->norm, alpha->a, alpha->b);
fmpz_factor(normFactors, alpha->norm);
//Set alphaCopy
fmpz_set(alphaCopy->a, alpha->a);fmpz_set(alphaCopy->b, alpha->b);fmpz_set(alphaCopy->norm, alpha->norm);
//printf("Norm: ");PrintFactors(normFactors);
for(slong i = 0; i < normFactors->num; i++)
{
fmpzE complexFactor = fmpzE_init();fmpzE complexConjugate = fmpzE_init();
fmpz_set(pE->a, &normFactors->p[i]);fmpz_set_ui(pE->b, 0);
if(fmpz_cmp_ui(&normFactors->p[i], 3) == 0)
{
//if p[i] = 3, 1 - w is a factor
fmpz_set_ui(complexFactor->a, 1);
fmpz_set_si(complexFactor->b, -1);
complexFactor->exponent = normFactors->exp[i];
for(int j = 0; j < normFactors->exp[i]; j++)
{
fmpzE_divide(divTemp, alphaCopy, complexFactor);
fmpz_set(alphaCopy->a,divTemp->a);fmpz_set(alphaCopy->b,divTemp->b);
}
}
else
{
int mod3 = fmpz_mod_ui(mod, &normFactors->p[i], 3);
assert(mod3 != 0);
if(mod3 == 2)
{
//if p[i] % 3 == 2, remove two factors of p[i]
fmpz_set(complexFactor->a, &normFactors->p[i]);
fmpz_set_ui(complexFactor->b, 0);
complexFactor->exponent = normFactors->exp[i] / 2;
for(int j = 0; j < normFactors->exp[i]; j += 2)
{
fmpzE_divide(divTemp, alphaCopy, complexFactor);
fmpz_set(alphaCopy->a,divTemp->a);fmpz_set(alphaCopy->b,divTemp->b);
}
}
else
{
int foundSquareRoot = fmpz_sqrtmod(neg3Root, neg3, &normFactors->p[i]);
assert(foundSquareRoot == 1);
fmpz_add_ui(gcdTemp->a, neg3Root, 1);fmpz_set_ui(gcdTemp->b, 2);
//Set complex factor to gcd result
fmpzE_gcd(complexFactor, gcdTemp, pE);
fmpzE_norm(complexFactor->norm,complexFactor->a,complexFactor->b);
//Find Complex conjugate
fmpzE_conjugate(complexConjugate, complexFactor);
for(int j = 0; j < normFactors->exp[i]; j++)
{
//Test if complex factor divides alphaCopy
fmpz_mod(mod, alphaCopy->norm,complexFactor->norm);
if(fmpz_cmp_ui(mod, 0) == 0)
{
fmpzE_divide(divTemp, alphaCopy, complexFactor);
fmpz_set(alphaCopy->a,divTemp->a);fmpz_set(alphaCopy->b,divTemp->b);
complexFactor->exponent += 1;
}
else
{
fmpzE_divide(divTemp, alphaCopy, complexConjugate);
fmpz_set(alphaCopy->a,divTemp->a);fmpz_set(alphaCopy->b,divTemp->b);
complexConjugate->exponent += 1;
}
}
}
}
if(complexConjugate->exponent > 0){arrput(primeFactors, complexConjugate);}else{fmpzE_clear(complexConjugate);}
if(complexFactor->exponent > 0){arrput(primeFactors, complexFactor);}else{fmpzE_clear(complexFactor);}
}
//Push alphaCopy as final factor if not unit
fmpzE_norm(alphaCopy->norm, alphaCopy->a, alphaCopy->b);
alphaCopy->exponent = 1;
arrput(primeFactors, alphaCopy);
fmpz_clear(mod);fmpz_clear(neg3);fmpz_clear(neg3Root);
fmpz_factor_clear(normFactors);
fmpzE_clear(pE);fmpzE_clear(gcdTemp);fmpzE_clear(divTemp);
return primeFactors;
}
bool fmpzE_FactorsNormCompare(fmpzE e, fmpzE *primeFactors)
{
//Find target norm
fmpzE_norm(e->norm,e->a,e->b);
for(size_t i = 0; i < arrlen(primeFactors); i++)
{
fmpzE_norm(primeFactors[i]->norm, primeFactors[i]->a, primeFactors[i]->b);
fmpz_pow_ui(primeFactors[i]->norm, primeFactors[i]->norm, primeFactors[i]->exponent);
fmpz_divexact(e->norm, e->norm, primeFactors[i]->norm);
}
return fmpz_cmp_ui(e->norm,1) == 0;
}
void TestFactors()
{
fmpzE test0 = fmpzE_init();
//x = 10 + 2w
fmpz_set_ui(test0->a, 10);fmpz_set_ui(test0->b, 2);
fmpzE_print(test0);
fmpzE *primeFactors0 = PrimeFactorization(test0);
fmpzE_printPrimeFactors(primeFactors0);
bool comp0 = fmpzE_FactorsNormCompare(test0,primeFactors0);
assert(comp0 == true);
//x = 5 + 16w
fmpz_set_ui(test0->a, 5);fmpz_set_ui(test0->b, 16);
fmpzE_print(test0);
fmpzE *primeFactors1 = PrimeFactorization(test0);
fmpzE_printPrimeFactors(primeFactors1);
bool comp1 = fmpzE_FactorsNormCompare(test0,primeFactors1);
assert(comp1 == true);
//X = 6-8w
fmpz_set_ui(test0->a, 6);fmpz_set_si(test0->b, -8);
fmpzE_print(test0);
fmpzE *primeFactors2 = PrimeFactorization(test0);
fmpzE_printPrimeFactors(primeFactors2);
bool comp2 = fmpzE_FactorsNormCompare(test0,primeFactors2);
assert(comp2 == true);
//X = -4+2w
fmpz_set_si(test0->a, -4);fmpz_set_si(test0->b, 2);
fmpzE_print(test0);
fmpzE *primeFactors3 = PrimeFactorization(test0);
fmpzE_printPrimeFactors(primeFactors3);
bool comp3 = fmpzE_FactorsNormCompare(test0,primeFactors3);
assert(comp3 == true);
fmpzE_freePrimeFactors(primeFactors0);
fmpzE_freePrimeFactors(primeFactors1);
fmpzE_freePrimeFactors(primeFactors2);
fmpzE_clear(test0);
}
void TestFMPZEArithmetic()
{
//Division test
fmpzE alpha = fmpzE_init();
fmpzE beta = fmpzE_init();
fmpzE result= fmpzE_init();
fmpz_set_si(alpha->a, -2);fmpz_set_ui(alpha->b, 6);
fmpz_set_ui(beta->a, 3);fmpz_set_ui(beta->b, 4);
fmpzE_divide(result, alpha, beta);
assert(fmpz_cmp_ui(result->a, 2) == 0);assert(fmpz_cmp_ui(result->b, 2) == 0);
printf("Division Test Pass: ");fmpzE_print(result);
//GCD Test
fmpz_set_si(alpha->a, 34);fmpz_set_si(alpha->b, -21);
fmpz_set_ui(beta->a, 17);fmpz_set_si(beta->b, 8);
fmpzE_gcd(result, alpha, beta);
assert(fmpz_cmp_si(result->a, -1) == 0);assert(fmpz_cmp_ui(result->b, 0) == 0);
printf("GCD Test0 Pass: ");fmpzE_print(result);
fmpz_set_si(alpha->a, 24);fmpz_set_si(alpha->b, 57);
fmpz_set_si(beta->a, -42);fmpz_set_si(beta->b, -12);
fmpzE_gcd(result, alpha, beta);
assert(fmpz_cmp_si(result->a, -6) == 0);assert(fmpz_cmp_ui(result->b, 15) == 0);
printf("GCD Test1 Pass: ");fmpzE_print(result);
fmpzE_clear(alpha);
fmpzE_clear(beta);
fmpzE_clear(result);
}
void TestBitcoin()
{
char *primitiveRootString10 = "70926458516305583538736597864263618897131839054605814437836363826034175889389";
char *primeNumberHexadecimal = "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F";
char *targetString = "3088993657925229173047110405354521151032325819440498983565";
fmpz_t z, prime, primitiveRoot,target, testLog, testResult;
fmpz_init(z);fmpz_init(prime);fmpz_init(primitiveRoot);fmpz_init(target);fmpz_init(testLog);fmpz_init(testResult);
fmpz_set_str(prime, primeNumberHexadecimal, 16);
fmpz_set_str(primitiveRoot, primitiveRootString10, 10);
System system = CreateSystem(primitiveRoot, prime);
PrintSystem(system);
DestroySystem(system);
fmpz_clear(z);fmpz_clear(prime);fmpz_clear(primitiveRoot);fmpz_clear(target);fmpz_clear(testLog);fmpz_clear(testResult);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment