A better way for me to recurse

[tclviii-dyalog]

in writing a function to generate gamma variates i find myself "double recursing" and feel that it is inefficient

generating gamma variates involves generating candidate variates and then running an acceptance/rejection test

so the "outer" function "gammavariates" recurses until there are enough variates

but there is also an inner function "foo" that is defined within "gammavariates" that recurses that tests iv variable "v" isnt ≤ 0, and recurses to get enough "v"(s)

my fear is that "foo" is being defined inside "gammavariates", so every time "gammavariates" recurses, a redefinition of "foo" occurs, leading to an inefficiency

i suppose i could have a function "testv" that does what "foo" does, but "testv" defined outside "gammavariates" ... but i'd rather keep it as one function.

any ideas about how not to redefine "foo" every time "gammavariates" recurses., or should i just bite the bullet and (SHUDDER) loop! instead of recurse



function below: (yes i know argument beta←1↓⍺ isnt being used, but baby steps)

gammavariates←{

⍝ gamma variates of shape parameter alpha←1↑⍺ , alpha must ≥ 1
⍝ and scale parameter beta←1↓⍺
⍝ ⍵ is number of variates

⍝ uses 'the final algorithm'
⍝ from 'a simple method for generating gamma variables' marsaglia and tsang
⍝ ACM transactions on Mathamatical Software Sept 2000
⍝ use cited paper's un-APLish notation so can be checked against paper

alpha beta←⍺
d←alpha-1÷3
c←1÷(9×d)*0.5

foo←{
⍝ test v (not ≤) 0, preserve x(s) that genned v, recurse till enough v(s)
x←NormRand ⍵
v←(1+c×x)*3
mask←~v≤0 ⋄ v←mask/v ⋄ x←mask/x
result←↑v x ⍝ we need to keep the x, and make result not nested

⍵>≢v:result←result,foo ⍵-≢v ⍝ recurse if not enough v(s)

↓result
}

v x←foo ⍵ ⍝ v(s) from paper and standard normal x(s) that genned em


u←(÷¯2+2*31)×?⍵⍴¯2+2*31 ⍝ u← ⍵ uniform variates twixt 0 and 1

accept1←u<1-0.331×x*4
accept2←(⍟u)<(0.5×x*2)+d×(1-v+⍟v)
accept←accept1∨accept2

result←accept/d×v

(⍵<≢result):result←result,⍺ gammavariates ⍵-≢result ⍝ recurse if not enough variates

result
}

[tclviii-dyalog]

oops, posted earlier function, the final "⍵< tally result" should have flipped "⍵>tally result"

(N.B., cant seem to type tally symbol from keyboard in browser, but can paste it)

[Phil Last]

Hey Tony

So far as I can see your ⍺ is passed back into gammavariates unchanged so c and d are constant throughout. These appear to be the only semi-globals that foo accesses.

If you create another function gv [say] at the same level as foo within gammavariates after creating c and d that does all the other stuff from

      v x←foo ⍵ ⍝ v(s) from paper and standard normal x(s) that genned em

onwards and call that recursively instead of gammavariates I think it'd work with only defining foo the once. gv wouldn't need ⍺ as c and d already exist.

      gammavariates←{
⍝ gamma variates of shape parameter alpha←1↑⍺ , alpha must ≥ 1
⍝ and scale parameter beta←1↓⍺
⍝ ⍵ is number of variates
⍝ uses 'the final algorithm'
⍝ from 'a simple method for generating gamma variables' marsaglia and tsang
⍝ ACM transactions on Mathamatical Software Sept 2000
⍝ use cited paper's un-APLish notation so can be checked against paper
     alpha beta←⍺
     d←alpha-1÷3
     c←1÷(9×d)*0.5
     foo←{
⍝ test v (not ≤) 0, preserve x(s) that genned v, recurse till enough v(s)
         x←NormRand ⍵
         v←(1+c×x)*3
         mask←~v≤0 ⋄ v←mask/v ⋄ x←mask/x
         result←↑v x ⍝ we need to keep the x, and make result not nested
         ⍵>≢v:result←result,foo ⍵-≢v ⍝ recurse if not enough v(s)
         ↓result
     }
     gv←{
         v x←foo ⍵ ⍝ v(s) from paper and standard normal x(s) that genned em
         u←(÷¯2+2*31)×?⍵⍴¯2+2*31 ⍝ u← ⍵ uniform variates twixt 0 and 1
         accept1←u<1-0.331×x*4
         accept2←(⍟u)<(0.5×x*2)+d×(1-v+⍟v)
         accept←accept1∨accept2
         result←accept/d×v
         (⍵<≢result):result←result,gv ⍵-≢result ⍝ recurse if not enough variates
         result
     }
 }

You may not know that recursion itself can be speeded up in a dfn if you replace the literal call to the same function name with a ∇ (del) so "foo" within foo and "gv" within gv could both benefit in that way.

      gammavariates←{
     [...]
     foo←{
         [...]
         ⍵>≢v:result←result,∇ ⍵-≢v ⍝ recurse if not enough v(s)
         [...]
     }
     gv←{
         [...]
         (⍵>≢result):result←result,∇ ⍵-≢result ⍝ recurse if not enough variates
         [...]
     }
 }

[Phil Last]

Sorry. Forgot to call it!
Add new last line

      gv ⍵

or if you prefer

      result←gv ⍵
      result

Then it might work.

[tclviii-dyalog]

Thanks Phil, splitting it up into two subfunctions and only recursing on the second is a great idea.

But now I'm curious as to why a del recurses "faster" then the name of the function itself. Does using the del guarantee a tail call regardless of how clumsily I write it the recursion? (That would be useful because I can never seem to assure myself that I'm not screwing up the tail call rules)

[paulmansour]

I don't know why a del runs faster, but I'm sure it does not guarantee a tail call. An easy way to verify that tail recursion is happening is to simply trace the function through two or more times and see )SI. If the stack does not get deeper, you are doing tail recursion.

[Phil Last]

I don't think it's quite that clever. I don't know the details but similarly to the tail call I'm fairly sure it has something to do with reusing "stack frames" of which I have only a very vague concept. But it certainly works:

      ∇f00←{⍺←⊢
[1]        ⍵<100:f00 ⍵+1
[2]        ⍵
[3]
[4]   ⍝
[5]   ⍝ Phil 2020-02-10 12:50 7.0.0+295
[6]    }
     ∇
      ∇f01∇
     ∇f01←{⍺←⊢
[1]        ⍵<100:∇ ⍵+1
[2]        ⍵
[3]
[4]   ⍝
[5]   ⍝ Phil 2020-02-10 12:50 7.0.0+295
[6]    }
     ∇
      )copy dfns cmpx
C:\Program Files\Dyalog\Dyalog APL-64 18.0 Unicode\ws\dfns.dws saved Thu Dec 19 09:58:56 2019
      cmpx 'f00 1' 'f01 1'
  f00 1 → 4.3E¯5 |   0% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕
  f01 1 → 2.0E¯5 | -55% ⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕⎕

[Phil Last]

It appears that I'm wrong. My example above is a tail call.

Changing lines[1] to read

      ⍵<100:z⊣z←f00 ⍵

and

      ⍵<100:z⊣z←∇ ⍵

respectively, wipes out the timing advantage.

So I guess the advantage is all in that you can call an anonymous function recursively and you can rename a named function without having to find all occurrencies of recursion. Nothing to do with timing.