$M=5000; # granularity
$N=2000;  # number of terms in the sum
$lambda=1; # must be <= 1
$alpha=1;

$beta=1/$alpha;

open(OUT,">solve.txt");

for ($k=0; $k<=$M; $k++) { $a[$k]=1; }

for ($iter=0; $iter<10; $iter++) {
  $sum=0;
  for ($k=0; $k<=$M; $k++) {
    $b[$k]=0;
    $x=$k/$M;
    for ($n=1; $n<=$N; $n++) {
      $y=($lambda/($x+$n))**$beta;
      $y=int(0.5+ $M*$y);  ## must be: 0 <= y < M
      if ($y<=$M) { 
        $b[$k]+=$a[$y]*$beta*($lambda**$beta)*(($x+$n)**(-1-$beta));
      } else {
        print "ERROR\n"; # ..... error? exit?
      }
    }
    $sum+=$b[$k];
    if ($k%20 ==0) {print OUT "$iter\t$x\t$b[$k]\n"; }
  }
  for ($k=0; $k<=$M; $k++) { $a[$k]=$b[$k]*$M/$sum; } 
  print "$iter | $M | $a[0] | $a[$M] | $sum\n";
}

close(OUT);