/* ---------------------------------------------------------------- */
/* Associative law problems:                                        */
/*    See the following paper:                                      */
/*    Jean-Francois Colonna,                                        */
/*    Kepler, von Neumann and God (More rounding-off error          */
/*         visualizations),                                         */
/*    The Visual Computer, Vol. 12 (1996), pp. 346-349.             */
/*                                                                  */
/* The formulas is X = (R + 1)*X - R*X^2, which has the following   */
/* mathematically equivalent variations:                            */
/*    (1)  X = (R + 1)*X - R*(X*X)                                  */
/*    (2)  X = (R + 1)*X - (R*X)*X                                  */
/*    (3)  X = ((R + 1) - R*X)*X                                    */
/*    (4)  X = R*X + (1 - R*X)*X                                    */
/*    (5)  X = X + R*(X - X*X)                                      */
/* where R > 2.57 causes the so-called "chaotic" hehavior.  In this */
/* program, R = 3.0 and the initial value of x is 0.5.              */
/* ---------------------------------------------------------------- */


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

typedef   double  REAL;

void  main(int argc, char *argv[])
{
     long  n, m, k;
     REAL  x1, x2, x3, x4, x5;
     const REAL  R = 3.0;
     const REAL  INIT = 0.5;

     if (argc != 3) {
          printf("Use %s #-of-iteration steps\n", argv[0]);
          exit(1);
     }
     sscanf(argv[1], "%ld", &n);
     sscanf(argv[2], "%ld", &m);

     x1 = x2 = x3 = x4 = x5 = INIT;
     printf("\n%5ld  %12g  %12g  %12g  %12g  %12g", 0, x1, x2, x3, x4, x5);
     for (k = 1; k <= n; k++) {
          if (k % m == 0)
               printf("\n%5d  %12g  %12g  %12g  %12g  %12g", 
                         k, x1, x2, x3, x4, x5);
          x1 = (R+1)*x1 - R*(x1*x1);
          x2 = (R+1)*x2 - (R*x2)*x2;
          x3 = ((R+1) - R*x3)*x3;
          x4 = R*x4 + (1-R*x4)*x4;
          x5 = x5 + R*(x5 - x5*x5);
     }
     printf("\n");
}