Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

nbody.perl-2a.perl: Beautify loop unroller with qv(). 0.1s slower.

  • Loading branch information...
commit c35bb85ed84941157eb01b7ca844d3b4472e0df3 1 parent 62b2167
Reini Urban authored October 06, 2012

Showing 1 changed file with 151 additions and 0 deletions. Show diff stats Hide diff stats

  1. 151  bench/nbody/nbody.perl-2a.perl
151  bench/nbody/nbody.perl-2a.perl
... ...
@@ -0,0 +1,151 @@
  1
+# The Computer Language Shootout
  2
+# http://shootout.alioth.debian.org/
  3
+#
  4
+# contributed by Christoph Bauer
  5
+# converted into Perl by Márton Papp
  6
+# fixed and cleaned up by Danny Sauer
  7
+# optimized by Jesse Millikan
  8
+# optimized by Reini Urban
  9
+
  10
+use constant PI            => 3.141592653589793;
  11
+use constant SOLAR_MASS    => (4 * PI * PI);
  12
+use constant DAYS_PER_YEAR => 365.24;
  13
+
  14
+sub energy;
  15
+sub advance($);
  16
+sub offset_momentum;
  17
+
  18
+my (@xs, @ys, @zs, @vxs, @vys, @vzs, @mass, $last);
  19
+my ($energy, $offset_momentum, $advance);
  20
+BEGIN {
  21
+# Global lexicals for arrays.
  22
+# Almost every iteration is a range, so I keep the last index rather than a count.
  23
+
  24
+# @ns = ( sun, jupiter, saturn, uranus, neptune )
  25
+@xs = (0, 4.84143144246472090e+00, 8.34336671824457987e+00, 1.28943695621391310e+01, 1.53796971148509165e+01);
  26
+@ys = (0, -1.16032004402742839e+00, 4.12479856412430479e+00, -1.51111514016986312e+01, -2.59193146099879641e+01);
  27
+@zs = (0, -1.03622044471123109e-01, -4.03523417114321381e-01, -2.23307578892655734e-01, 1.79258772950371181e-01);
  28
+@vxs = map {$_ * DAYS_PER_YEAR}
  29
+  (0, 1.66007664274403694e-03, -2.76742510726862411e-03, 2.96460137564761618e-03, 2.68067772490389322e-03);
  30
+@vys = map {$_ * DAYS_PER_YEAR}
  31
+  (0, 7.69901118419740425e-03, 4.99852801234917238e-03, 2.37847173959480950e-03, 1.62824170038242295e-03);
  32
+@vzs = map {$_ * DAYS_PER_YEAR}
  33
+  (0, -6.90460016972063023e-05, 2.30417297573763929e-05, -2.96589568540237556e-05, -9.51592254519715870e-05);
  34
+@mass = map {$_ * SOLAR_MASS}
  35
+  (1, 9.54791938424326609e-04, 2.85885980666130812e-04, 4.36624404335156298e-05, 5.15138902046611451e-05);
  36
+$last = $#xs;
  37
+
  38
+# Optimize array accesses: $a[const] are optimized to AELEMFAST, $a[$lexical] not.
  39
+# So unroll the loops in macro-like fashion (2x times faster). We do it in a BEGIN block,
  40
+# so perlcc can also benefit (again 2x faster).
  41
+sub qv {
  42
+  my $s = shift;
  43
+  my $env = shift;
  44
+  # expand our local loop vars
  45
+  $s =~ s/(\$\w+?)\b/exists($env->{$1})?$env->{$1}:$1/sge;
  46
+  $s
  47
+}
  48
+
  49
+$energy = '
  50
+sub energy
  51
+{
  52
+  my $e = 0.0;
  53
+  my ($dx, $dy, $dz, $distance);';
  54
+  for my $i (0 .. $last) {
  55
+    my $env = {'$i'=>$i,'$last'=>$last};
  56
+    $energy .= qv('
  57
+    # outer-loop $i..4
  58
+    $e += 0.5 * $mass[$i] *
  59
+          ($vxs[$i] * $vxs[$i] + $vys[$i] * $vys[$i] + $vzs[$i] * $vzs[$i]);', $env);
  60
+    for (my $j = $i + 1; $j < $last + 1; $j++) {
  61
+      $env->{'$j'} = $j;
  62
+      $energy .= qv('
  63
+      # inner-loop $j..4
  64
+      $dx = $xs[$i] - $xs[$j];
  65
+      $dy = $ys[$i] - $ys[$j];
  66
+      $dz = $zs[$i] - $zs[$j];
  67
+      $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
  68
+      $e -= ($mass[$i] * $mass[$j]) / $distance;', $env);
  69
+    }
  70
+  }
  71
+  $energy .= '
  72
+  return $e;
  73
+}';
  74
+eval $energy; die if $@;
  75
+
  76
+$advance = '
  77
+sub advance($)
  78
+{
  79
+  my $dt = $_[0];
  80
+  my ($mm, $mm2, $j, $dx, $dy, $dz, $distance, $mag);';
  81
+  for my $i (0..$last) {
  82
+    my $env = {'$i'=>$i};
  83
+    for (my $j = $i + 1; $j < $last + 1; $j++) {
  84
+      $env->{'$j'} = $j;
  85
+      $advance .= qv('
  86
+      # outer-loop $i..4
  87
+      # inner-loop $j..4
  88
+      $dx = $xs[$i] - $xs[$j];
  89
+      $dy = $ys[$i] - $ys[$j];
  90
+      $dz = $zs[$i] - $zs[$j];
  91
+      $distance = sqrt($dx * $dx + $dy * $dy + $dz * $dz);
  92
+      $mag = $dt / ($distance * $distance * $distance);
  93
+      $mm  = $mass[$i] * $mag;
  94
+      $mm2 = $mass[$j] * $mag;
  95
+      $vxs[$i] -= $dx * $mm2;
  96
+      $vxs[$j] += $dx * $mm;
  97
+      $vys[$i] -= $dy * $mm2;
  98
+      $vys[$j] += $dy * $mm;
  99
+      $vzs[$i] -= $dz * $mm2;
  100
+      $vzs[$j] += $dz * $mm;', $env);
  101
+    }
  102
+  }
  103
+  # We're done with planet $i at this point
  104
+  for my $i (0..$last) {
  105
+    my $env = {'$i'=>$i};
  106
+    $advance .= qv('
  107
+    $xs[$i] += $dt * $vxs[$i];
  108
+    $ys[$i] += $dt * $vys[$i];
  109
+    $zs[$i] += $dt * $vzs[$i];', $env);
  110
+  }
  111
+  $advance .= '
  112
+}';
  113
+eval $advance; die if $@;
  114
+
  115
+$offset_momentum = ';
  116
+sub offset_momentum
  117
+{
  118
+  my $px = 0.0; 
  119
+  my $py = 0.0;
  120
+  my $pz = 0.0;
  121
+  my $mass;
  122
+';
  123
+for my $i (0 .. $last) {
  124
+  my $env = {'$i'=>$i};
  125
+  $offset_momentum .= qv('
  126
+    $mass = $mass[$i];
  127
+    $px += $vxs[$i] * $mass;
  128
+    $py += $vys[$i] * $mass;
  129
+    $pz += $vzs[$i] * $mass;', $env);
  130
+}
  131
+$offset_momentum .= '
  132
+  $vxs[0] = - $px / SOLAR_MASS;
  133
+  $vys[0] = - $py / SOLAR_MASS;
  134
+  $vzs[0] = - $pz / SOLAR_MASS;
  135
+}';
  136
+eval $offset_momentum; die if $@;
  137
+
  138
+} #BEGIN
  139
+
  140
+offset_momentum();
  141
+printf ("%.9f\n", energy());
  142
+
  143
+my $n = $ARGV[0];
  144
+$n =~ s/[,_]//g; # allow 50_000_000 or 50,000,000
  145
+
  146
+# This does not, in fact, consume N*4 bytes of memory
  147
+for (1 .. $n) {
  148
+  advance(0.01);
  149
+}
  150
+
  151
+printf ("%.9f\n", energy());

0 notes on commit c35bb85

Please sign in to comment.
Something went wrong with that request. Please try again.