Return to Revisions

I investigated the sequence $Y_k=X_{3k}$, which has far more communal planes, and thus more useful to build a random generator. Of course, choosing large values for $b_1,b_2$ will further drastically improve the generator by adding a lot more planes. I suggest choosing values larger than (say) $2^{30}$ for $b_1,b_2$.

There are $M=b_2^3+3b_1b_2+b_1^3$ communal planes and they all have an equation of the form

$$b_1^3\cdot Y_k+b_2(b_2^2+3b_1) \cdot Y_{k+1}-Y_{k+2} =d.$$

The possible values for $d$ are $0,1,\cdots,M-1.$ Each plane (identified by $d$) contains a different proportion of triplets $(Y_k,Y_{k+1},Y_{k+2})$. The empirical distribution for these proportions is featured in the histogram below (corresponding to $b_1=5,b_2=3$), where the X-axis represents $d$, and the Y-axis the proportion of triplets lying in plan $d$.

Of course it is easy by looking at this chart to guess what the exact theoretical distribution is. To identify these planes, I used the program below and some experimental math.

# Compute equations of planes containing 3 random vectors
#   P(k) = (x[k], x[k+1], x[k+2])
#   P(l) = (x[l], x[l+1], x[l+2])
#   P(m) = (x[m], x[m+1], x[m+2])
# (k, l, m) are randomly selected (M triplets)
#
# Equation of planes is x + s*y + t*z = intercept
# For each (k,l,m) output the coefficients s, t, intercept 
#
# Goal: Find communal planes absorbing many (P(k), P(l), P(m))
# Once the planes are computed, sort them by s, t, intercept

$n=100000;

$b1=5; 
$b2=3;

# xx[] is the original sequence

$xx[0]=0.5;
$xx[1]=sqrt(2)/2;

for ($k=2; $k<$n; $k++) {
  $xx[$k]=$b2*$xx[$k-1]+$b1*$xx[$k-2]-int($b2*$xx[$k-1]+$b1*$xx[$k-2]); 
  if ($xx[$k]<0) { $xx[$k]=1+$xx[$k]; }
}

# we actually use 1 out of 3 consecutive terms from original sequence xx[]
# to see if it the new sequence x[] also has a small number of communal planes 

for ($k=0; $k< $n/3; $k++) {
  $x[$k]=$xx[3*$k];
}

$M=10000; # must be < n/3
open(OUT,">coplanes2.txt");

for ($iter=0; $iter<$M; $iter++) {

       $k=int($M*rand()); 
       $l=int($M*rand());
       $m=int($M*rand());

       # in case k=l or k=m or l=m, an ERROR message is reported

       $a=$x[$k]; $b=$x[$k+1]; $c=$x[$k+2];
       $d=$x[$l]; $e=$x[$l+1]; $f=$x[$l+2];
       $p=$x[$m]; $q=$x[$m+1]; $r=$x[$m+2];
       $u=($e-$b)*($r-$c)-($f-$c)*($q-$b);
       $v=-($d-$a)*($r-$c)+($f-$c)*($p-$a);
       $w=($d-$a)*($q-$b)-($e-$b)*($p-$a);

       if ($u != 0) {
         $s=$v/$u;
         $t=$w/$u;
         $intercept=($u*$a + $v*$b + $w*$c)/$u;

         print OUT "$k\t$l\t$m\t";
         print OUT "$s\t$t\t$intercept\n";

       } else {
         print OUT "$k\t$l\t$m\tERROR (u=0)\n";
       }
}  
close(OUT);