sequence fib = {1,1}

function zeckendorf(atom n)
-- Same as [[Zeckendorf_number_representation#Phix]]
atom r = 0
    while fib[$]<n do
        fib &= fib[$] + fib[$-1]
    end while
    integer k = length(fib)
    while k>2 and n<fib[k] do
        k -= 1
    end while
    for i=k to 2 by -1 do
        integer c = n>=fib[i]
        r += r+c
        n -= c*fib[i]
    end for
    return r
end function

function decimal(object z)
-- Convert Zeckendorf number(s) to decimal
atom dec = 0, bit = 2
    if sequence(z) then
        for i=1 to length(z) do
            z[i] = decimal(z[i])
        end for
        return z
    end if
    while z do
        if and_bits(z,1) then
            dec += fib[bit]
        end if
        bit += 1
        if bit>length(fib) then
            fib &= fib[$] + fib[$-1]
        end if
        z = floor(z/2)
    end while
    return dec
end function

function to_bits(integer x)
-- Simplified copy of int_to_bits(), but in reverse order,
-- and +ve only but (also only) as many bits as needed, and
-- ensures there are *two* trailing 0 (most significant)
    sequence bits = {}
    if x<0 then ?9/0 end if     -- sanity/avoid infinite loop
    while 1 do
        bits &= remainder(x,2)
        if x=0 then exit end if
        x = floor(x/2)
    end while
    bits &= 0 -- (since eg 101+101 -> 10000)
    return bits
end function

function to_bits2(integer a,b)
-- Apply to_bits() to a and b, and pad to the same length
    sequence sa = to_bits(a), sb = to_bits(b)
    integer diff = length(sa)-length(sb)
    if diff!=0 then
        if diff<0 then  sa &= repeat(0,-diff)
                  else  sb &= repeat(0,+diff)
        end if
    end if
    return {sa,sb}
end function

function to_int(sequence bits)
-- Copy of bits_to_int(), but in reverse order (lsb last)
    atom val = 0, p = 1
    for i=length(bits) to 1 by -1 do
        if bits[i] then
            val += p
        end if
        p += p
    end for
    return val
end function

function zstr(object z)
    if sequence(z) then
        for i=1 to length(z) do
            z[i] = zstr(z[i])
        end for
        return z
    end if
    return sprintf("%b",z)
end function

function rep(sequence res, integer ds, sequence was, wth)
-- helper for cleanup, validates replacements
    integer de = ds+length(was)-1
    if res[ds..de]!=was then ?9/0 end if
    if length(was)!=length(wth) then ?9/0 end if
    res[ds..de] = wth
    return res
end function

function zcleanup(sequence res)
-- (shared by zadd and zsub)
    integer l = length(res)
    -- first stage, left to right, {020x -> 100x', 030x -> 110x', 021x->110x, 012x->101x}
    for i=1 to l-3 do
        switch res[i..i+2]
            case {0,2,0}:   res[i..i+2] = {1,0,0}   res[i+3] += 1
            case {0,3,0}:   res[i..i+2] = {1,1,0}   res[i+3] += 1
            case {0,2,1}:   res[i..i+2] = {1,1,0}
            case {0,1,2}:   res[i..i+2] = {1,0,1}
        end switch
    end for
    -- first stage cleanup
    if l>1 then
        if res[l-1]=3 then      res = rep(res,l-2,{0,3,0},{1,1,1})      -- 030 -> 111
        elsif res[l-1]=2 then
            if res[l-2]=0 then  res = rep(res,l-2,{0,2,0},{1,0,1})      -- 020 -> 101
                          else  res = rep(res,l-3,{0,1,2,0},{1,0,1,0})  -- 0120 -> 1010
            end if
        end if
    end if
    if res[l]=3 then            res = rep(res,l-1,{0,3},{1,1})          -- 03 -> 11
    elsif res[l]=2 then
        if res[l-1]=0 then      res = rep(res,l-1,{0,2},{1,0})          -- 02 -> 10
                      else      res = rep(res,l-2,{0,1,2},{1,0,1})      -- 012 -> 101
        end if
    end if
    -- second stage, pass 1, right to left, 011 -> 100
    for i=length(res)-2 to 1 by -1 do
        if res[i..i+2]={0,1,1} then res[i..i+2] = {1,0,0} end if
    end for
    -- second stage, pass 2, left to right, 011 -> 100
    for i=1 to length(res)-2 do
        if res[i..i+2]={0,1,1} then res[i..i+2] = {1,0,0} end if
    end for
    return to_int(res)
end function

function zadd(integer a, b)
    sequence {sa,sb} = to_bits2(a,b)
    return zcleanup(reverse(sq_add(sa,sb)))
end function

function zinc(integer a)
    return zadd(a,0b1)
end function

function zsub(integer a, b)
    sequence {sa,sb} = to_bits2(a,b)
    sequence res = reverse(sq_sub(sa,sb))
    -- (/not/ combined with the first pass of the add routine!)
    for i=1 to length(res)-2 do
        switch res[i..i+2] do
            case {1, 0, 0}: res[i..i+2] = {0,1,1}
            case {1,-1, 0}: res[i..i+2] = {0,0,1}
            case {1,-1, 1}: res[i..i+2] = {0,0,2}
            case {1, 0,-1}: res[i..i+2] = {0,1,0}
            case {2, 0, 0}: res[i..i+2] = {1,1,1}
            case {2,-1, 0}: res[i..i+2] = {1,0,1}
            case {2,-1, 1}: res[i..i+2] = {1,0,2}
            case {2, 0,-1}: res[i..i+2] = {1,1,0}
        end switch
    end for
    -- copied from PicoLisp: {1,-1} -> {0,1} and {2,-1} -> {1,1}
    for i=1 to length(res)-1 do
        switch res[i..i+1] do
            case {1,-1}: res[i..i+1] = {0,1}
            case {2,-1}: res[i..i+1] = {1,1}
        end switch
    end for
    if find(-1,res) then ?9/0 end if -- sanity check
    return zcleanup(res)
end function

function zdec(integer a)
    return zsub(a,0b1)
end function

function zmul(integer a, b)
integer res = 0
    sequence mult = {a,zadd(a,a)}   -- (as per task desc)
    integer bits = 2
    while bits<b do
        mult = append(mult,zadd(mult[$],mult[$-1]))
        bits *= 2
    end while
    integer bit = 1
    while b do
        if and_bits(b,1) then
            res = zadd(res,mult[bit])
        end if
        b = floor(b/2)
        bit += 1
    end while
    return res
end function

function zdiv(integer a, b)
integer res = 0
    sequence mult = {b,zadd(b,b)}
    integer bits = 2
    while mult[$]<a do
        mult = append(mult,zadd(mult[$],mult[$-1]))
        bits *= 2
    end while
    for i=length(mult) to 1 by -1 do
        integer mi = mult[i]
        if mi<=a then
            res = zadd(res,bits)
            a = zsub(a,mi)
            if a=0 then exit end if
        end if
        bits = floor(bits/2)
    end for
    return {res,a} -- (a is the remainder)
end function

for i=0 to 20 do
    integer zi = zeckendorf(i)
    atom d = decimal(zi)
    printf(1,"%2d: %7b (%d)\n",{i,zi,d})
end for

procedure test(atom a, string op, atom b, object res, string expected)
    string zres = iff(atom(res)?zstr(res):join(zstr(res)," rem ")),
           dres = sprintf(iff(atom(res)?"%d":"%d rem %d"),decimal(res)),
           aka = sprintf("aka %d %s %d = %s",{decimal(a),op,decimal(b),dres}),
           ok = iff(zres=expected?"":" *** ERROR ***!!")
    printf(1,"%s %s %s = %s, %s %s\n",{zstr(a),op,zstr(b),zres,aka,ok})
end procedure

test(0b0,"+",0b0,zadd(0b0,0b0),"0")
test(0b101,"+",0b101,zadd(0b101,0b101),"10000")
test(0b10100,"-",0b1000,zsub(0b10100,0b1000),"1001")
test(0b100100,"-",0b1000,zsub(0b100100,0b1000),"10100")
test(0b1001,"*",0b101,zmul(0b1001,0b101),"1000100")
test(0b1000101,"/",0b101,zdiv(0b1000101,0b101),"1001 rem 1")

test(0b10,"+",0b10,zadd(0b10,0b10),"101")
test(0b101,"+",0b10,zadd(0b101,0b10),"1001")
test(0b1001,"+",0b1001,zadd(0b1001,0b1001),"10101")
test(0b10101,"+",0b1000,zadd(0b10101,0b1000),"100101")
test(0b100101,"+",0b10101,zadd(0b100101,0b10101),"1010000")
test(0b1000,"-",0b101,zsub(0b1000,0b101),"1")
test(0b10101010,"-",0b1010101,zsub(0b10101010,0b1010101),"1000000")
test(0b1001,"*",0b101,zmul(0b1001,0b101),"1000100")
test(0b101010,"+",0b101,zadd(0b101010,0b101),"1000100")

test(0b10100,"+",0b1010,zadd(0b10100,0b1010),"101000")
test(0b101000,"-",0b1010,zsub(0b101000,0b1010),"10100")

test(0b100010,"*",0b100101,zmul(0b100010,0b100101),"100001000001")
test(0b100001000001,"/",0b100,zdiv(0b100001000001,0b100),"101010001 rem 0")
test(0b101000101,"*",0b101001,zmul(0b101000101,0b101001),"101010000010101")
test(0b101010000010101,"/",0b100,zdiv(0b101010000010101,0b100),"1001010001001 rem 10")

test(0b10100010010100,"+",0b1001000001,zadd(0b10100010010100,0b1001000001),"100000000010101")
test(0b10100010010100,"-",0b1001000001,zsub(0b10100010010100,0b1001000001),"10010001000010")
test(0b10000,"*",0b1001000001,zmul(0b10000,0b1001000001),"10100010010100")
test(0b1010001010000001001,"/",0b100000000100000,zdiv(0b1010001010000001001,0b100000000100000),"10001 rem 10100001010101")

test(0b10100,"+",0b1010,zadd(0b10100,0b1010),"101000")
test(0b10100,"-",0b1010,zsub(0b10100,0b1010),"101")
test(0b10100,"*",0b1010,zmul(0b10100,0b1010),"101000001")
test(0b10100,"/",0b1010,zdiv(0b10100,0b1010),"1 rem 101")
integer m = zmul(0b10100,0b1010)
test(m,"/",0b1010,zdiv(m,0b1010),"10100 rem 0")