--- Basic arithmetic and calendar utilities.
-- Ported from "Calendrical Calculations" (4th edition)
-- by Nachum Dershowitz and Edward M. Reingold.
-- Original Lisp code (CALENDRICA 4.0) is Apache 2.0 licensed.
-- @module calendrica-basic
-- @release 0.1 2026-07-19


local M = {}

local begin, end_, to_radix, positions_in_range


-- === Constants ===

-- M.TRUE = true  (not used)
-- M.FALSE = false  (not used)

local BOGUS = "bogus"

local SUNDAY = 0

local MONDAY = 1

local TUESDAY = 2

local WEDNESDAY = 3

local THURSDAY = 4

local FRIDAY = 5

local SATURDAY = 6

local JD_EPOCH   = -1721424.5
local MJD_EPOCH  = 678576
local UNIX_EPOCH = 719163


-- === Arithmetic helpers ===

--- Whole part of `m`/`n`.
-- @tparam number m Dividend.
-- @tparam number n Divisor (non-zero).
-- @treturn number Integer quotient (floor division).
function M.quotient(m, n)
  return math.floor(m / n)
end
local quotient = M.quotient

--- The value of (`x` mod `y`) with `y` instead of 0.
-- @tparam number x Dividend.
-- @tparam number y Modulus (non-zero).
-- @treturn number Result in range [1..y].
function M.amod(x, y)
  return y + (x % (-y))
end

--- The value of `x` shifted into the range [`a`..`b`). Returns `x` if `a` = `b`.
-- @tparam number x Value to shift.
-- @tparam number a Lower bound.
-- @tparam number b Upper bound.
-- @treturn number
function M.mod3(x, a, b)
  if a == b then
    return x
  else
    return a + ((x - a) % (b - a))
  end
end
local mod3 = M.mod3

--- Sum `expression(index)` for `index` = `initial` and successive integers,
-- as long as `condition(index)` holds.
-- @tparam func expression Function mapping index to a number.
-- @tparam number initial Starting index.
-- @tparam func condition Predicate; iteration continues while true.
-- @treturn number
function M.sum(expression, initial, condition)
  local result = 0
  local i = initial
  while condition(i) do
    result = result + expression(i)
    i = i + 1
  end
  return result
end
local sum = M.sum

-- Product of expression(index) for index = initial and successive integers,
-- as long as condition(index) holds. (local helper, not exported)
local function prod(expression, initial, condition)
  local result = 1
  local i = initial
  while condition(i) do
    result = result * expression(i)
    i = i + 1
  end
  return result
end

--- First integer greater or equal to `initial` such that `condition` holds.
-- @tparam number initial Starting value.
-- @tparam func condition Predicate.
-- @treturn number
function M.next(initial, condition)
  local i = initial
  while not condition(i) do
    i = i + 1
  end
  return i
end

--- Last integer greater or equal to `initial` such that `condition` holds.
-- @tparam number initial Starting value.
-- @tparam func condition Predicate; stops when false.
-- @treturn number
function M.final(initial, condition)
  local i = initial
  while condition(i) do
    i = i + 1
  end
  return i - 1
end

--- Sum of `body(i1,...,in)` for indices running simultaneously through parallel lists.
-- `lists` is an array of arrays all of the same length; `body` receives one element
-- from each list on each iteration.
-- @tparam table lists Array of equal-length arrays.
-- @tparam func body Function applied to one element from each list.
-- @treturn number
function M.sigma(lists, body)
    local result = 0
    local n = #lists[1]
    for j = 1, n do
        local args = {}
        for i = 1, #lists do
            args[i] = lists[i][j]
        end
        result = result + body(table.unpack(args))
    end
    return result
end

--- Use bisection to find the inverse of angular function `f` at `y` within interval `r`.
-- @tparam func f Angular function to invert.
-- @tparam number y Target value.
-- @tparam table r Interval {lo, hi}.
-- @treturn number
function M.invert_angular(f, y, r)
  local varepsilon = 1 / 100000
  local l = begin(r)
  local u = end_(r)
  while (u - l) >= varepsilon do
    local x = (l + u) / 2
    if (f(x) - y) % 360 < 180 then
      u = x
    else
      l = x
    end
  end
  return (l + u) / 2
end

--- Bisection search for `x` in [`lo`..`hi`] such that `end_condition` holds.
-- `test` determines when to go left.
-- @tparam number lo Lower bound.
-- @tparam number hi Upper bound.
-- @tparam func test Predicate; go left (upper = x) when true.
-- @tparam func end_condition Predicate on (l, u); stop when true.
-- @treturn number
function M.binary_search(lo, hi, test, end_condition)
  local l = lo
  local u = hi
  local x
  while true do
    x = (l + u) / 2
    if end_condition(l, u) then
      return x
    end
    if test(x) then
      u = x
    else
      l = x
    end
  end
end

--- Evaluate polynomial in `x` with coefficients `a` (from order 0 up).
-- @tparam number x Argument.
-- @tparam table a Array of coefficients, a[1] is the constant term.
-- @treturn number
function M.poly(x, a)
  if #a == 0 then
    return 0
  end
  local result = 0
  for i = #a, 1, -1 do
    result = a[i] + x * result
  end
  return result
end

--- Sign of `y`: returns -1, 0, or 1.
-- @tparam number y
-- @treturn number -1, 0, or 1.
function M.sign(y)
  if y < 0 then return -1
  elseif y > 0 then return 1
  else return 0
  end
end

--- Round `x` to the nearest integer using banker's rounding (half-to-even).
-- Matches Common Lisp's `round`, used in `persian_from_fixed`.
-- @tparam number x Value to round.
-- @treturn number Integer result.
function M.round_half_to_even(x)
  local f    = math.floor(x)
  local diff = x - f          -- Fractional part, always in [0, 1).
  if diff < 0.5 then
    return f                  -- Closer to floor.
  elseif diff > 0.5 then
    return f + 1              -- Closer to ceiling.
  else
    -- Exactly halfway: return whichever of floor/ceiling is even.
    return f % 2 == 0 and f or f + 1
  end
end

-- === RD (Rata Die) ===

--- Identity function for fixed dates/moments (Rata Die).
-- If internal timekeeping is shifted, change `epoch` inside this function.
-- @tparam number tee Moment or fixed date.
-- @treturn number
function M.rd(tee)
  local epoch = 0
  return tee - epoch
end
local rd = M.rd


-- === Day of week ===

--- Day of the week of fixed `date` (0 = Sunday .. 6 = Saturday).
-- @tparam number date Fixed date.
-- @treturn number Day-of-week residue class.
function M.day_of_week_from_fixed(date)
  return (date - rd(0) - SUNDAY) % 7
end


-- === Standard date accessors ===

--- Month field of `date` = {year, month, day}.
-- @tparam table date Standard date tuple.
-- @treturn number Month.
function M.standard_month(date)
  return date[2]
end

--- Day field of `date` = {year, month, day}.
-- @tparam table date Standard date tuple.
-- @treturn number Day.
function M.standard_day(date)
  return date[3]
end

--- Year field of `date` = {year, month, day}.
-- @tparam table date Standard date tuple.
-- @treturn number Year.
function M.standard_year(date)
  return date[1]
end


-- === Clock time ===

--- TYPE (hour minute second) -> clock-time
--- Construct a clock time from hours, minutes, and seconds.
--- @param hour integer
--- @param minute integer
--- @param second real
--- @return clock-time
local function time_of_day(hour, minute, second) -- luacheck: ignore
  return {hour, minute, second}
end

--- TYPE clock-time -> hour
--- Hour field of `clock`.
--- @param clock clock-time
--- @return integer
local function hour(clock) -- luacheck: ignore
  return clock[1]
end

--- TYPE clock-time -> minute
--- Minute field of `clock`.
--- @param clock clock-time
--- @return integer
local function minute(clock) -- luacheck: ignore
  return clock[2]
end

--- TYPE clock-time -> second
--- Seconds field of `clock`.
--- @param clock clock-time
--- @return real
local function seconds(clock) -- luacheck: ignore
  return clock[3]
end


-- === Moment / time conversions ===

--- Fixed date (floor) from moment `tee`.
-- @tparam number tee Moment.
-- @treturn number Fixed date.
function M.fixed_from_moment(tee)
  return math.floor(tee)
end
local fixed_from_moment = M.fixed_from_moment

--- Fractional time-of-day from moment `tee`.
-- @tparam number tee Moment.
-- @treturn number Time in [0, 1).
function M.time_from_moment(tee)
  return tee % 1
end

--- TYPE (list-of-reals list-of-rationals list-of-rationals) -> real
--- The number corresponding to `a` in radix notation
--- with base `b` for whole part and `c` for fraction.
--- @param a real[]
--- @param b real[]
--- @param c real[]?
--- @return real
local function from_radix(a, b, c)
  c = c or {}
  local bc = {}
  for _, v in ipairs(b) do bc[#bc+1] = v end
  for _, v in ipairs(c) do bc[#bc+1] = v end
  local total = sum(
    function(i)
      return a[i+1] * prod(
        function(j) return bc[j+1] end,
        i, function(j) return j < #bc end
      )
    end,
    0, function(i) return i < #a end
  )
  local denom = prod(function(j) return c[j+1] end,
                       0, function(j) return j < #c end)
  if denom == 0 then return total end
  return total / denom
end

--- Convert `x` to mixed-radix notation with bases `b` (whole) and `c` (fraction).
-- @tparam number x Value to convert.
-- @tparam table b Array of bases for the integer part.
-- @tparam[opt] table c Array of bases for the fractional part.
-- @treturn table Array of digits.
function M.to_radix(x, b, c)
  c = c or {}
  if #c == 0 then
    if #b == 0 then
      return {x}
    else
      local last_b = b[#b]
      local prefix = {table.unpack(b, 1, #b - 1)}
      local head = to_radix(quotient(x, last_b), prefix)
      local tail = {x % last_b}
      for _, v in ipairs(tail) do head[#head+1] = v end
      return head
    end
  else
    local denom = prod(function(j) return c[j+1] end,
                         0, function(j) return j < #c end)
    local bc = {}
    for _, v in ipairs(b) do bc[#bc+1] = v end
    for _, v in ipairs(c) do bc[#bc+1] = v end
    return to_radix(x * denom, bc)
  end
end
to_radix = M.to_radix

--- TYPE moment -> clock-time
--- Clock time hour:minute:second from moment `tee`.
--- @param tee moment
--- @return clock-time
local function clock_from_moment(tee)
  local r = to_radix(tee, {}, {24, 60, 60})
  return {r[2], r[3], r[4]}  -- drop the day part (first element)
end

--- TYPE clock-time -> time
--- Time of day from `hms` = hour:minute:second.
--- @param hms clock-time
--- @return real
local function time_from_clock(hms)
  return from_radix(hms, {}, {24, 60, 60}) / 24
end


-- === Angles ===

--- TYPE (degree minute real) -> angle
--- Construct an angle from degrees, arcminutes, and arcseconds.
--- @param d integer
--- @param m integer
--- @param s real
--- @return angle
local function degrees_minutes_seconds(d, m, s) -- luacheck: ignore
  return {d, m, s}
end

--- TYPE angle -> list-of-reals
--- List of degrees-arcminutes-arcseconds from angle `alpha` in degrees.
--- @param alpha real
--- @return angle
local function angle_from_degrees(alpha) -- luacheck: ignore
  local dms = to_radix(math.abs(alpha), {}, {60, 60})
  if alpha >= 0 then
    return dms
  else
    return {-dms[1], -dms[2], -dms[3]}
  end
end


-- === Julian Day Numbers ===

--- TYPE julian-day-number -> moment
--- Moment of Julian Day Number `jd`.
--- @param jd real
--- @return moment
local function moment_from_jd(jd)
  return jd + JD_EPOCH
end

local function jd_from_moment(tee)
  return tee - JD_EPOCH
end

local function fixed_from_jd(jd)
  return math.floor(moment_from_jd(jd))
end

local function jd_from_fixed(date)
  return jd_from_moment(date)
end


-- === Modified Julian Day Numbers ===

--- TYPE julian-day-number -> fixed-date
--- Fixed date of Modified Julian Day Number `mjd`.
--- @param mjd real
--- @return fixed-date
local function fixed_from_mjd(mjd)
  return mjd + MJD_EPOCH
end

local function mjd_from_fixed(date)
  return date - MJD_EPOCH
end


-- === Unix time ===

--- TYPE second -> moment
--- Fixed date from Unix second count `s`.
--- @param s integer
--- @return moment
local function moment_from_unix(s)
  return UNIX_EPOCH + (s / 24 / 60 / 60)
end

local function unix_from_moment(tee)
  return 24 * 60 * 60 * (tee - UNIX_EPOCH)
end


-- === Intervals ===

--- Half-open interval [`t0`..`t1`).
-- @tparam number t0 Start moment.
-- @tparam number t1 End moment (excluded).
-- @treturn table
function M.interval(t0, t1)
  return {t0, t1}
end
local interval = M.interval

--- Closed interval [`t0`..`t1`].
-- @tparam number t0 Start moment.
-- @tparam number t1 End moment (included).
-- @treturn table
function M.interval_closed(t0, t1)
  return {t0, t1}
end

--- Start of interval `range`.
-- @tparam table range Interval {t0, t1}.
-- @treturn number t0.
function M.begin(range)
  return range[1]
end
begin = M.begin

--- End of interval `range`. Named `end_` to avoid the Lua keyword `end`.
-- @tparam table range Interval {t0, t1}.
-- @treturn number t1.
function M.end_(range)
  return range[2]
end
end_ = M.end_

--- True if `tee` is in half-open `range` [t0..t1).
-- @tparam number tee Moment to test.
-- @tparam table range Interval {t0, t1}.
-- @treturn boolean
function M.in_range(tee, range)
  return begin(range) <= tee and tee < end_(range)
end
local in_range = M.in_range

--- Those moments in list `ell` that occur in `range`.
-- @tparam table ell List of moments.
-- @tparam table range Interval {t0, t1}.
-- @treturn table Filtered list.
function M.list_range(ell, range)
  local result = {}
  for _, tee in ipairs(ell) do
    if in_range(tee, range) then
      result[#result + 1] = tee
    end
  end
  return result
end

-- List of fixed dates from a list of moments. (local helper)
local function list_of_fixed_from_moments(ell)
  local result = {}
  for _, tee in ipairs(ell) do
    result[#result + 1] = fixed_from_moment(tee)
  end
  return result
end

--- List of occurrences of moment `p` of `c`-day cycle within `range`.
-- `cap_delta` is the position in the cycle of RD moment 0.
-- @tparam number p Position within the cycle.
-- @tparam number c Cycle length in days.
-- @tparam number cap_delta Cycle offset at RD 0.
-- @tparam table range Interval {t0, t1}.
-- @treturn table List of moments.
function M.positions_in_range(p, c, cap_delta, range)
  local a = begin(range)
  local b = end_(range)
  local date = mod3(p - cap_delta, a, a + c)
  if date >= b then
    return {}
  end
  local result = {date}
  local rest = positions_in_range(p, c, cap_delta, interval(a + c, b))
  for _, v in ipairs(rest) do
    result[#result + 1] = v
  end
  return result
end
positions_in_range = M.positions_in_range


-- === JD / MJD / Unix conversions (public) ===

--- Moment from Julian Day Number `jd`.
-- @tparam number jd Julian Day Number.
-- @treturn number Moment.
function M.moment_from_jd(jd)
  return moment_from_jd(jd)
end

--- Julian Day Number from moment `tee`.
-- @tparam number tee Moment.
-- @treturn number Julian Day Number.
function M.jd_from_moment(tee)
  return jd_from_moment(tee)
end

--- Fixed date from Julian Day Number `jd`.
-- @tparam number jd Julian Day Number.
-- @treturn number Fixed date.
function M.fixed_from_jd(jd)
  return fixed_from_jd(jd)
end

--- Julian Day Number from fixed `date`.
-- @tparam number date Fixed date.
-- @treturn number Julian Day Number.
function M.jd_from_fixed(date)
  return jd_from_fixed(date)
end

--- Fixed date from Modified Julian Day Number `mjd`.
-- @tparam number mjd Modified Julian Day Number.
-- @treturn number Fixed date.
function M.fixed_from_mjd(mjd)
  return fixed_from_mjd(mjd)
end

--- Modified Julian Day Number from fixed `date`.
-- @tparam number date Fixed date.
-- @treturn number Modified Julian Day Number.
function M.mjd_from_fixed(date)
  return mjd_from_fixed(date)
end

--- Unix timestamp (seconds) from fixed `date`.
-- @tparam number date Fixed date.
-- @treturn number Unix timestamp.
function M.unix_from_fixed(date)
  return unix_from_moment(date)
end

--- Fixed date from Unix timestamp `s` (seconds since 1970-01-01).
-- @tparam number s Unix timestamp.
-- @treturn number Fixed date.
function M.fixed_from_unix(s)
  return math.floor(moment_from_unix(s))
end


-- === Exports ===

--- Sentinel value used to denote nonexistent dates.
M.BOGUS                  = BOGUS
--- Residue class for Friday.
M.FRIDAY                 = FRIDAY
--- Residue class for Monday.
M.MONDAY                 = MONDAY
--- Residue class for Saturday.
M.SATURDAY               = SATURDAY
--- Residue class for Sunday.
M.SUNDAY                 = SUNDAY
--- Residue class for Thursday.
M.THURSDAY               = THURSDAY
--- Residue class for Tuesday.
M.TUESDAY                = TUESDAY
--- Residue class for Wednesday.
M.WEDNESDAY              = WEDNESDAY

--- Product of `expression(i)` for `i` from `initial` while `condition(i)` holds.
-- @function prod
-- @tparam func   expression Function of i.
-- @tparam number initial    Starting value.
-- @tparam func   condition  Predicate.
-- @treturn number
M.prod = prod

--- Convert mixed-radix digits `a` with bases `b` (whole) and `c` (fraction) to a number.
-- Inverse of `to_radix`.
-- @function from_radix
-- @tparam table a Digits.
-- @tparam table b Bases for integer part.
-- @tparam[opt] table c Bases for fractional part.
-- @treturn number
M.from_radix = from_radix

--- Clock time {hour, minute, second} from moment `tee`.
-- @function clock_from_moment
-- @tparam number tee Moment.
-- @treturn table {hour, minute, second}
M.clock_from_moment = clock_from_moment

--- Fractional day from clock time `hms` = {hour, minute, second}.
-- @function time_from_clock
-- @tparam table hms Clock time {hour, minute, second}.
-- @treturn number Fractional day.
M.time_from_clock = time_from_clock

--- Moment from Unix timestamp `s` (seconds since 1970-01-01).
-- @function moment_from_unix
-- @tparam number s Unix timestamp.
-- @treturn number Moment.
M.moment_from_unix = moment_from_unix

--- Unix timestamp (seconds) from moment `tee`.
-- @tparam number tee Moment (fractional fixed date).
-- @treturn number Unix timestamp.
function M.unix_from_moment(tee)
  return unix_from_moment(tee)
end

--- List of fixed dates from a list of moments (floors each moment).
-- @function list_of_fixed_from_moments
-- @tparam table ell List of moments.
-- @treturn table List of fixed dates.
M.list_of_fixed_from_moments = list_of_fixed_from_moments

return M