use Math::Polynomial::Solve qw! poly_derivative poly_roots !;

# $rep is the right endpoint of the interval [0, $rep]

foreach my $a (3..$rep ) {
    foreach my $b (2..$a-1 ) {
        foreach my $c (1..$b-1 ) {
            
        # expanded form of x*(x - a)*(x - b)*(x - c)
        # coeffs are in an array
      
        my @quartic = (1, -$a - $b - $c, $a*$b + $a*$c + $b*$c, -$a*$b*$c, 0);        
        my @derivative = poly_derivative(@quartic);
        my @zeros = poly_roots(@derivative);        
      
        $haystack{"$a, $b, $c"} =  \@zeros;
        }
    }
}