ppoly := proc(pe, vb, n, x, nlist, max)
    local i, j, dcycle, clen, ob10, nb10, res, cyct, vs,
          vc, plist, k, dg, nsum, tmp, m, ct, tmax;
    vs := [];
    vc := [];
    plist := {};
    tmax := max;
    for i from 1 to n do
        vs := [vb[i], op(vs)];
        vc := [vb[i], op(vc)];
    od:
    res := 1;
    dg := 0;
    cyct := 0;
    if linalg[vectdim](pe) = 0 then
       res := res * x[1];
       dg := dg + 1;
    fi:
    if linalg[vectdim](pe) <> 0 then
       ob10 := convert([seq(vs[linalg[vectdim](vs) - ct + 1],
               ct = 1 .. linalg[vectdim](vs))], base, 2, 10);
       if linalg[vectdim](ob10) > 1 then
          m := linalg[vectdim](ob10);
          nsum := 0;
          for i from 1 to m do
              nsum := nsum + ob10[m-i+1]*10^(m-i);
          od:
          ob10 := subsop(1 = nsum, ob10);
       fi:
       if linalg[vectdim](ob10) = 0 then
          res := res*x[1];
          dg := dg+1;
          plist := plist union {0};
       else
          if (member(ob10[1], nlist) = false) and
             (linalg[vectdim](pe) <> 0)
          then
          plist := plist union {ob10[1]};
          nb10 := -1;
          cyct := 0;
          while nb10 <> ob10[1] do
              cyct := cyct + 1;
              for i from 1 to linalg[vectdim](pe) do
                  dcycle := pe[i];
                  clen := linalg[vectdim](dcycle);
                  for j from 1 to clen-1 do
                      vs := 
                      subsop(dcycle[j+1]= vc[dcycle[j]], vs);
                  od;
                  vs := 
                  subsop(dcycle[1] = vc[dcycle[clen]], vs);
                  for k from 1 to n do
                      vc := subsop(k = vs[k], vc);
                  od:
              od:
              plist := plist union {nb10};
              if linalg[vectdim](convert
                  ([seq(vs[linalg[vectdim](vs) - ct + 1],
                      ct = 1 .. linalg[vectdim](vs))],
                      base, 2, 10)) > 1
              then
                 nsum := 0;
                 tmp := convert
                     ([seq(vs[linalg[vectdim](vs) - ct + 1],
                         ct = 1 .. linalg[vectdim](vs))],
                         base, 2, 10);
                 m := linalg[vectdim](tmp);
                 for i from 1 to m do
                     nsum := nsum + tmp[m-i+1]*10^(m-i);
                 od:
                 nb10 := nsum;
              else
                 nb10 := convert
                     ([seq(vs[linalg[vectdim](vs) - ct + 1],
                         ct = 1 .. linalg[vectdim](vs))],
                         base, 2, 10)[1];
              fi:
          od;
          dg := dg + cyct;
          res := res*x[cyct];
          if cyct > tmax then
             tmax := cyct;
          fi:
          fi;
       fi;
    fi;
    RETURN(res, plist, dg, tmax);
end: