To Do

• -turbo
• using Stopwatch
• profilers
• benchmarks
• avoid when possible: slicing, dynamic, exceptions, i/o
• don't optimize if you don't need to, root of all evil

Example

From the discussion thread in the forums (Why so slow?), some modified versions of the code are available for comparison:

• With Python 3.6.0, running against a slightly modified version (for tracking time) with 500,000 loops yielded a run-time of 0.66819638 seconds.
• An initial port (which included substantial refactoring to identify the cause of the different final results) with 50,000 loops yielded a run-time of 98.284 seconds
• After numerous attempts at speed improvement, the cobra code was sped up to a run-time of 0.057 seconds (or 11x faster than python, or 1,724x faster than original implementation.

Some of the strategies attempting are listed below and whether they yielded much of a performance change:

• Made sure the the for loop variables had explicit types: yielded 33x improvement
• Used float instead of decimal: yielded 7x improvement
• Using multiplication instead of division: yielded 1.2x improvement
• Converted Math.pow to Math.sqrt for square roots: yielded 1.05x improvement
• Moved variable declarations into function-level scope: NO IMPROVEMENT
• Replaced properties with public fields: 1.6x improvement

Discussions

Here are some select discussions on performance:

•  Why so slow? 2010-12
•  Cobra vs Python - speed 2010-11

---

*Python Code:*

# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# originally by Kevin Carson
# modified by Tupteq, Fredrik Johansson, and Daniel Nanz
# modified by Maciej Fijalkowski
# 2to3

import sys

def combinations(l):
    result = []
    for x in range(len(l) - 1):
        ls = l[x+1:]
        for y in ls:
            result.append((l[x],y))
    return result

PI = 3.14159265358979323
SOLAR_MASS = 4 * PI * PI
DAYS_PER_YEAR = 365.24

BODIES = {
    'sun': ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], SOLAR_MASS),

    'jupiter': ([4.84143144246472090e+00,
                 -1.16032004402742839e+00,
                 -1.03622044471123109e-01],
                [1.66007664274403694e-03 * DAYS_PER_YEAR,
                 7.69901118419740425e-03 * DAYS_PER_YEAR,
                 -6.90460016972063023e-05 * DAYS_PER_YEAR],
                9.54791938424326609e-04 * SOLAR_MASS),

    'saturn': ([8.34336671824457987e+00,
                4.12479856412430479e+00,
                -4.03523417114321381e-01],
               [-2.76742510726862411e-03 * DAYS_PER_YEAR,
                4.99852801234917238e-03 * DAYS_PER_YEAR,
                2.30417297573763929e-05 * DAYS_PER_YEAR],
               2.85885980666130812e-04 * SOLAR_MASS),

    'uranus': ([1.28943695621391310e+01,
                -1.51111514016986312e+01,
                -2.23307578892655734e-01],
               [2.96460137564761618e-03 * DAYS_PER_YEAR,
                2.37847173959480950e-03 * DAYS_PER_YEAR,
                -2.96589568540237556e-05 * DAYS_PER_YEAR],
               4.36624404335156298e-05 * SOLAR_MASS),

    'neptune': ([1.53796971148509165e+01,
                 -2.59193146099879641e+01,
                 1.79258772950371181e-01],
                [2.68067772490389322e-03 * DAYS_PER_YEAR,
                 1.62824170038242295e-03 * DAYS_PER_YEAR,
                 -9.51592254519715870e-05 * DAYS_PER_YEAR],
                5.15138902046611451e-05 * SOLAR_MASS) }


SYSTEM = list(BODIES.values())
PAIRS = combinations(SYSTEM)


def advance(dt, n, bodies=SYSTEM, pairs=PAIRS):

    for i in range(n):
        for (([x1, y1, z1], v1, m1),
             ([x2, y2, z2], v2, m2)) in pairs:
            dx = x1 - x2
            dy = y1 - y2
            dz = z1 - z2
            mag = dt * ((dx * dx + dy * dy + dz * dz) ** (-1.5))
            
            b1m = m1 * mag
            b2m = m2 * mag
            v1[0] -= dx * b2m
            v1[1] -= dy * b2m
            v1[2] -= dz * b2m
            v2[0] += dx * b1m
            v2[1] += dy * b1m
            v2[2] += dz * b1m
        for (r, [vx, vy, vz], m) in bodies:
            r[0] += dt * vx
            r[1] += dt * vy
            r[2] += dt * vz


def report_energy(bodies=SYSTEM, pairs=PAIRS, e=0.0):
    #print ("Bodies Count:%s\nPairs Count:%s"%(len(bodies),len(pairs)))
    for (((x1, y1, z1), v1, m1),
         ((x2, y2, z2), v2, m2)) in pairs:
        dx = x1 - x2
        dy = y1 - y2
        dz = z1 - z2
        #print (m1*m2)
        e -= (m1 * m2) / ((dx * dx + dy * dy + dz * dz) ** 0.5)
    for (r, [vx, vy, vz], m) in bodies:
        e += m * (vx * vx + vy * vy + vz * vz) / 2.
    print("%.9f" % e)

def offset_momentum(ref, bodies=SYSTEM):
    px=0.0
    py=0.0
    pz=0.0
    
    for (r, [vx, vy, vz], m) in bodies:
        px -= vx * m
        py -= vy * m
        pz -= vz * m
    (r, v, m) = ref
    v[0] = px / m
    v[1] = py / m
    v[2] = pz / m

def main(n, ref='sun'):
    from time import clock
    start = clock()
    report_energy()
    offset_momentum(BODIES[ref])
    report_energy()
    advance(0.01, n)
    report_energy()
    duration = clock() - start

    print('duration: %s seconds \n\nDone...'%(duration))
if __name__ == '__main__':
    main(50000)

---

*Initial Cobra Port:*

# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# originally by Kevin Carson
# modified by Tupteq, Fredrik Johansson, and Daniel Nanz
# modified by Maciej Fijalkowski
# 2to3
#@number float
class Triple
    pro x from var as decimal 
    pro y from var as decimal 
    pro z from var as decimal
    get squared as decimal
        return .x*.x + .y*.y + .z*.z
    
    cue init(x as decimal, y as decimal, z as decimal)
        base.init
        .x = x
        .y = y
        .z = z
        
    def toString as String is override
        return "\[[.x], [.y], [.z]\]"

class Body
    pro coord from var as Triple
    pro velocity from var as Triple
    pro mass from var as decimal
    
    cue init(coord as Triple, velocity as Triple, mass as decimal)
        base.init
        .coord = coord
        .velocity = velocity
        .mass = mass  
        
class BodyPair
    pro body1 from var as Body
    pro body2 from var as Body  
    
    cue init(body1 as Body, body2 as Body)
        base.init
        .body1 = body1
        .body2 = body2

class Program
    def combinations(l) as List<of BodyPair>
        result = List<of BodyPair>()
        for x in l.count - 1
            ls = l[x+1:]
            for y in ls
                result.add(BodyPair(l[x],y))
        return result

    var bODIES
    cue init
        base.init
    
        .bODIES = {
            'sun': Body(Triple(0.0, 0.0, 0.0), Triple(0.0, 0.0, 0.0), .sOLAR_MASS),

            'jupiter': Body(Triple(4.84143144246472090,
                -1.16032004402742839,
                -0.103622044471123109),
                Triple(0.00166007664274403694 * .dAYS_PER_YEAR,
                0.00769901118419740425 * .dAYS_PER_YEAR,
                -0.0000690460016972063023 * .dAYS_PER_YEAR),
                0.000954791938424326609 * .sOLAR_MASS),

            'saturn': Body(Triple(8.34336671824457987,
                4.12479856412430479,
                -0.403523417114321381),
                Triple(-0.00276742510726862411 * .dAYS_PER_YEAR,
                0.00499852801234917238 * .dAYS_PER_YEAR,
                0.0000230417297573763929 * .dAYS_PER_YEAR),
                0.000285885980666130812 * .sOLAR_MASS),

            'uranus': Body(Triple(12.8943695621391310,
                -15.1111514016986312,
                -0.223307578892655734),
                Triple(0.00296460137564761618* .dAYS_PER_YEAR,
                0.00237847173959480950* .dAYS_PER_YEAR,
                -0.0000296589568540237556* .dAYS_PER_YEAR),
                0.0000436624404335156298* .sOLAR_MASS),

            'neptune': Body(Triple(15.3796971148509165,
                -25.9193146099879641,
                0.179258772950371181),
                Triple(0.00268067772490389322* .dAYS_PER_YEAR,
                0.00162824170038242295* .dAYS_PER_YEAR,
                -0.0000951592254519715870* .dAYS_PER_YEAR),
                0.0000515138902046611451* .sOLAR_MASS) 
        }
        .sYSTEM = [.bODIES['sun'], 
                    .bODIES['jupiter'],
                    .bODIES['saturn'],
                    .bODIES['uranus'],
                    .bODIES['neptune']
        ] 
        .pAIRS = .combinations(.sYSTEM)
    var sYSTEM 
    var _pI = 3.14159265358979323
    var sOLAR_MASS = 4 * 3.14159265358979323 *  3.14159265358979323
    var dAYS_PER_YEAR = 365.24

    var pAIRS


    def advance(dt as decimal, n as int)
        bodies=.sYSTEM
        pairs=.pAIRS

        for i in n
            for p in pairs
                x1, y1, z1 = p.body1.coord.x, p.body1.coord.y, p.body1.coord.z
                x2, y2, z2 = p.body2.coord.x, p.body2.coord.y, p.body2.coord.z
                m1, m2 = p.body1.mass, p.body2.mass
                v1, v2 =p.body1.velocity, p.body2.velocity
                dx = x1 - x2
                dy = y1 - y2
                dz = z1 - z2
                mag = dt * (Math.pow(Convert.toDouble(dx * dx + dy * dy + dz * dz), -1.5f) to decimal)
                b1m = m1 * mag
                b2m = m2 * mag
                v1.x = v1.x - dx * b2m
                v1.y = v1.y - dy * b2m
                v1.z = v1.z - dz * b2m
                v2.x = v2.x + dx * b1m
                v2.y = v2.y + dy * b1m
                v2.z = v2.z + dz * b1m
            for b in bodies
                r = b.coord
                vx,vy,vz = b.velocity.x, b.velocity.y, b.velocity.z
                r.x = r.x + dt * vx
                r.y = r.y + dt * vy
                r.z = r.z + dt * vz


    def report_energy
        bodies=.sYSTEM
        pairs=.pAIRS
        #print "Bodies Count:[bodies.count]\nPairs Count:[pairs.count]"
        e=0.0
        for p in pairs
            m = p.body1.mass * p.body2.mass
            #print m
            dx = p.body1.coord.x - p.body2.coord.x
            dy = p.body1.coord.y - p.body2.coord.y
            dz = p.body1.coord.z - p.body2.coord.z
            sqrs = Convert.toDouble(dx * dx + dy * dy + dz * dz) 
            e = e - (m / Math.pow(sqrs , 0.5f) to decimal)
        for b in bodies
            e = e + b.mass * b.velocity.squared / 2.0
        print e

    def offset_momentum(ref1 )
        bodies=.sYSTEM
        px=0.0 
        py=0.0
        pz=0.0

        for b in bodies
            vx ,vy ,vz  = b.velocity.x, b.velocity.y,b.velocity.z
            m = b.mass
            px = px - vx * m
            py =py - vy * m
            pz = pz - vz * m
        v, m = ref1.velocity, ref1.mass
        v.x = px / m
        v.y = py / m
        v.z = pz / m
        ref1.velocity = v
        print v

    def main
        sw = System.Diagnostics.Stopwatch()
        sw.start
        ref1 = 'sun'
        .report_energy
        .offset_momentum(  .bODIES[ref1])
        .report_energy
        .advance(0.01, 50000)
        .report_energy
        sw.stop
        print 'duration: [sw.elapsedMilliseconds] ms '

---

*Performance Improved version of Cobra code:*

# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# originally by Kevin Carson
# modified by Tupteq, Fredrik Johansson, and Daniel Nanz
# modified by Maciej Fijalkowski
# 2to3
@number float
class Triple
    var x as float 
    var y as float 
    var z as float
    #pro x from var as float 
    #pro y from var as float 
    #pro z from var as float
    get squared as float
        return .x*.x + .y*.y + .z*.z
    
    cue init(x as float, y as float, z as float)
        base.init
        .x = x
        .y = y
        .z = z
        
    def toString as String is override
        return "\[[.x], [.y], [.z]\]"

class Body
    var coord as Triple
    var velocity as Triple
    var mass as float
    #pro coord from var as Triple
    #pro velocity from var as Triple
    #pro mass from var as float
    
    cue init(coord as Triple, velocity as Triple, mass as float)
        base.init
        .coord = coord
        .velocity = velocity
        .mass = mass  

    def updateCoordinates(dt as float)
        .coord.x += dt * .velocity.x
        .coord.y += dt * .velocity.y
        .coord.z += dt * .velocity.z

        
class BodyPair
    var body1 as Body
    var body2 as Body  
    #pro body1 from var as Body
    #pro body2 from var as Body  
    
    cue init(body1 as Body, body2 as Body)
        base.init
        .body1 = body1
        .body2 = body2

class Program
    def combinations(l) as List<of BodyPair>
        result = List<of BodyPair>()
        for x in l.count - 1
            ls = l[x+1:]
            for y in ls
                result.add(BodyPair(l[x],y))
        return result

    cue init
        base.init
    
        bodies as Dictionary<of String, Body> = {
            'sun': Body(Triple(0.0, 0.0, 0.0), Triple(0.0, 0.0, 0.0), .sOLAR_MASS),

            'jupiter': Body(Triple(4.84143144246472090,
                -1.16032004402742839,
                -0.103622044471123109),
                Triple(0.00166007664274403694 * .dAYS_PER_YEAR,
                0.00769901118419740425 * .dAYS_PER_YEAR,
                -0.0000690460016972063023 * .dAYS_PER_YEAR),
                0.000954791938424326609 * .sOLAR_MASS),

            'saturn': Body(Triple(8.34336671824457987,
                4.12479856412430479,
                -0.403523417114321381),
                Triple(-0.00276742510726862411 * .dAYS_PER_YEAR,
                0.00499852801234917238 * .dAYS_PER_YEAR,
                0.0000230417297573763929 * .dAYS_PER_YEAR),
                0.000285885980666130812 * .sOLAR_MASS),

            'uranus': Body(Triple(12.8943695621391310,
                -15.1111514016986312,
                -0.223307578892655734),
                Triple(0.00296460137564761618* .dAYS_PER_YEAR,
                0.00237847173959480950* .dAYS_PER_YEAR,
                -0.0000296589568540237556* .dAYS_PER_YEAR),
                0.0000436624404335156298* .sOLAR_MASS),

            'neptune': Body(Triple(15.3796971148509165,
                -25.9193146099879641,
                0.179258772950371181),
                Triple(0.00268067772490389322* .dAYS_PER_YEAR,
                0.00162824170038242295* .dAYS_PER_YEAR,
                -0.0000951592254519715870* .dAYS_PER_YEAR),
                0.0000515138902046611451* .sOLAR_MASS) 
        }
        .sYSTEM = List<of Body>()
        .sYSTEM.add(bodies['sun'])
        .sYSTEM.add(bodies['jupiter'])
        .sYSTEM.add(bodies['saturn'])
        .sYSTEM.add(bodies['uranus'])
        .sYSTEM.add(bodies['neptune'])
        .pAIRS = .combinations(.sYSTEM)
    var sYSTEM as List<of Body>
    var sOLAR_MASS as float = 4 * 3.14159265358979323 *  3.14159265358979323
    var dAYS_PER_YEAR as float = 365.24

    var pAIRS as List<of BodyPair>


    def advance(dt as float, n as int)
        bodies=.sYSTEM
        pairs=.pAIRS
        b1 as Body = bodies[0]
        b2 as Body = bodies[0]
        v1 as Triple = b1.velocity
        v2 as Triple = b2.velocity
        dx as float
        dy as float
        dz as float
        mag as float

        for i in n
            for p as BodyPair in pairs
                b1 = p.body1
                b2 = p.body2
                v1, v2 =b1.velocity, b2.velocity
                dx = b1.coord.x - b2.coord.x
                dy = b1.coord.y - b2.coord.y
                dz = b1.coord.z - b2.coord.z
                mag = dt * (Math.pow(dx * dx + dy * dy + dz * dz, -1.5f) )
                b1m as float= b1.mass * mag
                b2m as float = b2.mass * mag
                v1.x = v1.x - dx * b2m
                v1.y = v1.y - dy * b2m
                v1.z = v1.z - dz * b2m
                v2.x = v2.x + dx * b1m
                v2.y = v2.y + dy * b1m
                v2.z = v2.z + dz * b1m
            for b as Body in bodies
                b.updateCoordinates(dt)


    def report_energy
        bodies=.sYSTEM
        pairs=.pAIRS
        #print "Bodies Count:[bodies.count]\nPairs Count:[pairs.count]"
        e as float=0.0
        for p as BodyPair in pairs
            b1 as Body = p.body1
            b2 as Body = p.body2
            m = b1.mass * b2.mass
            #print m
            dx = b1.coord.x - b2.coord.x
            dy = b1.coord.y - b2.coord.y
            dz = b1.coord.z - b2.coord.z
            sqrs = dx * dx + dy * dy + dz * dz 
            e -= m / Math.sqrt(sqrs ) 
        for b as Body in bodies
            e += b.mass * b.velocity.squared *0.5
        print e

    def offset_momentum(ref1 as Body)
        bodies=.sYSTEM
        px as float=0.0 
        py as float=0.0
        pz as float=0.0

        for b as Body in bodies
            m as float = b.mass
            px = px- b.velocity.x * m
            py = py-b.velocity.y * m
            pz = pz-b.velocity.z * m
        v, m = ref1.velocity, 1/ ref1.mass
        v.x = px * m
        v.y = py * m
        v.z = pz * m
        print v

    def main
        sw = System.Diagnostics.Stopwatch()
        sw.start
        .report_energy
        .offset_momentum( .sYSTEM[0]) #sun
        .report_energy
        .advance(0.01, 50000)
        .report_energy
        sw.stop
        print 'duration: [sw.elapsed] ms '