get "libhdr"

let calcsigmas(sig, n) be
$(  sig!0 := 0
    for i = 0 to n do sig!i := 0
    for i = 1 to n/2 do
    $(  let j = i
        while 0 < j <= n do
        $(  sig!j := sig!j + i
            j := j + i
        $)
    $)
$)

let gcd(m, n) = n=0 -> m, gcd(n, m rem n)

let duff(sig, n) = sig!n > n+1 & gcd(n, sig!n) = 1
let triple(sig, n) = duff(sig, n) & duff(sig, n+1) & duff(sig, n+2)

let first(sig, f, max, cb) be
$(  let n = 0
    for i = 1 to max
    $(  n := n+1 repeatuntil f(sig, n)
        cb(i, n)
    $)
$)

let start() be
$(  let showsingle(i, n) be
    $(  writef("%I4", n)
        if i rem 10=0 then wrch('*N')
    $)

    let showtriple(i, n) be writef("%I2: %I6 %I6 %I6*N", i, n, n+1, n+2)

    let sig = getvec(20000)
    calcsigmas(sig, 20000)

    writes("First 50 Duffinian numbers:*N")
    first(sig, duff, 50, showsingle)

    writes("*NFirst 15 Duffinian triples:*N")
    first(sig, triple, 15, showtriple)
    freevec(sig)
$)