#ifndef CPP2_CPP2TAYLOR_H
#define CPP2_CPP2TAYLOR_H

cpp2: namespace = {

taylor: <R, dim: int> type  = {
    public v : std::array<R, dim> = ();

    operator=:(out this) = {}
    operator=:(out this, d1: R) = {
        v = ();
        v[0] = d1;
    }
    operator=:(out this, that) = {}

    operator=:(out this, l: std::initializer_list<R>) = {
        (copy i := 1)
        for l do (cur) {
            set(i, cur);
        }
    }

    // C++ interface

    operator[]: (this, k: int) -> R = {
        assert(1 <= k <= dim);
        r: R = v[k - 1];
        (copy i := 2)
        while i <= k next i += 1 {
            r *= i;
        }
        return r;
    }

    set: (inout this, k: int, value: R) = {
        assert(1 <= k <= dim);
        v[k - 1] = value;

        (copy i := 2)
        while i <= k next i += 1 {
            v[k - 1] /= i;
        }
    }

    // C++2 interface / AD interface

    get: (this, i: int, v0: R) -> R = {
        if i == 0 {
            return v0;
        }
        else if i > dim {
            return 0.0;
        }

        return v[i - 1];
    }

    // Overload for reverse AD.
    operator+=: (inout this, o: taylor) -> forward_ref _ = {
        this = this + o;
        return this;
    }

    // Overload for reverse AD.
    operator-=: (inout this, o: taylor) -> forward_ref _ = {
        this = this - o;
        return this;
    }

    // Overload for simple handling of connected adds.
    operator+: (this, o: taylor) -> taylor = {
        return add(o, 0.0, 0.0); // Primal values are not required.
    }

    // Overload for simple handling of connected minuses.
    operator-: (this, o: taylor) -> taylor = {
        return sub(o, 0.0, 0.0); // Primal values are not required.
    }

    // Overload for simple handling of prefix +.
    operator+: (this) -> taylor = {
        return this;
    }

    // Overload for simple handling of prefix -.
    operator-: (this) -> taylor = {
        r: taylor = ();

        (copy k:= 1)
        while k <= dim next k += 1 {
            r.v[k - 1] = -v[k - 1];
        }

        return r;
    }

    add: (this, o: taylor, v0: R, o0: R) -> taylor = {
        r: taylor = ();

        (copy k:= 1)
        while k <= dim next k += 1 {
            r.v[k - 1] = get(k, v0) +  o.get(k, o0);
        }

        return r;
    }

    sub: (this, o: taylor, v0: R, o0: R) -> taylor = {
        r: taylor = ();

        (copy k:= 1)
        while k <= dim next k += 1 {
            r.v[k - 1] = get(k, v0) -  o.get(k, o0);
        }

        return r;
    }

    mul: <dim_o: int> (this, in_ref o: taylor<R, dim_o>, v0: R, o0: R) -> taylor<R, std::max(dim, dim_o)> = {
        dim_r : int == std::max(dim, dim_o);
        r: taylor<R, dim_r> = ();

        (copy k:= 1)
        while k <= dim_r next k += 1 {
            (copy j :=  0)
            while j <= k next j += 1 {
                r..v[k - 1] += get(j, v0) * o..get(k - j, o0);
            }            
        }
        return r;
    }

    div: (this, o: taylor, v0: R, o0: R) -> taylor = {
        r: taylor = ();
        r0: R = v0 / o0;

        factor : R = 1.0 / o0;

        (copy k:= 1)
        while k <= dim next k += 1 {
            r..v[k - 1] = get(k, v0);
            (copy j :=  0)
            while j < k next j += 1 {
                r..v[k - 1] -= r.get(j, r0) * o..get(k - j, o0);
            }            
            r..v[k - 1] *= factor;
        }
        return r;
    }

    sqrt: (this, v0: R) -> taylor = {
        r: taylor = ();
        r0: R = std::sqrt(v0);

        factor : R = 0.5 / r0;

        (copy k:= 1)
        while k <= dim next k += 1 {
            r..v[k - 1] = get(k, v0);
            (copy j :=  1)
            while j < k next j += 1 {
                r..v[k - 1] -= r..get(j, r0) * r..get(k - j, r0);
            }            
            r..v[k - 1] *= factor;
        }

        return r;
    }

    log: (this, v0: R) -> taylor = {
        r: taylor = ();
        r0: R = std::log(v0);

        factor : R = 1.0 / v0;

        (copy k:= 1)
        while k <= dim next k += 1 {
            r..v[k - 1] =k * get(k, v0);
            (copy j :=  1)
            while j < k next j += 1 {
                r..v[k - 1] -= j * get(k - j, v0) * r..get(j, r0);
            }            
            r..v[k - 1] *= factor / k;
        }

        return r;
    }

    exp: (this, v0: R) -> taylor = {
        r: taylor = ();
        r0: R = std::exp(v0);

        factor : R = 1.0 / v0;

        (copy k:= 1)
        while k <= dim next k += 1 {
            (copy j :=  1)
            while j <= k next j += 1 {
                r..v[k - 1] += j * r..get(k - j, r0) * get(j, v0);
            }
            r..v[k - 1] /= k;
        }

        return r;
    }


    comp_sin_cos: (inout s: taylor, inout c: taylor, u: taylor, u0: R) = {
        s0: R = std::sin(u0);
        c0: R = std::cos(u0);

        (copy k:= 1)
        while k <= dim next k += 1 {
            (copy j :=  1)
            while j <= k next j += 1 {
                s..v[k - 1] += j * u..get(j, u0) * c..get(k - j, c0);
                c..v[k - 1] -= j * u..get(j, u0) * s..get(k - j, s0);
            }
            s..v[k - 1] /= k;
            c..v[k - 1] /= k;
        }
    }


    sin: (this, v0: R) -> taylor = {
        t: taylor = ();
        r: taylor = ();
        
        comp_sin_cos(r, t, this, v0);
        _ = t;

        return r;
    }

    cos: (this, v0: R) -> taylor = {
        t: taylor = ();
        r: taylor = ();
        
        comp_sin_cos(t, r, this, v0);
        _ = t;

        return r;
    }
}

to_string: <R, order: int> (o: taylor<R, order>) -> std::string = {
    r : std::string = "(";
    (copy i := 1)
    while i <= order next i += 1 {
        r += " (o[i])$";
    }
    r += " )";

    return r;
}

} // cpp2 namespace

#endif // CPP2_CPP2TAYLOR_H