int p=3;  int n=3;  int d=5; int k=2;
  ring rp = p,(x(1..n)),dp;
  int s = size(maxideal(d));
  s;
→ 21
  // create a dense homogeneous ideal m, all generators of degree d, with
  // generic (random) coefficients:
  ideal m = maxideal(d)*random(p,s,n-2);
  m;
→ m[1]=x(1)^3*x(2)^2-x(1)*x(2)^4+x(1)^4*x(3)-x(1)^3*x(2)*x(3)+x(1)*x(2)^3*x\
   (3)+x(2)^4*x(3)+x(2)^3*x(3)^2+x(1)*x(2)*x(3)^3+x(1)*x(3)^4-x(3)^5
  // look for zeros on the torus by checking all points (with no component 0)
  // of the affine n-space over the field with p elements :
  ideal mt;
  int i(1..n);                    // initialize integers i(1),...,i(n)
  int l;
  s=0;
  for (i(1)=1;i(1)<p;i(1)=i(1)+1)
  {
    for (i(2)=1;i(2)<p;i(2)=i(2)+1)
    {
      for (i(3)=1;i(3)<p;i(3)=i(3)+1)
      {
        mt=m;
        for (l=1;l<=n;l=l+1)
        {
          mt=subst(mt,x(l),i(l));
        }
        if (size(mt)==0)
        {
          "solution:",i(1..n);
          s=s+1;
        }
      }
    }
  }
→ solution: 1 1 2
→ solution: 1 2 1
→ solution: 1 2 2
→ solution: 2 1 1
→ solution: 2 1 2
→ solution: 2 2 1
  "//",s,"solutions over GF("+string(p)+")";
→ // 6 solutions over GF(3)
  // Now go to the field with p^3 elements:
  // As long as there is no map from Z/p to the field with p^3 elements
  // implemented, use the following trick: convert the ideal to be mapped
  // to the new ring to a string and then execute this string in the
  // new ring
  string ms="ideal m="+string(m)+";";
  ms;
→ ideal m=x(1)^3*x(2)^2-x(1)*x(2)^4+x(1)^4*x(3)-x(1)^3*x(2)*x(3)+x(1)*x(2)^\
   3*x(3)+x(2)^4*x(3)+x(2)^3*x(3)^2+x(1)*x(2)*x(3)^3+x(1)*x(3)^4-x(3)^5;
  // define a ring rpk with p^k elements, call the primitive element z. Hence
  // 'solution exponent: 0 1 5' means that (z^0,z^1,z^5) is a solution
  ring rpk=(p^k,z),(x(1..n)),dp;
  rpk;
→ //   # ground field : 9
→ //   primitive element : z
→ //   minpoly        : 1*z^2+1*z^1+2*z^0
→ //   number of vars : 3
→ //        block   1 : ordering dp
→ //                  : names    x(1) x(2) x(3) 
→ //        block   2 : ordering C
  execute(ms);
  s=0;
  ideal mt;
  for (i(1)=0;i(1)<p^k-1;i(1)=i(1)+1)
  {
    for (i(2)=0;i(2)<p^k-1;i(2)=i(2)+1)
    {
      for (i(3)=0;i(3)<p^k-1;i(3)=i(3)+1)
      {
        mt=m;
        for (l=1;l<=n;l=l+1)
        {
          mt=subst(mt,x(l),z^i(l));
        }
        if (size(mt)==0)
        {
          "solution exponent:",i(1..n);
          s=s+1;
        }
      }
    }
  }
→ solution exponent: 0 0 2
→ solution exponent: 0 0 4
→ solution exponent: 0 0 6
→ solution exponent: 0 1 0
→ solution exponent: 0 3 0
→ solution exponent: 0 4 0
→ solution exponent: 0 4 4
→ solution exponent: 0 4 5
→ solution exponent: 0 4 7
→ solution exponent: 1 1 3
→ solution exponent: 1 1 5
→ solution exponent: 1 1 7
→ solution exponent: 1 2 1
→ solution exponent: 1 4 1
→ solution exponent: 1 5 0
→ solution exponent: 1 5 1
→ solution exponent: 1 5 5
→ solution exponent: 1 5 6
→ solution exponent: 2 2 0
→ solution exponent: 2 2 4
→ solution exponent: 2 2 6
→ solution exponent: 2 3 2
→ solution exponent: 2 5 2
→ solution exponent: 2 6 1
→ solution exponent: 2 6 2
→ solution exponent: 2 6 6
→ solution exponent: 2 6 7
→ solution exponent: 3 3 1
→ solution exponent: 3 3 5
→ solution exponent: 3 3 7
→ solution exponent: 3 4 3
→ solution exponent: 3 6 3
→ solution exponent: 3 7 0
→ solution exponent: 3 7 2
→ solution exponent: 3 7 3
→ solution exponent: 3 7 7
→ solution exponent: 4 0 0
→ solution exponent: 4 0 1
→ solution exponent: 4 0 3
→ solution exponent: 4 0 4
→ solution exponent: 4 4 0
→ solution exponent: 4 4 2
→ solution exponent: 4 4 6
→ solution exponent: 4 5 4
→ solution exponent: 4 7 4
→ solution exponent: 5 0 5
→ solution exponent: 5 1 1
→ solution exponent: 5 1 2
→ solution exponent: 5 1 4
→ solution exponent: 5 1 5
→ solution exponent: 5 5 1
→ solution exponent: 5 5 3
→ solution exponent: 5 5 7
→ solution exponent: 5 6 5
→ solution exponent: 6 1 6
→ solution exponent: 6 2 2
→ solution exponent: 6 2 3
→ solution exponent: 6 2 5
→ solution exponent: 6 2 6
→ solution exponent: 6 6 0
→ solution exponent: 6 6 2
→ solution exponent: 6 6 4
→ solution exponent: 6 7 6
→ solution exponent: 7 0 7
→ solution exponent: 7 2 7
→ solution exponent: 7 3 3
→ solution exponent: 7 3 4
→ solution exponent: 7 3 6
→ solution exponent: 7 3 7
→ solution exponent: 7 7 1
→ solution exponent: 7 7 3
→ solution exponent: 7 7 5
  "//",s,"solutions over GF("+string(p^k)+")";
→ // 72 solutions over GF(9)