map.md

map(a, b, x, y, t) maps [a, b] -> [x, y] wrt t

in a map function, most importantly we want exactness:

• map(a, b, x, y, a) = x
• map(a, b, x, y, b) = y

another nice guarantee to have is consistency:

• map(a, b, x, x, t) = x

because many algorithms assume going only forward or backward in time, monotonicity:

• d/dt map(a, b, x, y, t) is strictly >= 0 or <= 0

finally, if a or b and t are close to 0, we want to maintain full precision. For example:

• map(1, 0, 1, 0, 1e-300) = 1e-300

---

first, consider the lerp function

lerp(x, y, t) maps [0, 1] -> [x, y] wrt t

the following lerp algorithm has all the properties above:

function lerp(x, y, t)
    if t < 1/2 then
        return x + (y - x)*t
    else
        return y + (y - x)*(t - 1)
    end
end

guarantees exactness because:

• x + 0*(y - x) = x
• y + (1 - 1)*(y - x) = y + 0*(y - x) = y

guarantees consistency because:

• a + (a - a)*t = a + 0*t = a

The argument for monotonicity is harder. The following is my attempt to show the worst case by simple floating point error terms, however it fails to take into account that these calculations are not independent from each other. In reality, testing billions of pairs of randomly generated numbers has yielded no case where lerp(x, y, 1/2 - 2^-54) > lerp(x, y, 1/2)

under floating point, any expression a + b*t will be monotonic wrt t.
this means that the only time we need to worry about monotonicity is at the exchange point around t = 1/2.
specifically, we want to show that this is always true:
    x + (y - x)*(1/2 - 2^-54) <= y + (y - x)*(-1/2)

y + (y - x)/2:
    y - x worst case rounds to         (y - x)*(1 + 2^-54)
    (y - x)/2 worst case rounds to     (y - x)/2*(1 + 2^-54)
    y - (y - x)/2 worst case rounds to (y - (y - x)/2*(1 + 2^-54))*(1 - 2^-54)
    after expansion, this becomes      (x + y)/2 - 2^-54*y - 2^-109*x + 2^-109*y/2
    dropping inconsequential terms     (x + y)/2 - 2^-54*y

x + (y - x)*(1/2 - 2^-54):
    y - x worst case rounds to                     (y - x)*(1 + 2^-54)
    multiplication by (1/2 - 2^-54) guarantees division by two
        and the subtraction of the last bit with ulp^2 (none) error
    (y - x)*(1/2 - 2^-54) worst case rounds to     (y - x)*(1/2 - 2^-54)*(1 + 2^-54)
    x + (y - x)*(1/2 - 2^-54) worse case rounds to (x + (y - x)*(1/2 - 2^-54)*(1 + 2^-54))*(1 + 2^-54)
    after expansion, this becomes                  (x + y)/2 + 2^-54*x + 3*2^-109*(x - y) + 2^-162*(x - y)
    dropping inconsequential terms                 (x + y)/2 + 2^-54*x
    
the above does not take into account the fact that these caluculations are dependent on choice of x and y.
here is one such case:
    
when y/2 < x < y
    (y - x)/2     has no error
    y - (y - x)/2 also has no error
    x + (y - x)/2 also has no error

the following correctness test never seems to fail:

for i = 1, 10000000 do
    local x = math.random()*2^math.random(-1023, 1024)*(math.random() < 1/2 and -1 or 1)
    local y = math.random()*2^math.random(-1023, 1024)*(math.random() < 1/2 and -1 or 1)
    if y < x then x, y = y, x end

    local t0 = 1/2 - 2^-54
    local t1 = 1/2

    local m0 = lerp(x, y, t0)
    local m1 = lerp(x, y, t1)
    if m0 > m1 then
        print("interpolation is not monotonic")
        print("x, y", x, y)
        print("t:", t0, t1)
        print("m:", m0, m1)
        return
    end
end

---

the following map algorithm has all the properties above:

local function map(a, b, x, y, t)
    local s = (t - a)/(b - a)
    local r = (b - t)/(b - a)
    if s < 1/2 then
        return x + s*(y - x)
    elseif r < 1/2 then
        return y + r*(x - y)
    else
        return (x + y)/2
    end
end

guarantees exactness because

• lerp guarantees exactness
• (a - a)/(b - a) = 0
• (b - b)/(b - a) = 0

guarantees consistency because

• lerp guarantees consistency

guarantees monotonicity because

• lerp guarantees monotonicity
• (t - a)/(b - a) is monotonic
• we take care to make sure the interchange point is when s, i == 1/2 - 2^-54
• and give an intermediate value, (x + y)/2, which will always lie equal to or between the previous branches.

In the end, the following test for correctness never seems to fail:

for i = 1, 10000000 do
    local x = math.random()*2^math.random(-1023, 1024)*(math.random() < 1/2 and -1 or 1)
    local y = math.random()*2^math.random(-1023, 1024)*(math.random() < 1/2 and -1 or 1)
    local a = math.random()*2^math.random(-1023, 1024)*(math.random() < 1/2 and -1 or 1)
    local b = math.random()*2^math.random(-1023, 1024)*(math.random() < 1/2 and -1 or 1)
    if y < x then x, y = y, x end
    if b < a then a, b = b, a end

    local t0, t1, t2, t3, t4, t5
    if a + b > 0 then
        t0 = (a + b)/2*(1 - 2*2^-53)
        t1 = (a + b)/2*(1 - 1*2^-53)
        t2 = (a + b)/2*(1 - 0      )
        t3 = (a + b)/2*(1 + 1*2^-52)
        t4 = (a + b)/2*(1 + 2*2^-52)
    else
        t4 = (a + b)/2*(1 - 2*2^-53)
        t3 = (a + b)/2*(1 - 1*2^-53)
        t2 = (a + b)/2*(1 - 0      )
        t1 = (a + b)/2*(1 + 1*2^-52)
        t0 = (a + b)/2*(1 + 2*2^-52)
    end

    local m0 = map(a, b, x, y, t0)
    local m1 = map(a, b, x, y, t1)
    local m2 = map(a, b, x, y, t2)
    local m3 = map(a, b, x, y, t3)
    local m4 = map(a, b, x, y, t4)
    if m0 > m1 or m1 > m2 or m2 > m3 or m3 > m4 then
        print("interpolation is not monotonic")
        print("a, b, x, y", a, b, x, y)
        print("t:", t0, t1, t2, t3, t4)
        print("m:", m0, m1, m2, m3, m4)
        print("cases:", m0 > m1, m1 > m2, m2 > m3, m3 > m4)
        return
    end
end