vardef base(expr Y, M, d) = 
  save m, y; 
  if M < 3: m = M + 9; y = Y - 1;
  else: m = M - 3; y = Y; fi
  365/1024 y + (floor(y/4) - floor(y/100) + floor(y/400) 
    + floor((2+3m)/5) + 30m + d - 307) / 1024
enddef;