def breakpoints(next_horizon, first_horizon, horizon, debug=False):
    """ Return the list of restart point (breakpoints), if starting from ``first_horizon`` to ``horizon`` with growth function ``next_horizon``.

    - Also return the gap between the last guess for horizon and the true horizon. This gap should not be too large.
    - Nicely print all the values if ``debug=True``.

    - First examples:

    >>> first_horizon = 1000
    >>> horizon = 30000
    >>> breakpoints(next_horizon__arithmetic, first_horizon, horizon)  # doctest: +ELLIPSIS
    ([1000, 1200, 1400, ..., 29800, 30000], 0)
    >>> breakpoints(next_horizon__geometric, first_horizon, horizon)
    ([1000, 2000, 4000, 8000, 16000, 32000], 2000)
    >>> breakpoints(next_horizon__exponential, first_horizon, horizon)
    ([1000, 31622], 1622)
    >>> breakpoints(next_horizon__exponential_slow, first_horizon, horizon)
    ([1000, 1995, 4265, 9838, 24671, 67827], 37827)
    >>> breakpoints(next_horizon__exponential_fast, first_horizon, horizon)
    ([1000, 1000000], 970000)

    - Second examples:

    >>> first_horizon = 5000
    >>> horizon = 1000000
    >>> breakpoints(next_horizon__arithmetic, first_horizon, horizon)  # doctest: +ELLIPSIS
    ([5000, 5200, ..., 999600, 999800, 1000000], 0)
    >>> breakpoints(next_horizon__geometric, first_horizon, horizon)
    ([5000, 10000, 20000, 40000, 80000, 160000, 320000, 640000, 1280000], 280000)
    >>> breakpoints(next_horizon__exponential, first_horizon, horizon)
    ([5000, 353553, 210223755], 209223755)
    >>> breakpoints(next_horizon__exponential_slow, first_horizon, horizon)
    ([5000, 11718, 29904, 83811, 260394, 906137, 3572014], 2572014)
    >>> breakpoints(next_horizon__exponential_fast, first_horizon, horizon)
    ([5000, 25000000], 24000000)

    - Third examples:

    >>> first_horizon = 10
    >>> horizon = 1123456
    >>> breakpoints(next_horizon__arithmetic, first_horizon, horizon)  # doctest: +ELLIPSIS
    ([10, 210, 410, ..., 1123210, 1123410, 1123610], 154)
    >>> breakpoints(next_horizon__geometric, first_horizon, horizon)
    ([10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480, 40960, 81920, 163840, 327680, 655360, 1310720], 187264)
    >>> breakpoints(next_horizon__exponential, first_horizon, horizon)
    ([10, 31, 172, 2255, 107082, 35040856], 33917400)
    >>> breakpoints(next_horizon__exponential_slow, first_horizon, horizon)
    ([10, 12, 15, 19, 25, 34, 48, 70, 107, 170, 284, 499, 928, 1837, 3895, 8903, 22104, 60106, 180638, 606024, 2294768], 1171312)
    >>> breakpoints(next_horizon__exponential_fast, first_horizon, horizon)
    ([10, 100, 10000, 100000000], 98876544)
    """
    i = 0
    t = max(first_horizon, 2)
    times = [t]
    if debug:
        print(
            """

For the growth function {}, named '{}', first guess of the horizon = {} and true horizon = {} ...
 ==> The times will be:"""
            .format(next_horizon, getattr(next_horizon, '__latex_name__',
            '?'), first_horizon, horizon))
    while t < horizon:
        t = next_horizon(i, t)
        i += 1
        times.append(t)
        if debug:
            print('    The {}th breakpoint is {} ...'.format(i, t))
    assert horizon <= t, 'Error: the last guess for horizon = {} was found smaller than the true horizon = {}...'.format(
        t, horizon)
    gap = t - horizon
    if debug:
        print(
            'This last guess for horizon = {} gives a gap = {} against the true horizon {}. Relative difference = {:.3%}...'
            .format(t, gap, horizon, gap / float(horizon)))
    return times, gap
