#!/usr/bin/perl

# RockNRodl
# Matthew Mullin
# mdmullin@alumni.princeton.edu

# submission #4 : added a large affine plane (125.25.3) and induce
# coverings on it for larger $n

# named after V. Rodl who proved something that implies a solution
# that is asymptopic to N^3/1890
# the absolute lower bound is N^3/7350
# this program gives something like N^3/648; it gives N^3/721 for N=200
# a little bit better on submissions 3 and 4

# &cover generates the cover[][][], coverm[][][] and coverg[][][]
#  arrays according to the simple recursive algorithm
# &coverset actually generates a set of the above parameters
# &lottery computes a optimal lottery for small $n and a good? for any $n
# &lottery uses &lotteryaux to check for the best way to split $n
#  into 3 groups to minimize the solution size
# most other functions are auxillary functions

# GLOBAL VARIABLES:
# cover[][][] = size of cover given by simple recursion
# coverm[][][] =  recursion parameter for cover
# coverg[][][] = amount of 'garbage' - unneeded number
# coverset[][][] = sets already calculated
# ap15_7_3 and ap125_25_3 are affine planes used for coverings

sub fact {             # factorial function
  my ($x)=@_;
  my $res=1;

  while ($x>1) { $res *= $x; $x--; }
  return $res;
}
sub binom {            # binomial coefficient
  my ($n,$k)=@_;
  return &fact($n)/(&fact($k)*&fact($n-$k));
}
sub min {             # minimum of two numbers
  my ($x,$y)=@_;
  if ($x>$y) { return $y;}
  return $x;
}
sub max {             # maximum of two numbers
  my ($x,$y)=@_;
  if ($x>$y) { return $x;}
  return $y;
}
sub floor {           # to replace the broken int function
  my ($x)=@_;
  if ($x>=0) { return int $x; }
  if ($x==int $x) { return $x; }
  return int $x-1;
}
sub ceil {            # replace broken int function
  my ($x)=@_;
  if ($x<=0) { return int $x; }
  if ($x==int $x) { return $x; }
  return 1+int $x;
}
# takes two lists and returns their union; turns a multiset into a set
sub union {
  my ($lr1,$lr2)=@_;
  my @l1=@$lr1;
  my @l2=@$lr2;
  my %l=();

  grep( $l{$_}=1, @l1);
  grep( $l{$_}=1, @l2);
  return keys %l;
}
# takes two lists and returns their intersection; turns a multiset into a set
sub intersect {
  my ($lr1,$lr2)=@_;
  my @l1=@$lr1;
  my @l2=@$lr2;
  my %l1=();
  my %l2=();

  grep( $l1{$_}=1, @l1);
  grep(exists $l1{$_} && ($l2{$_}=1),@l2);
  return keys %l2;
}
# takes two lists and returns their difference; turns a multi into a set
sub diff {
  my ($lr1,$lr2)=@_;
  my @l1=@$lr1;
  my @l2=@$lr2;

  grep( $l{$_}=1, @l1);
  grep( delete $l{$_}, @l2);
  return keys %l;
}


# calculate a cover value for k,l,n and store it and the best m to use
# store the number of garbage elements in @coverg
# will be used for small values of $n -- left over from first draft
sub cover {
  my ($k,$l,$n)=@_;
  my ($m,$best,$bestm,$bestg,$thistry,$garbage);

  if (defined($cover[$k][$l][$n])) {
    return ;
  }
  if ($l==1) {
    $cover[$k][$l][$n]=&ceil($n/$k);
    $coverm[$k][$l][$n]=0;
    $coverg[$k][$l][$n]=$k*&ceil($n/$k)-$n;
  } elsif ($k==$n) {
    $cover[$k][$l][$n]=1;
    $coverm[$k][$l][$n]=0;
    $coverg[$k][$l][$n]=0;
  } elsif ($k==$l) {
    $cover[$k][$l][$n]=&binom($n,$l);
    $coverm[$k][$l][$n]=0;
    $coverg[$k][$l][$n]=0;
  } else {
    $best=1e10;
    $bestg=0;
    $bestm=0;
    $m=&min($k-$l+1,$n-$k);
    for(;$m>0;$m--) {
      &cover($k,$l,$n-$m);
      &cover($k-$m,$l-1,$n-$m);
      $thistry=$cover[$k][$l][$n-$m]+$cover[$k-$m][$l-1][$n-$m];
      $garbage=$coverg[$k][$l][$n-$m]+$coverg[$k-$m][$l-1][$n-$m];
      if ($thistry < $best) {     # first priority: least number
        $best=$thistry;
        $bestm=$m;
	$bestg=$garbage;
      } elsif (($thistry==$best)&&($garbage>$bestg)) {
                                  # second priority: more garbage
        $bestm=$m;
	$bestg=$garbage;
      }
    }
    $cover[$k][$l][$n]=$best;
    $coverm[$k][$l][$n]=$bestm;
    $coverg[$k][$l][$n]=$bestg;
  }
  return ;
}

# returns the cartesian product of two lists of lists
sub setprod {
  my ($x,$y)=@_;
  my ($a,$b,@tmp2);
  my @tmp=();

  foreach $a (@$x) {
    foreach $b (@$y) {
      @tmp2=@$a;
      push @tmp2,@$b;
      push @tmp,[@tmp2];
    }
  }
  return @tmp;
}

# takes k, n and returns a list of all k-subsets of n
sub allsubsets {
  my ($k,$n)=@_;
  my $x;
  my @tmp=();

  if ($k==$n) {
    return (([1..$n]));
  } elsif ($k==1) {
    for($x=1;$x<=$n;$x++) {
      push @tmp,[$x];
    }
    return @tmp;
  } else {
    @tmp=&allsubsets($k,$n-1);
    push @tmp,&setprod([&allsubsets($k-1,$n-1)],[[$n]]);
    return @tmp;
  }
}

# take k and a list and return all the subsets of size k of the list
sub listsubsets {
  my ($k,$lr)=@_;
  my @l=@$lr;

  return &listremap([&allsubsets($k,$#l+1)],[@l]);
}
  
# takes a list and adds a given number to each element, returning the new list
sub trans {
  my ($l,$n)=@_;

  return (map { $_+$n; } @$l);
}

# takes a list of list and trans's each list equally
sub listtrans {
  my ($ll,$n)=@_;

  return (map { [&trans($_,$n)]; } @$ll);
}

# take a mapping and return the inverse mapping
sub invertmap {
  my @l=@_;
  my ($i,@x);

  for($i=0;$i<=$#l;$i++) {
    $x[$l[$i]-1]=$i+1;      # fix off by 1 error
  }
  return @x;
}

# take a list and a mapping and return a new list according to the mapping
sub remap {
  my ($l,$m)=@_;

  return (map { $m->[$_-1]; } @$l);    # off by 1 error here
}

# take a list of lists and a mapping and remap all the lists
sub listremap {
  my ($ll,$m)=@_;

  return (map { [&remap($_,$m)]; } @$ll);
}

 
# afplane computes the affine planes that will be used for covers
# no arguments and sets global variables...
# AP(16,8,4),AP(15,7,3)
sub afplane {
  my ($a,$b,$c,$d,$e,$x,$y,$z,$w,@line,@ap,@ap1,@x); 
# AP(16.8.4), AP(15,7,3)
  @ap=();
  foreach $b (0..1) {
    foreach $c (0..1) {
      foreach $d (0..1) {
        foreach $e (0..1) {
          @line=();
          foreach $y (0..1) {
            foreach $z (0..1) {
              foreach $w (0..1) {
                push @line,(1+(($b*$y)^($c*$z)^($d*$w)^$e)+2*$y+4*$z+8*$w);
          } } }
          push @ap,[@line];
  } } } }
  foreach $c (0..1) {
    foreach $d (0..1) {
      foreach $e (0..1) {
        @line=();
        foreach $x (0..1) {
          foreach $z (0..1) {
            foreach $w (0..1) {
              push @line,(1+$x+2*(($c*$z)^($d*$w)^$e)+4*$z+8*$w);
        } } }
        push @ap,[@line];
  } } }
  foreach $d (0..1) {
    foreach $e (0..1) {
      @line=();
      foreach $x (0..1) {
        foreach $y (0..1) {
          foreach $w (0..1) {
            push @line,(1+$x+2*$y+4*(($d*$w)^$e)+8*$w);
      } } }
      push @ap,[@line];
  } }
  foreach $e (0..1) {
    @line=();
    foreach $x (0..1) {
      foreach $y (0..1) {
        foreach $z (0..1) {
          push @line,(1+$x+2*$y+4*$z+8*$e);
    } } }
    push @ap,[@line];
  }
  @ap1=();
  foreach $x (@ap) {
    @x=&intersect($x,[1..15]);
    if ($#x==6) {
      push @ap1,[@x];
    }
  }
  @ap16_8_4=@ap;
  @ap15_7_3=@ap1;
# set a cover for this
  foreach $a (13..15) {
    $cover[7][3][$a]=15;
    $coverg[7][3][$a]=0;
    @x=();
    foreach $x (@ap15_7_3) {
      push @x,[&intersect($x,[1..$a])];
    }
    $coverset[7][3][$a]=[@x];
  }
#AP(125,25,3)
  @ap=();
  foreach $a (0..4) {
    foreach $b (0..4) {
      foreach $d (0..4) {
        @line=();
        foreach $x (0..4) {
          foreach $y (0..4) {
            push @line,(1+$x+5*$y+25*(($d+$a*$x+$b*$y)%5));
        } }
        push @ap,[@line];
  } } }
  foreach $a (0..4) {
    foreach $d (0..4) {
      @line=();
      foreach $x (0..4) {
        foreach $z (0..4) {
          push @line,(1+$x+5*(($d+$a*$x)%5)+25*$z);
      } }
      push @ap,[@line];
  } }
  foreach $d (0..4) {
    @line=();
    foreach $y (0..4) {
      foreach $z (0..4) {
        push @line,(1+($d%5)+5*$y+25*$z);
    } }
    push @ap,[@line];
  }
  @ap125_25_3=@ap;
}

# randomset generates a random k-subset of 1..n
sub randomset {
  my ($k,$n)=@_;
  my (@set,$i);

  @set=();
  for($i=$n;$i>=1;$i--) {
    if ($i*rand()<$k) {
      unshift @set,$i;
      $k--;
    }
  }
  return @set;
}
# takes $n and pick a random set of $n elements from 1..125 and induce
# a covering on 1..$n from it
# try 5 times
sub induce {
  my ($n)=@_;
  my (@index,$x,@list,$i);
  my ($best,@bestlist);

  @index=&randomset($n,125);
  $best=1e10; @bestlist=();
  foreach $i (1..5) {
    foreach $x (@ap125_25_3) {
      @x=&intersect(\@index,$x);
      if ($#x<2) { next; }
      if ($#x<7) { push @list,[@x]; next; }
      push @list,&listremap([&coverset(7,3,$#x+1)],\@x);
    }
    if ($#list<$best) {
      $best=$#list;
      @bestlist=@list;
    }
  }
  @rindex=&invertmap(@index);
  @bestlist=&listremap([@bestlist],\@rindex);
  return @bestlist;
}
  
# takes k,l,n and returns a list of lists that make a k,l covering of n
# actually some subsets may have less than k elements in them; this is
# because of the 'garbage' elements; they will be added back in at a later time
sub coverset {
  my ($k,$l,$n)=@_;
  my @tmp=();
  my ($x,$m);

  if (defined($coverset[$k][$l][$n])) {
    return @{$coverset[$k][$l][$n]};
  }
  if (($k==7) && ($l==3) && ($n>25)) {
    @tmp=&induce($n);
    $coverset[7][3][$n]=[@tmp];
    return @tmp;
  }
  &cover($k,$l,$n);      # make sure the cover is calculated;
  $m=$coverm[$k][$l][$n];  # get our recursion parameter;
  if ($k==$l) {
    $coverset[$k][$l][$n]=[&allsubsets($k,$n)];
    return @{$coverset[$k][$l][$n]};
  } elsif ($k==$n) {
    $coverset[$k][$l][$n]=[[1..$n]];
    return (([1..$n]));
  } elsif ($l==1) {
    for($x=1;$x*$k<=$n;$x++) {
      push @tmp, [1+$k*($x-1)..$k*$x];
    }
    if ($n%$k!=0) {
      push @tmp, [1+$k*&floor($n/$k)..$n];
    }
    $coverset[$k][$l][$m]=[@tmp];
    return @tmp;
  } else {                 # general recursion step
    @tmp=&coverset($k,$l,$n-$m);
    push @tmp,&setprod([&coverset($k-$m,$l-1,$n-$m)],[[($n-$m+1)..$n]]);
    $coverset[$k][$l][$m]=[@tmp];
    return @tmp;
  }
}

# given $n, check split of $n into 3 approximately equal parts
# check with lower limit of 7 or $n/3+-3
# assume $n>=21
sub lotteryaux {
  my ($n)=@_;
  my ($lowerlimit,$upperlimit,@c,@cg);
  my ($best1,$best2,$best3,$best,$bestg);
  my ($i,$j,$k,$v,$g,@l);

  if ($n>=27) {
    $i=&floor($n/3); $j=&ceil($n/3); $k=$n-$i-$j;
    @l=&coverset(7,3,$i);
    push @l,&listtrans([&coverset(7,3,$j)],$i);
    push @l,&listtrans([&coverset(7,3,$k)],$i+$j);
    return @l;
  }
  $lowerlimit=&max(7,&floor($n/3)-1);
  $upperlimit=&ceil($n/3);
# $upperlimit is actually the upper limit for the first two parts
# the 3rd part is always $n-$part1-$part2 and could go to $n/3+6
# setup an initial best split
  &cover(7,3,$n-2*$lowerlimit);
  @c=@{$cover[7]->[3]};      # copy the interesting parts of the
  @cg=@{$coverg[7]->[3]};    # cover and coverg arrays
# setup the initial best score (lowest)
  $best1=$lowerlimit; $best2=$lowerlimit; $best3=$n-2*$lowerlimit;
  $best=$c[$best1]+$c[$best2]+$c[$best3];
  $bestg=$cg[$best1]+$cg[$best2]+$cg[$best3];
# now check all possibilities
  for($i=$lowerlimit;$i<=$upperlimit;$i++) {
    for($j=$lowerlimit;$j<=$upperlimit;$j++) {
      $k=$n-$i-$j;
      next if ($k<7);
      $v=$c[$i]+$c[$j]+$c[$k];
      next if ($best<$v);
      $g=$cg[$i]+$cg[$j]+$cg[$k];
      next if (($best==$v) && ($g<=$bestg));
# now we have a new best
      $best1=$i; $best2=$j; $best3=$k;
      $best=$v; $bestg=$g;
    }
  }
# now we have our best split - now make the lottery tickets
  @l=&coverset(7,3,$best1);
  push @l,&listtrans([&coverset(7,3,$best2)],$best1);
  push @l,&listtrans([&coverset(7,3,$best3)],$best1+$best2);
  return @l;
}

# take a list of lists and return a tabulation of the elements occuring in it
sub sumlist {
  my ($l,$len)=@_;
  my ($l1,$x,@sum);

  @sum=();
  foreach $l1 (@$l) {
    foreach $x (@$l1) {
      $sum[$x]+=($len+$#{@$l1}+2);
    }
  }
  return @sum;
}

# take a list of lists and renumber the elements so that the sum of all
# of them is as small as possible
sub minimize {
  my @x=@_;
  my (@sum,@y);

  @sum=&sumlist(\@x,$#x);     # total up what elements we use most
  @y=(1..$n);            # sort this wrt @sum to get best mapping
  @y=sort { $sum[$b] <=> $sum[$a] } @y;
  @y=&invertmap(@y);
  return &listremap([@x],[@y]);
}

# take a list of lists and a number for the minimum length of each list
# and add elements to each list to make it that long - add the minimum
# number not already there (mex=Minimum EXcluded)
sub addmex {
  my ($ll,$n)=@_;          # $n will be 7
  my @l=@$ll;              # dereference
  my ($lref,$x,@l2,@ll);
  my %l;

  @ll=();                  # reuse the variable for output
  foreach $lref (@l) {
    @l2=@$lref;
    %l=();
    foreach $_ (@l2) {
      $l{$_}=1;
    }
    for ($x=1;$#l2<$n-1;$x++) {
      next if (defined($l{$x}));
      push @l2,$x;
    }
    push @ll,[@l2];
  }
  return @ll;
}

# now assume $k=7,$l=3, given $n, do the small special cases
# then call lotteryaux; assume $n>=7
sub lottery {
  my ($n)=@_;
  my @x=();
  my (@sum,@y);

  if ($n<=11) {
    @x=([1..$n-4]);
  } elsif ($n<=16) {
    my $x=&floor(($n-2)/2);
    @x=([1..$x],[$x+1..$n-2]);
  } elsif ($n<=21) {
    my ($x,$y)=(&floor($n/3),&ceil($n/3));
    @x=([1..$x],[$x+1..$x+$y],[$x+$y+1..$n]);
  } else {
    @x=&lotteryaux($n);
  }
  @x=&minimize(@x);     # make the sum of elements as low as possible
  @x=&addmex([@x],7);   # make all tickets 7 long
  @x=&minimize(@x);     # minimize the sum of elements again
  return @x;
}

# print's out a list of lists
sub printset {
  my $x;

  foreach $x (@_) {
    print "@$x\n";
  }
}

# @cover holds the size of the covering
# @coverm holds the m-value for the covering for reconstruction
# @coverg holds the number of garbage elements

  srand();
  ($n)=$ARGV[0];
  &afplane;
  &cover(7,3,25);    # these lines set up cover patterns
  @x=&lottery($n);
  &printset(@x);