Finding the largest double/float (non-"infinite") number...

[Robert G. Brown]

Dear ListHumans,

OK, so I'm compulsive.  I worked this one out for the hell of it.  I
tried to post a tarball of the solution but got nailed by the dread 40K
message limit and the message has to await moderation, and the solution
isn't really that big.  The following subroutine(s), compiled with the
appropriate (ISO C99) define, will do a
not-terribly-efficient-but-adequate binary search for the biggest double
or float that isn't "infinity" on a given (ISO C99-supporting) system.
Using FPE signals turned out to be too bit a PITA and with all the new
ISO C99 calls, unnecessary.  Just learning about the ISO C99 macros in
math.h was worth the exercise -- it is now easy to trap Inf and NaN and
so forth in runtime code (all of these have plagued me at one time or
another in the past).

They are wrapped up for inclusion in arbitrary code below.  I actually
made the routines portable after sending out the tarball, so if you want
a full package that generates the numbers for both float and double
you'll have to ignore the email tarball and grab:

  http://www.phy.duke.edu/brahma/maxsize-0.1.1.tgz

   rgb

-- 
Robert G. Brown	                       http://www.phy.duke.edu/~rgb/
Duke University Dept. of Physics, Box 90305
Durham, N.C. 27708-0305
Phone: 1-919-660-2567  Fax: 919-660-2525     email:rgb at phy.duke.edu

%< --- Snip snip snip --------------------------------------------------

/*
 *========================================================================
 * $Id: maxvalue.c,v 1.1 2001/04/20 16:21:13 rgb Exp $
 *
 * See copyright in copyright.h and the accompanying file COPYING
 *  (a GPL v2b, for those that care:-)
 *========================================================================
 */

/*
 *========================================================================
 * Two routines to find the largest double or float that a system
 * considers not "infinite".  If verbose is defined globally and nonzero
 * you can watch it run.  I get values of:
 *   Found biggest double 8.9884656743115795e+307
 *   Found biggest float 1.7014118599577043e+38
 * on a PIII.  Eyeballing the convergence, 16 digits or thereabouts
 * are significant for a double, 8 or thereabouts are significant for a
 * float (the "8" of "118599..." is possibly significant, the rest
 * following is definitely not, ditto the "9" in the trailing "795" in
 * in the double result).
 *========================================================================
 */

#define _ISOC99_SOURCE
#include <math.h>
extern int verbose;

double maxdouble()
{

 double d,dold,dfrac,dnew;
 int i;

 d = 1.0;
 i = 0;

 /*
  * This routine finds the order of infinity (with a factor of two).
  * dold contains the last number that was NOT infinite.
  */
 while(fpclassify(d) != FP_INFINITE){
   if(verbose) printf("d[%d] = %21.16e\n",i,d);
   dold = d;
   d *= 2.0;
   i++;
 }

 /*
  * OK, dold contains the biggest number so far that isn't infinite.
  * We now embark on a binomial search for the boundary number by
  * subdividing d until the fraction added doesn't change the number
  * itself (and it still isn't infinite).
  */
 dfrac = dold;
 while(dnew != dold){
   dold = dnew;
   dfrac /= 2.0;
   d = dold + dfrac;
   /* If result isn't infinite, accept it and continue up */
   if(fpclassify(d) != FP_INFINITE) {
     dnew = d;
     i++;
     if(verbose) printf("d[%d] = %21.16e, dold = %21.16e\n",i,d,dold);
   }
 }
 if(verbose) printf("Biggest double that isn't infinite: %21.16e\n",d);

 return(d);

}

float maxfloat()
{

 float d,dold,dfrac,dnew;
 int i;

 d = 1.0;
 i = 0;

 /*
  * This routine finds the order of infinity (with a factor of two).
  * dold contains the last number that was NOT infinite.
  */
 while(fpclassify(d) != FP_INFINITE){
   if(verbose) printf("d[%d] = %21.16e\n",i,d);
   dold = d;
   d *= 2.0;
   i++;
 }

 /*
  * OK, dold contains the biggest number so far that isn't infinite.
  * We now embark on a binomial search for the boundary number by
  * subdividing d until the fraction added doesn't change the number
  * itself (and still isn't infinite).
  */
 dfrac = dold;
 while(dnew != dold){
   dold = dnew;
   dfrac /= 2.0;
   d = dold + dfrac;
   /* If result isn't infinite, accept it and continue up */
   if(fpclassify(d) != FP_INFINITE) {
     dnew = d;
     i++;
     if(verbose) printf("d[%d] = %21.16e, dold = %21.16e\n",i,d,dold);
   }
 }
 if(verbose) printf("Biggest float that isn't infinite: %21.16e\n",d);

 return(d);

}