Normalize the size of sqrt remainder

examples/multiprecision/Makefile

@@ -47,7 +47,7 @@ tests: extract check-gmp
 	./build/minitests
 
 bench-tests: extract
-	gcc $(CFLAGS) tests.c $(CFILES) -Ibench-include -Irandom -lgmp -o build/bench-tests
+	gcc $(CFLAGS) tests.c $(CFILES) -Iinclude -Ibench-include -Irandom -lgmp -o build/bench-tests
 
 
 build/why3%bench: extract check-gmp

examples/multiprecision/sqrtrem.mlw

@@ -859,7 +859,9 @@ module Sqrt
     raises   { StackOverflow -> true }
     ensures  { value np n = value sp (ceilhalf n) * value sp (ceilhalf n)
                             + value rp result }
+    ensures  { 0 <= result <= n }
     ensures  { value rp result <= 2 * value sp (ceilhalf n) }
+    ensures  { result > 0 -> (pelts rp)[offset rp + result - 1] > 0 }
     ensures  { max np = old max np }
     ensures  { min np = old min np }
     ensures  { plength np = old plength np }
@@ -1223,6 +1225,7 @@ module Sqrt
         end
       end;
       begin ensures { value rp rn = vn - vsp * vsp }
+            ensures { 0 < rn <= k+1 }
         if (not (c2 = 0))
         then begin
           let ghost b = wmpn_rshift rp tp tn c2 in
@@ -1268,6 +1271,24 @@ module Sqrt
                by [@case_split] (h = 0 \/ h = 1) };
       rn <- tn + to_int32 h;
     end;
+    let ghost orp = pure { rp } in
+    let ghost orn = pure { rn } in
+    assert { value np n = value sp k * value sp k + value orp orn };
+    assert { 1 <= rn <= n };
+    while C.get_ofs rp (rn - 1) = 0 do
+      invariant { value rp rn = value orp orn }
+      invariant { 1 <= rn <= orn }
+      variant   { rn }
+      value_tail rp (rn-1);
+      assert { value rp (rn - 1) = value rp rn };
+      rn <- rn - 1;
+      if rn = 0
+      then begin
+        assert { value orp orn = 0
+                 by 0 = value rp 0 = value rp 1 = value orp orn };
+        break
+      end
+    done;
     rn
 
 end
\ No newline at end of file

examples/multiprecision/tests.c

@@ -13,6 +13,7 @@
 #define TEST_TOOM
 #define TEST_DIV
 #define TEST_SQRT1
+#define TEST_SQRTREM
 #endif
 
 #ifdef TEST_MINIGMP
@@ -27,6 +28,7 @@
 #include "build/div.h"
 #include "build/toom.h"
 #include "build/sqrt1.h"
+#include "build/sqrt.h"
 #endif
 
 #include "mt19937-64.c"
@@ -56,7 +58,7 @@ void init_valid (mp_ptr ap, mp_ptr bp, mp_size_t an, mp_size_t bn) {
 
 int main () {
   mp_ptr ap, bp, rp, refp, rq, rr, refq, refr;
-  mp_size_t max_n, max_add, max_mul, max_toom, max_div, an, bn, rn;
+  mp_size_t max_n, max_add, max_mul, max_toom, max_div, max_sqrt, an, bn, rn;
 #ifdef BENCH
   int nb, nb_iter;
   struct timeval begin, end;
@@ -84,6 +86,7 @@ int main () {
   max_mul = 20;
   max_toom = 95;
   max_div = 20;
+  max_sqrt = 95;
   ap = TMP_ALLOC_LIMBS (max_n + 1);
   bp = TMP_ALLOC_LIMBS (max_n + 1);
   /* nap = TMP_ALLOC_LIMBS (max_n + 1); */
@@ -371,7 +374,7 @@ int main () {
 #endif
       an = bn = rn = 1;
       for (int iter = 0; iter != 500; ++iter) {
-        init_valid (ap, bp, 1, 1);
+        init_valid (bp, ap, 1, 1);
         a = *ap;
         if (a < 0x4000000000000000) continue;
 #ifdef BENCH
@@ -393,7 +396,7 @@ int main () {
         elapsed +=
           (end.tv_sec - begin.tv_sec)
           + ((end.tv_usec - begin.tv_usec)/1000000.0)
-        printf ("%f\n", an, bn, elapsed);
+        printf ("%f\n", elapsed);
         printf ("\n"); //for gnuplot
 #endif
 #ifdef COMPARE
@@ -415,7 +418,87 @@ int main () {
 #endif
 #endif
 
-          //TMP_FREE;
-          //tests_end ();
-          return 0;
+#ifdef TEST_SQRTREM
+#ifdef BENCH
+  printf ("#an t(µs)\n");
+#endif
+  bn=1;
+  for (an = 1; an <= max_sqrt; an += 1)
+    {
+      init_valid (bp, ap, 1, an);
+      
+#ifdef BENCH
+      elapsed = 0;
+      nb_iter = 1000;
+      for (int iter = 0; iter != nb_iter; ++iter) {
+        init_valid (ap, bp, an, 1);
+#ifdef TEST_MINIGMP
+        mpn_copyi(refr, ap, an);
+#endif
+        gettimeofday(&begin, NULL);
+        nb = 1500 / an;
+        for (int i = 0; i != nb; ++i)
+          {
+#endif
+#if defined(TEST_GMP) || defined(TEST_MINIGMP)
+            rn = mpn_sqrtrem(refq, refr, ap, an);
+#endif
+#ifdef TEST_WHY3
+            c = wmpn_sqrtrem(rq, rr, ap, an);
+#endif
+            
+#ifdef BENCH
+          }
+        gettimeofday(&end, NULL);
+        elapsed +=
+          (end.tv_sec - begin.tv_sec) * 1000000.0
+          + (end.tv_usec - begin.tv_usec);
+      }
+      elapsed = elapsed / (nb * nb_iter);
+      printf ("%d %f\n", an, elapsed);
+      printf ("\n"); //for gnuplot
+#endif
+#ifdef COMPARE
+      if (c != rn)
+        {
+          printf ("ERROR, an = %d, expected rn = %d, actual rn = %d\n",
+                  (int) an, (int) rn, (int) c);
+          printf ("a: "); mpn_dump (ap, an);
+          printf ("s: "); mpn_dump (rq, (an+1)/2);
+          printf ("refs: "); mpn_dump (refq, (an+1)/2);
+          printf ("r: "); mpn_dump (rr, c);
+          printf ("refr: "); mpn_dump (refr, rn);
+          abort ();
+        }
+      if (mpn_cmp (refr, rr, rn))
+        {
+          printf ("ERROR, an = %d, rn = %d\n",
+                  (int) an, (int) rn);
+          printf ("a: "); mpn_dump (ap, an);
+          printf ("s: "); mpn_dump (rq, (an+1)/2);
+          printf ("refs: "); mpn_dump (refq, (an+1)/2);
+          printf ("r: "); mpn_dump (rr, c);
+          printf ("refr: "); mpn_dump (refr, rn);
+          abort();
         }
+      if (mpn_cmp (refq, rq, (an+1)/2))
+        {
+          printf ("ERROR, an = %d, rn = %d\n",
+                  (int) an, (int) rn);
+          printf ("a: "); mpn_dump (ap, an);
+          printf ("s: "); mpn_dump (rq, (an+1)/2);
+          printf ("refs: "); mpn_dump (refq, (an+1)/2);
+          printf ("r: "); mpn_dump (rr, c);
+          printf ("refr: "); mpn_dump (refr, rn);
+          abort();
+        }
+#endif
+    }
+#ifdef COMPARE
+  printf ("sqrtrem ok\n");
+#endif
+#endif
+  //TMP_FREE;
+  //tests_end ();
+  return 0;
+}