Application of recursion coarsening to compute Obara-Saika recurrences

[jiahao]

Here is good thesis fodder, summarizing two separate discussions with @alanedelman yesterday and @stevengj today.

• Recursion over type parameters can be used for recursion coarsening.
  • Recursion is a natural way to express lots of scientific computations, particularly to implement divide-and-conquer algorithms.
  • Recursion incurs overhead from function calls.
  • However, if the base case is big enough, recursion is not that expensive in practice.
  • The trick is to don't recurse all the way down to the base case.
  • You can coarsen the base case (make it bigger) by unrolling the recursion.
  • The only paper @stevengj can find on recursion unrolling is \cite{RR00}.
  • Recursion unrolling is a generic optimization that can be implemented at the compiler level.
  • The problem is that one can have a combinatorial explosion of base cases. So this is a pain to do manually.
  • Inlining avoids function call overhead, and has the side benefit of facilitating instruction level parallelism.
  • Using recurrences at the level of type parameters, combining with inlining passes in the Julia compiler, allows the compiler to automatically generate many coarsened base cases.
• Application: computing $\int_{0}^{1} tan^n x dx$ using recurrences
• Application: Obara-Saika recurrences in quantum chemistry \cite{OS76}
  • This notebook implements Obara-Saika recurrences for the very simplest integrals - two center overlap integrals https://gist.github.com/jiahao/09f034be159e4959256b
  • The core computation deals with a large number of one-dimensional integrals.
  • The base case is very small, it is just computing a square root and an exponential
  • The recurrences can be applied in two different ways, and ordinarily one has to choose manually which way to apply them.
  • In this implementation we allow the compiler to decide how to inline the base cases
  • Has multiple dispatch because you are computing integrals over two functions, each of which has one parameter
  • Contrast with expert implementation like evaleev/libint using templated C++ functions. The relevant lines are preceded by lots of boilerplate.
  • Open question: why not C++ templates?
    • Julia is easier to use
    • Dynamic code generation - why is this better? Maybe flexibility in generating coarsened base cases?

@inproceedings{RR00,
    author = {Rugina, Radu and Rinard, Martin C.},
    title = {Recursion Unrolling for Divide and Conquer Programs},
    booktitle = {Proceedings of the 13th International Workshop on Languages
and Compilers for Parallel Computing-Revised Papers},
    series = {LCPC '00},
    year = {2001},
    isbn = {3-540-42862-3},
    pages = {34--48},
    numpages = {15},
    url = {http://dl.acm.org/citation.cfm?id=645678.663942},
    acmid = {663942},
    publisher = {Springer-Verlag},
    address = {London, UK},
    editors = {Samuel P. Midkiff and José E. Moreira and Manish Gupta and
        Siddhartha Chatterjee and Jeanne Ferrante and Jan Prins and
            William Pugh and Chau-Wen Tseng},
}

@article{OS86,
   author = {Obara, S. and Saika, A.},
   title = {Efficient recursive computation of molecular integrals over {Cartesian} {Gaussian} functions},
   journal = {Journal of Chemical Physics},
   year = 1986,
   volume = 84,
   number = 7, 
   pages = {3963--3974},
   doi = {10.1063/1.450106}, 
}