From d3bd187f75ac47a9a51d23da01e07b7e29a3a007 Mon Sep 17 00:00:00 2001
From: Julien Coolen <jtcoolen@pm.me>
Date: Tue, 29 Oct 2024 09:24:39 +0100
Subject: [PATCH] update snf

---
 examples/number_fields.ml |  7 +++++++
 src/pari.ml               | 21 ++++++++++++++++++++-
 src/pari.mli              | 30 +++++++++++++++++++++++++++---
 3 files changed, 54 insertions(+), 4 deletions(-)

diff --git a/examples/number_fields.ml b/examples/number_fields.ml
index 4542cdb..4171af3 100644
--- a/examples/number_fields.ml
+++ b/examples/number_fields.ml
@@ -51,6 +51,13 @@ let () =
   Printf.eprintf "%b\n"
     Number_field.(equal a (add gaussian_integers (mul gaussian_integers b q) r))
 
+(* Smith normal form *)
+let m = gp_read_str "[1, 0, x; 0, x, 0; 0,0,1+x]"
+let d, u, v = Number_field.smith_normal_form gaussian_integers m
+
+let _ =
+  Printf.eprintf "res=%s %s\n" (gentostr Matrix.(mul (mul u m) v)) (gentostr d)
+
 let _nf2 =
   Number_field.create (Option.get @@ Polynomial.of_string "x^4-2" |> inj_rat)
 
diff --git a/src/pari.ml b/src/pari.ml
index 41eee3f..8ed5348 100644
--- a/src/pari.ml
+++ b/src/pari.ml
@@ -17,6 +17,7 @@ type integer_mod = Integer_mod
 type finite_field = Finite_field
 type number_field = Number_field
 type 'a elliptic_curve = Elliptic_curve of 'a
+type 'a matrix = private Matrix of 'a
 
 let register_gc v =
   Gc.finalise_last (fun () -> pari_free Ctypes.(coerce gen (ptr void) v)) v
@@ -557,7 +558,23 @@ module Number_field = struct
   let nfmulmodpr = nfmulmodpr
   let nfpowmodpr = nfpowmodpr
   let nfreduce = nfreduce
-  let nfsnf = nfsnf
+
+  let smith_normal_form nf sipm =
+    (*square integral invertible pseudo-matrix sipm*)
+    let d, _ = Matrix.dimensions sipm in
+    let sipm =
+      Vector.of_array
+        [|
+          sipm;
+          Vector.init d ~f:(fun _ -> Integer.of_int 1);
+          Vector.init d ~f:(fun _ -> Integer.of_int 1);
+        |]
+    in
+    let res = nfsnf0 nf sipm Signed.Long.one in
+    ( Vector.(res.%[1]),
+      Vector.(res.%[2]) |> matbasistoalg nf,
+      Vector.(res.%[3]) |> matbasistoalg nf )
+
   let nfsnf0 = nfsnf0
   let nfsolvemodpr = nfsolvemodpr
   let nfsubfields = nfsubfields
@@ -914,6 +931,8 @@ module Elliptic_curve = struct
   let ell_is_inf = ell_is_inf
 end
 
+let isdiagonal e = isdiagonal e = 1
+
 let factor e =
   let m = factor e in
   let n_rows, n_cols = Matrix.dimensions m in
diff --git a/src/pari.mli b/src/pari.mli
index f82871a..bfad899 100644
--- a/src/pari.mli
+++ b/src/pari.mli
@@ -241,8 +241,10 @@ module Vector : sig
   end
 end
 
+type 'a matrix = private Matrix of 'a
+
 module Matrix : sig
-  type 'a t constraint 'a = 'b ty
+  type 'a t = 'a matrix ty constraint 'a = 'b ty
 
   val dimensions : 'a t -> int * int
   val id : int -> Integer.t t
@@ -742,7 +744,29 @@ module Number_field : sig
   val nfmulmodpr : 'a ty -> 'a ty -> 'a ty -> 'a ty -> 'a ty
   val nfpowmodpr : 'a ty -> 'a ty -> 'a ty -> 'a ty -> 'a ty
   val nfreduce : 'a ty -> 'a ty -> 'a ty -> 'a ty
-  val nfsnf : 'a ty -> 'a ty -> 'a ty
+
+  val smith_normal_form :
+    t -> elt Matrix.t -> elt Matrix.t * elt Matrix.t * elt Matrix.t
+  (**
+      {@ocaml[
+      # let gaussian_integers =
+        (* Gaussian integers: the ring Z[i] (here we work in the field Q(i)) *)
+        Number_field.create
+        (Polynomial.create [| Rational.of_int 1; Rational.of_int 0; Rational.of_int 1 |]);;
+      val gaussian_integers : Number_field.t = <abstr>
+      # let m = gp_read_str "[1, 0, x; 0, x, 0; 0,0,1+x]"
+      val m : '_weak1 ty = <abstr>
+      # let _, u, v = Number_field.smith_normal_form gaussian_integers m
+      val u : Number_field.elt Matrix.t = <abstr>
+      val v : Number_field.elt Matrix.t = <abstr>
+      # let d = Matrix.(mul (mul u m) v);;
+      val d : Number_field.elt Matrix.t = <abstr>
+      # gentostr d
+      - : string = "[1, 0, 0; 0, 1, 0; 0, 0, 1]"
+      # isdiagonal d 
+      - : bool = true
+      ]} *)
+
   val nfsnf0 : 'a ty -> 'a ty -> Signed.long -> 'a ty
   val nfsolvemodpr : 'a ty -> 'a ty -> 'a ty -> 'a ty -> 'a ty
   val nfsubfields : 'a ty -> Signed.long -> 'a ty
@@ -3639,7 +3663,7 @@ val rgv_zc_mul : 'a ty -> 'a ty -> 'a ty
 val rgv_zm_mul : 'a ty -> 'a ty -> 'a ty
 val rgx_rgm_eval : 'a ty -> 'a ty -> 'a ty
 val rgx_rgmv_eval : 'a ty -> 'a ty -> 'a ty
-val isdiagonal : 'a ty -> int
+val isdiagonal : 'a ty -> bool
 val scalarcol : 'a ty -> Signed.long -> 'a ty
 val scalarcol_shallow : 'a ty -> Signed.long -> 'a ty
 val scalarmat : 'a ty -> Signed.long -> 'a ty