Quelle  nicequotfinder.g   Sprache: unbekannt

LoadPackage("chop");
LoadPackage("recog");
DeclareInfoClass("InfoNiceQuot");
#agens := AtlasGenerators([ "HS", [ "HSG1-f2r56B0.m1", "HSG1-f2r56B0.m2" ], 1, 2 ]).generators;
#cgens := List(agens,CMat);
#s := AtlasStraightLineProgram([ "HS", "HSG1-max5W1", 1 ]).program;
#hgens := ResultOfStraightLineProgram(s,cgens);
#h := Group(hgens);
gens := AtlasGenerators([ "3.Fi22", [ "3F22G1-f4r27aB0.m1", "3F22G1-f4r27aB0.m2" ], 1, 4 ]).generators;
s := AtlasStraightLineProgram([ "Fi22", "F22G1-max11W1", 1 ]).program;
ggens := ResultOfStraightLineProgram(s,gens);
g := Group(ggens);

ProcessEstimatedGroup := function(estimate,pool,sizes)
  # This function should be called for every successfully estimated group
  local d, stabest, r;
  d := Length(estimate.orbests);
  if d >= 2 then
      if d = 2 then 
          stabest := estimate.stretch;
      else
          stabest := Product(estimate.orbests{[3..d]}) * estimate.stretch;
      fi;
      if stabest >= QuoInt(Minimum(sizes),10) and
         stabest <= Maximum(sizes)*10 then
          r := ShallowCopy(estimate);
          r.orbests := r.orbests{[3..Length(r.orbests)]};
          Add(pool,rec( gens := estimate.stabgens[2],
                        sizeest := stabest,
                        estimate := r,
                        cand := rec( points := estimate.cand.points{[3..d]},
                                     ops := estimate.cand.ops{[3..d]},
                                     used := 0 ) ));
      fi;
  fi;
  if d = 1 then
      stabest := estimate.stretch;
  else
      stabest := Product(estimate.orbests{[2..d]}) * estimate.stretch;
  fi;
  if stabest >= QuoInt(Minimum(sizes),10) and
     stabest <= Maximum(sizes)*10 then
      r := ShallowCopy(estimate);
      r.orbests := r.orbests{[2..Length(r.orbests)]};
      Add(pool,rec( gens := estimate.stabgens[1],
                    sizeest := stabest,
                    estimate := r,
                    cand := rec( points := estimate.cand.points{[2..d]},
                                 ops := estimate.cand.ops{[2..d]},
                                 used := 0 ) ));
  fi;
  return;
end;

LoadPackage("orb");
HappyBirthday := function( gens, pt, op, L, limit, timeout )
  local starttime, endtime, pr, ht, tries, coinc, grpcoinc, S, x, ptx, 
        elm, estimate, upper, lower;
  # Uses the birthday paradox to estimate the size of the orbit
  # of the group G generated by <gens> of the point <pt> with the action
  # given by the operation <op>.
  # Points of the orbit are produced with random elements of G until <L>
  # coincidences are found, or a maximum number of <limit> random elements
  # have been created, or we exceed the time limit of <timeout> ms. 
  # Note that <L> should be chosen at least 15 and <limit> should be 
  # sufficiently large, about sqrt(2*L*N) where N is the size of the orbit.
  # Output if <L> coincidences are found is a record with the components
  #  estimate: the estimated orbit size
  #  lowerbound, upperbound : the lower and upper bounds of the 5% confidence
  #            intervall of the estimate
  #  if <wantSLPs> is true then
  #    stabilizerSLPs : SLPs of stabilizer elements found via the coincidences
  # otherwise fail is output.
  # Addendum:
  # check for coincidences in the stored group elements: if there are very
  # few, i.e. next to none, then the group will be much larger than the orbit.
  # Thus we may assume it to be biggish.
  starttime := Runtime();
  endtime := starttime + timeout;
  pr := ProductReplacer( gens, rec( maxdepth:=400 ) );
  ht := HTCreate( pt, rec( hashlen := NextPrimeInt( Minimum( 100000, limit))) );
  tries := 0;
  coinc := 0;
  grpcoinc := 0;
  S := [];
  while tries <= 3 * limit and coinc < L do
    if ( tries mod 100 = 1 ) and 
       Runtime() >= endtime then
      # check every 100 random elements if we have exceeded the
      # maximum alloted time. If so, bail out.
      Info( InfoNiceQuot, 4, "Time limit exceeded." );
      return fail;
    fi;
      
    x := Next( pr );
    ptx := op( pt, x );
    elm := HTValue( ht, ptx );
    if elm = fail then
      # new
      HTAdd( ht, ptx, x );
    else
      # coincidence
      coinc := coinc + 1;
      # record stabilizer element
      if elm = x then
        # we have found a coincidence in the group
        grpcoinc := grpcoinc + 1;
      else 
        Add( S, elm * x^-1 );
      fi;
    fi;
    tries := tries + 1;
    if tries > limit and coinc < QuoInt( L, 3 ) then
      # if we don't find at least a third of the necessary coincidences
      # within the maximum number of tries, we bail out.
      # Otherwise we allow three times the limit, to find the remaining
      # coincidences.
      Info( InfoNiceQuot, 4, "Too few coincidences found within limit." );
      return fail;
    fi;
  od;
  if coinc = L then
    estimate := QuoInt( tries^2, 2*L );
    Info( InfoNiceQuot, 4, "Estimated orbit size is ", estimate, ".");
    upper := Int( 2*tries^2 / (-MACFLOAT_INT( 196 )/100 +
           SQRT_MACFLOAT( MACFLOAT_INT( 4*L - 1 ) ) )^2
    + MACFLOAT_INT(1)/2 );
    lower := Int( 2*tries^2 / ( MACFLOAT_INT( 196 )/ 100 + 
      SQRT_MACFLOAT( MACFLOAT_INT( 4*L - 1 ) ) )^2 );
    Info( InfoNiceQuot, 4, "Confidence interval with error 5% is [ ",
             lower, ", ", upper," ]." );
    #if grpcoinc = 0 then
    #  Info( InfoNiceQuot, 4, "Huge group with small orbit. Failing." );
    #  return fail;
    #else
      return rec( estimate := estimate, 
                  lowerbound := lower, 
                  upperbound := upper, 
                  Sgens := S,
                  grpcoinc := grpcoinc );
    #fi;
  else
    return fail;
  fi;
end;

BirthdayStatisticsOnVectors := function( gens, size, time )
  local starttime, endtime, v, tries, estimate, testrun;
  # Given generators <gens> of a matrix group G of known size <size>
  # the average estimated length of orbits of G on random vectors
  # is computed for the duration of <time> ms. For the estimate we
  # use 15 coincidences.
  # If no estimate was computable, it returns fail. Otherwise
  # a record with components
  #  estimate: the average of orbit length estimates
  #  tries:    the number of times an orbit length was estimated
  # is returned.
  starttime := Runtime();
  endtime := starttime + time;
  v := ShallowCopy( gens[1][1] );
  tries := 0;
  estimate := 0;
  while Runtime() < endtime do
    Randomize( v );
    testrun := HappyBirthday( gens, v, OnRight, 15, 
                   Maximum( 100, 
                      6*Int( SQRT_MACFLOAT( MACFLOAT_INT( size ) ) ) ),
                   time );
    # As the orbit size N is approx the sample size squared divided by
    # twice the number of coincidences, and we want 15 coincidences, the
    # sample size has to be at least sqrt(30N). Here we take the sample
    # size to be 6sqrt(N) for good measure.
    if testrun = fail then
      Info( InfoNiceQuot, 4, 
            "FAILED orbit size estimate. Too few coincidences." );
    else
      Info( InfoNiceQuot, 4,
            "Estimated orbit size to be ", testrun.estimate,"." );
      tries := tries + 1;
      estimate := estimate + testrun.estimate;
    fi;
  od;
  if tries > 0 then
    estimate := Int( estimate / tries );
    return rec( estimate := estimate, tries := tries );
  else
    return fail;
  fi;
end;

BirthdayStatisticsOnLines := function( gens, size, time )
  local starttime, endtime, v, tries, estimate, testrun;
  # Given generators <gens> of a matrix group G of known size <size>
  # the average estimated length of orbits of G on random 1-dimensional
  # subspaces (i.e. lines) is computed for the duration of <time> ms. 
  # For the estimate we use 15 coincidences.
  # If no estimate was computable, it returns fail. Otherwise
  # a record with components
  #  estimate: the average of orbit length estimates
  #  tries:    the number of times an orbit length was estimated
  # is returned.
  starttime := Runtime();
  endtime := starttime + time;
  v := ShallowCopy( gens[1][1] );
  tries := 0;
  estimate := 0;
  while Runtime() < endtime do
    Randomize( v );
    ORB_NormalizeVector( v );
    testrun := HappyBirthday( gens, v, OnLines, 15, 
                 Maximum( 100, 
                   6*Int( SQRT_MACFLOAT( MACFLOAT_INT( size ) ) ) ),
                   time );
    # As the orbit size N is approx the sample size squared divided by
    # twice the number of coincidences, and we want 15 coincidences, the
    # sample size has to be at least sqrt(30N). Here we take the sample
    # size to be 6sqrt(N) for good measure.
    if testrun = fail then
      Info( InfoNiceQuot, 4, 
            "FAILED orbit size estimate. Too few coincidences." );
    else
      Info( InfoNiceQuot, 4,
            "Estimated orbit size to be ", testrun.estimate,"." );
      tries := tries + 1;
      estimate := estimate + testrun.estimate;
    fi;
  od;
  if tries > 0 then
    estimate := Int( estimate / tries );
    return rec( estimate := estimate, tries := tries );
  else
    return fail;
  fi;
end;

SmallishGroupSize := function( gens, limit )
  local v, estimate, cand, H, S, stabsize, s, t;
  # <gens> generators
  # <limit> upper bound above which we no longer consider the group
  #        to be smallish
  # 1. take random vector and use birthday paradox to estimate orbit length
  # 2. if orbit length is short, do stabilizer chain using this orbit on
  #    first level
  #    if orbit length is long, take a new random vector and repeat step 1
  # 3. after 3 vectors, say, give up and consider the group to be biggish

  for t in [1..3] do
    v := ShallowCopy( gens[1][1] );
    Randomize( v );
    ORB_NormalizeVector( v );
    estimate := HappyBirthday( gens, v, OnLines, 15, 
                 Maximum( 100, 
                   6*Int( SQRT_MACFLOAT( MACFLOAT_INT( limit ) ) ) ),
                   5000 );
    if estimate = fail then
      Info(InfoNiceQuot,4,"SmallishGroupSize: Estimate failed.");
      continue;
    elif estimate.estimate <= limit then
      cand := rec( points := [v], ops := [OnLines], used := 0 );
      if estimate.grpcoinc > 0 then
        # orbit is essentially regular
        H := GroupWithGenerators(gens);
        S := StabilizerChain( H, rec( Cand := cand ) );
        return rec( size := Size(S), S := S );
      else
        Info( InfoNiceQuot,3,"Estimated orbit size is ",
              estimate.estimate,"; recursing down to stabilizer.");
        stabsize := SmallishGroupSize( estimate.Sgens, limit );
        if stabsize = fail then
          # Stabilizer itself is not smallish
          Info( InfoNiceQuot,3,"Group is most probably biggish." );
          return fail;
        else
          s := stabsize.S;
          while s <> false do
            Add( cand.points, s!.orb[1] );
            Add( cand.ops, s!.orb!.op );
            s := s!.stab;
          od;
          H := GroupWithGenerators(gens);
          S := StabilizerChain( H, rec( Cand := cand ) );
          return rec( size := Size(S), S := S );
        fi;
      fi;
    else
      continue;
    fi;
  od;
  Info( InfoNiceQuot,3,"Group is probably biggish." );
  return fail;
end;
      
EstimateStabilizerChain := function( gens, maxorblength )
  # With incredible ingenuity we estimate a stabilizer chain
  # for the group given by its generators <gens>. The nonnegative
  # integer <maxorblength> is an upper bound for the orbitlengths we 
  # allow on any level of the stabilizer chain.
  # On each level we allow three attempts to estimate the orbit length.
  # Output:
  # A record with components
  # orbests:  a list of estimated orbit lengths, one for each level
  # estimate: the estimated group size
  # cand:     a list of records giving candidates for base points
  # stabgens: a list of lists, giving generators for the stabilisers
  # stretch:  the stretch factor for the lowest orbit obtained from
  #           group coincidences
  # IsStabilierChainEstimateRecord: set to true to allow nice viewing
  #           method.
  local v, estimate, cand, orbests, stabgens, stabest, t;
  for t in [1..3] do
    v := ShallowCopy( gens[1][1] );
    Randomize( v );
    ORB_NormalizeVector( v );
    estimate := HappyBirthday( gens, v, OnLines, 15, 
                Maximum( 100, 
                  6*Int( SQRT_MACFLOAT( MACFLOAT_INT( maxorblength ) ) ) ),
                  5000 );
    # magic numbers: 
    # 15 is the number of coincidences look for in order to make use of
    # the birthday paradox. For this we need to produce at least about
    # sqrt(30 * orbit length) point of the orbit, but at least 100 points.
    # Additionally we allow a maximum of 5 seconds to be patient.
    if estimate = fail then
      continue;
    else
      if estimate.grpcoinc > 0 then
        # orbit is hopefully essentially regular
        cand := rec( points := [ v ], ops := [ OnLines ], used := 0 );
        return rec( orbests := [ estimate.estimate ],
                    estimate := QuoInt( 15 * estimate.estimate, 
                                               estimate.grpcoinc),
                    cand := cand,
                    stabgens := [estimate.Sgens],
                    stretch := QuoInt(15, estimate.grpcoinc),
                    IsStabilizerChainEstimateRecord:=true );
      else # no group coincidences, have stabilizer
        stabest := EstimateStabilizerChain( estimate.Sgens, maxorblength );
        Add( stabest.orbests, estimate.estimate, 1 );
        Add( stabest.cand.points, v, 1 );
        Add( stabest.cand.ops, OnLines, 1 );
        Add( stabest.stabgens, estimate.Sgens, 1 );
        if not stabest.estimate = fail then
          stabest.estimate := estimate.estimate * stabest.estimate;
        fi;
        return stabest;
      fi;
    fi;
  od;
  # rockbottom!
  return rec( estimate := fail,
              orbests := [], cand := rec( points:=[], ops := [], used := 0 ), 
              stabgens := [], 
              stretch := fail,
              IsStabilizerChainEstimateRecord := true );
end;
              
      

BlackJackMeth := function( r, pr, sizes )
  local max, min, hand, value, counter, i;
  max := Maximum( sizes );
  min := Minimum( sizes );
  for i in [1..3] do
    # pick up new hand
    hand := [ Next( pr ), Next( pr ) ];
    #value := SmallishGroupSize( hand, max );
    value := EstimateStabilizerChain( hand, max );
    if value.estimate <> fail then
      Info(InfoNiceQuot,2,"BlackJack: Found group of size ",value.estimate);
          ProcessEstimatedGroup(value,r.grppool,sizes);
      if value.estimate < 10 * max then
        r.S := StabilizerChain(Group(hand),rec(Cand := value.cand));
        value.estimate := Size(r.S); #if we know size, use it
        if Size(r.S) in sizes then
          r.hsize := Size(r.S);
          r.hgens := hand;
          return true;
        fi;
      else
        r.S := fail;
      fi;
      
      if value.estimate < min then
        counter := 1;
        repeat
          Info( InfoNiceQuot, 3, "Hit me!" );
          Add( hand, Next(pr) );
          #value := SmallishGroupSize( hand, max );
          value := EstimateStabilizerChain( hand, max );
          if value.estimate = fail then
            # cannot estimate, discard card
            hand := hand{[1..Length(hand)-1]};
          elif value.estimate < min then
            # if we are still too small, try again
            value.estimate := fail;
          fi;
          counter := counter + 1;
        until value.estimate <> fail or counter = 3;
      fi;
      if value.estimate = fail then
        # tried 3 times to enlarge without being able to estimate
        # make completely new attempt
        Info( InfoNiceQuot,3,"I fold." );
        continue;
      else
        # have enlarged and estimated
        ProcessEstimatedGroup( value, r.grppool, sizes );
        if value.estimate <= 10 * max then
          r.S := StabilizerChain( Group(hand), rec( Cand := value.cand ) );
          value.estimate := Size(r.S); # if we know size, use it.
          if Size(r.S) in sizes then
            r.hgens := hand;
            r.hsize := value.estimate;
            #r.S := value.S;
            return true;
          fi;
        fi;
      fi;
    else
    Info(InfoNiceQuot,3,
         "BlackJack: Could not sensibly estimate group size");
    fi;
  od;
  # nothing worked
  Info( InfoNiceQuot, 2, "BlackJack: BUSTED!" );
  return fail;
end;

UpAndDownMeth := function( r, pr, sizes )
  local max, min, x, cent, tol, guck, y, o, p, pr2, counter, xx;
  # Going Up and Down: Use involution centralizer method to go down
  # if we end up too small use black jack to go up, if we are to large
  # go down again by taking an involution centralizer in this grp etc.
  max := Maximum(sizes);
  min := Minimum(sizes);
  x := RECOG.InvolutionSearcher(pr,Order,20);
  if x = fail then 
      Info(InfoNiceQuot,2,"UpAndDown: Did not find involution.");
      return false; 
  fi;
  # Now x is an involution in G
  cent := RECOG.InvolutionCentraliser(pr,Order,x,5);
  tol := 10;
  while tol > 0 do
    tol := tol - 1;
    #guck := SmallishGroupSize(cent,max);
    guck := EstimateStabilizerChain( cent, max );
    if guck.estimate <> fail then
      Info(InfoNiceQuot,2,"UpAndDown: Current size=",guck.estimate);
      ProcessEstimatedGroup( guck, r.grppool, sizes );      
      if guck.estimate < 10 * max then
         r.S := StabilizerChain(Group(cent),rec(Cand := guck.cand));
         guck.estimate := Size(r.S); #if we know the group size, we use it
         if Size(r.S) in sizes then
           Info(InfoNiceQuot,3,"UpAndDown: Success!");
           r.hsize := Size(r.S);
           r.hgens := cent;
           # r.S := guck.stabgens;
           return true;
         fi;
      fi;
      # we might be too small
      if guck.estimate < min then
        Info(InfoNiceQuot,3,"UpAndDown: Going up...");
        # black jack method
        repeat
          y := Next(pr);
          o := Order(y);
        until o > 1;
        if IsEvenInt(o) then
          y := y^(o/2);
        elif o mod 3 = 0 then
          y := y^(o/3);
        else
          p := Factors(o)[1];
          y := y^(o/p);
        fi;
        Add( cent, y );
        continue;
      fi;
    else
      Info(InfoNiceQuot,2,"UpAndDown: Could not estimate group size");
    fi;
    # we are too large! Have to go down again.
    # If we get here, the involution centraliser seems to be too big
    Info(InfoNiceQuot,3,"UpAndDown: Going down...");
    pr2 := ProductReplacer(cent,
                           rec( scramble := 10,scramblefactor := 1,
                                maxdepth := 400 ));
    x := RECOG.InvolutionSearcher(pr2,Order,20);
    if x = fail then 
        Info(InfoNiceQuot,2,"UpAndDown: Did not find involution.");
        return false; 
    fi;
    if ForAll(cent,y->x*y=y*x) then
      if ForAll(cent,a->ForAll(cent,b->a*b=b*a)) then
        Info(InfoNiceQuot,3,"UpAndDown: Ran into abelian group!");
        return fail;
      fi;
      counter := 0;
      repeat
        counter := counter + 1;
        xx := RECOG.InvolutionSearcher(pr2,Order,20);
        if xx = fail then 
          Info(InfoNiceQuot,3,
               "UpAndDown: Could not find involution in group");
          return fail; 
        fi;
        if InfoLevel(InfoNiceQuot) >= 3 then Print(".\c"); fi;
      until not(ForAll(cent,y->xx*y=y*xx)) or counter = 20;
      x := xx;
    fi;
    cent := RECOG.InvolutionCentraliser(pr2,Order,x,5);
  od;
  Info(InfoNiceQuot,2,"UpAndDown: Lost patience");
  return fail;
end;



FindPartSets := function(a,sums)
  local p, l, max, DoRecurse, res;
  # a a list of integers, sums a set of integers
  # this function finds all subsets of the indices of l such that
  # the sum of the numbers in these positions is in the set sums
  a := ShallowCopy(a);
  p := Sortex(a);
  a := Reversed(a);  # now ascending
  l := Length(a);
  max := Maximum(sums);

  DoRecurse := function(res,sub,d,m,s)
    local x, i;
    d := d + 1;
    for i in [m..l] do
      if a[i] > 0 then   # ignore zeros
        x := a[i];
        sub[d] := i;
        s := s + x;
        if s <= max then
            if s in sums then
                Add(res,ShallowCopy(sub));
            fi;
            DoRecurse(res,sub,d,i+1,s);
        fi;
        s := s - x;
      fi;
    od;
    Unbind(sub[d]);
  end;
  res := [];
  DoRecurse(res,[],0,1,0);
  if Length(res) = 0 then return fail; fi;
  return List(res,v->OnTuples(List(v,x->l+1-x),p^-1));
end;
  

EvaluateCandidateGroup := function(Hgens,Hsize,codims,avorbmin)
  local f, q, d, m, r, l, dims, di, pos1, pos2, gens, c, qgens, tol, 
        ti, avorbest, rad, laydims, laydimsums, min, presum, FindDim, 
        localdims, res, basind, bas, basi, i, re, count;
  # Hgens must generate a matrix group of size Hsize in GL(d,q)
  # This function tries to find a submodule of the natural module
  # such that its codimension c has c in codims and such
  # that the average orbit length in the quotient are estimated to be
  # bigger than avorbmin.

  f := BaseDomain(Hgens[1]);
  q := Size(f);
  d := Length(Hgens[1]);
  
  m := Module(Hgens);
  r := Chop(m,rec(compbasis := true));
  l := Length(r.acs);
  dims := 0*[1..l];
  di := 0;
  for i in [1..l] do
      di := di + Dimension(r.db[r.acs[i]]);
      dims[i] := di;
  od;
  Info(InfoNiceQuot,2,"Dimensions in composition series: ",dims);
  pos1 := First([1..l],i->d-dims[i] in codims);
  pos2 := First([l,l-1..1],i->d-dims[i] in codims);
  if pos1 <> fail then
      r.basis := Reversed(r.basis);
      r.ibasis := r.basis^-1;
      gens := List(Hgens,x->r.basis * x * r.ibasis);
      for i in [pos2,pos2-1..pos1] do
          c := d-dims[i];
          Info(InfoNiceQuot,2,"Trying codimension ",c);
          qgens := List(gens,x->ExtractSubMatrix(x,[1..c],[1..c]));
          tol := 5; ti := 1000;
          repeat
              avorbest := BirthdayStatisticsOnVectors(qgens,Hsize,ti);
              if avorbest = fail or avorbest.tries < 5 then
                  Info(InfoNiceQuot,3,"...doubling time limit to ",ti);
                  ti := ti * 2;
              fi;
              tol := tol - 1;
          until (avorbest <> fail and avorbest.tries >= 5)
                or tol = 0;
          if avorbest <> fail and 
             avorbest.tries >= 5 then
              Info(InfoNiceQuot,2,"found average orbit length estimate of ",
                   avorbest.estimate," (",avorbest.tries," tries)");
              if avorbest.estimate >= avorbmin then
                  return rec( avorbest := avorbest.estimate, codim := c,
                              basis := r.basis, ibasis := r.ibasis,
                              choprec := r, compseriesmember := i,
                              quotgens := qgens );
              fi;
          fi;
          Info(InfoNiceQuot,2,"Do not like codimension ",c);
      od;
      Info(InfoNiceQuot,2,
           "Giving up on initial composition series, now trying radical");
  else
      Info(InfoNiceQuot,2,
        "No suitable quotients in original comp series, now trying radical.");
  fi;
  rad := RadicalSeries(m,r.db);
  l := Length(rad.isotypes);
  dims := List(rad.isotypes,v->List(v,j->Dimension(rad.db[j])));
  Info(InfoNiceQuot,2,"Dimensions in radical series: ",dims);
  laydims := List(dims,Sum);
  laydimsums := List([1..l],i->Sum(laydims{[1..i]}));
  pos1 := 1;
  min := Minimum(codims);
  while laydimsums[pos1] < min and pos1 <= l do pos1 := pos1 + 1; od;
  if pos1 > l then
      Info(InfoNiceQuot,2,
           "No suitable quotients in radical series, giving up.");
      return fail;
  fi;
  if pos1 = 1 then
      presum := 0;
  else
      presum := laydimsums[pos1-1];
  fi;

  FindDim := function(m)
    if ForAll(RepresentingMatrices(m),IsOne) then 
      return 0;
    else
      return Dimension(m);
    fi;
  end;

  # Ignore trivial summands by setting their dimension to 0:
  localdims := List(rad.isotypes[pos1],i->FindDim(rad.db[i]));
  res := FindPartSets(localdims,List(codims,x->x-presum));
  if res = fail then
      Info(InfoNiceQuot,2,
           "No suitable quotients in layer of radical series, giving up.");
      return fail;
  fi;
  count := 0;
  for re in res do
      count := count + 1;
      if count > 30 then return fail; fi;
      Info(InfoNiceQuot,2, "Trying radical layer cut ",localdims{re});
      basind := [];
      for i in [1..pos1-1] do
          for r in rad.cfposs[i] do
              Append(basind,r);
          od;
      od;
      for r in re do
          Append(basind,rad.cfposs[pos1][r]);
      od;
      c := Length(basind);
      for r in Difference([1..Length(rad.cfposs[pos1])],re) do;
          Append(basind,rad.cfposs[pos1][r]);
      od;
      for i in [pos1+1..l] do
          for r in rad.cfposs[i] do
              Append(basind,r);
          od;
      od;
      
      bas := rad.basis{basind};
      basi := bas^-1;
      gens := List(Hgens,x->bas * x * basi);

      qgens := List(gens,x->ExtractSubMatrix(x,[1..c],[1..c]));
      tol := 5; ti := 1000;
      repeat
          avorbest := BirthdayStatisticsOnVectors(qgens,Hsize,ti);
          if avorbest = fail or avorbest.tries < 5 then
              Info(InfoNiceQuot,3,"...doubling time limit to ",ti);
              ti := ti * 2;
          fi;
          tol := tol - 1;
      until (avorbest <> fail and avorbest.tries >= 5)
            or tol = 0;

      if avorbest <> fail and 
         avorbest.tries >= 5 then
          Info(InfoNiceQuot,2,"found average orbit length estimate of ",
               avorbest.estimate," (",avorbest.tries," tries)");
          if avorbest.estimate >= avorbmin then
              return rec( avorbest := avorbest.estimate, codim := c,
                          basis := bas, ibasis := basi,
                          choprec := r, rad := rad, basind := basind,
                          quotgens := qgens );
          fi;
      fi;
  od;

  Info(InfoNiceQuot,2,"Radical series did not work as well, giving up.");
  return fail;
end;

TwoRandElsMeth := function(r,pr,sizes)
  local max, hgens, guck, i;
  max := Maximum(sizes);
  for i in [1..10] do
      hgens := [Next(pr),Next(pr)];
      guck := EstimateStabilizerChain(hgens,1000000);
      if guck.estimate <> fail then
          Info(InfoNiceQuot,2,"TwoRandEls: Found group of estimated size ",
               guck.estimate);
          ProcessEstimatedGroup(guck,r.grppool,sizes);
          if guck.estimate < 10 * max then
              r.S := StabilizerChain(Group(hgens),rec(Cand := guck.cand));
              if Size(r.S) in sizes then
                  r.hsize := Size(r.S);
                  r.hgens := hgens;
                  return true;
              fi;
          fi;
      else
          Info(InfoNiceQuot,3,
               "TwoRandEls: Could not sensibly estimate group size");
      fi;
  od;
  return fail;
end;

DihedralMeth := function(r,pr,sizes)
  local tol, x, y, z, o, d;
  tol := 0;
  repeat
      x := RECOG.InvolutionSearcher(pr,Order,20);
      if x = fail then
          Info(InfoNiceQuot,2,"Dihedral: Did not find involution");
      fi;
      y := RECOG.InvolutionSearcher(pr,Order,20);
      if y = fail then
          z := Next(pr);
          y := x^z;
      fi;
      z := x*y;
      o := Order(z);
      if 2*o in sizes then
          Info(InfoNiceQuot,2,"Dihedral: Success, order ",2*o);
          r.hgens := [x,y];
          r.hsize := 2*o;
          d := Group(x,y);
          SetSize(d,2*o);
          r.S := StabilizerChain(d);
          return true;
      fi;
      Info(InfoNiceQuot,2,"Dihedral: Failure, order ",2*o);
      tol := tol + 1;
  until tol > 10;
  return fail;
end;

InstallMethod( ViewObj, "for stabilizer chain estimate records",
  [ IsRecord ],
  function(r)
    if not IsBound(r.IsStabilizerChainEstimateRecord) then
       TryNextMethod();
    fi;
    Print("<StabilizerChainEstimateRecord estimate:",r.estimate,
          "\n    orbit estimates:",r.orbests,
          "\n    stretch:",r.stretch,">");
  end );

InvCentMeth := function(r,pr,sizes)
  local max, x, pr2, tol, cent, guck, S, y, o, p, count, xx;
  max := Maximum(sizes);
  x := RECOG.InvolutionSearcher(pr,Order,20);
  if x = fail then 
      Info(InfoNiceQuot,3,"InvCent: Did not find involution.");
      return false; 
  fi;
  # Now x is an involution in G
  pr2 := pr;
  tol := 10;
  while tol > 0 do
      tol := tol - 1;
      Info(InfoNiceQuot,2,"InvCent: Computing an involution centraliser...");
      cent := RECOG.InvolutionCentraliser(pr2,Order,x,5);
      guck := EstimateStabilizerChain(cent,max);
      if guck.estimate <> fail then
          Info(InfoNiceQuot,2,"InvCent: Estimated size: ",guck.estimate);
          ProcessEstimatedGroup(guck,r.grppool,sizes);
          if guck.estimate <= 100 * max then
              S := StabilizerChain(Group(cent),rec(Cand := guck.cand));
              Info(InfoNiceQuot,2,"InvCent: Found size ",Size(S));
              if Size(S) in sizes then
                  Info(InfoNiceQuot,2,"InvCent: Success: ",Size(S));
                  r.hgens := cent;
                  r.hsize := Size(S);
                  r.S := S;
                  return true;
              fi;
          else
              S := fail;
          fi;
          if S <> fail and Size(S) <= max then
              Info(InfoNiceQuot,2,"InvCent: Involution centraliser too small: ",
                   Size(S)," adding random element.");
              repeat
                y := Next(pr);
                o := Order(y);
              until o > 1;
              if IsEvenInt(o) then
                y := y^(o/2);
              elif o mod 3 = 0 then
                y := y^(o/3);
              else
                p := Factors(o)[1];
                y := y^(o/p);
              fi;
              Add( cent, y );
          fi;
      else
          Info(InfoNiceQuot,2,"InvCent: Could not estimate inv centraliser");
      fi;
      # If we get here, the involution centraliser seems to be too big
      pr2 := ProductReplacer(cent,rec( scramble := 10,scramblefactor := 1,
                                       maxdepth := 400 ));
      x := RECOG.InvolutionSearcher(pr2,Order,20);
      if x = fail then 
          Info(InfoNiceQuot,2,"InvCent: Did not find involution.");
          return false; 
      fi;
      if ForAll(cent,y->x*y=y*x) then
          if ForAll(cent,a->ForAll(cent,b->a*b=b*a)) then
              Info(InfoNiceQuot,2,"InvCent: Ran into abelian group, baeh!");
              return fail;
          fi;
          count := 0;
          repeat
              #xx := RECOG.InvolutionJumper(pr2,Order,x,100,true);
              xx := RECOG.InvolutionSearcher(pr2,Order,20);
              if xx = fail then 
                  Info(InfoNiceQuot,3,
                       "InvCent: Could not find involution in invcent");
                  return fail; 
              fi;
              count := count + 1;
          until not(ForAll(cent,y->xx*y=y*xx)) or count > 20;
          if count > 20 then 
              Info(InfoNiceQuot,2,
                   "InvCent: Did not find noncentral involution");
              return fail;
          fi;
          x := xx;
      fi;
  od;
  Info(InfoNiceQuot,2,"InvCent: Lost patience");
  return fail;
end;

# For solvable groups:
DerivedMeth := function(r,pr,sizes)
  local max, min, der, tol, guck, S, y, o, p, pr2, i;

  max := Maximum(sizes);
  min := Minimum(sizes);

  # Go down one step in derived series:
  der := [];
  for i in [1..5] do
    Add(der,Comm(Next(pr),Next(pr)));
  od;
  if ForAll(der,IsOne) then
      Info(InfoNiceQuot,2,"Derived: Group is abelian, giving up.");
      return false;
  fi;

  tol := 10;
  while tol > 0 do
    tol := tol - 1;
    guck := EstimateStabilizerChain(der,max);
    if guck.estimate <> fail then
      Info(InfoNiceQuot,2,"Derived: Estimated size: ",guck.estimate);
      ProcessEstimatedGroup(guck,r.grppool,sizes);
      if guck.estimate <= 100 * max then
          S := StabilizerChain(Group(der),rec(Cand := guck.cand));
          Info(InfoNiceQuot,2,"Derived: Found size ",Size(S));
          if Size(S) in sizes then
              Info(InfoNiceQuot,3,"Derived: Success!");
              r.hgens := der;
              r.hsize := Size(S);
              r.S := S;
              return true;
          fi;
      else
          S := fail;
      fi;
      # we might be too small
      if S <> fail and Size(S) < min then
        Info(InfoNiceQuot,2,"Derived: Going up...");
        repeat
          y := Next(pr);
          o := Order(y);
        until o > 1;
        if IsEvenInt(o) then
          y := y^(o/2);
        elif o mod 3 = 0 then
          y := y^(o/3);
        else
          p := Factors(o)[1];
          y := y^(o/p);
        fi;
        Add( der, y );
        continue;
      fi;
    else
      Info(InfoNiceQuot,2,"Derived: Could not estimate group size");
    fi;
    # we are too large! Have to go down again.
    # If we get here, the group seems to be too big
    Info(InfoNiceQuot,3,"Derived: Going down in derived series...");
    pr2 := ProductReplacer(der,
                           rec( scramble := 10,scramblefactor := 1,
                                maxdepth := 400 ));
    # Go down one step in derived series:
    der := [];
    for i in [1..5] do
      Add(der,Comm(Next(pr),Next(pr)));
    od;
    if ForAll(der,IsOne) then
        Info(InfoNiceQuot,2,"Derived: Group is abelian, giving up.");
        return fail;
    fi;

  od;
  Info(InfoNiceQuot,2,"Derived: Lost patience");
  return fail;
end;


FindSmallishSubgroupMethodsDB := [];

AddMethod( FindSmallishSubgroupMethodsDB,
           DihedralMeth, 130, "Dihedral",
           "try some involutions to find a dihedral subgroup" );
AddMethod( FindSmallishSubgroupMethodsDB,
           UpAndDownMeth, 120, "UpAndDown",
           "try to go up and down and up and down" );
AddMethod( FindSmallishSubgroupMethodsDB,
           InvCentMeth, 110, "InvCent",
           "try repeated involution centralisers" );
AddMethod( FindSmallishSubgroupMethodsDB,
           DerivedMeth, 105, "Derived",
           "try repeated derived subgroups" );
AddMethod( FindSmallishSubgroupMethodsDB,
           BlackJackMeth, 100, "BlackJack",
           "play a round of black jack with group elements" );
AddMethod( FindSmallishSubgroupMethodsDB,
           TwoRandElsMeth, 90, "TwoRandEls",
           "try group generated by two random elements" );

FindSmallishSubgroup := function(G,sizes,timeout)
  local starttime, gensm, pr, r, res;
  starttime := Runtime();

  gensm := GeneratorsWithMemory(GeneratorsOfGroup(G));
  pr := ProductReplacer(gensm,rec( maxdepth := 400, scramblefactor := 1,
                                   scramble := 10 ));
  while true do
      r := rec(grppool := []);
      res := CallMethods(FindSmallishSubgroupMethodsDB,3,
                         r,pr,sizes);
      if res.result = true then
          r.methsel := res;
          return r;
      fi;
      if Runtime() - starttime > timeout then
          return fail;
      fi;
  od;
end;

RingelPiez := function(G,sizes,timeout)
  local starttime, gensm, pr, i, pool, grps, r, res;
  starttime := Runtime();

  gensm := GeneratorsWithMemory(GeneratorsOfGroup(G));
  pr := ProductReplacer(gensm,rec( maxdepth := 400, scramblefactor := 1,
                                   scramble := 10 ));
  i := 1;
  pool := [];
  grps := [];
  while true do
      r := rec(grppool := []);
      Info(InfoNiceQuot,1,"RingelPiez: Doing ",
           FindSmallishSubgroupMethodsDB[i].stamp);
      res := FindSmallishSubgroupMethodsDB[i].method(r,pr,sizes);
      if res = true then
          Append(pool,r.grppool);
          Unbind(r.grppool);
          Add(grps,r);
      fi;
      i := i + 1;
      if i > Length(FindSmallishSubgroupMethodsDB) then
          i := 1;
      fi;
      if Runtime() - starttime > timeout then
          return rec( groups := grps, pool := pool, G := G, sizes := sizes,
                      timeout := timeout );
      fi;
  od;
end;
  
AnalysiereRingelPiezResult := function(r)
  local q, d, h, cgens, codims, i;
  r.spar := [];
  r.sparfactor := [];
  q := Size(FieldOfMatrixGroup(r.G));
  d := DimensionOfMatrixGroup(r.G);
  for i in [1..Length(r.groups)] do
      h := r.groups[i];
      cgens := List(h.hgens,x->CMat(StripMemory(x)));
      codims := [First([1..d],i->q^i > h.hsize)..
                 First([d,d-1..1],i->q^i < 10000000)];
      if Length(codims) > 0 then
        r.spar[i] := EvaluateCandidateGroup(cgens,h.hsize,codims,
                                            QuoInt(h.hsize,100));
        Info(InfoNiceQuot,1,"Allowing codimensions ",codims);
      else
        Info(InfoNiceQuot,1,"Bad luck, no possible codimensions!");
        r.spar[i] := fail;
      fi;
      if r.spar[i] = fail then
          Info(InfoNiceQuot,1,"Wir sparen nix!");
          r.sparfactor[i] := fail;
      else
          Info(InfoNiceQuot,1,"Wir sparen ",r.spar[i].avorbest,"!");
          r.sparfactor[i] := r.spar[i].avorbest;
      fi;
  od;
end;

G := Group(());

GuckManyGroups := function(dir,timepergroup)
  local l, ll, i, r, ff, fff, f, ra;
  l := Set(IO_ListDir(dir));
  ll := ShallowCopy(l);
  RemoveSet(l,".");
  RemoveSet(l,"..");
  i := 1;
  while i <= Length(l) do
      if l[i]{[Length(l[i])-1..Length(l[i])]} <> ".g" then
          Remove(l,i);
      else
          i := i + 1;
      fi;
  od;
  for f in l do
    ff := Concatenation(f,".result");
    if ff in ll then 
        Print("Have already done ",f,"\n");
        continue; 
    fi;
    ff := Concatenation(dir,"/",ff);
    Print("Doing ",f,"\n");
    Read(Concatenation(dir,"/",f));
    r := RingelPiez(G,[1000..1000000],timepergroup);
    AnalysiereRingelPiezResult(r);
    ra := Filtered([1..Length(r.groups)],i->r.sparfactor[i] <> fail);
    PrintTo(ff,"Group sizes: ",List(r.groups{ra},x->x.hsize),"\n",
               "Saving factors: ",r.sparfactor{ra},"\n");
    Print("P\c");
    ff := Concatenation(ff,".gp");
    fff := IO_File(ff,"w");
    IO_Pickle(fff,List(r.groups{ra},x->SLPOfElms(x.hgens)));
    IO_Pickle(fff,List(r.groups{ra},x->x.hsize));
    IO_Pickle(fff,r.sparfactor{ra});
    IO_Close(fff);
    Print("\rSaving factors: ",r.sparfactor,"\n");
  od;
end;

Unbind(G);