linked big integer arithmetic to algebra routines

Author: Alan Szepieniec (laptop 'Akaros')

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2016 aszepieniec
+Copyright (c) 2017 Alan Szepieniec
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
@@ -19,3 +19,4 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.
+

README.md

@@ -9,8 +9,12 @@ Ansi C library for basic algebra routines for prime finite fields, chiefly targe
 
 # TODO
  * polynomial and extension field arithmetic
- * link big integer routines with finite field arithmetic
+ * make primality test faster
 
 # How To
- * compile: `gcc -o test_algebra test_algebra.c bi.c csprng.c Keccak-readable-and-compact-c89.c -ansi -Wpedantic`
+ * compile:
+   * small field unit tests: `gcc -o test_algebra test_algebra.c algebra.c gfp.c csprng.c Keccak-readable-and-compact-c89.c -ansi -Wpedantic`
+   * big field unit tests: `gcc -o test_algebra test_algebra.c algebra.c bi.c gfbi.c Keccak-readable-and-compact-c89.c csprng.c -DBIG -ansi -Wpedantic`
+   * big number arithmetic: `gcc -o test_bi test_bi.c bi.c csprng.c Keccak-readable-and-compact-c89.c -ansi -Wpedantic`
+   * prime sampling benchmark: `gcc -o benchmark_primality benchmark_primality.c bi.c csprng.c Keccak-readable-and-compact-c89.c -ansi -Wpedantic -O3`
 

algebra.c

@@ -48,10 +48,18 @@ int xgcd( int a, int b, int* x, int* y, int* gcd )
  */
 gfpmatrix gfpm( unsigned  int height, unsigned  int width, gfp_element* pdata )
 {
+    int i, j;
     gfpmatrix mat;
     mat.height = height;
     mat.width = width;
     mat.data = pdata;
+    for( i = 0 ; i < height ; ++i )
+    {
+        for( j = 0 ; j < width ; ++j )
+        {
+            mat.data[i*width + j] = gfp_init(0);
+        }
+    }
     return mat;
 }
 
@@ -89,9 +97,14 @@ gfpmatrix gfpm_init( unsigned  int height, unsigned  int width )
  */
 int gfpm_destroy( gfpmatrix fm )
 {
-    /*
-    printf("destroying gfpm object with data member set to memory address %#010x\n", fm.data);
-    */
+    int i, j;
+    for( i = 0 ; i < fm.height ; ++i )
+    {
+        for( j = 0 ; j < fm.width ; ++j )
+        {
+            gfp_destroy(fm.data[i*fm.width + j]);
+        }
+    }
     free(fm.data);
     fm.width = 0;
     fm.height = 0;
@@ -817,7 +830,6 @@ int gfpm_redech( gfpmatrix mat )
         {
             continue;
         }
-
         gfp_inverse(&inv, mat.data[row*mat.width + col]);
 
         if( gfp_is_one(inv) != 1 )
@@ -1226,15 +1238,27 @@ int gfpm_inverse( gfpmatrix inv, gfpmatrix mat )
 int gfpm_print( gfpmatrix mat )
 {
     unsigned int i, j;
+    printf("[");
     for( i = 0 ; i < mat.height ; ++i )
     {
+        printf("[");
         for( j = 0 ; j < mat.width ; ++j )
         {
             gfp_print(mat.data[i*mat.width+j]);
-            printf(" ");
+            if( j < mat.width - 1 )
+            {
+                printf(",");
+                printf(" ");
+            }
+        }
+        printf("]");
+        if( i < mat.height - 1 )
+        {
+            printf(",");
+            printf("\n");
         }
-        printf("\n");
     }
+    printf("]\n");
 }
 
 /**
@@ -1280,9 +1304,7 @@ hqsystem hqs_init( unsigned  int n, unsigned  int m )
     hqs.quadratic_forms = malloc(sizeof(gfpmatrix) * m);
     for( i = 0 ; i < m ; ++i )
     {
-        hqs.quadratic_forms[i].data = malloc(n*n*sizeof(gfp_element));
-        hqs.quadratic_forms[i].width = n;
-        hqs.quadratic_forms[i].height = n;
+        hqs.quadratic_forms[i] = gfpm_init(n,n);
     }
     return hqs;
 }
@@ -1373,7 +1395,7 @@ int hqs_random( hqsystem sys, unsigned char * randomness )
             for( j = 0 ; j < sys.n ; ++j )
             {
                 gfp_random(&sys.quadratic_forms[k].data[i*sys.n + j], &randomness[l]);
-                l = l + GFP_NUMBYTES + 1;
+                l = l + (GFP_NUMBITS + sizeof(unsigned long int) * 8 - 1) / (sizeof(unsigned long int) * 8);
             }
         }
     }

benchmark_primality.c

@@ -53,6 +53,7 @@ int main( int argc, char ** argv )
             csprng_generate(&rng, sizeof(unsigned long int) * num_limbs, (unsigned char*)randomness);
             bi_random(&p, bitsize, (unsigned char *)randomness);
             bi_setbit(&p, 0, 1);
+            printf("testing candidate number "); bi_print(p); printf(" for primality ...\n");
         }
         while( bi_is_prime(p, random_ints, certainty) == 0 );
     }

bi.c

@@ -564,6 +564,13 @@ int bi_subtract_ignoresign( bi * res, bi lhs, bi rhs )
 int bi_subtract( bi * res, bi lhs, bi rhs )
 {
     int cmp;
+
+    if( bi_is_zero(rhs) == 1 )
+    {
+        bi_copy(res, lhs);
+        return 1;
+    }
+
     cmp = bi_compare_absolute(lhs, rhs);
 
     if( lhs.sign == 1 && rhs.sign == 1 )
@@ -726,7 +733,6 @@ int bi_multiply( bi * res, bi lhs, bi rhs )
     unsigned int i;
     bi product;
 
-    //printf("multiplying lhs %i/%i with rhs %i/%i\n", bi_bitsize(lhs), lhs.num_limbs, bi_bitsize(rhs), rhs.num_limbs);
     if( rhs.num_limbs > lhs.num_limbs )
     {
         return bi_multiply(res, rhs, lhs);
@@ -1275,12 +1281,12 @@ int bi_random( bi * dest, unsigned int num_bits, unsigned char * randomness )
         }
     }
 
-    rem = num_bits % (sizeof(unsigned long int) * 8);
+    rem = num_bits % 8;
     if( rem == 0 )
     {
-        rem = rem + sizeof(unsigned long int) * 8;
+        rem = rem + 8;
     }
-    dest->data[num_limbs-1] = dest->data[num_limbs-1] >> (sizeof(unsigned long int) * 8 - rem);
+    dest->data[num_limbs-1] = dest->data[num_limbs-1] >> (8 - rem);
 
     return 1;
 }
@@ -1339,14 +1345,14 @@ int bi_xgcd( bi * x, bi * y, bi * g, bi a, bi b )
 
     while( !bi_is_zero(r) )
     {
-        //printf("r : %i / %i; s: %i / %i; t: %i / %i\n", bi_bitsize(r), r.num_limbs, bi_bitsize(s), s.num_limbs, bi_bitsize(t), t.num_limbs);
         bi_divide(&quotient, &remainder, old_r, r);
 
         /*
         printf("numerator: "); bi_print_bitstring(old_r); printf("\n");
         printf("divisor: "); bi_print_bitstring(r); printf("\n");
         printf("quotient: "); bi_print_bitstring(quotient); printf("\n");
         printf("remainder: "); bi_print_bitstring(remainder); printf("\n");
+        getchar();
         */
 
         /**/
@@ -1814,7 +1820,7 @@ int bi_getbit( bi integer, int bit_index )
 /**
  * bi_print
  * Send the decimal representation of this integer to stdout.
- * TODO: debig
+ * TODO: debug
  */
 int bi_print( bi integer )
 {
@@ -1824,25 +1830,38 @@ int bi_print( bi integer )
     int numdigits;
     int bitsize;
 
+    if( bi_is_zero(integer) == 1 )
+    {
+        printf("0");
+        return 1;
+    }
+
+    if( integer.sign == -1 )
+    {
+        printf("-");
+        rem1 = bi_init(0);
+        bi_copy(&rem1, integer);
+        bi_negate(&rem1);
+        bi_print(rem1);
+        bi_destroy(rem1);
+        return 1;
+    }
+
     bitsize = bi_bitsize(integer);
     expansion = malloc(bi_bitsize(integer));
     numdigits = 0;
     quo = bi_init(0);
     rem1 = bi_init(0);
     rem2 = bi_init(0);
     ten = bi_cast(10);
-    bi_copy(&rem2, integer);
-
-    if( integer.sign == -1 )
-    {
-        printf("-");
-    }
 
+    bi_copy(&rem2, integer);
 
-    for( i = 0 ; i < bitsize && bi_is_zero(rem2) == 0 ; ++i )
+    bi_one(&quo);
+    for( i = 0 ; i < bitsize && bi_is_zero(quo) == 0 ; ++i )
     {
         bi_divide(&quo, &rem1, rem2, ten);
-        bi_copy(&rem2, rem1);
+        bi_copy(&rem2, quo);
         expansion[numdigits++] = rem1.data[0];
     }
 

gfbi.c

@@ -115,7 +115,7 @@ int gfp_inverse( gfp_element * res, gfp_element elm )
     x = bi_init(0);
     y = bi_init(0);
     g = bi_init(0);
-    bi_xgcd(&g, &x, &y, elm, prime_modulus);
+    bi_xgcd(&x, &y, &g, elm, prime_modulus);
     bi_modulo(res, x, prime_modulus);
     bi_destroy(x);
     bi_destroy(y);

gfp.c

@@ -78,6 +78,24 @@ int gfp_compare( gfp_element lhs, gfp_element rhs )
     return 0;
 }
 
+int gfp_is_one( gfp_element elm )
+{
+    if( elm == 1 )
+    {
+        return 1;
+    }
+    return 0;
+}
+
+int gfp_is_zero( gfp_element elm )
+{
+    if( elm == 0 )
+    {
+        return 1;
+    }
+    return 0;
+}
+
 int gfp_add( gfp_element * res, gfp_element lhs, gfp_element rhs )
 {
     *res = ((int)(lhs) + (int)(rhs)) % GF_PRIME_MODULUS;