magma check works

magma_utils/magma_utils.mgm

 // In this file, we provide some magma functions to test that the output of rrspace is correct.
 
-RRSPACE_PATH="/users/pspaenle/repos/rrspace"
+RRSPACE_PATH        :="/users/pspaenle/repos/rrspace";
+RRSPACE_TMP_INFILE  :="/tmp/tmpin.rrspace";
+RRSPACE_TMP_OUTFILE :="/tmp/tmpout.rrspace";
 
 // Given a function f and a finite set B of functions on a curve, check if the k-vector space generated by B
 // contains f
@@ -17,11 +19,7 @@ function CheckBelongsSpan(f, B)
     lM := [];
     lv := [];
     for i in [1..nbtests] do
-	repeat
-	    P := Random(C);
-	until P ne Infinity() and
-	      &and[Evaluate(g, P) ne Infinity() : g in B] and
-	      Evaluate(f, P) ne Infinity();
+	P := RandomPlace(C, 1);
 	lM cat:= [Evaluate(g, P) : g in B];
 	lv cat:= [Evaluate(f, P)];
     end for;
@@ -36,22 +34,154 @@ function CheckEqualSpan(B1, B2)
 end function;
 
 // Set variable order y > x > z
+// Ambient space of the curve should be the projective space
 function Check_rrspace(C, D)
     // first compute the input representation for rrspace
     K := BaseRing(C);
     assert Type(K) eq FldFin;
-    assert AmbientSpace(C) eq ProjectiveSpace(K, 2);
     Dp, Dm := SignDecomposition(D);
     Ip := Ideal(Dp);
     Im := Ideal(Dm);
-    R := Universe(Generators(Ip));
-    Gp := GroebnerBasis(ChangeOrder(Ip+Ideal(R.3-1), "lex"));
-    Gm := GroebnerBasis(ChangeOrder(Im+Ideal(R.3-1), "lex"));
+    R := PolynomialRing(K, 3);
+    eqnC  := R!Equation(C);
+    eqnC /:= LeadingCoefficient(eqnC);
+    assert LeadingMonomial(eqnC) eq R.1^Degree(C);
+    Is := Ideal([Derivative(eqnC, 1), Derivative(eqnC, 2), Derivative(eqnC, 3), R.3-1]); // There should not be any node at infinity
+    Gp := GroebnerBasis(Ideal([R!g : g in Generators(Ip)] cat [R.3-1]));
+    Gm := GroebnerBasis(Ideal([R!g : g in Generators(Im)] cat [R.3-1]));
+    Gs := GroebnerBasis(Is);
     assert(Degree(Gp[2]) eq Degree(Ip));
-    assert(Degree(Gm[2]) eq Degree(Im));
+    assert(Gm eq [1] or Degree(Gm[2]) eq Degree(Im));
+    assert(Gs eq [1] or Degree(Gs[2]) eq Dimension(R/Is));
+    if Gm eq [1] then
+	Gm := [R!0, R!1];
+    end if;
+    if Gs eq [1] then
+	Gs := [R!0, R!1];
+    end if;
     Runiv := PolynomialRing(K);
     fp := Evaluate( Gp[2], [0, Runiv.1, 0]);
     gp := Evaluate(-Gp[1], [0, Runiv.1, 0]);
     fm := Evaluate( Gm[2], [0, Runiv.1, 0]);
     gm := Evaluate(-Gm[1], [0, Runiv.1, 0]);
+    fs := Evaluate( Gs[2], [0, Runiv.1, 0]);
+    gs := Evaluate(-Gs[1], [0, Runiv.1, 0]);
+
+    // Write input file for rrspace
+    SetOutputFile(RRSPACE_TMP_INFILE : Overwrite := true);
+    printf "%o\n", #K;
+    printf "[ %o ", Degree(C);
+    tmp_eqn := eqnC;
+    for i in [0..Degree(C)] do
+	tmp_c := Evaluate(tmp_eqn, [0, Runiv.1, 1])/Factorial(i);
+	tmp_eqn := Derivative(tmp_eqn, 1);
+	printf "[";
+	lc := Eltseq(tmp_c);
+	for j in [1..#lc] do
+	    printf "%o", lc[j];
+	    if j ne #lc then
+		printf " ";
+	    end if;
+	end for;
+	printf "] ";
+    end for;
+    printf " ]\n";
+    printf "< [";
+
+    lc := Eltseq(fs);
+    for j in [1..#lc] do
+	printf "%o", lc[j];
+	if j ne #lc then
+	    printf " ";
+	end if;
+    end for;
+    printf "] ";
+    printf "[";
+    lc := Eltseq(gs);
+    for j in [1..#lc] do
+	printf "%o", lc[j];
+	if j ne #lc then
+	    printf " ";
+	end if;
+    end for;
+    printf "] ";
+    printf ">\n";
+    
+    printf "{ < [";
+    lc := Eltseq(fp);
+    for j in [1..#lc] do
+	printf "%o", lc[j];
+	if j ne #lc then
+	    printf " ";
+	end if;
+    end for;
+    printf "] ";
+    printf "[";
+    lc := Eltseq(gp);
+    for j in [1..#lc] do
+	printf "%o", lc[j];
+	if j ne #lc then
+	    printf " ";
+	else
+	    printf "] ";
+	end if;
+    end for;
+    printf "> < [";
+    lc := Eltseq(fm);
+    for j in [1..#lc] do
+	printf "%o", lc[j];
+	if j ne #lc then
+	    printf " ";
+	end if;
+    end for;
+    printf "] ";
+    printf "[";
+    lc := Eltseq(gm);
+    for j in [1..#lc] do
+	printf "%o", lc[j];
+	if j ne #lc then
+	    printf " ";
+	end if;
+    end for;
+    printf "] ";
+    printf "> }\n";
     
+    UnsetOutputFile();
+
+    // Compute with rrspace
+
+    System(RRSPACE_PATH cat
+	   "/rrspace -c < " cat
+	   RRSPACE_TMP_INFILE cat
+	   " > " cat
+	   RRSPACE_TMP_OUTFILE);
+
+    FF<y,x> := FunctionField(C);
+    basis_rrspace := eval Read(RRSPACE_TMP_OUTFILE);
+
+    vs_rr, phi := RiemannRochSpace(D);
+
+    print Dimension(vs_rr), #basis_rrspace;
+//    print basis_rrspace;
+    
+    return Dimension(vs_rr) eq #basis_rrspace and
+	   CheckEqualSpan([phi(b) : b in Basis(vs_rr)], basis_rrspace);
+end function;
+
+K := GF(1009);
+R<Y,X,Z> := PolynomialRing(K, 3);
+
+Q := &+[Random(K) * m : m in MonomialsOfDegree(R, 4)];
+C := Curve(ProjectiveSpace(K, 2), Q);
+
+D := Divisor(RandomPlace(C, 5));
+Check_rrspace(C, D);
+
+K := GF(1009);
+R<Y,X,Z> := PolynomialRing(K, 3);
+
+Q := - Y^2*Z^2 + X^2*Z^2 + Y^4 -X^3*Z;
+C := Curve(ProjectiveSpace(K, 2), Q);
+
+D := Divisor(RandomPlace(C, 5));
+Check_rrspace(C, D);

rrspace.cc

@@ -21,11 +21,12 @@ Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 #include <iostream>
 #include <sstream>
 #include "algos.h"
+#include <unistd.h>
 
 using namespace NTL;
 using namespace std;
 
-int main() {
+int main(int argc, char** argv) {
   long p;
   BivPol eq_curve;
   cin >> p;
@@ -41,13 +42,40 @@ int main() {
   cin >> D;
 
   RRspace basisRR = RiemannRochBasis(D, E);
+
+  bool MAGMACHECK_FLAG = false;
+
+  int c;
+  while ((c = getopt(argc, argv, "chvo:s:")) != -1)
+    switch (c) {
+      case 'c':
+        MAGMACHECK_FLAG = true;
+        break;
+      default:
+        std::cerr << "option not recognized, aborting." << std::endl;
+        abort();
+    }
   
-  cout << "Dimension: " << basisRR.num_basis.size() << endl;
-  cout << "Denominator: " << endl;
-  PrintMagma(std::cout, basisRR.denom) << endl;
-  cout << "Numerators: " << endl;
-  for (size_t i = 0; i < basisRR.num_basis.size(); ++i)
-    PrintMagma(std::cout, basisRR.num_basis[i]) << endl;
-
-  return 0;
+  if (!MAGMACHECK_FLAG) {
+    cout << "Dimension: " << basisRR.num_basis.size() << endl;
+    cout << "Denominator: " << endl;
+    PrintMagma(std::cout, basisRR.denom) << endl;
+    cout << "Numerators: " << endl;
+    for (size_t i = 0; i < basisRR.num_basis.size(); ++i)
+      PrintMagma(std::cout, basisRR.num_basis[i]) << endl;
+  } else {
+    cout << "[" << endl;
+    for (size_t i = 0; i < basisRR.num_basis.size(); ++i) {
+      cout << "(";
+      PrintMagma(std::cout, basisRR.num_basis[i]);
+      cout << ")/(";
+      PrintMagma(std::cout, basisRR.denom);
+      cout << ")";
+      if (i != basisRR.num_basis.size()-1)
+        cout << ",";
+      cout << endl;
+    }
+    cout << "]" << endl;
+  }
+  return EXIT_SUCCESS;
 }