congruences.c in congruence_solver-0.5.0

- old
+ new

@@ -1,35 +1,36 @@
 #include "arith_utils.h"
 #include "prime_gen.h"
 #include <stdlib.h>
 #include <stdio.h>
 
-static int * adjust_coeffs_to_mod(int degree, int * coeffs, int mod);
-static int * solve_prime_power_congruence(int degree, int coeffs[], int prime, int power);
-static int * solve_system_of_order_1_congruence_sets(int numOfSets, int * lengthsOfSets, int ** sets, int mods[]);
+static long * adjust_coeffs_to_mod(int degree, long * coeffs, long mod);
+static long * solve_prime_power_congruence(int degree, long coeffs[], long prime, int power);
+static long * solve_system_of_order_1_congruence_sets(int numOfSets, int * lengthsOfSets, long ** sets, long mods[]);
 
-int chinese_remainder_solution(int numberOfEquations, int scals[], int mods[]){
+long chinese_remainder_solution(int numberOfEquations, long scals[], long mods[]){
   int i;
-  int x = 0;
-  int m = mods[0];
-  int modCoeff;
+  long long x = 0;
+  long m = mods[0];
+  long modCoeff;
 
   for(i=1; i<numberOfEquations; i++){
     m *= mods[i];
   }
 
   for(i=0; i<numberOfEquations; i++){
     modCoeff = m/mods[i];
-    x += modCoeff*mod_inv(modCoeff % mods[i], mods[i])*scals[i];
+    x += (long long) modCoeff*mod_inv(modCoeff % mods[i], mods[i])*scals[i] % m;
+    x %= m;
   }
 
   return x % m;
 }
 
 
-int * adjust_coeffs_to_mod(int degree, int * coeffs, int mod){
-  int * adjustedCoeffs = calloc(degree+1, sizeof(int));
+long * adjust_coeffs_to_mod(int degree, long * coeffs, long mod){
+  long * adjustedCoeffs = calloc(degree+1, sizeof(long));
   int i;
 
   for(i = 0; i <= degree; i++){
     adjustedCoeffs[i] = coeffs[i] % mod;
     if(adjustedCoeffs[i] < 0){
@@ -39,23 +40,23 @@
 
   return adjustedCoeffs;
 }
 
 
-int * brute_force_congruence(int degree, int coeffs[], int primeMod){
+long * brute_force_congruence(int degree, long coeffs[], long primeMod){
   //assumes a prime modulus. split congruences of composite modulus into systems of congrueneces
   //of prime modulus and/or apply the lifting theorem to make use of this function
   //solve a0x^n + a1x^n-1... = 0 (mod mod) where n is the order a0, a1, ... are coeffieicients
   //also assumes positive representation of coeffs
-  int * adjustedCoeffs = adjust_coeffs_to_mod(degree, coeffs, primeMod);
-  int * solutionList = calloc(degree+1, sizeof(int));
-  int * solutions = solutionList+1;
+  long * adjustedCoeffs = adjust_coeffs_to_mod(degree, coeffs, primeMod);
+  long * solutionList = calloc(degree+1, sizeof(long));
+  long * solutions = solutionList+1;
   int numberOfSolutions = 0;
-  int x;
+  long x;
 
-
-  for(x = 0; x < primeMod && numberOfSolutions <= degree; x++){
+  #pragma omp parallel for
+  for(x = 0; x < primeMod; x++){
     if(mod_eval_polynomial(degree, adjustedCoeffs, primeMod, x) == 0){
       solutions[numberOfSolutions++] = x;
     }
   }
 
@@ -65,45 +66,42 @@
 
   return solutionList;
 }
 
 
-static int * solve_prime_power_congruence(int funcDegree, int funcCoeffs[], int prime, int power){
-
-  int * adjustedCoeffs;
-
-  int * baseSolutionList;
+static long * solve_prime_power_congruence(int funcDegree, long funcCoeffs[], long prime, int power){
+  long * baseSolutionList;
   int numOfBaseSolutions;
-  int * baseSolutions;
+  long * baseSolutions;
 
-  int * liftedSolutions;
+  long * liftedSolutions;
   int numOfLiftedSolutions;
 
-  int coeff;
+  long coeff;
 
   int derivDegree;
-  int * derivCoeffs;
-  int deriv;
-  long int divFunc;
+  long * derivCoeffs;
+  long deriv;
+  long func;
 
   int i, j, t;
-  int currentMod;
+  long currentMod;
 
   if(power == 1){
     baseSolutions = brute_force_congruence(funcDegree, funcCoeffs, prime);
     return baseSolutions;
   }
 
   baseSolutionList = solve_prime_power_congruence(funcDegree, funcCoeffs, prime, power-1);
   numOfBaseSolutions = *baseSolutionList;
   baseSolutions = baseSolutionList+1;
 
-  liftedSolutions = calloc(prime*numOfBaseSolutions+1, sizeof(int));
+  liftedSolutions = calloc(prime*numOfBaseSolutions+1, sizeof(long));
   numOfLiftedSolutions = 0;
 
   derivDegree = funcDegree-1;
-  derivCoeffs = calloc(derivDegree+1, sizeof(int));
+  derivCoeffs = calloc(derivDegree+1, sizeof(long));
 
   currentMod = prime;
   for(j = 1; j < power; j++){
     currentMod *= prime;
   }
@@ -118,18 +116,18 @@
 
 
   for(j = 0; j < numOfBaseSolutions; j++){
 
     deriv = mod_eval_polynomial(derivDegree, derivCoeffs, prime, baseSolutions[j]);
-    divFunc = (eval_polynomial(funcDegree, funcCoeffs, baseSolutions[j]) / (currentMod/prime)) % prime;
+    func = mod_eval_polynomial(funcDegree, funcCoeffs, currentMod, baseSolutions[j]);
 
     if(deriv % prime != 0){
-      t = (-divFunc*mod_inv(deriv, prime) % prime) + prime;
+      t = ((-func / (currentMod/prime))*mod_inv(deriv, currentMod) % prime) + prime;
       liftedSolutions[++numOfLiftedSolutions] = baseSolutions[j] + t*prime;
     }
 
-    else if(divFunc % prime == 0){
+    else if(func == 0){
       for(t = 1; t <= prime; t++){
         liftedSolutions[++numOfLiftedSolutions] = baseSolutions[j] + t*(currentMod/prime);
       }
     }
   }
@@ -142,26 +140,26 @@
 
   return liftedSolutions;
 }
 
 
-static int * solve_system_of_order_1_congruence_sets(int numOfSets, int * setLengths, int * * sets, int * mods){
+static long * solve_system_of_order_1_congruence_sets(int numOfSets, int * setLengths, long * * sets, long * mods){
   //allocate perumtation array
-  int * divAry = calloc(numOfSets, sizeof(int));
-  int * scalAry = calloc(numOfSets, sizeof(int));
+  long * divAry = calloc(numOfSets, sizeof(long));
+  long * scalAry = calloc(numOfSets, sizeof(long));
   int i, j;
   int numOfSolutions;
-  int * solutionAry;
-  int * dest;
+  long * solutionAry;
+  long * dest;
   int idx;
 
   for(i = 0, numOfSolutions = 1; i < numOfSets; i++){
     divAry[i] = numOfSolutions;
     numOfSolutions *= setLengths[i];
   }
 
-  solutionAry = calloc(numOfSolutions+1, sizeof(int));
+  solutionAry = calloc(numOfSolutions+1, sizeof(long));
   solutionAry[0] = numOfSolutions;
   dest = solutionAry+1;
 
   for(i = 0; i < numOfSolutions; i++){
     for(j = 0; j < numOfSets; j++){
@@ -173,28 +171,28 @@
   }
 
   return solutionAry;
 }
 
-int * solve_congruence(int funcDegree, int funcCoeffs[], int mod){
-  int * solutionList;
+long * solve_congruence(int funcDegree, long funcCoeffs[], long mod){
+  long * solutionList;
 
-  int * modFactorList = prime_factors(mod);
-  int numOfModFactors = *modFactorList;
-  int * modFactors = modFactorList+1;
+  long * modFactorList = prime_factors(mod);
+  long numOfModFactors = *modFactorList;
+  long * modFactors = modFactorList+1;
 
-  int * * primePowerSolutions = calloc(numOfModFactors, sizeof(int *));
-  int * primePowers = calloc(numOfModFactors, sizeof(int));
+  long * * primePowerSolutions = calloc(numOfModFactors, sizeof(long *));
+  long * primePowers = calloc(numOfModFactors, sizeof(long));
   int * primePowerSolutionLengths = calloc(numOfModFactors, sizeof(int *));
 
   int power;
-  int i;
+  int i,j,k;
 
   for(i = 0; i < numOfModFactors; i++){
-    primePowers[i] = modFactors[i]; 
+    primePowers[i] = modFactors[i];
     power = 1;
-    
+
     while(mod % (primePowers[i]*modFactors[i]) == 0){
       primePowers[i] *= modFactors[i];
       power++;
     }
 
@@ -212,16 +210,5 @@
   free(primePowers);
   free(modFactorList);
 
   return solutionList;
 }
-
-/*
-int * solve_system_of_congruences(int numOfFuncs, int * funcDegrees, int ** funcCoeffs, int * mods){
-  int i;
-  int * * funcSolutionSets = calloc(numOfFuncs, sizeof(int *));
-  for(i=0; i<numOfFuncs; i++){
-    funcSolutionSets[i] = solve_congruence(funcDegrees[i], funcCoeffs[i], mods[i]);
-  }
-  return solve_system_of_congruence_sets(numOfFuncs, funcSolutionSets, mods);
-} 
-*/
\ No newline at end of file