get the internal state from a PHP lcg in 3 outputs

fast_reverse.c

/* 
 *
 */

#include <stdio.h>
#include <math.h>

#define MODMULT(a, b, c, m, s) q = s/a;s=b*(s-a*q)-c*q;if(s<0)s+=m;
/* MODMULT computes s*b mod m, provided that m=a*b+c and 0<=c<m */

#define EPSILON 0.00001

double lcg_forward(int *s1, int *s2)
{
  double z;
  int q;
  MODMULT(53668, 40014, 12211, 2147483563L, *s1);
  MODMULT(52774, 40692, 3791, 2147483399L, *s2);
  z = *s1 - *s2;
  if (z < 1){
    z += 2147483562;
  }
    z *= 4.656613e-10;
    return z;
}


int main(int argc, char** argv)
{
  double y[10];
  int num_outputs;

  if (argc < 3)
  {
      printf("usage: %s y1 y2 ... (where y1, y2, are successive outputs from RNG)\n\n", argv[0]);
      return -1;
  }

  int i;
  int possible_candidate;

  // read in outputs from RNG
  for (i=1; i < argc; ++i) {
    sscanf( argv[i], "%lf", &y[i-1] );
    // don't overflow the buffer
    if (i == sizeof(y)/sizeof(double) - 1)
      break;
  }
  num_outputs = argc-1;

  printf("\nOutputs of RNG: ");
  for (i=0; i < num_outputs; ++i)
    printf("%lf ",y[i]);
  printf("\n\n");

  int q;
  int s1AfterFirstRound;
  int s2AfterFirstRound;
  double z;
  double zFirstRound;

  zFirstRound = y[0]/4.656613e-10;

  // try all possibilities for s1AfterFirstRound
  for (s1AfterFirstRound = 0; s1AfterFirstRound < (int) ((1L<<31)-1); ++s1AfterFirstRound) {
    // If  s1AfterFirstRound  is right, then
    // s2AfterFirstRound  is either  s1AfterFirstRound - zFirstRound  or  s1AfterFirstRound - zFirstRound + 2147483563 .
    int s1AfterNextRound;
    int s2AfterNextRound;

    s2AfterFirstRound = (int) s1AfterFirstRound - zFirstRound + 0.5;
    // s2AfterFirstRound is not allowed to be negative
    if (s2AfterFirstRound < 0) {
      // If  s2AfterFirstRound  was negative then we had an overflow, so we must correct it.
      s2AfterFirstRound = (int) s1AfterFirstRound - zFirstRound + 2147483562 + 0.5;
    }

    possible_candidate = 1;
    s1AfterNextRound = s1AfterFirstRound;
    s2AfterNextRound = s2AfterFirstRound;
    int j;
    for (j=0; j < num_outputs-1; ++j) {
      MODMULT(53668, 40014, 12211, 2147483563L, s1AfterNextRound);
      MODMULT(52774, 40692, 3791, 2147483399L, s2AfterNextRound);
      z = s1AfterNextRound - s2AfterNextRound;
      if (z < 1)
        z += 2147483562;
      z *= 4.656613e-10;
      if (fabs(z-y[j+1]) > EPSILON) {
        possible_candidate = 0;
        break;
      }
    }
    if (possible_candidate) {
      printf("===>CANDIDATE: s1 = %d, s2 = %d\n", s1AfterFirstRound, s2AfterFirstRound);
      printf("next is %f\n", lcg_forward(&s1AfterFirstRound, &s2AfterFirstRound));
      printf("next is %f\n", lcg_forward(&s1AfterFirstRound, &s2AfterFirstRound));
      printf("next is %f\n", lcg_forward(&s1AfterFirstRound, &s2AfterFirstRound));
      printf("next is %f\n", lcg_forward(&s1AfterFirstRound, &s2AfterFirstRound));
    }
  
  }

  printf("\n");

  return 0;
}