Re: [Bug-gsl] Re: Wrong orthogonalization in gsl_rng_gfsr4

[Andreas Schneider]

Brian Gough wrote:
> Andreas Schneider writes:
>  > I should have searched the bug-gsl list first. Now I found that
>  > this is a known bug:
>  > http://mail.gnu.org/archive/html/bug-gsl/2003-07/msg00001.html
>  > 
>  > But the proposed fix for 2.0 has a bug too: The array elements
>  > containing the orthogonalization will be destroyed before ever
>  > being used.
>
> Can you point out where that happens exactly. Thanks.

The orthogonalization affects indices up to 220.
state->nd is initialized with M (==16383).

The first call to gfsr4_get() will set state->nd to 0 and does

state->ra[0]=
  state->ra[15913] ^ state->ra[14798] ^ state->ra[9396] ^ state->ra[6695];

Next calls do

state->ra[1]=
  state->ra[15914] ^ state->ra[14799] ^ state->ra[9397] ^ state->ra[6696];

state->ra[2]=
  state->ra[15915] ^ state->ra[14800] ^ state->ra[9398] ^ state->ra[6697];

...

The initial lower 6695 indices do never participate in calculations of 
subsequent random numbers.

A possible fix would be to replace

  for (i = 0; i <= M; i++)

with

  for (i = 0; i < D; i++)

in gfsr4_set(). Then, the first call to gfsr4_get() yields

state->ra[9689]=
  state->ra[9218] ^ state->ra[8103] ^ state->ra[2701] ^ state->ra[0];

Andreas Schneider