rem.py

import numpy as np

from math import sin, pi
from scipy.integrate import quad


def remainder_map(a, b, xs, xe):
    period = xe - xs
    interval = b - a
    n_periods, left = divmod(interval, period)

    a = xs + (a - xs) % period
    b = a + left

    if b <= xe:
        ranges = [(a, b)]
    else:
        ranges = [(xs, xs + left + a - xe), (a, xe)]

    return n_periods, ranges


def doit(a, b):
    rper = 1/3.
    
    def f(x):
        return sin(pi*x/rper)**2

    ig0 = quad(f, a, b)[0]

    n_periods, ranges = remainder_map(a, b, 0, rper)
    print(a, b, n_periods, ranges)

    ig1 = quad(f, 0, rper)[0] * n_periods
    for xa, xb in ranges:
        ig1 += quad(f, xa, xb)[0]

    msg = repr((a, b, n_periods, ranges))        
    np.testing.assert_allclose(ig0, ig1, err_msg=msg)


def main():
    # Check all fp numbers around b=1
    b = 1
    for k in range(100):
        b = np.nextafter(b, 0)

    for k in range(200):
        b = np.nextafter(b, np.inf)
        doit(0, b)


if __name__ == "__main__":
    main()