2010-08-05

Reviving an Old Tiger

An old machine that proved its mettle on the line can be a handy development server, on a pinch. Especially if you're starting up by pulling your boot straps, as I currently am. I have this Blue & White G3 that served me faithfully for the last, say, 10 years. Bitrot, alas, is not a hardware problem: software version numbers climb and software authors drop support for older operating systems and configurations, and we misers are left struggling and scavenging for information to make things compile and install. So, for my own reference, here's what I did this time. This applies to Mac OS X 10.4.11, AKA Tiger:

Sybase 12
Get the latest EBF for Mac OS, EBF17473. It works beautifully.
XCode 2.5
I had to root quite a bit around the newly redesigned Apple Developer site. The image name is xcode25_8m2558_developerdvd.dmg
MySQL 5.1.49
A challenge for GCC 4.0.1, it kernel-panicked the machine once. I followed these fine instructions. I set it up to start up automatically on boot, via the following launchd script:

<!DOCTYPE plist PUBLIC "-//Apple Computer//DTD PLIST 1.0//EN" "http://www.apple.com/DTDs/PropertyList-1.0.dtd">
<plist version="1.0">
<dict>
 <key>KeepAlive</key>
 <true/>
 <key>Label</key>
 <string>com.mysql.mysqld</string>
 <key>Program</key>
 <string>/usr/local/mysql/bin/mysqld_safe</string>
 <key>RunAtLoad</key>
 <true/>
 <key>UserName</key>
 <string>mysql</string>
 <key>WorkingDirectory</key>
 <string>/usr/local/mysql</string>
</dict>
</plist>

Put it in /Library/LaunchDaemons/com.mysql.mysqld.plist.

Apache httpd 2.2.16
I followed these instructions. The layout I used is this:
<Layout DarwinLocal>
    prefix:        /usr/local/apache2
    exec_prefix:   ${prefix}
    bindir:        ${exec_prefix}/bin
    sbindir:       ${exec_prefix}/sbin
    libdir:        ${exec_prefix}/lib
    libexecdir:    ${exec_prefix}/modules
    mandir:        ${prefix}/man
    datadir:       /Library/Documents
    sysconfdir:    /etc/httpd
    installbuilddir: ${datadir}/build
    errordir:      ${datadir}/error
    iconsdir:      ${datadir}/icons
    htdocsdir:     /Library/WebServer/Documents
    manualdir:     ${datadir}/manual
    cgidir:        /Library/WebServer/CGI-Executables
    includedir:    ${prefix}/include+
    localstatedir: /var
    runtimedir:    ${localstatedir}/run
    logfiledir:    ${localstatedir}/log/httpd
    proxycachedir: ${runtimedir}/proxy
</Layout>

and the configuration is this (culled from here):

CFLAGS="-arch ppc -isysroot /Developer/SDKs/MacOSX10.4u.sdk" ./configure \
--enable-layout=DarwinLocal \
--enable-access \
--enable-actions --enable-alias --enable-asis --enable-auth \
--enable-auth_dbm --enable-auth_digest --enable-autoindex \
--enable-cache --enable-cgi --enable-dav --enable-dav_fs \
--enable-deflate --enable-dir --enable-disk_cache \
--enable-dumpio --enable-env --enable-expires --enable-fastcgi \
--enable-file_cache --enable-headers --enable-imap \
--enable-include --enable-info --enable-log_config \
--enable-log_forensic --enable-logio --enable-mem_cache --enable-mime \
--enable-mime_magic --enable-negotiation --enable-perl \
--enable-rewrite --enable-setenvif --enable-speling --enable-ssl \
--enable-status --enable-suexec --enable-unique_id --enable-userdir \
--enable-usertrack --enable-version --enable-vhost_alias --enable-so \
--enable-module=all --enable-shared=max

These instructions allow automatic launching of Apache at startup.

zlib 1.2.5
libpng 1.4.3
Beware, as the DLL gets installed with the wrong name: it should be libpng14.14.dylib.
jpegsrc v8b
freetype 2.4.1
gd-2.0.35
I configured it with the usual caveats.
libmcrypt 2.5.8
mhash 0.9.9.9
mcrypt-2.6.8
php 5.3.3
It took me quite a long while sort out the correct incantation to configure PHP. This is what worked for me:
export CFLAGS="-DSQLITE_ENABLE_LOCKING_STYLE=0 -DBIND_8_COMPAT"
./configure \
-prefix=/usr/local/php5 \
-with-config-file-path=/etc \
-with-zlib=/usr/local \
-with-xml \
-with-xsl \
-enable-ctype \
-enable-dom \
-enable-exif \
-enable-filter \
-enable-ftp \
-enable-gd-native-ttf \
-enable-libxml \
-enable-mbregex \
-enable-mbstring \
-enable-pcntl \
-enable-pdo \
-enable-posix \
-enable-shmop \
-enable-simplexml \
-enable-soap \
-enable-xml \
-enable-bcmath \
-with-bz2=/usr/bin \
-with-curl \
-with-freetype-dir=/usr/local \
-with-gd=/usr/local \
-with-iconv \
-with-jpeg-dir=/usr/local \
-with-ldap \
-with-mcrypt=/usr/local \
-with-mysql=/usr/local/mysql \
--with-mysqli=/usr/local/mysql/bin/mysql_config \
-with-openssl \
-with-pdo-mysql=/usr/local/mysql \
-with-png-dir=/usr/local \
-with-zlib=/usr/local \
-enable-dbx \
-enable-sockets \
-with-apxs2=/usr/local/apache2/sbin/apxs \
-with-kerberos=/usr \
-with-mysqlsock=/var/mysql/mysql.sock \
-without-pear

Don't forget to apply the patch described in this post. Be sure to make test!

Next, configure Apache to run phpMyAdmin, Drupal and Mantis (I won't have anything to do with either Python or Ruby, sorry). I expect this old Tiger to bend under the pressure but not break.

2010-07-02

Sukhotin's Algorithm

There is an obscure algorithm devised by one Boris Viktorovich Sukhotin, a Russian… linguist? mathematician? cryptographer? let's say researcher, in the early '70s, usually called Sukhotin's Vowel Identification Algorithm, whose initial aim was to quickly attack an alphabetic cyphertext encoded with an unknown substitution cypher. Jacques Guy, a linguist and programmer well-known for his interest in the Voynich manuscript, translated the algorithm from Russian and simplified it but his publishing in Cryptologia didn't help it being more widely known. Besides, his explanation is not very illuminating, so I'll show the algorithm and then try to explain it from first principles.

This classification algorithm, given a text over an alphabet Σ of size N = |Σ|, partitions Σ into two disjoint subsets V and K (for vowels and consonants). It takes as input the symmetric digraph frequencies F[xy] = F[yx], that is, the frequency with which letter xΣ appears next to letter yΣ in the text. The algorithm is predicated on the assumption that in natural languages vowels tend to appear next to consonants rather than to other vowels. On input, the array d records the frequencies F[xy], subject to the condition:

⟨∀i, j : 0 ≤ i, j < N : d.i.j = d.j.i ≥ 0⟩

On exit, the stack v will contain the letters identified as vowels by the algorithm. I use Dijkstra's guarded command language, except that I omit trivial skip branches in single-branched ifs.

begin
   var v: stack of int;
   var r: array [0 … N) of int;
   var m, i, j: int;
   m, i := 0, 0;
   do iNr.i, j := 0, 0;
      do jNr.i := r.i + d.i.j;
         j := j + 1
      od;
      if r.m < r.im := i fi;
      i := i + 1
   od;
   { P }
   do r.m > 0 →
      v:push.m;
      i, j := 0, 0;
      do iNif imr.i :=  r.i − 2⋅d.i.mi = mr.i := −r.i
         fi;
         if r.j < r.ij := i fi;
         i := i + 1
      od;
      m := j
      { P }
   od
   { v = V }
end

The algorithm initially assumes that every letter is a consonant, that is, initially K = Σ. The first loop computes in r.i the cumulative frequency of character i, recording in m the character with maximum frequency, thus establishing the initial validity of the loop invariant P :

P ≡ ⟨∀ i : 0 ≤ i < N : r.i = ⟨∑j : jK : F[ij]⟩ ∧ r.mr.i

The second loop tries to maximize:

Q = ⟨∑i : iV : ⟨∑j : jK : F[ij]⟩⟩

under invariance of P. This measure is motivated by the original assumption that vowels are more likely to be preceded or followed by consonants than by other vowels, and that a partition of Σ should maximize this likelihood. The second loop proceeds to hill-climb to this maximum by assuming that the most frequent letter is a vowel, adding each r.m to V and removing it from K under the invariance of P. Its inner loop restores it by discounting twice the frequency for m: once for d.m.i and then for d.i.m, except for d.m.m, since after reduction it is precisely equal to r.m. It also keeps track of the new maximum for the next iteration. This process continues as long as Q can increase, which is implied by the loop condition.

References

  1. B. V. Sukhotin. DECIPHERING METHODS AS A MEANS OF LINGUISTIC RESEARCH, COLING Vol.1 1973. Online.
  2. ———. OPTIMIZATION ALGORITHMS OF DECIPHERING AS THE ELEMENTS OF A LINGUISTIC THEORY, COLING Budapest Vol.1 1988. Online.
  3. Jacques B. M. Guy. VOWEL IDENTIFICATION: AN OLD (BUT GOOD) ALGORITHM, CRYPTOLOGIA Vol.XV No.3 1991. Online.
  4. David M. W. Powers. Unsupervised learning of linguistic structure. IJCL Vol.2 No.1 1997. Online.

2010-06-11

Diagonalize, Now!

It is easy to append two lists, and the method is well-known:

let rec append xs ys = match xs with
| []      -> ys
| x :: xs -> x :: append xs ys

It is also easy to merge two lists, so that elements from the first are interleaved with those of the second:

let rec merge xs ys = match xs with
| []      -> ys
| x :: xs -> x :: merge ys xs

A simple matter of swapping two arguments! The similarities between both functions hide deep differences. In a strict setting, the equation append xs ys = xs has just a solution ys = [] if xs is finite; but in a lazy setting this equation has infinite solutions for ys if xs is infinite, because AωA* ≡ Aω. By contrast, merge xs ys = xs has as its only solution ys = [] independently of whether xs is finite or not (in a lazy language).

What this all means is that merge is fairer than append, in this sense: merge xs ys includes all the elements of both xs and ys. This is, in transfinite terms, one half of the constructive proof that ω + ω > ω but 2⋅ω = ω, the other being an effective procedure for un-merge-ing an infinite list (edit: thanks Andrej Bauer for your thorough comments!). In practice, it is the reason why you cannot enumerate all the integers by first counting the nonnegatives 0, 1, 2, 3… and then the negatives -1, -2, -3… but you can if you interleave them 0, -1, 1, -2, 2, -3, 3…

The downside to this is that append is associative but merge is not, so it doesn't form a monoid over lists (even though [] is its left and right identity) and thus neither endows the lists with monadic structure.

Can we generalize to lists of lists? Sure we can:

let rec concat = function
| []        -> []
| xs :: xss -> append xs (concat xss)

(a simple consequence of lists being a monoid), or:

let concat xss = List.fold_right (fun xs ys -> append xs ys) xss []

If we try to do the same with merge:

let rec diagonalize = function
| []        -> []
| xs :: xss -> merge xs (diagonalize xss)

again, a simple change (I omit the folded version for brevity). The problem is that, since merge is not associative, the result of diagonalize [xs0; xs1; xs2; …] is not very fair: it includes every other element of xs0 followed by every fourth element of xs1 followed by every eighth element of xs2…, so that you have to wait 2n + 1 elements before finding the next element in xsn! The result does include every element in every list even if there are an infinite number of infinite lists, because ½ + ¼ + ⅛ + … = 1, but you have to wait exponentially long for them. Georg Cantor devised a clever way out of this: read the (possibly doubly infinite) matrix:

Matrix of elements

(thank you codecogs for the wonderful LATEX rendering service!) along ascending diagonals, as if it were rotated 45° to the right: x00, x01, x10, x02, x11, x20, x03, x12, x21, x30, … You get the idea. Now you only have to wait quadratically long before you see an element. However, how should we proceed in order to flatten a list of lists in this order? By induction, of course!

Suppose every column in the matrix is a list, and that we can "half-transpose" such a list of lists so that each list in the result corresponds to a diagonal in the original matrix: [[x00]; [x01; x10]; [x02; x11; x20]; [x03; x12; x21; x30]; …]. Call this operation rotate. Then, by elementwise consing the head list with each list in the rotation of the rest we diagonalize the entire list of lists. The code for pairing up these lists is:

let rec pairup xs yss = match xs, yss with
| [], yss -> yss
| xs, []  -> List.map (fun x -> [x]) xs
| x :: xs, ys :: yss -> (x :: ys) :: pairup xs yss

(being careful not to lose any element!). The code to perform the rotation is:

let rec rotate = function
| []        -> []
| xs :: xss -> pairup xs ([] :: rotate xss)

Note that the corner element gets to be by itself. Then the diagonalization is simply the flattening of this rotation:

let diagonalize xss = concat (diag xss)

The question now is, how to generalize this to higher dimensions (lists of lists of lists of …)? Since nowhere to be seen is the nice symmetry between appending and merging here I don't know what the guiding principle could be.

2010-06-10

You are Here

I am a terrible person: I read code, I think to myself "I can do better" and I set out thus, out of pride. I did just that with an OCaml implementation of the geohash geographical encoding. Even if the code I read was perfectly workable, I found it unidiomatic and somewhat bloated. I wrote two versions, the first one compact but impenetrable, the second one from first principles: what would be an executable description of the algorithm? To wit:

let decode =
   (narrow (-90., 90.) *** narrow (-180., 180.)) % swap % split
  % List.concat % List.map (to_bits 5 % decode_base32) % explode

Because of (or despite!) the point-free, combinator-heavy style I intended this to be read as a pipeline or as a recipe, from right to left (it is an applicative-order pipeline), much as the algorithm description does:

  1. Explode the code string into a list of characters
  2. For each character, decode it as a base-32 digit decomposing it into a list of 5 bits
  3. Flatten the list of lists of bits into a big list of bits
  4. Split the list of bits into a pair of lists for the even bits and the odd bits
  5. Swap both components of the pair, since the coding uses the even bits for the longitude
  6. Narrow each list independently into a coordinate starting from the corresponding interval

The result is a (latitude, longitude) pair. Supposing that every step is invertible, the algorithm for encoding a coordinate into a geohash is:

let encode nchars =
  let bits = 5 * nchars in
  let lo = bits / 2 in
  let hi = bits - lo in
    implode % List.map (encode_base32 % of_bits) % group 5 % join
  % swap % (expand lo (-90., 90.) *** expand hi (-180., 180.))

Again, this is easy to follow as a recipe or as a succint description:

  1. Given the desired hash length, compute the number of bits for the latitude and for the longitude
  2. Expand each coordinate independently into a list of bits describing its position in the corresponding interval
  3. Swap both components of the pair
  4. Join both lists by alternating even and odd elements
  5. Group the bits in the list into sublists by taking 5 bits at a time
  6. For each group, convert the list of 5 bits into an integer and encode it as a base-32 digit
  7. Implode the list of characters into a code string

Now I have a number of inverse or near-inverse pairs of functions to write, namely:

  • explode and implode
  • decode_base32 and encode_base32
  • to_bits and of_bits
  • List.concat and group
  • split and join
  • narrow and expand

Let's start with the (well-known to Haskellers) combinators (though they pronounce % as .):

let ( $ ) f x = f x
let ( % ) f g x = f (g x)

The first is the application operator and the second is the composition operator. Next are some operators on pairs borrowed from Control.Arrow:

let first f (x, y) = (f x, y)
let second f (x, y) = (x, f y)
let ( *** ) f g (x, y) = (f x, g y)
let swap (x, y) = (y, x)

They allow applying a function to one or both members of a pair, and to manipulate its components. Note that swapswapid. This completes the functional scaffolding. The first and more imperative pair of functions is explode and implode:

let explode s =
  let res = ref [] in
  for i = String.length s - 1 downto 0 do
    res := String.unsafe_get s i :: !res
  done;
  !res

let implode l =
  let b = Buffer.create 10 in
  List.iter (Buffer.add_char b) l;
  Buffer.contents b

The first traverses the string from right to left to accumulate each character on the head of the resulting list; the second uses a string buffer to build the resulting string a character at a time. I use unsafe_get because, well, these are as low-level functions as they come, so squeezing a bit extra from them doesn't seem inappropriate. Each is the other's inverse, or in more abstract terms the pair (explode, implode) witnesses the monoid isomorphism between OCaml strings and lists of characters (that's why Haskell identifies both datatypes).

Aside: I won't attempt any proofs to avoid making a long post even longer. Some are easy, most are tedious, and the preceding is probably difficult. I've verified a few of the inverse pairs but not all of them, so I don't consider this code verified by any stretch of the word. I am, however, confident that the code works.

Now decode_base32and encode_base32 are also fairly imperative, but just because I chose to use an array as the map between characters and their base-32 encodings:

let codetable = [|
'0'; '1'; '2'; '3'; '4'; '5'; '6'; '7';
'8'; '9'; 'b'; 'c'; 'd'; 'e'; 'f'; 'g';
'h'; 'j'; 'k'; 'm'; 'n'; 'p'; 'q'; 'r';
's'; 't'; 'u'; 'v'; 'w'; 'x'; 'y'; 'z'
|]

This means that for encoding a 5-bit number indexing into the array is enough, but decoding a base-32 character requires a lookup. Since the array is sorted by character, a binary search is a good choice:

let binsearch ?(compare=Pervasives.compare) v e =
  let i = ref 0
  and j = ref (Array.length v) in
  if !j = 0 then raise Not_found else
  while !i <> !j - 1 do
    let m = (!i + !j) / 2 in
    if compare v.(m) e <= 0
      then i := m
      else j := m
  done;
  if v.(!i) = e then !i else raise Not_found

let encode_base32 i = codetable.(i land 31)
and decode_base32 c = binsearch codetable (Char.lowercase c)

(This binary search is completely general). Now decode_base32encode_base32id (for the restricted domain of 5-bit integers) if and only if binsearch v v.(i) = i for i in range of v. Similarly, encode_base32decode_base32id (again, for the restricted domain of base-32 characters) if and only if v.(binsearch v x) = x. The conditions on binsearch amount to it being correct (and would make for excelent unit or random testing).

The next inverse pair is to_bits and of_bits:

let iota =
  let rec go l i =
    if i <= 0 then l else let i = i - 1 in go (i :: l) i
  in go []

let to_bits k n = List.rev_map (fun i -> (n land (1 lsl i)) != 0) $ iota k

let of_bits = List.fold_left (fun n b -> 2 * n + if b then 1 else 0) 0

The first maps a list of bit positions (given by iota, the initial natural interval) into the corresponding bits by direct testing, most-significant (highest) bit first; the second computes the binary value of the list of bits by a Horner iteration.

The following functions are very general functions on lists. Any function f : α listα list list such that List.concatfid is called a partition. One such function is group:

let group k l =
  let rec go gs gss i = function
  | []           -> List.rev (if gs = [] then gss else List.rev gs :: gss)
  | l when i = 0 -> go [] (List.rev gs :: gss) k l
  | x :: xs      -> go (x :: gs) gss (i - 1) xs
  in go [] [] k l

It collects up to k elements in the stack gs, pushing complete groups into the stack gss. The last group might be incomplete and it is added to the top of the stack before returning it. A function that I will need later is the the converse of List.fold_right, appropriately called unfold:

let rec unfold f e n =
  if n <= 0 then [] else
  let (x, e') = f e in
  x :: unfold f e' (n - 1)

unfold generates a list of the given length n from a seed e and a function f that computes the head and the next seed. Now split and join are another very functional couple:

let cons x xs = x :: xs

let split l =
  let rec even = function
  | [] -> [], []
  | x :: xs -> first (cons x) $ odd xs
  and odd = function
  | [] -> [], []
  | x :: xs -> second (cons x) $ even xs
  in even l

let rec join (xs, ys) = match xs with
| [] -> ys
| x :: xs -> x :: join (ys, xs)

In general, splitjoin is not the identity on pairs of lists (think of ([], l) for nonempty l), but indeed joinsplitid, and a proof could proceed by induction. Now all these functions are completely general (except perhaps for the base-32 encoding and decoding) and have little if anything to do with geocaching. The last pair, narrow and expand, are the meat of the algorithm:

let avg (x, y) = (x +. y) /. 2.

let bisect (lo, hi as interval) b =
  let md = avg interval in
  if b then (md, hi) else (lo, md)

let narrow interval = avg % List.fold_left bisect interval

narrow repeatedly bisects the interval, keeping the upper or lower half depending on the next bit in the input list, finally returning the midpoint of the final interval. Its inverse is expand:

let divide x (lo, hi as interval) =
  let md = avg interval in
  if x > md
    then true , (md, hi)
    else false, (lo, md)

let expand bits interval x = unfold (divide x) interval bits

It unfolds each bit in turn, by divideing the interval and finding which half contains x. These functions are inverses in the sense that, for any interval i = (min, max) and list l, expand (List.length l) i (narrow i l) = l and conversely, for any n ≥ 0 and minxmax:

|narrow i (expand n i x) - x| ≤ (max - min)2-n-1

These are the crucial properties that make encode the true inverse of decode. The astute reader might have noticed the resemblance of geohashing with arithmetic coding: in a sense, a geohash is the base-32 encoding of the coordinates lossly compressed with an uniform probability model. One way to improve the precision of the encoding would be by allocating a predefined weight to each hemispherical quadrant, by taking into account that more landmass is present in the north-eastern quadrant of the globe. This would require splitting the intervals not along the middle but, say, the first third, thus representing more precise locations with the same number of bits in continental Europe and the Far East, then in North America, then in Africa and Oceania, and lastly in South America.

2010-05-28

Gematric Search

Given an array A.(0 ≤ i < N) of nonnegative integers, that is:

(0)   ⟨∀ i : 0 ≤ i < N : A.i ≥ 0⟩

and a positive integer s > 0, we want to find a segment of A that sums to s, if it exists. Define:

      S.i.j = ⟨∑ k : 0 ≤ ik < jN : A.k

The gematric search for s in A is to find, if possible, a pair (p, q) such that S.p.q = s (it is hoped that conoisseurs would readily make the connection). Note a number of things:

  • S.i.j is a total function on i, j
  • By definition, the sum on an empty range is zero:
    (1)   S.i.i = 0
    

    for all integer i.

  • For arguments in range of A, S.i.j is monotone on its second argument and antitone on its first. That is, given 0 ≤ i′i < jj′N:
    (2)   S.i.j ≤ S.i.j′
    (3)   S.i.j ≤ S.i′.j
    

What should the result be if no such segment is present in A? One possibility would be setting a flag found, such that the program ensures:

      ¬found ∨ S.p.q = s

An alternative would be to guarantee some condition on (p, q) equivalent by design to found. Given that s ≠ 0, one such condition would be satisfied by S.p.q = 0 as a suitable sentinel. By (1), this is implied by p = q, and we choose the weaker postcondition:

(4)   p = q ∨ S.p.q = s

Now (4) is trivially satisfied by:

{ A.(0 ≤ i < N) ∧ s > 0 }
gemsearch ≡
   p := 0; q := 0; sum := 0;
   { reduce |sum - s| under invariance of (4) }
{ p = q ∨ S.p.q = s }

where sum = S.p.q at each step. How should we meet the refinement obligation in the program? We clearly have three cases:

  • sum = s

    We have found our segment, and we can stop looking.

  • sum < s

    By (2) we can extend the segment on the right, provided we can, that is, if qN.

  • sum > s

    By (3) we can shrink the segment on the left, provided we can, that is, if pq. But since s > 0, this is implied by pN, a weaker condition symmetric to the previous one and thus preferable. It is also suitable for the repetition guard, since its negation implies the postcondition.

We have:

{ A.(0 ≤ i < N) ∧ s > 0 }
gemsearch ≡
   p := 0; q := 0; sum := 0;
   do pNsumsif sum < sqNsum := sum + A.q; q := q + 1
      [] sum > ssum := sum - A.p; p := p + 1
      [] … → ?
      fi
   od
{ p = q ∨ S.p.q = s }

The conditional is not exhaustive. Let's calculate what the obligations are for third case:

   ¬((sum < sqN) ∨ sum > s)
≡ { DeMorgan, twice }
   (sumsq = N) ∧ sums
≡ { Distribute }
   (sumssums) ∨ (q = Nsums)
≡ { Tertium non datur }
   sum = s ∨ (q = Nsums)
⇐ { Loop condition }
   q = Nsum < s

This condition expresses a deficient suffix of A which is, by (3), as large as it can be, hence we can quit the search. Left-factoring the common condition on s, we arrive at our final program:

{ A.(0 ≤ i < N) ∧ s > 0 }
gemsearch ≡
   p := 0; q := 0; sum := 0;
   do pNsumsif sum < sif qNsum := sum + A.q; q := q + 1
         [] q = Np := N
         fi
      [] sum > ssum := sum - A.p; p := p + 1
      fi
   od
{ p = q ∨ S.p.q = s }

The algorithm is of complexity obviously linear on N, since no element of A is accessed more than twice.