Mandelbrot set code for the Xerox Alto, written in BCPL.

mandel.bcpl

get "streams.d"
external
[
Ws;
Wns;
MulFull;
DoubleAdd;
keys;
Gets;
]

let Main() be
[
let v = vec 30705 // Pixel buffer
v = (v + 1) & -2 // Data needs to be 32-bit aligned
let dcb = vec 5 // Display control block: defines display region
dcb = (dcb + 1) & -2
dcb!0 = 0 // End of display list
dcb!1 = 38 // 38 words per line
dcb!2 = v // Data pointer
dcb!3 = 404 // # lines / 2
let lvdas = #420 // Address holds pointer to display control block
lvdas!0 = dcb

for i = 0 to 30703 do v!i = 0 // Clear display

// Values are represented as fixed point with 12 bits to right of decimal point.
let x0 = (-2) lshift 12 // Left boundary: x = -2
let x1 = 1 lshift 12 // Right boundary: x = 1
let y0 = (-1) lshift 12 // Top boundary: y = -1
let y1 = 1 lshift 12 // Bottom boundary: y = 1
let xstep = (x1 - x0) / 600 // Render 600 pixels horizontally
let ystep = (y1 - y0) / 400 // Render 400 pixels vertically
let x2 = vec 2 // double word to hold x^2
let y2 = vec 2 // double word to hold y^2
let xy = vec 2 // double word to hold x * y

let cy = y0 // Constant value, y part. I.e. the z value for this pixel
for ypos = 0 to 400 do // line count. Note "for" limits are inclusive.
   [
   let cx = x0 // Constant value, x part.
   for h = 0 to 37 do // horizontal word count
      [
      for b = 0 to 15 do // horizontal bit count
         [
         let x = cx // The complex z value is represented as x + i*y
         let y = cy
         for n = 0 to 999 do // Will bail out long before 999
            [
            MulFull(y, y, y2) // y2 = y*y. Integer multiplication, y2 is double word.
            MulFull(x, x, x2) // x2 = x*x
            if x2!0 + y2!0 ge 1024 then break // Quit if x^2 + y^2 > 4

            if n eq 20 then // Last iteration. Still inside set, so set pixel.
               [
               let adr = (200 + ypos) * 38 + h // 200 blank pixels at top, 38 words per line
               v!adr = v!adr % (1 lshift (15-b)) 
               break
               ]

            // Convert to single precision by dropping 12 bits.
            // rshift 12 = lshift top word 4
            let x2sp = (x2!0 lshift 4) % (x2!1 rshift 12) // x2sp = x^2, single precision
            let y2sp = (y2!0 lshift 4) % (y2!1 rshift 12) // y2sp = y^2, single precision

            MulFull(x, y, xy) // xy = x*y
            let xysp = (xy!0 lshift 4) % (xy!1 rshift 12) // xysp = x*y, single precision

            // z = z^2 + c (complex arithmetic, z = x+iy)
            // i.e. y = 2xy + cy
            //      x = x^2 - y2 + cx
            y = xysp + xysp + cy
            x = x2sp - y2sp + cx
            ]

         cx = cx + xstep // Move to next cx value
         ]
      ]
   cy = cy + ystep // Move to next cy value
   ]
Gets(keys) // Get a key, i.e. wait for a keypress
]