module ProjectionMethods.Infrastructure

let getNextJacobiPolynomial previousPolynomial ancientPolynomial k n =
    let float_n = n |> float
    let float_k = k |> float

    fun x ->
        (-)
            ((/)
                ((float_n + float_k + 2.0)
                 * (2.0 * float_n + 2.0 * float_k + 3.0)
                 * x
                 * (x |> previousPolynomial))
                ((float_n + 2.0 * float_k + 2.0) * (float_n + 2.0)))
            ((/)
                ((float_n + float_k + 2.0) * (float_n + float_k + 1.0) * (x |> ancientPolynomial))
                ((float_n + 2.0 * float_k + 2.0) * (float_n + 2.0)))

let getJacobiPolynomials k n =
    let polynomial_0 = (fun _ -> 1.0)

    if n < 0 then
        raise (invalidArg "n" "n should be more than 0")
    else if n = 0 then
        [| polynomial_0 |]
    else
        let polynomial_1 = ((*) (k + 1.0))

        let rec loop i result =
            if i = n then
                result |> List.rev |> List.toArray
            else
                let nextPolynomial = getNextJacobiPolynomial result[0] result[1] (i + 1) k
                loop (i + 1) (nextPolynomial :: result)

        loop 1 [ polynomial_1; polynomial_0 ]
