github.com/coyove/nj@v0.0.0-20221110084952-c7f8db1065c3/tests/bench/n-body.pl (about)

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