// Direct algebraic simulation for: "Playing pool with pi", G. Galperin,
//  Regular and Chaotic Dynamics, v.8, no.4, 2003.
// See also: "Playing Pool with |Psi>: from Bouncing Billiards to Quantum Search",
//  Adam R. Brown, 30 Oct. 2020, https://arxiv.org/abs/1912.02207
// And: "Bouncing Balls and Quantum Computing", Don Monroe, CACM, Oct. 2020.
//
// R. Perry, June 2021
//
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define V (-1) // initial velocity of M, value doesn't really matter if U=0
#define U (0)  // initial velocity of m

void disp( int table, int c, int w, double t, double tmin, double tmax,
 double x, double y, double u, double v)
{
  if( table)
    printf( "%6i %2c %12.6g %12.6g %12.6g %12.6g %12.6g\n", c, w, t, x, y, u, v);
  else if( t >= tmin && t <= tmax)
    printf( "%g %g %g\n", t, x, y);
}

int main( int argc, char *argv[]) // args: -p -g M tmin tmax
{
 // block m (m=1) starts at position x=0.5 with velocity u=U
 // block M starts at position y=1 with velocity v=V
 // E is 2*energy, P is momentum
 //
  int c = 0, w = ' ', plot = 0, table = 1, gflag = 0;
  double N = 0, M = 1, M1, u = U, v = V, P;
  double dt, t = 0, tmin = 0, tmax = HUGE_VAL, x = 0.5, y = 1;
  long double pi = 4*atanl(1);

 // table (default) or plot, mutually exclusive
 //
  if( argc > 1 && strcmp( argv[1], "-p")  == 0) { plot = 1; --argc; ++argv; }
  if( argc > 1 && strcmp( argv[1], "-g") == 0) { gflag = 1; --argc; ++argv; }
  if( argc > 1) { M = atof( argv[1]); N = log(M)/log(100); } // M = pow(100,N)
  if( argc > 2) tmin = atof( argv[2]);
  if( argc > 3) tmax = atof( argv[3]);

  table = !plot; M1 = M+1; P = u + M*v; // double E = u*u + M*v*v;

  if( gflag) u = v = -1/sqrt(M1); // for analogy with Grover search

  long double Q = pi*powl(10,N); unsigned long long q = Q;

  if( table) printf( "N = %g, M = %g, q = %llu ", N, M, q); // printf( "E = %g, P = %g\n", E, P);

  if( q+1 < Q) { printf( "overflow\n"); return 1; } else if( table) printf( "\n");

  if( table)
    printf( "%6s %2s %12s %12s %12s %12s %12s\n", "coll", "BW", "t", "x", "y", "u", "v");

  disp( table, c, w, t, tmin, tmax, x, y, u, v);

 // collisions
 //
  while( u < 0 || u > v)
  {
    ++c;

    if( u < 0) // m-wall collision: u and P changed, v unchanged
    {
      w = 'W'; // Wall collision: 0 = x + u*dt
      dt = -x/u; t += dt; x = 0; y += v*dt;
      u = -u; P = u + M*v;
    }
    else // m-M collision: u and v changed, P unchanged
    {
      w = 'B'; // Blocks collision: x + u*dt = y + v*dt
      dt = (y-x)/(u-v); t += dt; x = y = x + u*dt;
      v = (P + u - v)/M1; u = P - M*v;
    }

    disp( table, c, w, t, tmin, tmax, x, y, u, v);

    // printf( " P = %g\n", P);
    // double EE = u*u + M*v*v;
    // printf( " EE = %g, diff = %g\n", EE, EE-E);
  }

  w = ' '; dt = 0.5; t += dt; x += u*dt; y += v*dt;

  disp( table, c, w, t, tmin, tmax, x, y, u, v);
}