/***************************************************\
 * Author: Steve Palmer (remlaps@ccil.org)          *
 * Date: 3/2/8                                      *
 *                                                  *
 * Revision history:                                *
 * 3/16/8 - Changed int to long & float to double.  *
 *                                                  *
 * Purpose: Collect first digits of representations *
 * in various bases.  Build squares with base along *
 * the X axis, 2 through argv[1] and numbers along  *
 * the Y axis... also 2 through argv[1].            *
 *                                                  *
 * This program was inspired by "Benford's Law".    *
\***************************************************/
#include <iostream>
#include <cstdlib>
#include <iomanip>

using namespace std;
unsigned _stklen = 1048576 * 512;  /* 512 MB Stack */

int main(int argc, char* argv[])
{
  int MAX=11;                  /* Default if no args */
  if (argc == 2 )
    MAX=atoi(argv[1])+1;      /* Add 1 to 1st arg */

  int FD;

  /* If N > 30, don't bother printing a table.  It's too big. */
    if ( MAX < 31 )
    {
      // Print headings for small-enough N.
      cout << endl;
      cout << "    First Digit of Representation" << endl;
      cout << "          **** Base ****";
      cout << endl << "   ";
      for ( int i=2; i < MAX; i++) {
	cout << i << " ";
        if ( i < 10 )
	  cout << " ";
      }
      cout << endl << "N: ";
      for ( int i=2; i < MAX; i++)
	cout << "---";
      cout << endl;
    } // End if MAX < 31 *** print headings ***

    //Initialize count-array. (first digit can't be 0)
    long count[MAX-2], total=0;
    for (int lcv=0; lcv<MAX-1; lcv++)
      count[lcv]=0;

    // Loop through numbers (X axis)
    for ( int y=2; y < MAX; y++ )
    {
      // Print base on x-axis
      if ( MAX < 31 )
	cout << y << ": ";
       // Loop through bases (Y axis)
       for (int x=2; x < MAX; x++ )
       {
	 int DIV=y/x;       // Integer division, number / base.
	 int MOD=y%x;       // Remainder, number % base.

	 if ( DIV == 0 ) {   //  0 through base-1
	   FD=MOD;
	 } else {          // N = base and above.
	   while ( DIV >= x )   // Get down to first digit.
             DIV = DIV/x;
	   FD=DIV;
	 } // End if DIV == 0
	 if ( MAX < 31 )
	   cout << FD << ", ";

         count[FD-1] += 1;            // Increment first digit count.
	 total += 1;                  // Increment total values checked.
       }  // End loop through N (Y Axis)

       /* Print a line return for small N's */
       if ( MAX < 31 )
          cout << "\n";
    } // End loop through base (X axis)

    /* Display top-ten first digit counts */
    cout << endl << "First Digit Counts:\n";
    for ( int lcv=0; lcv < 10 && lcv < MAX - 1; lcv++ )
      cout << lcv << ": " << count[lcv] << " "
           << count[lcv] / total << endl;

    /* Display top-ten ratios */
    cout << endl;
    for ( int lcv=0; lcv < 10 && lcv < MAX-1; lcv++ ) {
      double ratio = (double) count[lcv] / (double) count[lcv+1];
      cout << setprecision(12) << lcv+1 << " / " << lcv+2 << ": "
           << ratio << endl;
    }
    cout << endl << "Total: " << total << endl;
} // End main()