Perl Guru Forums: Perl Programming Help: Intermediate: C to Perl conversion

CGI/Perl Guide | Learning Center | Forums | Advertise | Login

Site Search: in


	MAIN INDEX	SEARCH POSTS	WHO'S ONLINE	LOG IN

Home: Perl Programming Help: Intermediate:
C to Perl conversion

katomich
New User

May 26, 2009, 11:12 PM

Post #1 of 4 (387 views)

C to Perl conversion	Can't Post

Does anybody know of a working version of ran2.c or ran2.f from Numerical Recipes in Perl? It seems that, whatever I do, I cannot simulate the pointers in C to Perl. Here is my attempt below. Thanks.

#!/usr/bin/env perl

use constant IM1 => 2147483563;
use constant IM2 => 2147483399;
use constant AM => (1.0/IM1);
use constant IMM1 => (IM1-1);
use constant IA1 => 40014;
use constant IA2 => 40692;
use constant IQ1 => 53668;
use constant IQ2 => 52774;
use constant IR1 => 12211;
use constant IR2 => 3791;
use constant NTAB => 32;
use constant NDIV => (1+IMM1/NTAB);
use constant EPS => 1.2e-7;
use constant RNMX => (1.0-&EPS);

sub ran2 {

my $idum = $_[0];
$idum2 = 123456789;
$iy = 0;
my @iv = ();

if ($idum <= 0) {
if (-($idum) < 1) {$idum=1;}
else {$idum = -($idum);}
$idum2=($idum);
for ($j=NTAB+7; $j>=0; $j--) {
$k=($idum)/IQ1;
$idum=IA1*($idum-$k*IQ1)-$k*IR1;
if ($idum < 0) {$idum += IM1;}
if ($j < NTAB) {$iv[$j] = $idum;}
}
$iy=$iv[0];
}

$k=$idum/IQ1;
$idum=IA1*($idum-$k*IQ1)-$k*IR1;
if ($idum < 0) {$idum += IM1;}
$k=$idum2/IQ2;
$idum2=IA2*($idum2-$k*IQ2)-$k*IR2;
if ($idum2 < 0) {$idum2 += IM2;}
$j=$iy/NDIV;
$iy=$iv[$j]-$idum2;
$iv[$j] = $idum;
if ($iy < 1) {$iy += IMM1;}
if (($temp=AM*$iy) > RNMX) {return RNMX;}
else {return $temp;}
}

my $idum=-1;
for($m=1; $m<=10; $m++) {
$x = ran2($idum);
print "idum=" . $idum . ", " . "output=" . $x . "\n";
}

**spacing not the issue

1arryb
User

May 27, 2009, 8:10 AM

Post #2 of 4 (372 views)

Re: [katomich] C to Perl conversion [In reply to]	Can't Post

Hi katomich,

My first reaction to your post was, "what, like we need another pseudo-random number generator in Perl?"
A quick search on CPAN yields 10 or 20 random number generation packages. Surely one of them will work for you.

But maybe you are just curious or there's something special about this function that I'm too ignorant about numerical
methods to appreciate, so we carry with the analysis.

First of all, your question about pointers is just a distraction. It has nothing to do with why your program doesn't work.
If you are curious about how variables and references work in perl, read the perlvar and perlref perldocs.
But before you jump down that rabbit hole, fix the bugs in your script:

1. You are always calling ran2() with the same value (the loop in your main program is on $m but you call ran2($idum).

2. Even after I fixed this problem, ran2() always returns the same value, no matter what the argument is.
I stepped through the function using the perl debugger (perl -d).
Here are a couple of problems that jumped out at me:

Code

main::ran2(tmp.pl:55):      $j=$iy/NDIV;   # $iy is zero here. 
main::ran2(tmp.pl:56):      $iy=$iv[$j]-$idum2; # @iv is an empty array here. Also, $j has no value.

Note: I reformatted your script so I could figure out what was going on. Line numbers above refer to my version, not yours.

Code

#!/usr/bin/perl 
 
use constant IM1 => 2147483563; 
use constant IM2 => 2147483399; 
use constant AM => (1.0/IM1); 
use constant IMM1 => (IM1-1); 
use constant IA1 => 40014; 
use constant IA2 => 40692; 
use constant IQ1 => 53668; 
use constant IQ2 => 52774; 
use constant IR1 => 12211; 
use constant IR2 => 3791; 
use constant NTAB => 32; 
use constant NDIV => (1+IMM1/NTAB); 
use constant EPS => 1.2e-7; 
use constant RNMX => (1.0-&EPS); 
 
sub ran2 { 
 
    my $idum = $_[0]; 
    $idum2 = 123456789; 
    $iy = 0; 
    my @iv = (); 
 
    if ($idum <= 0) { 
        if (-($idum) < 1) { 
            $idum=1; # How will -($idum) ever be < 1 here? 
        } else { 
            $idum = -($idum); 
        } 
        $idum2=($idum); 
        for ($j=NTAB+7; $j>=0; $j--) { 
            $k=($idum)/IQ1; 
            $idum=IA1*($idum-$k*IQ1)-$k*IR1; 
            if ($idum < 0) { 
                $idum += IM1; 
            } 
            if ($j < NTAB) { 
                $iv[$j] = $idum; 
            } 
        } 
        $iy=$iv[0]; 
    } 
 
    $k=$idum/IQ1; 
    $idum=IA1*($idum-$k*IQ1)-$k*IR1; 
    if ($idum < 0) { 
        $idum += IM1; 
    } 
    $k=$idum2/IQ2; 
    $idum2=IA2*($idum2-$k*IQ2)-$k*IR2; 
    if ($idum2 < 0) { 
        $idum2 += IM2; 
    } 
    $j=$iy/NDIV; 
    $iy=$iv[$j]-$idum2; 
    $iv[$j] = $idum; 
    if ($iy < 1) { 
        $iy += IMM1; 
    } 
    if (($temp=AM*$iy) > RNMX) { 
        return RNMX; 
    } 
    else {return $temp;} 
} 
 
#my $idum=-1; 
for($m=1; $m<=10; $m++) { 
    $x = ran2($m); 
    print "idum=" . $m . ", " . "output=" . $x . "\n"; 
}

Cheers,

Larry

katomich
New User

May 27, 2009, 8:56 AM

Post #3 of 4 (369 views)

Re: [1arryb] C to Perl conversion [In reply to]	Can't Post

First, thanks for the reply. I require a better random number generator than CPAN can give me. If you refer to Numerical Recipes: The Art of Scientific Computing by Teukolsky, Vetterling, and Flannery, there is a portable routine named ran2.(either C or Fortran version) which has stood under rigorous statistical tests. In this routine, one should only initialize the seed -- idum in this case -- once. (If you think about it, using $m as the seed in the loop no longer makes the output random.) ran2 then uses this value and generates new idums to compute the other random uniform deviates outside the if block. You can see my predicament outside the subroutine in main. I initialize idum to -1, but somehow get the same values. So basically, the idea is to use one idum to get 10 random values. That is why the C version uses a pointer *idum in place of my $idum scalar variable. Once a new idum is computed inside the subroutine, memory is allocated to this new value in idum.

1arryb
User

Jun 2, 2009, 11:56 AM

Post #4 of 4 (349 views)

Re: [katomich] C to Perl conversion [In reply to]	Can't Post

Hi katomich,

You explicitly initialize $idum to the first argument of ran() each time it's called. Therefore, $idum can only have 10 values: 1..10.
If you want $idum's value to persist between ran() calls, declare it outside of ran() in your main program. However, if you do that,
then I'm clueless about what is the purpose of $m except as a loop counter.

If you don't want ran() to rely on a global $idum, you can pass it as a reference. This is the closest you'll get to a pointer in Perl:

Code

... 
sub ran { 
    my ($m, $idum_ref) = @_; 
    ... 
    $$idum_ref = <whatever>; # '$$' is how you dereference a scalar reference. 
    ... 
} 
... 
my $idum = -1; 
for ( my $m = 0; $m < 10; $m++) { 
    my $x = ran($m, \$idun); # The backslash means "take a reference to $idum". 
    print "ran($m) = $x\n"; 
}

Cheers,

Larry

Search for (options)

Powered by Gossamer Forum v.1.2.0

Web Applications & Managed Hosting Powered by Gossamer Threads
Visit our Mailing List Archives