Fast Square Root using BIGNUMs.

After doing some research, turned out that probably the best algorithm for
computing the square root of a BIGNUM is dijkstra algorithm, only using binay
shifts. The algorithm has complexity O(lg n), which is probabily best than
Newton's Algorithm (search the zero for f(x) = x² - N), which complexity is
instead O(M(n)), where M(n) is the complexity for a single multiplication.

src/questions/test/test_wiener.c

@@ -5,11 +5,11 @@
 #include <math.h>
 #include <string.h>
 
-#include <openssl/x509.h>
-#include <openssl/pem.h>
 
+#include <openssl/pem.h>
 #include <openssl/ssl.h>
 #include <openssl/x509.h>
+#include <openssl/err.h>
 
 #include "questions.h"
 #include "qwiener.h"
@@ -120,6 +120,7 @@ void test_BN_sqrtmod(void)
   BIGNUM *mayzero;
   BN_CTX *ctx;
 
+  a = b = expected = NULL;
   root = BN_new();
   rem = BN_new();
   mayzero = BN_new();
@@ -161,11 +162,11 @@ void test_BN_sqrtmod(void)
 void test_wiener(void)
 {
   X509 *crt;
-  FILE *fp = fopen("test/wiener_test.crt", "r");
+  FILE *fp = fopen("questions/test/wiener_test.crt", "r");
 
+  if (!fp) exit(EXIT_FAILURE);
   crt = PEM_read_X509(fp, NULL, 0, NULL);
   if (!crt) {
-    ERR_print_errors();
     exit(EXIT_FAILURE);
   }
 

src/questions/wiener.c

@@ -150,44 +150,43 @@ static void BN_int2bn(BIGNUM** a, short int i)
 /**
  * \brief Square Root for bignums.
  *
+ * An implementation of Dijkstra's Square Root Algorithm.
+ * A Discipline of Programming, page 61 - Fifth Exercise.
+ *
  */
 int BN_sqrtmod(BIGNUM* dv, BIGNUM* rem, BIGNUM* a, BN_CTX* ctx)
 {
-  char *abn2dec;
-  int g[100];
-  long al;
-  long x = 0, r = 0;
-  int i, j;
-  int d;
-  long y, yn;
-
-  abn2dec = BN_bn2dec(a);
-  sscanf(abn2dec, "%ld", &al);
-
-  r = 0;
-  x = 0;
-  for (i=0; al > 0; i++) {
-    g[i] = al%100;
-    al /= 100;
-  }
-
-  for (j=i-1; j>=0; j--) {
-    r = r*100 + g[j];
-    y = 0;
-    for (d=1; d!=10; d++) {
-      yn = d*(20*x + d);
-      if (yn <= r) y = yn; else break;
+  BIGNUM *shift;
+  BIGNUM *adj;
+
+  shift = BN_new();
+  adj = BN_new();
+  BN_zero(dv);
+  BN_copy(rem, a);
+
+  /* hacking into internal sequence to skip some cycles. */
+  /* for  (BN_one(shift);     original */
+  for (bn_wexpand(shift, a->top+1), shift->top=a->top, shift->d[shift->top-1] = 1;
+       BN_ucmp(shift, rem) != 1;
+       /* BN_rshift(shift, shift, 2); */
+       BN_lshift1(shift, shift), BN_lshift1(shift, shift));
+
+
+  while (!BN_is_one(shift)) {
+    /* BN_rshift(shift, shift, 2); */
+    BN_rshift1(shift, shift);
+    BN_rshift1(shift, shift);
+
+    BN_uadd(adj, dv, shift);
+    BN_rshift1(dv, dv);
+    if (BN_ucmp(rem, adj) != -1) {
+      BN_uadd(dv, dv, shift);
+      BN_usub(rem, rem, adj);
     }
-    r -= y;
-    x = 10*x + d -1;
   }
 
-  sprintf(abn2dec, "%ld", r);
-  BN_dec2bn(&rem, abn2dec);
-  sprintf(abn2dec, "%ld", x);
-  BN_dec2bn(&dv, abn2dec);
-  OPENSSL_free(abn2dec);
-
+  BN_free(shift);
+  BN_free(adj);
   return BN_is_zero(rem);
 }