daynames: ["Sunday", "Monday", "Tuesday",
    "Wednesday", "Thursday", "Friday", "Saturday"]$

isleapyear(y) := mod(y, 4)=0 and (mod(y, 100)#0 or mod(y, 400)=0)$
anchors(y) := [if isleapyear(y) then 4 else 3,
    if isleapyear(y) then 1 else 7,
    7, 4, 2, 6, 4, 1, 5, 3, 7, 5]$

doomsday(y) := mod(2+5*mod(y, 4)+4*mod(y, 100)+6*mod(y, 400), 7)$
wd(y, m, d) := mod(doomsday(y)+d-anchors(y)[m], 7)$

test_dates: [[1800, 1, 6], [1875, 3, 29], [1915, 12, 7], [1970, 12, 23], [2043, 5, 14], [2077, 2, 12], [2101, 4, 2]]$

raw_wds: map(lambda([l], apply(wd, l)), test_dates)$
for i in raw_wds do(sprint(daynames[i+1]))$