/*
 * From: David Goldberg <goldberg@parc.xerox.com>
 * To: numeric-interest@validgh.com
 * Subject: Re: 7 digit decimal number
 * Message-Id: <99Apr12.091046pdt.177505@arcturia.parc.xerox.com>
 * Date: Mon, 12 Apr 1999 09:10:41 PDT
 * Sender: owner-numeric-interest@validgh.com
 * Precedence: bulk
 * Reply-To: David Goldberg <goldberg@parc.xerox.com>
 *
 * Jerome beat me to the punch.  My answer was slightly smaller,
 * 9.0e9.  Here is the program with comments that prints out
 * a lot of examples (include 9.0e9):
 *
 * -david
 *
 */
/*
 * Find 7 digit decimal number d so that d rounded to 24 bits
 * (single prec) and back to 7 digits does not recover d.
 *
 *  Use fact that 2^33 =  8589934592
 *           and 10^10 = 10000000000
 *
 *   are close, so that there are more decimal numbers
 *   with 7 digits in the range (2^33, 10^10) than binary numbers
 *   with 24 bits.  Computing how many numbers in that range:
 *
 *    7 digit decimal numbers in (2^33, 10^10) are of form
 *                  i*10^3,  2^33/10^3 < i < 10^7
 *    24 bit binary  numbers   in (2^33, 10^10) are of form
 *                  j*2^10,       2^23 < j < 10^10/2^10
 *
 *    10^7 - 2^33/10^3 =  1410065 decimal numbers
 *    10^10/2^10 - 2^23 = 1377017 binary numbers
 *
 *   Other pairs besides (2^33, 10^10) are (2^43, 10^13),
 *     (2^53,10^16), (2^63,10^19), etc.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include "args.h"

#if STDC
double anint(double x)			/* missing function: added by NHFB */
#else
double 
anint(x)				/* missing function: added by NHFB */
double x;
#endif
{
    if (x != x)
	return (double)INT_MIN;
    else if ((x - 0.5) < (double)INT_MIN)
	return (double)INT_MIN;
    else if ((x + 0.5) > (double)INT_MAX)
	return (double)INT_MAX;
    else if (x < 0.0)
	return (double)(long)(x - 0.5);
    else
	return (double)(long)(x + 0.5);
}

int
main(VOID_ARG)
{
    double dlow, dhi, i, x, m;
    double d;  /* the 7 decimal number */
    double b;  /* rounded to 24 bit binary */

    dlow = ceil(scalb(1.0, 33.0)/1000.0); /* 2^33/10^3 */
    dhi = 10000000.0;                     /* 10^7 */
    m = 536870913.0;                      /* 2^k + 1 = 2^29 + 1 */

    for (i = dlow; i < dhi; i++) {
	d = i*1000.0;
	/* round to 24 bits, k=29 p-k=24 */
	b = x = (m*d) - (m*d - d);

	/* round back to 7 decimal digits */
	x = anint(x/1000.0)*1000.0;
	if (x != d) {
	    (void)printf("%9.0f, %.0f x 2^10,  %9.0f\n",
			 d, b/1024.0, x);
	}
    }
    return (EXIT_SUCCESS);
}

/*
 *  First answer is 8589973000, 8388646 x 2^10,  8589974000
 *
 *  Checking:
 *     8589973000 is 100000000000000000100101 1000001000
 *   which rounds to 100000000000000000100110 0000000000
 *   which in decimal is 8589973 504
 *   which rounds to 8589974 000
 */