mandel.bcpl

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

let Main() be
[
Ws("Hello World!*N")
let v = vec 30705
v = (v + 1) & -2 // Data needs to be 32-bit aligned (?)
let dcb = vec 5
dcb = (dcb + 1) & -2
dcb!0 = 0 // End of disply list
dcb!1 = 38 // 38 words per line
dcb!2 = v // Data pointer
dcb!3 = 404 // # lines / 2
let lvdas = #420
lvdas!0 = dcb

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

let x0 = (-2) lshift 12
let x1 = 1 lshift 12
let y0 = (-1) lshift 12
let y1 = 1 lshift 12
let xstep = (x1 - x0) / 600
let ystep = (y1 - y0) / 400
let x2 = vec 2
let y2 = vec 2
let xy = vec 2

let cy = y0
for ypos = 0 to 400 do // line count
   [
   let cx = x0
   for h = 0 to 37 do // horizontal word count
      [
      for b = 0 to 15 do // horizontal bit count
         [
         let x = cx
         let y = cy
         for n = 0 to 999 do
            [
            MulFull(y, y, y2)
            MulFull(x, x, x2)
            if x2!0 + y2!0 ge 1024 then break // Quit if x^2 + y^2 > 4

            if n eq 20 then // Last iteration. Inside set.
               [
               let adr = (200 + ypos) * 38 + h
               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)
            let y2sp = (y2!0 lshift 4) % (y2!1 rshift 12)

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

            // 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
         ]
      ]
   cy = cy + ystep
   ]
Gets(keys)
]