anchor.diff:          t=1.2*sq*(1.+y**2)*dexp(oldg-gammln(em+1.) 
        +-gammln(en-em+1.)+em*plog+(en-em)*pclog) 
   C       Reject. This happens about 1.5 times per deviate, on average.
  -        call random_number(rn)
  +!         call random_number(rn)
  +          rn = rand()
           if (rn.gt.t) goto 2 
           x=em 
           endif 
  
positive.code:         sq=sqrt(2.*am*pc) 
  C       The following code should by now seem familiar: rejection 
  C       method with a Lorentzian comparison function.
          call random_number(rn)
      2   y=tan(PI*rn) 
          em=sq*y+am 
  C       Reject.
  
negative.code:       DOUBLE PRECISION X, Y, S, ONE, TWO
        DATA ONE /1.0/, TWO /2.0/
  C
  !     1 call random_number(x)
      1 X = rand()
  !       call random_number(y)
        Y = rand()
        X = TWO * X - ONE
        Y = TWO * Y - ONE
        S = X * X + Y * Y
  