module LinearMultiStep {
    record empty {}
    record cons {
        param weight;
        type tail;
    }

    proc initial(type x : empty) param return 0;
    proc initial(type x : cons(?w, ?t)) param return 1 + initial(t);
    proc cff(param x : int, type ct : cons(?w, ?t)) param {
        if x == 1 {
            return w;
        } else {
            return cff(x-1, t);
        }
    }

    proc runMethod(type method, step : real, count : int, start : real,
                   n : real ... initial(method)): real {
        param coeffCount = initial(method);
        // Repeat the methods as many times as requested
        for i in 1..count {
            // We're computing by adding h*b_j*f(...) to y_n.
            // Set total to y_n.
            var total = n(coeffCount - 1);
            for param j in 1..coeffCount do
                // For each coefficient b_j given by cff(j, method)
                // increment the total by h*bj*f(...)
                total += step * cff(j, method) *
                    f(start + step*(i-1+coeffCount-j), n(coeffCount-j));
            // Shift each y_i over by one, and set y_{n+s} to the
            // newly computed total.
            for j in 0..< coeffCount - 1 do
                n(j) = n(j+1);
            n(coeffCount - 1) = total;
        }
        // return final y_{n+s}
        return n(coeffCount - 1);
    }
}

proc f(t: real, y: real) return y;

use LinearMultiStep;

type euler = cons(1.0, empty);
type adamsBashforth = cons(3.0/2.0, cons(-0.5, empty));
type someThirdMethod = cons(23.0/12.0, cons(-16.0/12.0, cons(5.0/12.0, empty)));

// prints 5.0625 (correct)
writeln(runMethod(euler, step=0.5, count=4, start=0, 1));

// For Adams-Bashforth, pick second initial point from Euler's method
// returns 6.0234 (correct)
writeln(runMethod(adamsBashforth, step=0.5, count=3, start=0, 1,
    runMethod(euler, step=0.5, count=1, start=0, 1)));

writeln(runMethod(someThirdMethod, step=0.5, count=2, start=0,
    1,
    runMethod(euler, step=0.5, count=1, start=0, 1),
    runMethod(adamsBashforth, step=0.5, count=1, start=0, 1, runMethod(euler, step=0.5, count=1, start=0, 1))));