#!perl -w
# vim: ts=4 sw=4 et
$|=1;

# about volumes:
#  http://mathforum.org/library/drmath/view/63013.html
# about areas:
#  http://mathforum.org/library/drmath/view/52053.html
#  http://www.btinternet.com/~se16/hgb/triangle.htm
#
# platonic solids:
#   http://davidf.faricy.net/polyhedra/Platonic_Solids.html
#   http://davidf.faricy.net/polyhedra/Polytopes.html
#
# overview of generation tactics:
#
#        0   1   2   3   4   5   6  ..  n-1  n <-- subdim
# tet :  h   ggggggggggggggggggggggggggggg   1
# cub :  h   ggggggggggggggggggggggggggggg   1
# oct :  h   g   b   b   xxxxxxxxxxxxxxxxx   1
# dod :  h   b   h   1   -----------------
# ico :  h   b   b   1   -----------------
# c24 :  h   b   b   x   1   -------------
# c12 :  h   b   x   x   1   -------------
# c6c :  h   b   b   x   1   -------------
# 
# b : bruteforced  ( 1:floatequal, 2:'testlines', ... )
# h : hardcoded
# x : not yet implemented
# g : generated

use strict;
use Math::Combinatorics;

our @subnames=qw(point line face cell);
sub subname { 
    return $subnames[$_[0]] if $_[0]<@subnames;
    return "obj$_[0]";
}
sub floatequal {
    return abs($_[0]-$_[1]) < 0.00001;
}
sub signbit {
    return $_[0] ? -1 : 1;
}

our @even4perms=map { [ split //, $_] } qw(0123 0231 0312 1032 1203 1320 2013 2130 2301 3021 3102 3210);
sub even4permute {
    my $i= shift;
    return map { $_[$even4perms[$i][$_]] } (0..3)
}
our @perm6=map { [ split //, $_] } qw(0123 0213 0231 2013 2031 2301);
sub permute6 {
    my $i= shift;
    return map { $_[$perm6[$i][$_]] } (0..3)
}
sub allsigns {
    my $i= shift;
    return ( signbit($i&1)*$_[0], signbit($i&2)*$_[1], signbit($i&4)*$_[2], signbit($i&8)*$_[3] );
}
# generate m out of n
# use one of these modules:
#   Algorithm-FastPermute
#   Algorithm-Permute
#   Algorithm-Combinatorics
#   List-Permutor
#   List-Permutor-LOL
#   Math-Combinatorics  <---

sub generate_combinations {
    my ($m, $n)= @_;
    return map { my $val= 0; $val |= $_ for @$_;  $val } combine($m, map { 1<<$_ } 0..$n-1);
}
sub binom {
    my $n=shift;
    my $val=1;
    for (my $m=1 ; $m<=$_[0] ; $m++)
    {
        $val *= $n--;
        $val /= $m;
    }
    return $val;
}

# common interface between 'Point' and 'Object':
#    containsObject
#    distance
#    as_string
#
# 0-dim objects contain 1 point object.
package Point;
use strict;

my %di;
sub new {
    my ($class, $name, @coords)= @_;
    return bless {
        dim=>0,
        name=>$name,
        coords=>\@coords,
    }, $class;
}
sub zero {
    my ($class, $name, $spacedim)= @_;
    my @coords;
    push @coords, 0 for 1..$spacedim;
    return $class->new($name, @coords);
}
sub add {
    my ($a, $b)= @_;
    $a->{coords}[$_] += $b->{coords}[$_] for 0..$#{$a->{coords}};
}
sub scalarmul {
    my ($a, $factor)= @_;
    $a->{coords}[$_] /= $factor for 0..$#{$a->{coords}};
}

sub distance {
    my ($a, $b) = @_;
    if (exists $di{$a}{$b}) {
        return $di{$a}{$b};
    }

    my $ip=0;
    $ip += ($a->{coords}[$_] - $b->{coords}[$_])**2 for (0..$#{$a->{coords}});

    return $di{$a}{$b}=sqrt($ip);
}
sub equal {
    my ($a, $b)= @_;
    if ($#{$a->{coords}} != $#{$b->{coords}}) {
        return 0;
    }
    for (0..$#{$a->{coords}}) {
        if ($a->{coords}[$_]!=$b->{coords}[$_]) {
            return 0;
        }
    }
    return 1;
}
sub containsObject {
    my ($self, $p)= @_;
    if (ref $p ne ref $self) {
        return 0;
    }
    return $self->equal($p);
}
sub as_string {
    my ($self)= @_;
    return sprintf("(%s)", join(",", @{$self->{coords}}));
}
package Object;
use strict;

my %eq;
my %ct;
# this represents a <dim> dimensional object
# an object is bounded by <dim-1> dimensional objects.
#
# 0-D objects contain 1 Point object
# n-D objects contain 'Object's with {dim}==n-1

sub new {
    my ($class, $dim, $index, @parts) = @_;
    return bless {
        dim=>$dim,
        index=>$index,
        name=>'{'.join("", map { $_->{name} } @parts).'}',
        parts=>\@parts,
    }, $class;
}

sub fromprevbyidx {
    my ($class, $dim, $index, $prev, @parts)= @_;

    return $class->new($dim, $index, map { $prev->[$_] } @parts);
}
sub fromprevbypt {
    my ($class, $dim, $index, $prev, @points)= @_;

    #printf("fromprevbypt{%d points) : %s ", $#points+1, join(", ", map { $_->as_string } @points));
    my @parts;
    for my $o (@$prev) {
        if (scalar @{$o->{parts}}==scalar grep { $o->containsObject($_) } @points) {
            push @parts, $o;
        }
    }
    #printf(" -> %d parts\n", $#parts+1);
    if (@parts==0) {
        printf("WARN: empty object created\n");
    }
    return $class->new($dim, $index, @parts);
}
# this currently only works for 1 dimensional objects
sub length {
    my ($self)= @_;
    return $self->{parts}[0]->distance($self->{parts}[1]);
}

sub determinant {
    my @x= @{$_[0]{parts}[0]{coords}};
    my @y= @{$_[1]{parts}[0]{coords}};
    my @z= @{$_[2]{parts}[0]{coords}};
    return ($z[0]+$z[1]+$z[2])* ( $x[0]*($y[1]-$y[2]) + $x[1]*($y[2]-$y[0]) + $x[2]*($y[0]-$y[1]) );
}
# determine if point is inside polygon:
#
# polygon: points { x[i], i=0 .. n-1 }
#          line segments { [x[i],x[(i+1)%n]], i=0..n-1 }
# point p
# for i=0..n-1
#     if (polyY[i]<y<=polyY[i-1]  ||  polyY[i-1]<y<=polyY[i]) {
#          if ((y-polyY[i])*(polyX[i-1]-polyX[i])<(x-polyX[i])*(polyY[i-1]-polyY[i])) {
#              oddNodes=!oddNodes; 
#          }
#     }
#  return oddNodes;

#  http://forums.wolfram.com/mathgroup/archive/2000/Sep/msg00387.html
#
sub isinside {
}
# this only works for 3 d objects
sub volume {

# het volume van een veelvlak begrensd door driehoeken :

# shape = set of triangles ( x[i], y[i], z[i],  i=1 .. n )
# where x,y,z are vectors (the order is important:  counter clock wise, as
# seen from the outside of the shape)

# -> volume = 1/6 * sum (  (z0+z1+z2) * ( x0*(y1-y2) + x1*(y2-y0) + x2*(y0-y1) ) for all (x,y,z) )

# ook te schrijven als 1/6 * som (  det (x[i],y[i],z[i]) )

# het oppervlak van een polygon =  1/2 * som ( det(x[i],y[i]) )


    my ($self)= @_;

    my $center;
    my $n=0;
    for my $point ($self->listof2(0)) {
        if (!defined $center) {
            $center= Point->zero('x', scalar @{$point->{parts}[0]{coords}});
        }
        $center->add($point->{parts}[0]);
        $n++;
    }
    $center->scalarmul(1/$n) if ($n);

    my $volume=0;
    for my $triangle ($self->listof2(2)) {
        # calc determinant
        # todo: figure out if all triangles have the same direction.
        # .... question: what is 'direction' in an n-D space?
        $volume += determinant($triangle->listof2(0))
    }
    return $volume/6;
}

sub surface {
    my ($self)= @_;
    # calcs circumfence for 2 d objects
    # calcs surface area for 3 d objects
    # calcs volume for 4 d objects

    # todo
}
sub distance {
    my ($a, $b)= @_;
    return $a->{parts}[0]->distance($b->{parts}[0]);
}
# assuming each object is allocated only once. not really comparing the parts.
# TODO: do a 'real' compare
sub equals {
    my ($a, $b)= @_;
    if (exists $eq{$a}{$b}) {
        return $eq{$a}{$b};
    }
    my %x;
    $x{$_}++ for @{$a->{parts}};
    $x{$_}-- for @{$b->{parts}};
    return $eq{$a}{$b} = !grep { $_ } values %x;
}
sub containsObject {
    my ($self, $p)= @_;
    if (exists $ct{$self}{$p}) {
        return $ct{$self}{$p};
    }
    if ($self->{dim} < $p->{dim}) {
        return $ct{$self}{$p}=0;
    }
    if ($self->{dim} == $p->{dim}) {
        return $ct{$self}{$p}=$self->equals($p);
    }
    for my $o (@{$self->{parts}}) {
        if ($o->containsObject($p)) {
            return $ct{$self}{$p}=1;
        }
    }
    return $ct{$self}{$p}=0;
}
sub as_string {
    my ($self)= @_;
    return sprintf("(%s)", join(",", map { $_->as_string() } @{$self->{parts}}));
}
sub findobjbyparts {
    my $list= shift;
    for (my $i=0 ; $i<@$list ; $i++) {

        # the set of items not containing a part should be empty.
        # -> true when all parts are contained.
        if (!grep { !$list->[$i]->containsObject($_) } @_) {
            #printf("line %d\n", $i);
            return $i;
        }
    }
}
sub containedIn {
    my $self= shift;
    # for every part there is a '$_[$i]'
    for (my $i=0 ; $i<@{$self->{parts}} ; $i++)
    {
        if (!grep { $self->{parts}[$i]->containsObject($_) } @_)
        {
            return undef;
        }
    }
    return 1;
}
sub listof {
    my ($self, $dim)= @_;
    if ($self->{dim}==$dim) {
        return $self->{index};
    }
    else {
        my %x;
        $x{$_}++ for map { $_->listof($dim) } @{$self->{parts}};
        return sort { $a<=>$b } keys %x;
    }
}
sub listof2 {
    my ($self, $dim)= @_;
    if ($self->{dim}==$dim) {
        return $self;
    }
    else {
        my %x;
        $x{$_->{index}}=$_ for map { $_->listof2($dim) } @{$self->{parts}};
        return values %x;
    }
}
sub testlines {
    my ($edgelen, $i,$j, $dim, $index, $prevpart, $newpart)= @_;
    #printf("test %d %d\n", $i, $j);
    for (my $zj=0 ; $zj<@{$prevpart->[$j]{parts}} ; $zj++) {
        for (my $zi=0 ; $zi<@{$prevpart->[$i]{parts}} ; $zi++) {
            if (::floatequal($prevpart->[$i]{parts}[$zi]->distance($prevpart->[$j]{parts}[$zj]), 0)
               && ::floatequal($prevpart->[$i]{parts}[1-$zi]->distance($prevpart->[$j]{parts}[1-$zj]), $edgelen)) {
                my $k= Object::findobjbyparts($prevpart, $prevpart->[$i]{parts}[1-$zi], $prevpart->[$j]{parts}[1-$zj]);
                if ($i<$k) {
                    # NOTE: $_[4] = the index by reference
                    push @$newpart, Object->fromprevbyidx($dim, $_[4]++, $prevpart, $i, $j, $k);
                    printf("%2d: %s(%d,%d,%d)  [%s]\n", $#$newpart, ::subname($dim), $i, $j, $k, join(" ", $newpart->[-1]->listof(0)));
                }

                #printf("%2d:%d=%s  %2d:%d=%s   %f   %f\n\n", $i, 1-$zi, $prevpart->[$i]->as_string(), $j, 1-$zj, $prevpart->[$j]->as_string(), $prevpart->[$i]{parts}[$zi]->distance($prevpart->[$j]{parts}[$zj]), $prevpart->[$i]{parts}[1-$zi]->distance($prevpart->[$j]{parts}[1-$zj]));
            }
        }
    }
}

##########################################################################
#
# shapes:
#
# N(n,m) = the number of m-dimensional parts of this n-dimensional shape
#
#  dim | N(n,n-4)| N(n,n-3)| N(n,n-2)| N(n,n-1) | (n-1)-shape | common name
#   3  |        -|       20|       30|       12 | pentagon    | dodecaeder
#   3  |        -|       12|       30|       20 | triangle    | icosaeder
#   4  |       24|       96|       96|       24 | octaeder    | 24-cel
#   4  |      600|     1200|      720|      120 | dodecaeder  | 120-cel
#   4  |      120|      720|     1200|      600 | tetraeder   | 600-cel

#   n  |         |         |         |          |  S[n-1]     | n-simplex
#   n  |         |         |         |          |  C[n-1]     | n-cube
#   n  |         |         |         |          |  S[n-1]     | n-orthoplex


#   3  |        -|        4|        6|         4|  triangle   | tetraeder
#   3  |        -|        8|       12|         6|  square     | cube
#   3  |        -|        6|       12|         8|  triangle   | octaeder

#   4  |        5|       10|       10|         5|  tetraeder  | simplex
#   4  |       16|       32|       24|         8|  cube       | hypercube
#   4  |        8|       24|       32|        16|  tetraeder  | orthoplex

# 3d: faces+vertices-edges == 2
# 4d : vertices-edges+faces-cells == 0
# n  : sum((-1)^i * N[i], i=0..n-1) == 1 - (-1)^n
#    N[i] = nr of i-dimensional shapes.
#############################################
#
#
# truncated forms:
#
#
# tetraeder -> octaeder
# octaeder -> X
# cube -> X
# icosaeder -> dodecaeder
# dodecaeder -> Y


package Tetraeder;
use strict;
# 3d
# punten | ribben | vlakken | vorm     | rib/hoek | vlak/rib | vlak/hoek
#      4 |      6 |       4 | driehoek |        3 |        2 |         3

# number of 'a' per 'b' ( dim=3)
#
#    b\0  1  2  3  <-a
#    0|-  3  3  1
#    1|2  -  2  1
#    2|3  3  -  1
#    3|4  6  4  -
#

# number of 'a' per 'b' (dim=4)
#
#    b\0  1  2  3  4 <-a
#    0|-  4  6  4  1
#    1|2  -  3  3  1
#    2|3  3  -  2  1
#    3|4  6  4  -  1
#    4|5 10 10  5  -

# a n-dimensional tetraeder
#
# N(n,m) = binom(n+1, m+1)

# mn0  1  2  3  4  5
# 0 1  2  3  4  5  6    n+1
# 1 0  1  3  6 10 15    (n+1)*n/2
# 2 0  0  1  4 10 20    (n+1)*n*(n-1)/6
# 3 0  0  0  1  5 15
# 4 0  0  0  0  1  6
# 5 0  0  0  0  0  1

# dim  name
#  0    point
#  1    line
#  2    triangle
#  3    tetraeder or simplex
#  4    pentachoron or 5-cell or 4-simplex

# volume: 3d-tetraeder: sqrt(2)/12 * a^3, a=rib-length
#  -> in ourcase, the volume would be 1/3

# volume: 4d-tetraeder: sqrt(5)/96 * a^4, a=rib-length

use constant EDGELEN=>sqrt(2);

our $name= "tetraeder";
sub new {
    my ($class, $dim)= @_;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub makePoint {
    my ($self, $i)= @_;
    my @coords;
    if ($i==$self->{dim}) {
        # N-dim tetraeder:
        # (N-1) points at axis
        #   ->  distance((1,0,...) , (0,1,...)) = sqrt(2)
        #
        # Nth point such that
        #       distance((a,a,...), (1,0,...)) = sqrt(2)
        #
        # -> (a-1)^2+(N-1)*a^2 = N*a^2-2*a+1 == 2
        # -> N*a^2-2*a-1 == 0
        #
        #   ( -b +- sqrt(b^2-4*a*c) ) / (2*a)
        #
        #  ( 2 +- sqrt(4+4*N) ) / (2*N) = (1 +- sqrt(N+1))/N
        #
        # dim :  possible values of 'a'
        # 1 : 2.414214  -0.414214
        # 2 : 1.366025  -0.366025
        # 3 : 1.000000  -0.333333
        # 4 : 0.809017  -0.309017
        # 5 : 0.689898  -0.289898
        push @coords, (1 + sqrt($self->{dim}+1))/$self->{dim} for (1..$self->{dim});
    }
    elsif ($i<$self->{dim}) {
        push @coords, (($i == $_)?1:0) for (0..$self->{dim}-1);
    }
    else {
        die "$name: invalid point $i\n";
    }
    return new Point(chr($i+48+($i>9?7:0)), @coords);
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..$self->{dim}) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # the lines in a tetraeder are all possible combinations of 2 points
        #
        # 0 1 2 3
        # * *      0
        # *   *    1
        # *     *  2
        #   *   *  3
        #     * *  4
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
            }
        }
    }
    elsif ($dim==2) {
        # done: find better way than bruteforcing tetraeder faces
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                Object::testlines(EDGELEN, $i,$j, $dim, $index, $prevpart, $newpart);
            }
        }
    }
    elsif ($dim<$self->{dim}-1) {
        # n-tetraeder faces .. N-2     ( for n >= 5 )

        #done: check if tetraeders dim==1, dim==2 and dim==n-1
        #      can also be created with this method -> YES

        #                     0  1  2  3  4
        # points : ($dim+1, 1)  1  2  3  4  5
        # lines  : ($dim+1, 2)  -  1  3  6 10
        # faces  : ($dim+1, 3)  -  -  1  4 10
        # cells  : ($dim+1, 4)  -  -  -  1  5
        #
        for my $ptsel (::generate_combinations($dim+1, $self->{dim}+1)) {
            # ptsel is a bitmask containing the points of this obj.
            my @pts= map { $self->{parts}[0][$_] } grep { $ptsel&(1<<$_) } (0..31);
            my @new;
            for (my $i=0 ; $i<@$prevpart ; $i++)
            {
                if ($prevpart->[$i]->containedIn(@pts)) {
                    push @new, $prevpart->[$i];
                }
            }
            push @$newpart, new Object($dim, $index++, @new);
        }
    }
    elsif ($dim==$self->{dim}-1) {
        # the n-1 objects consist of all n-2 objects not containing one specific point
        #
        # 0 1 2 3 4 5
        # * * *        0
        # *     * *    1
        #   *     * *  2
        #     * *   *  3

        for my $p (@{$self->{parts}[0]}) {
            push @$newpart, new Object($dim, $index++, grep { !$_->containsObject($p) } @$prevpart);
        }
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    if ($dim==3) {
        printf("tet volume[%s]= %f\n", join(",",$newpart->[-1]->listof(0)), $newpart->[-1]->volume());
    }
    printf("done: found %d objects\n", scalar @$newpart);
}

package Cube;
use strict;
# 3d
# punten | ribben | vlakken | vorm     | rib/hoek | vlak/rib | vlak/hoek
#      8 |     12 |       6 | vierkant |        3 |        2 |         3

# number of 'a' per 'b'
#
#    b\0  1  2  3  <-a
#    0|-  3  3  1
#    1|2  -  2  1
#    2|4  4  -  1
#    3|8 12  6  -

# number of 'a' per 'b'  for dim=4 ( hypercube )
#
#    b\ 0  1  2  3  4 <-a
#    0| -  4  6  4  1
#    1| 2  -  3  3  1
#    2| 4  4  -  2  1
#    3| 8 12  6  -  1
#    4|16 32 24  8  -


# N(n,m) = 2^(n-m) * binom(n,m)
#
# mn0  1  2  3  4  5
# 0 1  2  4  8 16 32     2^n
# 1 0  1  4 12 32 80     2^(n-1)*n
# 2 0  0  1  6 24 80     2^(n-3)*n*(n-1)
# 3 0  0  0  1  8 40     2^(n-4)*n*(n-1)*(n-2)/3
# 4 0  0  0  0  1 10     2^(n-4)*n!/4!/(n-4)!
# 5 0  0  0  0  0  1

# dim  name
#  0    point
#  1    line
#  2    square
#  3    cube
#  4    hypercube / tesseract
# ...
#  n    hypercube
#
# 3-d cube: volume = a^3
# 4-d cube: volume = a^4
#
use constant EDGELEN=>1;
our $name= "cube";
sub new {
    my ($class, $dim)= @_;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub makePoint {
    my ($self, $i)= @_;
    my @coords;
    for (my $bit=1 ; $bit < (1<<$self->{dim}) ; $bit<<=1)
    {
        push @coords, ($i&$bit)?1:0;
    }
    return new Point(chr($i+48+($i>9?7:0)), @coords);
}
sub countones {
    my ($a)= @_;
    my $count=0;
    while ($a) {
        $count++;
        $a &= $a-1;
    }
    return $count;
}
sub bitdistance {
    my ($a, $b)= @_;

    return countones($a^$b);
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1] if ($dim);
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..(1<<$self->{dim})-1) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
#   elsif (0 && $dim==1) {
#       #done: find less bruteforce method of finding cube lines

#       # the lines are all points which have bitdistance == 1
#       #
#       # 0 1 2 3 4 5 6 7
#       # * *               0 : 01
#       # *   *             1 : 02
#       # *       *         2 : 04
#       #   *   *           3 : 13
#       #   *       *       4 : 15
#       #     * *           5 : 23
#       #     *       *     6 : 26
#       #       *       *   7 : 37
#       #         * *       8 : 45
#       #         *   *     9 : 46
#       #           *   *   a : 57
#       #             * *   b : 67
#       for (my $i=1 ; $i<@$prevpart ; $i++)
#       {
#           for (my $j=0 ; $j<$i ; $j++)
#           {
#               if (bitdistance($i,$j)==1) {
#                   push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
#                   printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
#               }
#           }
#       }
#   }
    elsif ($dim<$self->{dim}) {
        #done: check if cube dim==1
        #      can also be created with this method -> YES

        # the m dim subobject consists of all combinations of m out of n bits
        # -> there are 2^(m-n) * binom(n,m)   m dim objects


        my $points= $self->{parts}[0];
        for my $bitsel (::generate_combinations($dim, $self->{dim})) {
           for my $xx (0..(1<<($self->{dim}-$dim))-1) {
               my $fixed=setbits($bitsel^((1<<$self->{dim})-1), $xx);
               #printf("bitsel=%03b  fixed=%03b  ~sel=%03b ", $bitsel, $fixed, $bitsel^((1<<$self->{dim})-1));
               #printf("  using points: [%s]\n", join(", ", map { $fixed|setbits($bitsel, $_) } (0..(1<<$dim)-1) ));

               # TODO: fix bug, when dim>=5  this returns empty objs.
               push @$newpart, Object->fromprevbypt($dim, $index++, $prevpart, map { $points->[$fixed|setbits($bitsel, $_)] } (0..(1<<$dim)-1));
               #printf("%2d: %s\n", $#$newpart, $newpart->[-1]->as_string());
           }
        }
        for my $p (@$newpart) {
            printf("%2d: %s(%s) [%s]\n", $p->{index}, ::subname($dim), join(",",$p->listof($dim-1)), join(",",$p->listof(0)));
        }
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    printf("done: found %d objects\n", scalar @$newpart);
}

sub setbits {
    my ($mask, $n)= @_;
    my $value=0;
    my $bit=1;
    while ($n) {
        if ($mask&$bit) {
            $value |= $bit if ($n&1);
            $n>>=1;
        }
        $bit<<=1;
    }
    return $value;
}

package Octaeder;
use strict;
use Set::Scalar;
# 3d
# punten | ribben | vlakken | vorm     | rib/hoek | vlak/rib | vlak/hoek
#      6 |     12 |       8 | driehoek |        4 |        2 |         4

# number of 'a' per 'b'  for dim=3
#
#    b\0  1  2  3  <-a
#    0|-  4  4  1
#    1|2  -  2  1
#    2|3  3  -  1
#    3|6 12  8  -

# number of 'a' per 'b'  for dim=4 ( 16-cell )
#
#    b\ 0  1  2  3  4 <-a
#    0| -  6 12  8  1
#    1| 2  -  4  4  1
#    2| 3  3  -  2  1      ... face shape = triangle
#    3| 4  6  4  -  1      ... cell shape = tetraeder
#    4| 8 24 32 16  -


# N(n,m) = 2^(m+1) * binom(n, m+1)
#
# mn0  1  2  3  4  5
# 0 1  2  4  6  8 10   2^(n+1)*n
# 1 0  1  4 12 24 40   2^(n)*n*(n-1)/2
# 2 0  0  1  8 32 80   2^(n)*n*(n-1)*(n-2)/3
# 3 0  0  0  1 16 80   2^(n-2)*n*(n-1)*(n-2)*(n-3)/3
# 4 0  0  0  0  1 32   2^(n+1)*n!/(n-5)!/120
# 5 0  0  0  0  0  1

# dim name
#  2   square
#  3   octaeder
#  4   16-cell or orthoplex or hexadecachoron
# ...
#  n   'cross-polytopes'

# possibly my mistake here is that the 16-cells are tetrahedra.

# volume: 3d-octaeder: sqrt(2)/3 * a^3, a=rib-length
#    -> in our case: 4/3
# volume: 4d-octaeder: 1/6 * a^4, a=rib-length

use constant EDGELEN=>sqrt(2);
our $name= "octaeder";
sub new {
    my ($class, $dim)= @_;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub makePoint {
    my ($self, $i)= @_;
    my @coords;
    for (0 .. $self->{dim}-1)
    {
        push @coords, 0;
    }
    $coords[$i/2]= ::signbit($i&1);
    return new Point(chr($i+48+($i>9?7:0)), @coords);
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..2*$self->{dim}-1) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # 0 1 2 3 4 5
        # *   *        0
        # *     *      1
        # *       *    2
        # *         *  3
        #   * *        4
        #   *   *      5
        #   *     *    6
        #   *       *  7
        #     *   *    8
        #     *     *  9
        #       * *    a
        #       *   *  b
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                if (($i^$j)!=1) {
                    push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                    printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
                }
            }
        }
    }
    elsif ($dim==2) {
        # todo: find better way than bruteforcing octaeder faces
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                Object::testlines(EDGELEN, $i,$j, $dim, $index, $prevpart, $newpart);
            }
        }
    }
    elsif ($dim==3 && $dim<$self->{dim}) {
        # n-octaeder cells .. N-1  ( for octaeder n >= 5 )
        my @faces;
        push @faces, new Set::Scalar( $_->listof(0) ) for @$prevpart;

        my @f;
        for my $i (0..$#faces) {
            push @{$f[$i]}, new Set::Scalar()  for 0..3;
            for my $j (0..$#faces) {
                my $ni= scalar($faces[$i]->intersection($faces[$j])->members);
                $f[$i][$ni]->insert($j);
            }
        }
        #todo: find better way than bruteforcing octaeder cells
        # bruteforcing it
        for my $i (0..$#faces) {
            for my $j ($f[$i][2]->elements) {
                next if $i>=$j;
                for my $k ($f[$i][2]->elements) {
                    next if $j>=$k;
                    for my $l ($f[$i][2]->elements) {
                        next if $k>=$l;
                        if (is_cell($faces[$i], $faces[$j], $faces[$k], $faces[$l])) {
                            my $pts= new Set::Scalar();
                            $pts->insert($faces[$_]->elements) for ($i, $j, $k, $l);
                            printf("cell: faces(%s)  points%s : %s\n", join(',',$i, $j, $k, $l), $pts, join("  ", map { $self->{parts}[0][$_]->as_string() } $pts->members));
                            push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j, $k, $l);
                        }
                    }
                }
            }
        }
    }
    elsif ($dim<$self->{dim}) {
        # todo octaeder dim=4 .. n-1
        printf("TODO: %s dim=%d - tetraeders, nprev=%d\n", $name, $dim, scalar @$prevpart);
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    if ($dim==3) {
        printf("oct volume[%s]= %f\n", join(",",$newpart->[-1]->listof(0)), $newpart->[-1]->volume());
    }
    printf("done: found %d objects\n", scalar @$newpart);
}
sub is_cell {
    my $c= new Set::Scalar;
    $c->insert($_->elements) for @_;
    #printf("cellsize=%d: %s\n", $c->size(), $c);
    return $c->size()==4;
}
package Dodecaeder;
use strict;
# 3d
# punten | ribben | vlakken | vorm     | rib/hoek | vlak/rib | vlak/hoek
#     20 |     30 |      12 | vijfhoek |        3 |        2 |         3

# number of 'a' per 'b'
#
#    b\ 0  1  2  3 <-a
#    0| -  3  3  1
#    1| 2  -  2  1
#    2| 5  5  -  1
#    3|20 30 12  -

# riblength: 1.236068
# volume: 1/4 * (15+7*sqrt(5))*a^3
#  -> in our case: 14.4721

our $name= "dodecaeder";
use constant PHI => (1+sqrt(5))/2;
use constant EDGELEN => 2/PHI;
sub new {
    my ($class, $dim)= @_;
    die "$name exists only in 3d" unless $dim==3;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub pt111 {
    my ($i)= @_;

    return ( ::signbit($i&4), ::signbit($i&2), ::signbit($i&1) );
}
sub pt011 {
    my ($i)= @_;
    return ( 0, ::signbit($i&1)/PHI, ::signbit($i&2)*PHI );
}
sub pt110 {
    my ($i)= @_;
    return ( ::signbit($i&1)/PHI, ::signbit($i&2)*PHI, 0 );
}
sub pt101 {
    my ($i)= @_;
    return ( ::signbit($i&2)*PHI, 0, ::signbit($i&1)/PHI );
}
sub makePoint {
    my ($self, $i)= @_;
    my $chr= chr($i+48+($i>9?7:0));

    #  (+-1, +-1, +-1)  -> these form the corners of a cube
    #  (0, +-1/F, +-F)
    #  (+-1/F, +-F, 0)
    #  (+-F, 0, +-1/F)
    return new Point($chr, pt111($i)) if ($i<8);
    $i-=8;
    return new Point($chr, pt011($i)) if ($i<4);
    $i-=4;
    return new Point($chr, pt110($i)) if ($i<4);
    $i-=4;
    return new Point($chr, pt101($i)) if ($i<4);
    $i-=4;
    die "$name: invalid point $i\n";
}

sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..19) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # todo: find better way than bruteforcing dodecaeder lines
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                if (::floatequal($prevpart->[$i]->distance($prevpart->[$j]), EDGELEN)) {
                    push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                    printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
                }
            }
        }
    }
    elsif ($dim==2) {
        # todo: find better way than handcoding dodecaeder faces
        my $points= $self->{parts}[0];
        push @$newpart, Object->fromprevbypt($dim, 0, $prevpart, map { $points->[$_] } 8, 9, 2, 16, 0);
        push @$newpart, Object->fromprevbypt($dim, 1, $prevpart, map { $points->[$_] } 8, 9, 6, 18, 4);
        push @$newpart, Object->fromprevbypt($dim, 2, $prevpart, map { $points->[$_] }10,11, 3, 17, 1);
        push @$newpart, Object->fromprevbypt($dim, 3, $prevpart, map { $points->[$_] }10,11, 7, 19, 5);
        push @$newpart, Object->fromprevbypt($dim, 4, $prevpart, map { $points->[$_] }12,13, 4,  8, 0);
        push @$newpart, Object->fromprevbypt($dim, 5, $prevpart, map { $points->[$_] }12,13, 5, 10, 1);
        push @$newpart, Object->fromprevbypt($dim, 6, $prevpart, map { $points->[$_] }14,15, 6,  9, 2);
        push @$newpart, Object->fromprevbypt($dim, 7, $prevpart, map { $points->[$_] }14,15, 7, 11, 3);
        push @$newpart, Object->fromprevbypt($dim, 8, $prevpart, map { $points->[$_] }16,17, 1, 12, 0);
        push @$newpart, Object->fromprevbypt($dim, 9, $prevpart, map { $points->[$_] }16,17, 3, 14, 2);
        push @$newpart, Object->fromprevbypt($dim,10, $prevpart, map { $points->[$_] }18,19, 5, 13, 4);
        push @$newpart, Object->fromprevbypt($dim,11, $prevpart, map { $points->[$_] }18,19, 7, 15, 6);

        for my $p (@$newpart) {
            printf("%2d: %s(%s) [%s]\n", $p->{index}, ::subname($dim), join(",",$p->listof(1)), join(",",$p->listof(0)));
        }
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    printf("done: found %d objects\n", scalar @$newpart);
}

package Icosaeder;
use strict;
# 3d
# punten | ribben | vlakken | vorm     | rib/hoek | vlak/rib | vlak/hoek
#     12 |     30 |      20 | driehoek |        5 |        2 |         5

# number of 'a' per 'b'
#
#    b\ 0  1  2  3  <-a
#    0| -  5  5  1
#    1| 2  -  2  1
#    2| 3  3  -  1
#    3|12 30 20  -

# riblength = 2
# volume: 5/12*(3+sqrt(5))*a^3
#  -> in our case: 17.45

our $name= "icosaeder";

use constant PHI => (1+sqrt(5))/2;
use constant EDGELEN => 2;
sub new {
    my ($class, $dim)= @_;
    die "$name exists only in 3d" unless $dim==3;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub pt011 {
    my ($i)= @_;
    return ( 0,            ::signbit($i&1), ::signbit($i&2)*PHI );
}
sub pt110 {
    my ($i)= @_;
    return ( ::signbit($i&1), ::signbit($i&2)*PHI,           0 );
}
sub pt101 {
    my ($i)= @_;
    return ( ::signbit($i&1)*PHI,             0, ::signbit($i&2) );
}

sub makePoint {
    my ($self, $i)= @_;
    my $chr= chr($i+48+($i>9?7:0));
    #  (0, +-1, +-F)
    #  (+-1, +-F, 0)
    #  (+-F, 0, +-1)
    return new Point($chr, pt011($i)) if ($i<4);
    $i-=4;
    return new Point($chr, pt110($i)) if ($i<4);
    $i-=4;
    return new Point($chr, pt101($i)) if ($i<4);
    $i-=4;
    die "$name: invalid point $i\n";
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..11) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # todo: find better way than bruteforcing icosaeder lines
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                if (::floatequal($prevpart->[$i]->distance($prevpart->[$j]), EDGELEN)) {
                    push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                    printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
                }
            }
        }
    }
    elsif ($dim==2) {
        # todo: find better way than bruteforcing icosaeder faces
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                Object::testlines(EDGELEN, $i,$j, $dim, $index, $prevpart, $newpart);
            }
        }
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    if ($dim==3) {
        printf("ico volume[%s]= %f\n", join(",",$newpart->[-1]->listof(0)), $newpart->[-1]->volume());
    }
    printf("done: found %d objects\n", scalar @$newpart);
}



package Cell24;
use strict;
# punten | ribben | vlakken | lichamen | vorm       | rib/hoek | vlak/rib | vlakhoek | lichaam/vlak | lichaam/rib | lichaam/hoek
#     24 |     96 |      96 |       24 | octaeder   |        8 |       ?3 |      12? |            2 |           3 |            6

# number of 'a' per 'b'
#
#    b\ 0  1  2  3  4<-a
#    0| -  8 12? 6  1
#    1| 2  -  3? 3  1
#    2| 3  3  -  2  1
#    3| 6 12  8  -  1
#    4|24 96 96 24  -
#

# 24-cell or icositetrachoron or octaplex

# volume : 2*a^4

our $name= "24-cell";
use constant EDGELEN=>1;

sub new {
    my ($class, $dim)= @_;
    die "$name exists only in 4d" unless $dim==4;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub pt_unit {
    my ($i)= @_;

    my @coords;
    for (0 .. 3)
    {
        push @coords, int($i/2)==$_ ? ::signbit($i&1) : 0;
    }

    return @coords;
}
sub pt_half {
    my ($i)= @_;

    return ( ::signbit($i&8)/2, ::signbit($i&4)/2, ::signbit($i&2)/2, ::signbit($i&1)/2 );
}
sub makePoint {
    my ($self, $i)= @_;
    my $chr= chr($i+48+($i>9?7:0));
    return new Point($chr, pt_unit($i)) if ($i<8);
    $i-=8;
    return new Point($chr, pt_half($i)) if ($i<16);
    $i-=16;
    die "$name: invalid point $i\n";
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..23) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            #printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # todo: find better way than bruteforcing 24-cell lines
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                if (::floatequal($prevpart->[$i]->distance($prevpart->[$j]), EDGELEN)) {
                    push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                    #printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
                }
            }
        }
    }
    elsif ($dim==2) {
        # todo: find better way than bruteforcing 24-cell faces
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                Object::testlines(EDGELEN, $i,$j, $dim, $index, $prevpart, $newpart);
            }
        }
    }
    elsif ($dim<$self->{dim}) {
        # todo  24-cell cells : octaeder
        printf("TODO: %s dim=%d - octaeders, nprev=%d\n", $name, $dim, scalar @$prevpart);
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    if ($dim==3) {
        printf("c24 volume[%s]= %f\n", join(",",$newpart->[-1]->listof(0)), $newpart->[-1]->volume());
    }
    printf("done: found %d objects\n", scalar @$newpart);
}

package Cell120;
use strict;
# punten | ribben | vlakken | lichamen | vorm       | rib/hoek | vlak/rib | vlakhoek | lichaam/vlak | lichaam/rib | lichaam/hoek
#    600 |   1200 |     720 |      120 | dodecaeder |        4 |        3 |        6 |            2 |           3 |            4

#  1200*3==5*720  -> each rib is in 3 faces
#  720*2==12*120  -> each face in in 2 cells

# number of 'a' per 'b'
#
#    b\   0    1    2    3    4 <-a
#    0|   -    4    6    4    1
#    1|   2    -    3    3    1
#    2|   5    5    -    2    1
#    3|  20   30   12    -    1
#    4| 600 1200  720  120    -

# 120-cell or hecatonicosachoron or dodecaplex
#  ... like 4d version of dodecaeder

# volume: 15/4*(47*sqrt(5)+105) * a^4

our $name= "120-cell";
use constant PHI => (1+sqrt(5))/2;
use constant EDGELEN => 2/PHI/PHI;
sub new {
    my ($class, $dim)= @_;
    die "$name exists only in 4d" unless $dim==4;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub pt0022 {
    my ($i)= @_;
    return (0, 0, ::signbit($i&1)*2, ::signbit($i&2)*2);
}
sub pt_sqrt5 {
    my ($i)= @_;
    my @coords= (1,1,1,1);
    $coords[$i]= sqrt(5);
    return @coords;
}
sub pt_F_2 {
    my ($i)= @_;
    my @coords= (PHI, PHI, PHI, PHI);
    $coords[$i]= 1/PHI/PHI;
    return @coords;
}
sub pt_F_F {
    my ($i)= @_;
    my @coords= (1/PHI, 1/PHI, 1/PHI, 1/PHI);
    $coords[$i]= PHI*PHI;
    return @coords;
}
sub pt_F0 {
    my ($i)= @_;
    return ( 0, ::signbit($i&1)/PHI/PHI, ::signbit($i&2), ::signbit($i&4)*PHI*PHI );
}
sub pt_F1 {
    my ($i)= @_;
    return ( 0, ::signbit($i&1)/PHI, ::signbit($i&2)*PHI, ::signbit($i&4)*sqrt(5) );
}
sub pt_F2 {
    my ($i)= @_;
    return ( ::signbit($i&8)/PHI, ::signbit($i&1), ::signbit($i&2)*PHI, ::signbit($i&4)*2 );
}
sub makePoint {
    my ($self, $i)= @_;
    my $chr= chr($i+48+($i>9?7:0));

    # (0, 0, ±2, ±2)
    return new Point($chr, ::permute6($i>>2, pt0022($i&3))) if ($i<4*6);
    $i-=4*6;
    # (±1, ±1, ±1, ±sqrt(5))
    return new Point($chr, ::allsigns($i>>2, pt_sqrt5($i&3))) if ($i<64);
    $i-=64;
    # (±f^-2, ±f, ±f, ±f)
    return new Point($chr, ::allsigns($i>>2, pt_F_2($i&3))) if ($i<64);
    $i-=64;
    # (±f^-1, ±f^-1, ±f^-1, ±f^2)
    return new Point($chr, ::allsigns($i>>2, pt_F_F($i&3))) if ($i<64);
    $i-=64;

    # and all even permutations of

    # (0, ±f^-2, ±1, ±f^2)
    return new Point($chr, ::even4permute($i>>3, pt_F0($i&7))) if ($i<8*12);
    $i-=8*12;
    # (0, ±f^-1, ±f, ±sqrt(5))
    return new Point($chr, ::even4permute($i>>3, pt_F1($i&7))) if ($i<8*12);
    $i-=8*12;
    # (±f^-1, ±1, ±f, ±2) 
    return new Point($chr, ::even4permute($i>>4, pt_F2($i&15))) if ($i<16*12);
    $i-=16*12;

    die "$name: invalid point $i\n";
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..599) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            #printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # todo: find better way than bruteforcing 120-cell lines
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                if (::floatequal($prevpart->[$i]->distance($prevpart->[$j]), EDGELEN)) {
                    push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                    #printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
                }
            }
        }
    }
    elsif ($dim==2) {
        # todo 120-cell faces pentagon
        printf("TODO: %s dim=%d - pentagons, nprev=%d\n", $name, $dim, scalar @$prevpart);
    }
    elsif ($dim<$self->{dim}) {
        # todo 120-cell cells dodecaeder
        printf("TODO: %s dim=%d - dodecaeders, nprev=%d\n", $name, $dim, scalar @$prevpart);
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    printf("done: found %d objects\n", scalar @$newpart);
}


package Cell600;
use strict;
# punten | ribben | vlakken | lichamen | vorm       | rib/hoek | vlak/rib | vlakhoek | lichaam/vlak | lichaam/rib | lichaam/hoek
#    120 |    720 |    1200 |      600 | tetraeder  |      12? |       ?5 |      30? |            2 |           5 |           20

# number of 'a' per 'b'
#
#    b\   0    1    2    3    4 <-a
#    0|   -   12?  30?  20    1
#    1|   2    -    5?   5    1
#    2|   3    3    -    2    1
#    3|   4    6    4    -    1
#    4| 120  720 1200  600    -


# 600-cell or hexacosichoron or tetraplex
#  ... like 4d version of icosaeder

# volume: 25/4*(sqrt(5)+2)

our $name= "600-cell";

use constant PHI => (1+sqrt(5))/2;
use constant EDGELEN => 1/PHI;
sub new {
    my ($class, $dim)= @_;
    die "$name exists only in 4d" unless $dim==4;
    my $self= bless { dim=>$dim }, $class;

    for my $subdim (0..$dim) {
        $self->createSubdim($subdim);
    }
    return $self;
}
sub pt_other {
    my ($i)= @_;
    return ( ::signbit($i&4)/2, ::signbit($i&2)*PHI/2, ::signbit($i&1)/PHI/2, 0 );
}
sub makePoint {
    my ($self, $i)= @_;
    my $chr= chr($i+48+($i>9?7:0));
    my @coords;
    for (0 .. $self->{dim}-1)
    {
        push @coords, 0;
    }
    # (±½,±½,±½,±½)  if $i<16
    return new Point($chr, Cell24::pt_half($i)) if ($i<16);
    $i-=16;
    # (0,0,0,±1)     if $i<16+8       # all perms
    return new Point($chr, Cell24::pt_unit($i)) if ($i<8);
    $i-=8;
    # ½(±1,±f,±1/f,0)  if $i<16+8+96  # even perms
    return new Point($chr, ::even4permute($i>>3, pt_other($i&7)));
    $i-=12*8;

    #   f = (1+v5)/2
    #   edge len = 1/f
    die "$name: invalid point $i\n";
}
sub createSubdim {
    my ($self, $dim)= @_;
    my $index= 0;
    my $newpart= $self->{parts}[$dim]= [];
    my $prevpart= $self->{parts}[$dim-1];
    printf("creating %dD-$name sub%d\n", $self->{dim}, $dim);
    if ($dim==0) {
        for my $i (0..119) {
            push @$newpart, new Object($dim, $index++, $self->makePoint($i));
            #printf("%2d: Point(%s)\n", $i, $newpart->[-1]->as_string());
        }
    }
    elsif ($dim==1) {
        # todo: find better way than bruteforcing 600 cell lines
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                if (::floatequal($prevpart->[$i]->distance($prevpart->[$j]), EDGELEN)) {
                    push @$newpart, Object->fromprevbyidx($dim, $index++, $prevpart, $i, $j);
                    #printf("%2d: line(%d..%d)  l=%f\n", $#$newpart, $i, $j, $newpart->[-1]->length());
                }
            }
        }
    }
    elsif ($dim==2) {
        # todo: find better way than bruteforcing 600 cell faces
        for (my $i=1 ; $i<@$prevpart ; $i++)
        {
            for (my $j=0 ; $j<$i ; $j++)
            {
                Object::testlines(EDGELEN, $i,$j, $dim, $index, $prevpart, $newpart);
            }
        }
    }
    elsif ($dim<$self->{dim}) {
        # todo 600-cell cells : tetraeder
        printf("TODO: %s dim=%d - tetraeders, nprev=%d\n", $name, $dim, scalar @$prevpart);
    }
    elsif ($dim==$self->{dim}) {
        push @$newpart, new Object($dim, $index++, @$prevpart);
    }
    if ($dim==3) {
        printf("c600 volume[%s]= %f\n", join(",",$newpart->[-1]->listof(0)), $newpart->[-1]->volume());
    }
    printf("done: found %d objects\n", scalar @$newpart);
}

package main;
use strict;
use Dumpvalue;
my $dd= new Dumpvalue;

my $dim= shift || 3;
my $shape= shift;

if ($shape) {
    my $s= $shape->new($dim);
}
else {
    my @shapes= qw(Tetraeder Cube Octaeder);
    push @shapes, qw(Dodecaeder Icosaeder) if ($dim==3);
    push @shapes, qw(Cell24 Cell120 Cell600) if ($dim==4);
    for my $shape (@shapes) {
        $shape->new($dim);
    }
}