Minimise rounding error for Real input to Periodic

[bencottier]

For TimeType data, Periodic transform uses the _periodic function to compute the periodic function for each instant in time.

The purpose of the _periodic function is to find an "origin" in time that minimises rounding error. This is particularly important when working with dates, because the calculations are done on milliseconds, and the number of milliseconds since "year 0" is enormous.

julia> sin(2π)
-2.4492935982947064e-16

julia> sin(2π*100)
3.928773447456944e-15

julia> x = Millisecond(DateTime(2021, 02, 09) - DateTime(0)).value
63780048000000

julia> sin(2π*x)
-0.020117088546809787

So _periodic finds the most recent time, period_begin, such that period_begin / period is a whole number, and in turn sin(2π * period_begin / period) ≈ 0 with the highest precision.

But this consideration could just as well apply to Real numbers, albeit less often in practice. So we should accomodate Real input in _periodic or define a new equivalent method, e.g.

function _periodic(f, instant, period, phase_shift=0)
    return f(2π * (mod(instant, period) - phase_shift) / period)
end

Examples:

julia> _periodic(sin, 7/4, 7.)
1.0

julia> _periodic(sin, 0, 7.)
0.0

julia> _periodic(sin, 7e14, 7.)
0.0

julia> _periodic(sin, 7e14 + 7/4, 7.)
1.0

[bencottier]

Note that even in my example implementation, performance degrades for large enough input:

julia> _periodic(sin, 7e16 + 7/4, 7)
0.0

This might be a more fundamental issue that we can't reasonably avoid.