bench/shootout/ackermann.ruby-5.ruby

   1 #!/usr/bin/ruby
   2 # The Computer Language Benchmark Shootout
   3 # http://shootout.alioth.debian.org
   4 #
   5 # original code by Martin DeMello
   6 # modified by Jabari Zakiya 3/20/05
   7 # modified by Glenn Parker 3/28/05 (format results)
   8
   9 include Math
  10
  11 SOLAR_MASS = 4*PI*PI
  12 DAYS_PER_YEAR = 365.24
  13
  14 Vector3D = Struct.new("Vector3D", :x, :y, :z)
  15
  16 class Vector3D
  17
  18   def *(val)
  19     Vector3D.new(*self.map {|i| i * val})
  20   end
  21
  22   def /(val)
  23     Vector3D.new(*self.map {|i| i / val})
  24   end
  25
  26   #in-place add with scale
  27   # a.adds(b, s) -> a += b*s
  28
  29   def adds(other, scale)
  30     self[0] += other[0]*scale; self[1] += other[1]*scale
  31     self[2] += other[2]*scale
  32   end
  33
  34   def subs(other, scale)
  35     self[0] -= other[0]*scale; self[1] -= other[1]*scale
  36     self[2] -= other[2]*scale
  37   end
  38
  39   def magnitude
  40     x=self[0]; y=self[1]; z=self[2]
  41     sqrt(x*x + y*y + z*z)
  42   end
  43
  44   # |self - other|
  45   def dmag(other)
  46     x=self[0]-other[0]; y=self[1]-other[1]; z=self[2]-other[2]
  47     sqrt(x*x + y*y + z*z)
  48   end
  49 end
  50
  51 class Planet
  52   attr_accessor :pos, :v, :mass
  53
  54   def initialize(x, y, z, vx, vy, vz, mass)
  55     @pos = Vector3D.new(x, y, z)
  56     @v = Vector3D.new(vx, vy, vz) * DAYS_PER_YEAR
  57     @mass = mass * SOLAR_MASS
  58   end
  59
  60   def distance(other)
  61     self.pos.dmag(other.pos)
  62   end
  63 end
  64
  65 jupiter = Planet.new(
  66    4.84143144246472090e+00,
  67    -1.16032004402742839e+00,
  68    -1.03622044471123109e-01,
  69    1.66007664274403694e-03,
  70    7.69901118419740425e-03,
  71    -6.90460016972063023e-05,
  72    9.54791938424326609e-04)
  73
  74 saturn = Planet.new(
  75    8.34336671824457987e+00,
  76    4.12479856412430479e+00,
  77    -4.03523417114321381e-01,
  78    -2.76742510726862411e-03,
  79    4.99852801234917238e-03,
  80    2.30417297573763929e-05,
  81    2.85885980666130812e-04)
  82
  83 uranus = Planet.new(
  84    1.28943695621391310e+01,
  85    -1.51111514016986312e+01,
  86    -2.23307578892655734e-01,
  87    2.96460137564761618e-03,
  88    2.37847173959480950e-03,
  89    -2.96589568540237556e-05,
  90    4.36624404335156298e-05)
  91
  92 neptune = Planet.new(
  93    1.53796971148509165e+01,
  94    -2.59193146099879641e+01,
  95    1.79258772950371181e-01,
  96    2.68067772490389322e-03,
  97    1.62824170038242295e-03,
  98    -9.51592254519715870e-05,
  99    5.15138902046611451e-05)
 100
 101 sun = Planet.new(0, 0, 0, 0, 0, 0, 1)
 102
 103 class Array
 104   def each_pair
 105     a = []
 106     each_index {|i|
 107       ((i+1)...length).each {|j|
 108         yield at(i), at(j)
 109       }
 110     }
 111   end
 112 end
 113
 114 bodies = [sun, jupiter, saturn, uranus, neptune]
 115
 116 class << bodies
 117   def advance(dt)
 118     mag = m1 = m2 = nil
 119     each_pair {|b1, b2|
 120       d = b1.distance(b2)
 121       mag = dt/(d*d*d)
 122
 123       m1 = b1.mass * mag
 124       m2 = b2.mass * mag
 125
 126       b1.v.adds(b2.pos, m2)
 127       b1.v.subs(b1.pos, m2)
 128       b2.v.adds(b1.pos, m1)
 129       b2.v.subs(b2.pos, m1)
 130     }
 131
 132     each {|b| b.pos.adds(b.v, dt)}
 133   end
 134
 135   def energy
 136     e = 0
 137     each {|b| e += 0.5 * b.mass * (b.v.magnitude ** 2) }
 138     each_pair {|b1, b2| e -= (b1.mass * b2.mass) / b1.distance(b2) }
 139     e
 140   end
 141
 142   def offset_momentum
 143     p = Vector3D.new(0,0,0)
 144     sun = self[0]
 145     each {|b| p.adds(b.v, b.mass) }
 146     sun.v.subs(p, 1.0/sun.mass)
 147   end
 148 end
 149
 150 bodies.offset_momentum
 151 puts "%.9f" % bodies.energy
 152 Integer(ARGV[0]).times { bodies.advance(0.01) }
 153 puts "%.9f" % bodies.energy