function rel(p, a, m, r, b)  /* the key relations */
    return (1 + a^(p-m) - b^(p-m))*
           (1 - b^(p+m-r))*(1 - (a-1)^(p+m-r));
end function;

function reltab(p, r, R, n) /* matrix of key relations */
  assert Type(p) eq RngIntElt and IsPrime(p) and IsOdd(p);
  assert Type(r) eq RngIntElt and IsEven(r);

  nc := (p-3) div 2; 
  nr := (p-1) div 2;

  M  := RMatrixSpace(R,nr,nc)!0;
  for i in [1..nc] do
     for j in [1..nr] do
        a := (R!(2*(i+1)))^(p^n);
        b := (R!2)^p;
        M[j, i] := rel(p, a, 2*j-1, r, b);
     end for;
  end for;
  return M;
end function;

function antientry(p, r, x, y) /* antisymmetry relations */
  if ((x-y) mod (p-1) eq 0) or ((x+y-r) mod (p-1) eq 0) then
     return 1;
  else
     return 0;
  end if;
end function;

function antitab(p, r, R) /* matrix of antisymmetry relations */
  assert Type(p) eq RngIntElt and IsPrime(p) and IsOdd(p);
  assert Type(r) eq RngIntElt and IsEven(r);

  nc := 1 + Floor((p-1)/4); 
  nr := (p-1) div 2;
  c  := (r div 2) - ((r div 2) mod 2);

  M  := RMatrixSpace(R,nr,nc)!0;
  for i in [1..nc] do
     for j in [1..nr] do
        M[j, i] := antientry(p, r, 2*i+c-1, 2*j-1);
     end for;
  end for;
  return M;
end function;

function kerp(p,r: n:=1) /* this calculates the nullspace for (p,r) */
  assert Type(p) eq RngIntElt and IsPrime(p) and IsOdd(p);
  assert Type(r) eq RngIntElt and IsEven(r);

  nranti := 1 + Floor((p-1)/4); 
  nrrel := (p-3) div 2;
  nc := (p-1) div 2;
  c  := (r div 2) - ((r div 2) mod 2);
  
  if (n eq 1) then R := GF(p);
  else R := IntegerRing(p^n);
  end if;
  V := RSpace(R,nc);
  W := sub<V|V>;
  W := NullSpace(Hom(W,RSpace(R,nrrel))!reltab(p,r,R,n-1));
  W := NullSpace(Hom(W,RSpace(R,nranti))!antitab(p,r,R));
  return W;
end function;

function bern(p,r)  /* this calculates B_r mod p (for p not really small) */
  assert Type(p) eq RngIntElt and IsPrime(p);
  assert Type(r) eq RngIntElt and IsEven(r);
 
  k := r div 2;
  m := GF(p)!k;
  fac1 := m/(2*m+4);
  fac2 := GF(p)!1;
  term := fac1;
  
  for j in [2..k] do
    i := GF(p)!j;
    fac1 *:= -(i-1)*(m-i+1)/((i+1)*(m+i+1));
    fac2 +:= i^r;
    term +:= fac1*fac2;
  end for;
  return term;  
end function;

function powsum(p,a,b,k)
  sum := GF(p)!0;
  for i in [a+1..b] do
    sum +:= (GF(p)!i)^k;
  end for;
  return sum;
end function;

function irreg(p,r) /* this will test if p divides B_r:
                       first, it uses classical congruences for
                       Bernoulli numbers to try to determine it quickly.
		       if these cannot be used, it calculates B_r mod p */
  assert Type(p) eq RngIntElt and IsPrime(p);
  assert Type(r) eq RngIntElt and IsEven(r);

  j := (GF(p)!2)^(r-1);
  k := (GF(p)!3)^(r-1);
  if not ((k-1)*j^2 eq k-j) then
    a := (p-1) div 6;
    b := (p-1) div 4;
    return IsZero(powsum(p,a,b,r-1));
  else
    if not ((j eq -1) or (2*j eq 1)) then
      b := (p-1) div 6;
      return IsZero(powsum(p,0,b,r-1));
    else
      if not (3*k eq 1) then
        b := (p-1) div 3;
        return IsZero(powsum(p,0,b,r-1));
        else 
          if not (2*j eq 1) then
             b := (p-1) div 2;
             return IsZero(powsum(p,0,b,r-1));
          else return IsZero(bern(p,r));
        end if;      
      end if;
    end if;
  end if; 
end function; 

function isreg(p)  /* this will test if a prime is regular */
  S := { r : r in [2..p-3] | IsEven(r) };
  T := { r : r in S | irreg(p,r) };

  return IsEmpty(T);
end function;



procedure allpair(p1,p2) /* this computes the nullspace for all pairs 
                            (p,r) with p1 <= p <= p2 */
 
  p := NextPrime(p1-1);
  while p le p2  do
    k2 := (p-3) div 2;
    for k in [1..k2] do
      r := 2*k;
      V := kerp(p,r);  
      printf "p = %o, r = %o, dim = %o,\n%o\n", p, r,
        Dimension(V), Basis(V);        
    end for;
    p := NextPrime(p);
  end while;
end procedure;

procedure pair(p1,p2) /* this computes the nullspace for irregular pairs */

  p := NextPrime(p1-1);
  while p le p2  do
    k2 := (p-3) div 2;
    for k in [1..k2] do
      r := 2*k;
      if irreg(p,r) then
        V := kerp(p,r);  
        printf "p = %o, r = %o, dim = %o,\n%o\n", p, r,
          Dimension(V), Basis(V);   
      end if;
    end for;
    p := NextPrime(p);
  end while;
end procedure;

procedure pairneq2(p1,p2) /* this checks if the nullspace mod p^2 
                             has order p, assuming it does mod p */
  p := NextPrime(p1-1);
  while p le p2  do
    k2 := (p-3) div 2;
    for k in [1..k2] do
      r := 2*k;
      if irreg(p,r) then
        test := Basis(kerp(p,r: n := 2))[1]*p;  
        printf "p = %o, r = %o, order p = %o\n", p, r, (test eq 0);   
      end if;
    end for;
    p := NextPrime(p);
  end while;
end procedure;