Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Computing recurrence relations in Haskell

Greetings, StackOverflow.

Let's say I have two following recurrence relations for computing S(i,j)

S_{i,j+1} = X_{PA}S_{i,j} + \frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1}) \\ S_{i+1,j} = X_{PB}S_{i,j} + \frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1})

I would like to compute values S(0,0), S(0,1), S(1,0), S(2,0), etc... in asymptotically optimal way. Few minutes with pencil and paper reveal that it unfolds into treelike structure which can be transversed in several ways. Now, it's unlikely tree will be useful later on, so for now I'm looking to produce nested list like [[S(00)],[S(10),S(01)],[S(20),S(21),S(12),S(02)],...]. I have created a function to produce a flat list of S(i,0) (or S(0,j), depending on first argument):

osrr xpa p predexp = os00 : os00 * (xpa + rp) : zipWith3 osrr' [1..] (tail osrr) osrr
  where
    osrr' n a b = xpa * a + rp * n * b
    os00  = sqrt (pi/p) * predexp
    rp    = recip (2*p)

I am, however, at loss as how to proceed further.

like image 342
Victor Zvodin Avatar asked Jan 21 '23 04:01

Victor Zvodin


1 Answers

I would suggest writing it in a direct recursive style and using memoization to create your traversal:

import qualified Data.MemoCombinators as Memo

osrr p = memoed
    where
    memoed = Memo.memo2 Memo.integral Memo.integral osrr'
    osrr' a b = ...  -- recursive calls to memoed (not osrr or osrr')

The library will create an infinite table to store values you have already computed. Because the memo constructors are under the p parameter, the table exists for the scope of p; i.e. osrr 1 2 3 will create a table for the purpose of computing A(2,3), and then clean it up. You can reuse the table for a specific p by partially applying:

osrr1 = osrr p

Now osrr1 will share the table between all its calls (which, depending on your situation, may or may not be what you want).

like image 118
luqui Avatar answered Jan 30 '23 22:01

luqui