extrapolate Function

    pure function extrapolate(ivals) result(iextrap)
    !! extrapolate the next value in the sequence using
    !! a difference table. Straightfoward implemention:
    !! create the full table, and then evaluate it.
    integer(ip),dimension(:),intent(in) :: ivals
    integer(ip) :: iextrap
    integer :: i, n
    type(int64_vec),dimension(:),allocatable :: diff_table !! difference table (vector of vectors)

    ! create the difference table:
    diff_table = [int64_vec(ivals)]; n = 1
    do
        n = n + 1
        diff_table = [diff_table, int64_vec(diff(diff_table(n-1)%vals))] ! next line is diff of previous line
        if (all(diff_table(n)%vals==0)) exit
    end do

    ! extrapolate
    iextrap = 0
    do i = n-1, 1, -1
        associate( ilast => diff_table(i)%vals(size(diff_table(i)%vals)) )
            iextrap = iextrap + ilast
        end associate
    end do

    end function extrapolate