Folds/iterators for primitive arrays

[palmskog]

Although Coq has primitive arrays, it currently provides no folds/iterators for such arrays.

One possible application of folds over primitive arrays is in refinement of abstract matrices to arrays of arrays (https://github.com/coq-community/coqeal/issues/72).

@maximedenes @strub

[zhengpushi]

It's an interesting problem. I have studied several formal matrix models ago, and focused on the properties of the matrix operations such as "mapping, fold, addition, subtraction, scalar multiplication, transposition, multiplication". Meanwhile, all of them could be used to compute with concrete types such as Z, Q, Qc, R.

The effecient in Coq (and extracted code in OCaml) is very different in these models. For example, matrix multiplication with the element of R type (extracted to float type in OCaml), and matrix size from 20x20 to 211x211, the running time (second) in OCaml is approximately below:

size	DL	DP	DR	NF	FF
20	0.00	0.00	0.00	0.00	3.50
50	0.26	0.02	0.01	0.37	?
60	0.52	0.04	0.02	0.76	?
72	1.04	0.06	0.02	1.58	?
86	2.06	0.10	0.04	3.22	?
103	4.18	0.18	0.07	6.52	?
123	8.42	0.30	0.12	13.11	?
147	17.26	0.51	0.19	26.63	?
176	34.46	0.86	0.33	54.43	?
211	72.90	1.48	0.55	111.90	?

DL : dependent list, (Coq.Vectors.Vector.t) DP : dependent pair, (a fixpoint function) DR : dependent record, (inner is a list of list) NF : function with natrual indexing (function style) FF : function with fin indexing (from MathComp)

---

From above, the DR model is the fastest. Multiply two 211x211 float matrices, need 0.55 seconds. But, this is still too slow than primitive array.

Conclusion: I think the fasted way to compute matrix operation in Coq directly is primitive array, and worthwhile to study the properties of matrix operation in this model.

The OCaml source file at https://github.com/zhengpushi/coq-matrix/blob/main/MatrixComparison/src/Matrix/OCamlExtractionComparison/test.ml (Note that, you may need to compile these project to get the extracted ocaml file from Coq. and the compile is easy with make.) The detail could be checked at https://github.com/zhengpushi/coq-matrix/blob/main/MatrixComparison/doc/ocaml-extraction-comparison.md The proejct could be checked at https://github.com/zhengpushi/coq-matrix/ The comparison of these models could be found in https://link.springer.com/chapter/10.1007/978-3-031-21213-0_11 or the draft here: https://github.com/zhengpushi/coq-matrix/blob/main/MatrixComparison/doc/matrix-coq-setta2022.pdf

Note that, this project is written with Chinese mixed with English, not very friendly for English users. I will correct it gradually

[SkySkimmer]

What kind of API would this be? Is it enough to have interruptible finite loops over ints ie

Parameter for_loop : forall {A B}, int -> int -> (A -> B) -> (int -> A -> A + B) -> A -> B

which in ocaml would be

let rec loop start end_ finish f x =
  if start = end_ then finish x
  else match f start x with
    | Inl x -> loop (start+1) end_ finish f x
    | Inr x -> x

let for_loop start end_ finish f x = if start <= end then loop start end_ finish f x else finish x

? Then eg

Definition array_fold_left {A B} (f : A -> B -> A) (acc : A) (a : array B) : A
  := for_loop 0 (Array.length a) id (fun i acc => inl (f acc (Array.get a i))) acc.

and

Definition array_exists {A} (p : A -> bool) (a : array A) : bool
  := for_loop 0 (Array.length a) (fun _ => false) (fun i _ => if p (Array.get a i) then inr true else inl tt) tt.

??

[zhengpushi]

This demo program maybe work well in OCaml. @SkySkimmer . Looks very well.

I have writen some code of PArray.array, which can be evaluated in Coq directly. But it is not so efficient, which use nat type to write a fixpoint function, and with the help of "n2i" to translate "int" to "nat" type.

(*
  remark    :
  1. A loop need to trace the index of int type, but int is axiomized, 
     havn't a inductive structure, hard to write fixpoint function.
     (1). we can use nat, the structure is easy, but max-size-of-nat is 5000
     (2). we can use Z, it is a bit complicate, and need to solve negative.
     (3). we can use positive, but also complicate.
 *)

Require Import Coq.Numbers.Cyclic.Int63.Uint63. (* native int *)
Require Import PArray.          (* native array *)

Module parray_loop_with_nat.

  (* int <-> nat *)
  Definition i2n : int -> nat := fun i => BinInt.Z.to_nat (to_Z i).
  Definition n2i : nat -> int := fun n : nat => of_Z (BinInt.Z.of_nat n).
  Compute i2n 2. (* = 2 : nat *)
  Compute n2i 2. (* = 2%uint63 : int *)

  (** loop : init + a[i] + a[i+1] + ... *)
  Fixpoint loop {A:Type} (arr : array A) (cnt : nat) (i : int) (f : A -> A -> A)
    (init : A) : A :=
    match cnt with
    | O => init
    | S cnt' => loop arr cnt' (i + 1) f (f init arr.[i])
    end.

  Definition data1 := make 50 1.
  Compute data1.
  (* = [| 
     1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 
     1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 
     1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 
     1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 
     1; 1; 1; 1; 1; 1; 1; 1; 1; 1 | 1 : nat |]
     : array nat *)

  Variable f : nat -> nat -> nat.
  Compute loop data1 10 0 f 0. (*
     = f (f (f (f (f (f (f (f (f (f 0 1) 1) 1) 1) 1) 1) 1) 1) 1) 1
     : nat *)
  Compute loop data1 10 0 (Nat.add) 0. (*  = 10 
                                           : nat *)
End parray_loop_with_nat.