def stream_tail(stream):
    return tail(stream)()
def stream_ref(s, n):
    return (head(s)
            if n == 0
            else stream_ref(stream_tail(s), n - 1))

def stream_map(f, s):
    return (None
            if is_none(s)
            else pair(f(head(s)),
                      lambda: stream_map(f, stream_tail(s))))

def stream_for_each(fun, s):
    if is_none(s):
        return True
    else:
        fun(head(s))
        return stream_for_each(fun, stream_tail(s))
def memo(fun):
    already_run = False
    result = None
    def memoized():
        nonlocal already_run, result
        if not already_run:
            result = fun()
            already_run = True
            return result
        else:
            return result
    return memoized

# note the use of the memoization optimization
def stream_combine(f, s1, s2):
    return (None
            if is_none(s1) and is_none(s2)
            else error("unexpected argument -- stream_combine", None)
            if is_none(s1) or is_none(s2)
            else pair(f(head(s1), head(s2)),
                      memo(lambda: stream_combine(f, stream_tail(s1),
                                                  stream_tail(s2)))))

def add_streams(s1, s2):
    return stream_combine(lambda x1, x2: x1 + x2, s1, s2)

def scale_stream(stream, factor):
    return stream_map(lambda x: x * factor,
                      stream)

# note the use of the memoization optimization
def integral(delayed_integrand, initial_value, dt):
    integ = pair(initial_value,
                 memo(lambda: add_streams(
                                  scale_stream(delayed_integrand(), dt),
                                  integ)))
    return integ

def solve(f, y0, dt):
    y = integral(lambda: dy, y0, dt)
    dy = stream_map(f, y)
    return y


print(stream_ref(solve(lambda y: y, 1, 0.001), 1000))
