slow loops with GVSV indices

nbody-perl-3.perl

# The Computer Language Shootout
# http://shootout.alioth.debian.org/
#
# contributed by Christoph Bauer
# converted into Perl by Márton Papp
# fixed and cleaned up by Danny Sauer
# optimized by Jesse Millikan
# de-optimized by Reini Urban

use constant PI            => 3.141592653589793;
use constant SOLAR_MASS    => (4 * PI * PI);
use constant DAYS_PER_YEAR => 365.24;

#  Globals for arrays... Oh well.
#  Almost every iteration is a range, so I keep the last index rather than a count.
my (@xs, @ys, @zs, @vxs, @vys, @vzs, @mass, $last);

sub energy
{
  my ($e);

  $e = 0.0;
  for my $i (0 .. $last) {
    $e += 0.5 * $mass[$i] *
          ($vxs[$i] * $vxs[$i] + $vys[$i] * $vys[$i] + $vzs[$i] * $vzs[$i]);
    for my $j ($i + 1..$last) {
      my $dx = $xs[$i] - $xs[$j];
      my $dy = $ys[$i] - $ys[$j];
      my $dz = $zs[$i] - $zs[$j];
      my $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
      $e -= ($mass[$i] * $mass[$j]) / $distance;
    }
  }
  return $e;
}

sub offset_momentum
{
  my ($px, $py, $pz) = (0.0, 0.0, 0.0);

  for my $i (0 .. $last) {
    $px += $vxs[$i] * $mass[$i];
    $py += $vys[$i] * $mass[$i];
    $pz += $vzs[$i] * $mass[$i];
  }
  $vxs[0] = - $px / SOLAR_MASS;
  $vys[0] = - $py / SOLAR_MASS;
  $vzs[0] = - $pz / SOLAR_MASS;
}

# @ns = ( sun, jupiter, saturn, uranus, neptune )
@xs = (0, 4.84143144246472090e+00, 8.34336671824457987e+00, 1.28943695621391310e+01, 1.53796971148509165e+01);
@ys = (0, -1.16032004402742839e+00, 4.12479856412430479e+00, -1.51111514016986312e+01, -2.59193146099879641e+01);
@zs = (0, -1.03622044471123109e-01, -4.03523417114321381e-01, -2.23307578892655734e-01, 1.79258772950371181e-01);
@vxs = map {$_ * DAYS_PER_YEAR}
  (0, 1.66007664274403694e-03, -2.76742510726862411e-03, 2.96460137564761618e-03, 2.68067772490389322e-03);
@vys = map {$_ * DAYS_PER_YEAR}
  (0, 7.69901118419740425e-03, 4.99852801234917238e-03, 2.37847173959480950e-03, 1.62824170038242295e-03);
@vzs = map {$_ * DAYS_PER_YEAR}
  (0, -6.90460016972063023e-05, 2.30417297573763929e-05, -2.96589568540237556e-05, -9.51592254519715870e-05);
@mass = map {$_ * SOLAR_MASS}
  (1, 9.54791938424326609e-04, 2.85885980666130812e-04, 4.36624404335156298e-05, 5.15138902046611451e-05);

$last = @xs - 1;

offset_momentum();
printf ("%.9f\n", energy());

my $n = $ARGV[0];
$n =~ s/[,_]//g; # allow 50_000_000 or 50,000,000

# create faster local scalars, instead of aelem
my ($xs_0, $xs_1, $xs_2, $xs_3, $xs_4);
my ($ys_0, $ys_1, $ys_2, $ys_3, $ys_4);
my ($zs_0, $zs_1, $zs_2, $zs_3, $zs_4);
my ($vxs_0, $vxs_1, $vxs_2, $vxs_3, $vxs_4);
my ($vys_0, $vys_1, $vys_2, $vys_3, $vys_4);
my ($vxs_0, $zs_1, $zs_2, $zs_3, $zs_4);
my ($mass_0, $mass_1, $mass_2, $mass_3, $mass_4);

for my $i (0 .. $last) {
  ${"xs_$i"} = $xs[$i];
  ${"ys_$i"} = $ys[$i];
  ${"zs_$i"} = $zs[$i];
  ${"vxs_$i"} = $vxs[$i];
  ${"vys_$i"} = $vys[$i];
  ${"vzs_$i"} = $vzs[$i];
  ${"mass_$i"} = $mass[$i];
}

# This does not, in fact, consume N*4 bytes of memory
for (1 .. $n) {
  my $dt = 0.01;
  for my $i (0 .. $last) {
    #  But not in the inner loop.
    for (my $j = $i + 1; $j <= $last; $j++) {
      my $dx = ${"xs_$i"} - ${"xs_$j"};
      my $dy = ${"ys_$i"} - ${"ys_$j"};
      my $dz = ${"zs_$i"} - ${"zs_$j"};
      my $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
      my $mag = $dt / ($distance * $distance * $distance);
      my $mm  = ${"mass_$i"} * $mag;
      my $mm2 = ${"mass_$j"} * $mag;
      ${"vxs_$i"} -= $dx * $mm2;
      ${"vxs_$j"} += $dx * $mm;
      ${"vys_$i"} -= $dy * $mm2;
      ${"vys_$j"} += $dy * $mm;
      ${"vzs_$i"} -= $dz * $mm2;
      ${"vzs_$j"} += $dz * $mm;
    }

    # We're done with planet $i at this point
    # This could be done in a seperate loop, but it's slower
    ${"xs_$i"} += $dt * ${"vxs_$i"};
    ${"xs_$i"} += $dt * ${"vys_$i"};
    ${"xs_$i"} += $dt * ${"vzs_$i"};
  }
}

for my $i (0 .. $last) {
  $xs[$i] = ${"xs_$i"};
  $ys[$i] = ${"ys_$i"};
  $zs[$i] = ${"zs_$i"};
  $vxs[$i] = ${"vxs_$i"};
  $vys[$i] = ${"vys_$i"};
  $vzs[$i] = ${"vzs_$i"};
  $mass[$i] = ${"mass_$i"};
}

printf ("%.9f\n", energy());