/* -*- Mode: C -*- */

/* problem-loop.c */

#ifndef _PROBLEM_LOOP_I
#define _PROBLEM_LOOP_I

#include <stdio.h>

#define KK 100

void problem_loop(long seed) {

  register int j;
  
  double ulp;
  double ss;

  printf("step (0): seed = %ld, (seed & 0x3fffffff) = %ld, seed = %ld\n",
	 seed, (seed & 0x3fffffff), seed);

  ulp = (1.0 / (1L << 30)) / (1L << 22);	/* 2 to the -52 */
  ss  = 2.0 * ulp * ((seed & 0x3fffffff) + 2);

  printf("step (1):\n\t ulp  = %.20f\n\t ss   = %.20f\n\n", ulp, ss);

  for (j = 0; j < KK; j++) {
    
    printf("\t j = %d; ss = %.20f\n", j, ss);
    
    ss += ss;
    if (ss >= 1.0) ss -= 1.0 - 2 * ulp;	/* cyclic shift of 51 bits */
  }

  printf("step (2):\n\t ulp  = %.20f\n\t ss   = %.20f\n\n", ulp, ss);
}


/* main --
 * Testing...
 */

int main() {
  problem_loop(310952);
}

#endif /* _PROBLEM_LOOP_I */

/* end of file -- problem-loop.c */