#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <stdbool.h>
#include <inttypes.h>
#include <math.h>
// Define constants for the size of lookup tables
#define TABLE_SIZE_103 103
#define TABLE_SIZE_139 139
#define TABLE_SIZE_163 163
#define TABLE_SIZE_19 19
#define TABLE_SIZE_241 241
#define TABLE_SIZE_37 37
#define TABLE_SIZE_61 61
#define TABLE_SIZE_67 67
#define TABLE_SIZE_7 7
#define TABLE_SIZE_79 79
#define TABLE_SIZE_97 97
#define TABLE_SIZE_9 9
#define TABLE_SIZE_13 13
#define TABLE_SIZE_337 337
#define TABLE_SIZE_379 379
#define TABLE_SIZE_199 199
#define TABLE_SIZE_271 271
#define TABLE_SIZE_421 421
#define TABLE_SIZE_409 409
// Define the size of the matrix
#define MATRIX_SIZE 3
// Structure to represent a Matrix
typedef struct {
uint64_t data[MATRIX_SIZE][MATRIX_SIZE];
} Matrix;
// Structure to represent a Vector
typedef struct {
uint64_t data[MATRIX_SIZE];
} Vector;
// Lookup tables (Initialize these based on your specific logic)
int tab103[TABLE_SIZE_103]={1,1,0,1,0,0,0,0,1,1,1,0,0,1,1,0,0,0,0,0,0,0,1,1,1,0,0,1,0,0,1,1,0,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0,1,0,0,1,1,0,0,1,0,0,1,1,1,0,0,0,0,0,0,0,1,1,0,0,1,1,1,0,0,0,0,1,0,1};
int tab139[TABLE_SIZE_139]={1,1,0,0,0,0,1,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,0,1,1,0,0,1,0,0,0,1,0,0,1,0,1,0,1,1,0,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,0,1,1,0,1,0,1,0,0,1,0,0,0,1,0,0,1,1,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,1,0,0,0,0,1};
int tab163[TABLE_SIZE_163]={1,1,0,0,0,1,1,0,1,0,0,0,0,1,0,0,0,1,0,0,0,1,1,1,0,1,0,1,1,0,1,1,0,0,0,0,1,1,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,1,0,1,1,1,0,0,0,0,1,1,0,1,1,0,1,0,1,1,1,0,0,0,1,0,0,0,1,0,0,0,0,1,0,1,1,0,0,0,1};
int tab67[TABLE_SIZE_67]={1,1,0,1,0,1,0,0,1,1,0,0,0,0,1,1,0,0,0,0,0,0,1,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,1,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,0,1,0,1,0,1};
int tab241[TABLE_SIZE_241]={1,1,0,0,0,1,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,1,1,1,1,0,1,0,0,1,0,0,1,0,0,0,1,1,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,0,0,1,1,1,0,1,1,0,0,0,0,1,0,0,0,1,1,1,0,0,0,0,0,0,1,1,1,0,0,0,1,0,0,0,0,1,1,0,1,1,1,0,0,1,0,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,0,0,0,1,0,0,1,0,0,1,0,1,1,1,1,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,1};
int tab37[TABLE_SIZE_37]={1,1,0,0,0,0,1,0,1,0,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0,1,0,1,0,0,0,0,1};
int tab61[TABLE_SIZE_61]={1,1,0,1,0,0,0,0,1,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,1,1,0,0,0,0,1,1,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,0,1,0,1};
int tab79[TABLE_SIZE_79]={1,1,0,0,0,0,0,0,1,0,1,0,1,0,1,1,0,1,1,0,0,1,1,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,1,1,0,0,1,1,0,1,1,0,1,0,1,0,1,0,0,0,0,0,0,1};
int tab97[TABLE_SIZE_97]={1,1,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,1,1,1,0,1,0,0,0,0,1,1,0,1,0,0,1,1,0,0,0,0,0,0,0,1,0,0,1,1,1,0,0,1,1,1,0,0,1,0,0,0,0,0,0,0,1,1,0,0,1,0,1,1,0,0,0,0,1,0,1,1,1,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,1};
// Add other tables as needed
int tab337[TABLE_SIZE_337]={1,1,0,0,0,1,1,1,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,1,1,0,0,1,1,0,1,1,0,0,0,1,1,1,0,0,1,0,0,1,1,1,1,1,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,1,0,1,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,1,1,1,0,1,0,1,0,0,0,0,0,0,0,1,1,1,0,0,0,0,1,0,0,0,1,0,1,0,1,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0,0,1,0,1,0,1,0,0,0,1,0,0,0,0,1,1,1,0,0,0,0,0,0,0,1,0,1,0,1,1,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,1,0,0,1,0,1,0,1,0,0,1,1,1,1,1,0,0,1,0,0,1,1,1,0,0,0,1,1,0,1,1,0,0,1,1,0,0,0,0,1,0,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1,1,1,1,0,0,0,1};
int tab379[TABLE_SIZE_379]={1,1,0,0,0,1,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,1,0,0,1,0,0,1,1,0,1,1,1,0,0,1,0,0,0,1,0,0,1,1,0,0,0,0,1,0,1,0,0,0,1,1,0,0,1,1,0,1,0,0,1,0,0,1,1,0,0,0,0,0,1,1,0,1,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,1,0,1,1,1,0,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,1,0,1,0,0,1,1,0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,1,1,1,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,0,1,1,1,0,0,0,0,1,0,0,0,0,0,0,1,0,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0,1,0,0,1,0,1,0,0,1,0,0,1,1,1,0,1,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,0,0,0,1,0,1,1,0,0,0,0,0,1,1,0,0,1,0,0,1,0,1,1,0,0,1,1,0,0,0,1,0,1,0,0,0,0,1,1,0,0,1,0,0,0,1,0,0,1,1,1,0,1,1,0,0,1,0,0,1,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,1};
int tab199[TABLE_SIZE_199]={1,1,0,0,0,1,0,0,1,0,0,1,1,0,0,0,0,1,1,0,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,1,1,1,1,1,0,0,1,0,0,0,0,0,0,1,0,1,0,1,0,0,0,1,1,0,1,0,0,1,0,1,0,1,1,0,0,1,0,1,0,0,1,0,1,0,0,1,1,0,1,0,1,0,0,1,0,1,1,0,0,0,1,0,1,0,1,0,0,0,0,0,0,1,0,0,1,1,1,1,1,1,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,0,1,0,0,1,0,0,0,1};
int tab271[TABLE_SIZE_271]={1,1,0,1,0,0,0,0,1,1,1,0,0,1,0,0,0,0,0,1,0,0,0,1,1,0,0,1,1,1,1,1,0,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,1,1,0,0,1,0,1,1,0,0,1,0,0,1,0,0,0,0,1,0,1,0,1,0,1,1,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,1,1,1,0,1,0,1,0,1,0,0,0,0,1,0,0,1,0,0,1,1,0,1,0,0,1,1,1,0,0,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,1,1,0,0,1,1,1,1,1,0,0,1,1,0,0,0,1,0,0,0,0,0,1,0,0,1,1,1,0,0,0,0,1,0,1};
int tab421[TABLE_SIZE_421]={1,1,0,0,0,0,1,1,1,0,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,1,1,0,0,0,0,1,0,1,1,0,1,1,1,0,1,0,0,0,1,1,0,0,1,1,1,1,0,1,0,0,1,1,1,1,0,0,0,0,1,0,0,1,0,1,0,0,0,0,1,0,0,0,1,0,1,1,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,1,0,1,0,1,0,1,1,1,0,1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0,1,1,0,0,0,1,0,1,0,0,0,0,1,0,0,1,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,0,0,1,1,0,1,1,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,0,1,0,1,0,0,0,1,1,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,1,1,1,0,1,0,1,0,1,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,1,1,0,1,0,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,0,1,1,1,1,0,0,1,0,1,1,1,1,0,0,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,0,1,1,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,0,1};
int tab409[TABLE_SIZE_409]={1,1,0,0,0,1,1,0,1,0,0,1,0,1,1,0,0,0,0,1,0,0,0,0,0,1,0,1,0,0,1,1,0,0,0,0,1,0,0,0,1,0,0,1,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,1,0,1,0,1,1,1,1,0,1,1,1,0,0,0,0,0,0,0,1,1,0,0,1,1,1,0,0,0,1,1,0,0,1,0,1,1,0,0,0,0,0,0,0,1,1,0,1,0,0,0,0,1,1,0,1,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,0,0,1,1,1,0,1,1,1,0,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,0,0,0,1,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,1,1,0,0,0,0,0,1,1,0,1,1,1,0,1,1,1,0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,1,0,1,1,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,1,0,0,0,1,1,1,0,0,1,1,0,0,0,0,0,0,0,1,1,1,0,1,1,1,1,0,1,0,1,1,0,0,1,0,0,0,1,0,0,1,0,0,0,0,1,0,0,1,0,0,0,1,0,0,0,0,1,1,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,0,0,1};
/********************************************************************************
// Initialize lookup tables with example values (Replace with actual values)
void initialize_lookup_tables() {
// Example initialization: set all to 0
for(int i = 0; i < TABLE_SIZE_103; i++) tab103[i] = 0;
for(int i = 0; i < TABLE_SIZE_139; i++) tab139[i] = 0;
for(int i = 0; i < TABLE_SIZE_163; i++) tab163[i] = 0;
for(int i = 0; i < TABLE_SIZE_67; i++) tab67[i] = 0;
for(int i = 0; i < TABLE_SIZE_241; i++) tab241[i] = 0;
for(int i = 0; i < TABLE_SIZE_37; i++) tab37[i] = 0;
for(int i = 0; i < TABLE_SIZE_61; i++) tab61[i] = 0;
for(int i = 0; i < TABLE_SIZE_79; i++) tab79[i] = 0;
for(int i = 0; i < TABLE_SIZE_97; i++) tab97[i] = 0;
// Populate tables with actual data as needed
}
*****************************************************************************/
// Function Prototypes for simpleXXX functions
bool simple103(uint64_t X);
bool simple139(uint64_t X);
bool simple163(uint64_t X);
bool simple19(uint64_t X);
bool simple241(uint64_t X);
bool simple37(uint64_t X);
bool simple61(uint64_t X);
bool simple67(uint64_t X);
bool simple7(uint64_t X);
bool simple79(uint64_t X);
bool simple97(uint64_t X);
bool simple9(uint64_t X);
bool simple13(uint64_t X);
// Implementation of simpleXXX functions
bool simple103(uint64_t X) {
uint64_t rem = X % 103;
if(rem == 0) return true;
if(rem > 0) return tab103[rem];
return false;
}
bool simple139(uint64_t X) {
uint64_t rem = X % 139;
if(rem == 0) return true;
if(rem > 0) return tab139[rem];
return false;
}
bool simple163(uint64_t X) {
uint64_t rem = X % 163;
if(rem == 0) return false;
if(rem > 0) return tab163[rem];
return false;
}
bool simple19(uint64_t X) {
uint64_t rem = X % 19;
if(rem == 0) return false;
if(rem > 0) return (rem == 1 || rem == 7 || rem == 8 || rem == 11 || rem == 12 || rem == 18);
return false;
}
bool simple241(uint64_t X) {
uint64_t rem = X % 241;
if(rem == 0) return true;
if(rem > 0) return tab241[rem];
return false;
}
bool simple337(uint64_t X) {
uint64_t rem = X % 337;
if(rem == 0) return true;
if(rem > 0) return tab337[rem];
return false;
}
bool simple421(uint64_t X) {
uint64_t rem = X % 421;
if(rem == 0) return true;
if(rem > 0) return tab421[rem];
return false;
}
bool simple379(uint64_t X) {
uint64_t rem = X % 379;
if(rem == 0) return true;
if(rem > 0) return tab379[rem];
return false;
}
bool simple199(uint64_t X) {
uint64_t rem = X % 199;
if(rem == 0) return true;
if(rem > 0) return tab199[rem];
return false;
}
bool simple271(uint64_t X) {
uint64_t rem = X % 271;
if(rem == 0) return true;
if(rem > 0) return tab271[rem];
return false;
}
bool simple409(uint64_t X) {
uint64_t rem = X % 409;
if(rem == 0) return true;
if(rem > 0) return tab409[rem];
return false;
}
bool simple37(uint64_t X) {
uint64_t rem = X % 37;
return (rem == 1 || rem == 6 || rem == 8 || rem == 11 || rem == 14 ||
rem == 23 || rem == 26 || rem == 29 || rem == 31 || rem == 36 ||
rem == 10 || rem == 27);
}
bool simple61(uint64_t X) {
uint64_t rem = X % 61;
return (rem == 1 || rem == 3 || rem == 8 || rem == 9 || rem == 11 ||
rem == 23 || rem == 27 || rem == 28 || rem == 33 || rem == 34 ||
rem == 38 || rem == 50 || rem == 52 || rem == 53 || rem == 58 ||
rem == 60 || rem == 24 || rem == 37 || rem == 20 || rem == 41 ||
rem == 0);
}
bool simple67(uint64_t X) {
uint64_t rem = X % 67;
if(rem == 0) return true;
if(rem > 0) return tab67[rem];
return false;
}
bool simple7(uint64_t X) {
uint64_t rem = X % 7;
return (rem == 1 || rem == 6 || rem == 0);
}
bool simple79(uint64_t X) {
uint64_t rem = X % 79;
if(rem == 0) return true;
if(rem > 0) return tab79[rem];
return false;
}
bool simple97(uint64_t X) {
uint64_t rem = X % 97;
if(rem == 0) return false;
if(rem > 0) return tab97[rem];
return false;
}
bool simple9(uint64_t X) {
uint64_t rem = X % 9;
if(rem == 0) return true;
if(rem > 0) return (rem == 1 || rem == 8);
return false;
}
bool simple13(uint64_t X) {
uint64_t rem = X % 13;
return (rem == 1 || rem == 5 || rem == 8 || rem == 12);
}
// Discriminant Functions
bool disc160757041(uint64_t X) {
return !simple409(X);
}
// disc160757041 =
// (X)->val=!simple409(X);return(val)
bool disc63984001(uint64_t X) {
return !simple421(X);
}
bool disc2839225(uint64_t X) {
return !simple337(X);
}
bool disc120802081(uint64_t X) {
return !simple379(X);
}
bool disc4791721(uint64_t X) {
return !simple199(X);
}
bool disc61763881(uint64_t X) {
return !simple271(X);
}
bool disc10220809(uint64_t X) {
// printf("This is function disc10220809\n");
return !simple139(X);
}
bool disc112225(uint64_t X) {
// printf("This is function disc112225\n");
return !simple67(X);
}
bool disc165649(uint64_t X) {
// printf("This is function disc165449\n");
return !simple37(X);
}
bool disc16605625(uint64_t X) {
// printf("This is function disc16605625\n");
return !simple163(X);
}
bool disc16785409(uint64_t X) {
// printf("This is function disc16785409\n");
return !simple241(X);
}
bool disc17689(uint64_t X) {
// printf("This is function disc17689\n");
return !simple19(X);
}
bool disc1792921(uint64_t X) {
// printf("This is function disc1792921\n");
return !simple103(X);
}
bool disc1803649(uint64_t X) {
// printf("This is function disc1803649\n");
return !simple79(X);
}
bool disc3396649(uint64_t X) {
// printf("This is function disc3396649\n");
return !simple97(X);
}
bool disc3721(uint64_t X) {
// printf("This is function disc3721\n");
return !simple61(X);
}
bool disc4225(uint64_t X) {
// printf("This is function disc4225\n");
return !simple13(X);
}
bool disc49(uint64_t X) {
// printf("This is function disc49\n");
return !simple7(X);
}
bool disc81(uint64_t X) {
// printf("This is function disc81\n");
return !simple9(X);
}
// Function to initialize matrix M and vector v3
void init3(uint64_t p, uint64_t r, uint64_t s, Matrix *M, Vector *v3) {
// Initialize matrix M
M->data[0][0] = 0;
M->data[0][1] = r % p;
M->data[0][2] = s % p;
M->data[1][0] = 1 % p;
M->data[1][1] = 0;
M->data[1][2] = 0;
M->data[2][0] = 0;
M->data[2][1] = 1 % p;
M->data[2][2] = 0;
// Initialize vector v3
v3->data[0] = (2 * r) % p;
v3->data[1] = 0;
v3->data[2] = 3 % p;
}
// Function to multiply two matrices: result = a * b mod p with Overflow Protection
void multiply_matrices(Matrix *a, Matrix *b, Matrix *result, uint64_t p) {
// printf("Multiplying two matrices modulo %" PRIu64 "\n", p);
for(int i = 0; i < MATRIX_SIZE; i++) {
for(int j = 0; j < MATRIX_SIZE; j++) {
__int128 temp = 0; // Use 128-bit integer to prevent overflow
for(int k = 0; k < MATRIX_SIZE; k++) {
// Cast to __int128 before multiplication to prevent overflow
__int128 product = (__int128)a->data[i][k] * (__int128)b->data[k][j];
temp += product;
temp %= p; // Apply modulo at each step to keep temp manageable
}
result->data[i][j] = (uint64_t)(temp % p);
// printf("Result[%d][%d] = %" PRIu64 "\n", i, j, result->data[i][j]);
}
}
}
/***************************************************************************************************
// Function to multiply two matrices: result = a * b mod p with Overflow Protection
void multiply_matrices(uint64_t a[3][3], uint64_t b[3][3], uint64_t result[3][3], uint64_t p) {
// printf("Multiplying two matrices modulo %" PRIu64 "\n", p);
for(int i = 0; i < 3; i++) {
for(int j = 0; j < 3; j++) {
__int128 temp = 0; // Use 128-bit integer to prevent overflow
for(int k = 0; k < 3; k++) {
__int128 product = (__int128)a[i][k] * (__int128)b[k][j];
temp += product;
temp %= p; // Apply modulo at each step to keep temp manageable
}
result[i][j] = (uint64_t)(temp % p);
// printf("Result[%d][%d] = %" PRIu64 "\n", i, j, result[i][j]);
}
}
}
***************************************************************************************/
/*******************************************************************
// Function to multiply two matrices: result = a * b mod p
void multiply_matrices(Matrix *a, Matrix *b, Matrix *result, uint64_t p) {
for(int i = 0; i < MATRIX_SIZE; i++) {
for(int j = 0; j < MATRIX_SIZE; j++) {
result->data[i][j] = 0;
for(int k = 0; k < MATRIX_SIZE; k++) {
result->data[i][j] += (a->data[i][k] * b->data[k][j]) % p;
result->data[i][j] %= p;
}
}
}
}
**************************************************************/
// Function to compute matrix exponentiation: result = a^n mod p
void matrix_power_func(Matrix *a, uint64_t n, Matrix *result, uint64_t p) {
// Initialize result as identity matrix
Matrix temp;
for(int i = 0; i < MATRIX_SIZE; i++) {
for(int j = 0; j < MATRIX_SIZE; j++) {
temp.data[i][j] = (i == j) ? 1 : 0;
}
}
Matrix base = *a;
while(n > 0) {
if(n % 2 == 1) {
multiply_matrices(&temp, &base, &temp, p);
}
multiply_matrices(&base, &base, &base, p);
n /= 2;
}
*result = temp;
}
/************************************************************
// Function to multiply matrix and vector: result = M * v mod p
void multiply_matrix_vector(Matrix *M, Vector *v, Vector *result, uint64_t p) {
for(int i = 0; i < MATRIX_SIZE; i++) {
result->data[i] = 0;
for(int j = 0; j < MATRIX_SIZE; j++) {
result->data[i] += (M->data[i][j] * v->data[j]) % p;
result->data[i] %= p;
}
}
}
********************************************************************/
// Function to multiply a matrix by a vector: result = M * v mod p with Overflow Protection
void multiply_matrix_vector(Matrix *M, Vector *v, Vector *result, uint64_t p) {
for(int i = 0; i < MATRIX_SIZE; i++) {
result->data[i] = 0;
for(int j = 0; j < MATRIX_SIZE; j++) {
// Cast to __int128 before multiplication to prevent overflow
__int128 product = (__int128)M->data[i][j] * (__int128)v->data[j];
// Compute product modulo p
uint64_t mod_product = (uint64_t)(product % p);
// Accumulate the result modulo p
result->data[i] += mod_product;
result->data[i] %= p;
}
}
}
// Function to compute power of a matrix: a^n
Matrix power(Matrix *a, uint64_t n, uint64_t p) {
Matrix result;
// Initialize result as identity matrix
for(int i = 0; i < MATRIX_SIZE; i++) {
for(int j = 0; j < MATRIX_SIZE; j++) {
result.data[i][j] = (i == j) ? 1 : 0;
}
}
Matrix base = *a;
while(n > 0) {
if(n % 2 == 1) {
Matrix temp;
multiply_matrices(&result, &base, &temp, p);
result = temp;
}
Matrix temp;
multiply_matrices(&base, &base, &temp, p);
base = temp;
n /= 2;
}
return result;
}
// Function to perform Frobenius test (isfrob3select)
int isfrob3select(uint64_t p) {
bool found = false;
uint64_t r, s;
// Handle special cases
if(p == 2 || p == 3) return 1;
if((p % 2) == 1 && (p % 3) > 0) {
if(disc49(p)) { r = 7; s = 7; found = true; }
if(!found && disc81(p)) { r = 3; s = 1; found = true; }
if(!found && disc165649(p)) { r = 37; s = 37; found = true; }
if(!found && disc4225(p)) { r = 13; s = 13; found = true; }
if(!found && disc3721(p)) { r = 61; s = 183; found = true; }
if(!found && disc3396649(p)) { r = 97; s = 97; found = true; }
if(!found && disc1803649(p)) { r = 79; s = 79; found = true; }
if(!found && disc17689(p)) { r = 19; s = 19; found = true; }
if(!found && disc10220809(p)) { r = 139; s = 139; found = true; }
if(!found && disc16605625(p)) { r = 163; s = 163; found = true; }
if(!found && disc112225(p)) { r = 67; s = 201; found = true; }
if(!found && disc1792921(p)) { r = 103; s = 309; found = true; }
if(!found && disc16785409(p)) { r = 241; s = 1205; found = true; }
if(!found && disc61763881(p)) { r = 271; s = 813; found = true; }
if(!found && disc4791721(p)) { r = 199; s = 995; found = true; }
if(!found && disc120802081(p)) { r = 379; s = 1895; found = true; }
if(!found && disc2839225(p)) { r = 337; s = 2359; found = true; }
if(!found && disc63984001(p)) { r = 421; s = 2947; found = true; }
if(!found && disc160757041(p)) { r = 409; s = 2045; found = true; }
// if((!found)&&disc63984001(p),r=421;s=2947;found=1);
if(!found) return -1;
// Initialize matrix M and vector v3
Matrix M;
Vector v3;
init3(p, r, s, &M, &v3);
// Compute M^p mod p
Matrix Mpp = power(&M, p, p);
// Multiply Mpp with v3
Vector output;
multiply_matrix_vector(&Mpp, &v3, &output, p);
// Perform the test
bool test = (output.data[1] == ((2*p - r) % p)) && (output.data[2] == 0);
return test ? 1 : 0;
}
return 0;
}
// Function to check if a number is prime (simple trial division, can be optimized)
bool is_prime(uint64_t n) {
if(n < 2) return false;
if(n == 2 || n == 3) return true;
if(n % 2 == 0 || n % 3 == 0) return false;
for(uint64_t i = 5; i * i <= n; i += 6) {
if(n % i == 0 || n % (i + 2) == 0) return false;
}
return true;
}
int main() {
// Define the range for X
uint64_t start = 15000000000;
uint64_t end = 100000000000;
uint64_t pcount = 0;
// Iterate through the range
for(uint64_t X = start; X <= end; X++) {
if((X % 2) == 1) { // Check if X is odd
int c1 = isfrob3select(X);
// int c2 = is_prime(X) ? 1 : 0; // Convert boolean to integer (1 for prime, 0 for not prime)
if(c1 == 1) {
pcount++;
}
// if((c1==-1)&&(is_prime(X)))
// {
// printf("%" PRIu64 " passses isprime but not isfrob3select\n", X);
// }
if(99999999==(X%100000000))
{
printf("X= %" PRIu64 " , pcount=%" PRIu64 "\n", X+1, pcount);
}
}
}
printf("%" PRIu64 " probable primes found\n", pcount);
return 0;
}
Related
