leonardo - April 5th, 2008

Tags:	d language, programming, psyco, python
Subject:	Euler problem 14
Time:	05:23 pm

The Euler problem n.14:

The following iterative sequence is defined for the set of positive integers:
n -> n/2 (n is even)
n -> 3n + 1 (n is odd)
Using the rule above and starting with 13, we generate the following sequence of 10 terms:
13 -> 40 -> 20 -> 10 -> 5 -> 16 -> 8 -> 4 -> 2 -> 1
Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1. Which starting number, under one million, produces the longest chain of Collatz numbers?

My Python solution, that uses memoization of results to speed up computation:

import array

def main(m):
  lens = array.array('H', [0]) * m
  max_steps = max_steps_pos = 1
  for i in xrange(1, m):
    n = i
    n_steps = 1
    while n != 1:
      if n < i:
        n_steps += lens[n] - 1
        break
      n = 3 * n + 1 if n & 1 else n // 2
      n_steps += 1
    if n_steps > max_steps:
      max_steps = n_steps
      max_steps_pos = i
    lens[i] = n_steps
  return max_steps, max_steps_pos

import psyco; psyco.bind(main)
print "max_steps, max_steps_pos:", main(1000000)

Psyco is often able to manage array.array very efficiently.
'H' represents unsigned shorts [0 .. 2^16-1], because the chains are never too much long. This allows to reduce memory used and improve CPU cache efficiency (mostly in the D version). In Python you can use 'l' that is safer and doesn't slow down the code much, but uses twice the memory (so you can manage smaller problems).

The D version is essentially the same, the only significant difference is that the array isn't pre-initialized:

import std.stdio, std.c.stdlib;

void main() {
  const uint m = 1_000_000;
  auto lens = (cast(ushort*)malloc(m * ushort.sizeof))[0 .. m];
  uint max_steps = 1;
  uint max_steps_pos = 1;

  for (uint i = 1; i < m; i++) {
    uint n = i;
    uint n_steps = 1;

    while (n != 1) {
      if (n < i) {
        n_steps += lens[n] - 1;
        break;
      }
      n = n & 1 ? 3 * n + 1 : n / 2;
      n_steps++;
    }
    if (n_steps > max_steps) {
      max_steps = n_steps;
      max_steps_pos = i;
    }
    lens[i] = n_steps;
  }
  writefln("max_steps, max_steps_pos: ", max_steps, " ", max_steps_pos);
}


Timings:
m = 1_000_000:
D/C:    0.29 s  ( 1   X)
Psyco:  1.60 s  ( 5.5 X) with array 'H'
Python: 19.20 s (66   X) with list

A C version is very similar to that D version, and the running time (with GCC 4.2.1) is the same, but the executable is almost 20 times smaller.

The D code scans m = 100_000_000 in 25 s finding a chain of length 758 on n= 83_706_505.

An Haskell version:
http://www.haskell.org/haskellwiki/Euler_problems/11_to_20#Problem_14

To improve cache efficiency I have tried to store half of the values, but the code is slower (it allows to scan twice the values with the same max memory).

comments: Leave a comment

Tags:	d language, programming, python
Subject:	Self-recursive generators in D
Time:	07:45 pm

I have (slightly) adapted the Python-like generators for D written by Witold Baryluk to my libs. They allow a simpler syntax that can be used for single iterations (not parallel scan allowed, you have to use some threads or light threads). So I have used them to translate to D some self-recursive generators from Python code by David Eppstein:

http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/221457

As an example this yields Fibbinary numbers (no consecutive ones in binary representation):
http://www.research.att.com/projects/OEIS?Anum=A003714

struct A003714 {
  void generator() {
    yield(1);
    foreach(x; A003714()) {
      yield(2 * x);
      if (!(x & 1))
        yield(2 * x + 1);
    }
  }
  mixin Generator!(int);
}

The original Python version:

def A003714():
  yield 1
  for x in A003714():
    yield 2 * x
    if not (x & 1):
      yield 2 * x + 1

The syntax isn't as clean as Python one, but it's simple enough anyway. They are 2-3 times slower than handwritten opApply struct methods, so you may want to use them only where code size/readability matters more.

( Click here to see more examples )

comments: Leave a comment

April 5th, 2008

April 5th, 2008

Advertisement