Add better computation of carry chain

We port the computation of the carry chain from generate_parameters.py
to Coq, for unsaturated solinas.

Note that while we now bounds-check p448, we do not yet support
goldilocks nor karatsuba.

However, there is still an issue with the synthesized p448 code, which
is that on 64-bit, it tries to use 256-bit and 512-bit integers.  I'm
not sure what's up with that.

Partial progress towards #507

After    | File Name                                    | Before   || Change    | % Change
------------------------------------------------------------------------------------------
8m51.64s | Total                                        | 8m34.60s || +0m17.04s | +3.31%
------------------------------------------------------------------------------------------
0m15.16s | p448_solinas_64.c                            |   N/A    || +0m15.16s | ∞
3m09.12s | p384_32.c                                    | 3m09.36s || -0m00.24s | -0.12%
0m44.99s | ExtractionHaskell/word_by_word_montgomery    | 0m44.91s || +0m00.08s | +0.17%
0m39.58s | p521_32.c                                    | 0m39.22s || +0m00.35s | +0.91%
0m32.54s | p521_64.c                                    | 0m32.49s || +0m00.04s | +0.15%
0m30.87s | ExtractionHaskell/unsaturated_solinas        | 0m31.04s || -0m00.16s | -0.54%
0m24.31s | ExtractionHaskell/saturated_solinas          | 0m24.32s || -0m00.01s | -0.04%
0m18.62s | PushButtonSynthesis/UnsaturatedSolinas.vo    | 0m17.90s || +0m00.72s | +4.02%
0m17.53s | ExtractionOCaml/word_by_word_montgomery      | 0m17.44s || +0m00.08s | +0.51%
0m13.36s | secp256k1_32.c                               | 0m13.58s || -0m00.22s | -1.62%
0m13.21s | p256_32.c                                    | 0m13.15s || +0m00.06s | +0.45%
0m11.47s | p484_64.c                                    | 0m11.39s || +0m00.08s | +0.70%
0m11.27s | ExtractionOCaml/unsaturated_solinas          | 0m10.71s || +0m00.55s | +5.22%
0m10.48s | ExtractionOCaml/word_by_word_montgomery.ml   | 0m10.34s || +0m00.14s | +1.35%
0m07.97s | ExtractionOCaml/saturated_solinas            | 0m07.98s || -0m00.01s | -0.12%
0m07.05s | ExtractionOCaml/unsaturated_solinas.ml       | 0m06.98s || +0m00.06s | +1.00%
0m06.58s | ExtractionHaskell/word_by_word_montgomery.hs | 0m06.48s || +0m00.09s | +1.54%
0m06.09s | p224_32.c                                    | 0m06.04s || +0m00.04s | +0.82%
0m05.24s | p384_64.c                                    | 0m05.34s || -0m00.09s | -1.87%
0m05.13s | ExtractionOCaml/saturated_solinas.ml         | 0m05.19s || -0m00.06s | -1.15%
0m05.00s | ExtractionHaskell/unsaturated_solinas.hs     | 0m04.98s || +0m00.01s | +0.40%
0m04.14s | ExtractionHaskell/saturated_solinas.hs       | 0m04.04s || +0m00.09s | +2.47%
0m02.22s | curve25519_32.c                              | 0m02.22s || +0m00.00s | +0.00%
0m01.49s | curve25519_64.c                              | 0m01.53s || -0m00.04s | -2.61%
0m01.46s | CLI.vo                                       | 0m01.44s || +0m00.02s | +1.38%
0m01.29s | SlowPrimeSynthesisExamples.vo                | 0m01.24s || +0m00.05s | +4.03%
0m01.08s | p256_64.c                                    | 0m01.00s || +0m00.08s | +8.00%
0m01.06s | StandaloneOCamlMain.vo                       | 0m00.96s || +0m00.10s | +10.41%
0m01.06s | secp256k1_64.c                               | 0m01.17s || -0m00.10s | -9.40%
0m01.02s | p224_64.c                                    | 0m01.08s || -0m00.06s | -5.55%
0m00.99s | StandaloneHaskellMain.vo                     | 0m01.08s || -0m00.09s | -8.33%
0m00.27s | TAPSort.vo                                   |   N/A    || +0m00.27s | ∞

2 files changed, 48 insertions, 1 deletions

diff --git a/src/PushButtonSynthesis/UnsaturatedSolinas.v b/src/PushButtonSynthesis/UnsaturatedSolinas.v
index 4f22070c5..cd70771d2 100644
--- a/src/PushButtonSynthesis/UnsaturatedSolinas.v
+++ b/src/PushButtonSynthesis/UnsaturatedSolinas.v
@@ -6,6 +6,7 @@ Require Import Coq.MSets.MSetPositive.
 Require Import Coq.Lists.List.
 Require Import Coq.QArith.QArith_base Coq.QArith.Qround.
 Require Import Coq.derive.Derive.
+Require Crypto.TAPSort.
 Require Import Crypto.Util.ErrorT.
 Require Import Crypto.Util.ListUtil.
 Require Import Crypto.Util.ListUtil.FoldBool.
@@ -72,7 +73,35 @@ Section __.
           (machine_wordsize : Z).
 
   Let limbwidth := (Z.log2_up (s - Associational.eval c) / Z.of_nat n)%Q.
-  Let idxs := (List.seq 0 n ++ [0; 1])%list%nat.
+  (** Translating from https://github.com/mit-plv/fiat-crypto/blob/c60b1d2556a72c37f4bc7444204e9ddc0791ce4f/src/Specific/solinas64_2e448m2e224m1_8limbs/CurveParameters.v#L11-L35
+<<
+if len(p) > 2:
+    # do interleaved carry chains, starting at where the taps are
+    starts = [(int(t[1] / (num_bits(p) / sz)) - 1) % sz for t in p[1:]]
+    chain2 = []
+    for n in range(1,sz):
+        for j in starts:
+            chain2.append((j + n) % sz)
+    chain2 = remove_duplicates(chain2)
+    chain3 = list(map(lambda x:(x+1)%sz,starts))
+    carry_chains = [starts,chain2,chain3]
+else:
+    carry_chains = "default"
+>> *)
+  (* p is [(value, weight)]; c is [(weight, value)] *)
+  Let p := [(s / 2^Z.log2 s, Z.log2 s)] ++ TAPSort.sort (List.map (fun '(w, v) => (-v, Z.log2 w)) c).
+  Let idxs : list nat
+    := if (2 <? List.length p)%nat
+       then (* do interleaved carry chains, starting at where the taps are *)
+         let starts := List.map (fun '((v, w) : Z * Z) => (Qfloor (w / limbwidth) - 1) mod n) (tl p) in
+         let chain2 := flat_map
+                         (fun n' : nat
+                          => List.map (fun j => (j + n') mod n) starts)
+                         (List.seq 1 (pred n)) in
+         let chain2 := remove_duplicates Z.eqb chain2 in
+         let chain3 := List.map (fun x => (x + 1) mod n) starts in
+         List.map Z.to_nat (starts ++ chain2 ++ chain3)
+       else (List.seq 0 n ++ [0; 1])%list%nat.
   Let coef := 2.
   Let n_bytes := bytes_n (Qnum limbwidth) (Qden limbwidth) n.
   Let prime_upperbound_list : list Z
diff --git a/src/TAPSort.v b/src/TAPSort.v
new file mode 100644
index 000000000..a9882c882
--- /dev/null
+++ b/src/TAPSort.v
@@ -0,0 +1,18 @@
+Require Import Coq.ZArith.ZArith.
+Require Import Coq.Sorting.Mergesort.
+
+
+Module TAPOrder <: Orders.TotalLeBool.
+  Local Coercion is_true : bool >-> Sortclass.
+  Definition t := (Z * Z)%type.
+  Definition leb (x y : t) := Z.geb (snd x) (snd y).
+  Infix "<=?" := leb.
+  Theorem leb_total : forall a1 a2, a1 <=? a2 \/ a2 <=? a1.
+  Proof.
+    intros; cbv [is_true leb]; rewrite !Z.geb_le; omega.
+  Qed.
+End TAPOrder.
+
+Module TAPSort := Sort TAPOrder.
+
+Notation sort := TAPSort.sort.