[comp.lang.fortran] rounding a real to a whole number

u714092@eagle.larc.nasa.gov (prichard devon ) (11/02/90)

what I want to do is to have a subroutine which, given a range of real
numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
outside that range. obviously the first application is axis values for
x-y plots.  optionally, picking a 10-multiple of whole numbers ought to
be included; for  -1.45567 -> 3.34112 , it should return -1.4 and 3.4 ...

while this seems so simple, I don't see a way to do this reasonably with
Fortran intrinsics.  I just can't get it to work correctly with truncations.

anyone done this sort of thing listening out there?

--
 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 | Devon Prichard             making the world safe for helicopters ... |
 | u714092@eagle.larc.nasa.gov, prichard@ias.larc.nasa.gov              |
 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

buckland@cheddar.ucs.ubc.ca (Tony Buckland) (11/02/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:
|
|what I want to do is to have a subroutine which, given a range of real
|numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
|outside that range. obviously the first application is axis values for
|x-y plots.  optionally, picking a 10-multiple of whole numbers ought to
|be included; for  -1.45567 -> 3.34112 , it should return -1.4 and 3.4 ...
|
|while this seems so simple, I don't see a way to do this reasonably with
|Fortran intrinsics.  I just can't get it to work correctly with truncations.
|
|anyone done this sort of thing listening out there?

 I imagine you meant -1.5, not -1.4.
 That quibble aside, a traditional way of rounding up is to add
 a floating-point quantity before truncating to integer.  Watch
 what happens in this example:
 
    f = 3.34112                    unrounded quantity
    f = f + 0.09999                add some variation on "just
                                   less than one"
    f = 10.0 * f                   we want to keep two digits
                                   f is now 34.4111
    i = f                          i is now 34
    f = i                          f is now 34.0000
    f = f / 10.0                   f is now 3.4
 
 The negative case is left as an exercise for the student.
 The problem with inexact representation of floating-point
 machine representation of most non-integers is ignored.

merritt@buphy.BU.EDU (Sean Merritt) (11/02/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:
>
>what I want to do is to have a subroutine which, given a range of real
>numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
>outside that range. obviously the first application is axis values for
>x-y plots.  optionally, picking a 10-multiple of whole numbers ought to
>be included; for  -1.45567 -> 3.34112 , it should return -1.4 and 3.4 ...
                    ^^^^^^^
		     -1.5		
>Fortran intrinsics.  I just can't get it to work correctly with truncations.
>
>anyone done this sort of thing listening out there?
>
>--
> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
> | Devon Prichard             making the world safe for helicopters ... |
> | u714092@eagle.larc.nasa.gov, prichard@ias.larc.nasa.gov              |
> ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||


 This is how I've always done it:


	A = number to be rounded

	B = number of decimal places rounded to

	C = rounded number

	Simply include this line in a subroutine or function :

	C = INT(A*10**B+0.5)/10**B


	-sjm
_________________________________________________________________________
Sean J. Merritt			  |	 
Dept. of Physics,Boston University|	SURSUM CORDA !
e-mail to:merritt@buphy.bu.edu	  |	

ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe) (11/02/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov>, u714092@eagle.larc.nasa.gov (prichard devon ) writes:
> what I want to do is to have a subroutine which, given a range of real
> numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
> outside that range. obviously the first application is axis values for
> x-y plots.

Have you tried the ACM algorithms collection?  I'm sure there's something
there.  Look through old issues of Computer Journal.  There's a sampling
of algorithms from the journal Applied Statistics has been published as a
book.  Most of my library is in another city, but there's a subroutine for
finding "nice numbers" in
	Applications, Basics, and Computing of
	Exploratory Data Analysis;
	Velleman & Hoaglin
	Duxbury Press
The subroutine NPOSW on page 316 is good for picking axis intervals.
Many books providing statistical algorithms will offer something similar.

-- 
The problem about real life is that moving one's knight to QB3
may always be replied to with a lob across the net.  --Alasdair Macintyre.

ok@goanna.cs.rmit.oz.au (Richard A. O'Keefe) (11/02/90)

In article <67611@bu.edu.bu.edu>, merritt@buphy.BU.EDU (Sean Merritt) writes:
> In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:
> >what I want to do is to have a subroutine which, given a range of real
> >numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
> >outside that range.

> 	C = INT(A*10**B+0.5)/10**B

C The original poster wanted something that rounds AWAY from zero.
      SUBROUTINE IRANGE(XLO, XHI, ILO, IHI, IERROR)
C Imagine we are given
	 REAL XLO, XHI
C and are required to find
	 INTEGER ILO, IHI
C such that
C	 ILO.LE.XLO .AND. XHI.LE.IHI

C Here's how.
	 IF (XLO .GT. XHI) THEN
	    IERROR = 1
	 ELSE IF (XHI .GT. MAXINT) THEN
C	    Where do you get MAXINT?  That's _your_ business.
	    IERROR = 2
	 ELSE IF (XLO .LT. MININT) THEN
C	    Where do you get MININT?  Try MACHAR.
	    IERROR = 3
	 ELSE
	    IHI = INT(XHI)
	    IF (IHI .LT. XHI) IHI = IHI+1
	    ILO = INT(XLO)
	    IF (ILO .GT. XLO) ILO = ILO-1
	    IERROR = 0
	 END IF
      END

Crude.  Totally lacking in imagination.  But it _does_ "choose whole
numbers _outside_ the range" as requested.

-- 
The problem about real life is that moving one's knight to QB3
may always be replied to with a lob across the net.  --Alasdair Macintyre.

dodson@convex.COM (Dave Dodson) (11/03/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:
>what I want to do is to have a subroutine which, given a range of real
>numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
>outside that range. obviously the first application is axis values for
>x-y plots.  optionally, picking a 10-multiple of whole numbers ought to
>be included; for  -1.45567 -> 3.34112 , it should return -1.4 and 3.4 ...

If you are going to plot data, you may want the axis divisions to be "round"
numbers.  I wrote a subroutine many years ago that determined such an axis.
The following is how I remember it working.  Caution: not checked out...

	subroutine axis (n,x,ndiv,xmin,delx,xmax)
c  input:
	integer n	! length of x
	integer ndiv	! number of axis divisions
	real x(n)	! data array
c  output:
	real xmin	! axis minimum, multiple of delx
	real delx	! axis increment per division
	real xmax	! axis maximum, multiple of delx
c  local:
	real table(4)
	data table / 1.0, 2.0, 4.0, 5.0 /
	
	xmini = x(1)
	xmaxi = x(1)
	do 10 i = 2, n
	    xmini = min ( xmini, x(i)
	    xmaxi = max ( xmaxi, x(i)
   10	continue
	if ( xmini .eq. xmaxi ) then
	    < I don't remember what I did when all x(i) are equal >
	end if

	delta = ( xmaxi - xmini ) / ndiv
	ratio = 0.0
	do 40 i = 1, 4
	    trial = table(i)
   20	    if ( trial .lt. delta ) then
		trial = 10.0 * trial
		go to 20
	    end if
   30	    if ( trial / 10.0 .gt. delta ) then
		trial = trial / 10.0
		go to 30
	    end if
	    xmint = trial * aint ( xmini / trial )
	    if ( xmint .gt. xmini ) xmint = xmint - trial
	    xmaxt = trial * aint ( xmaxi / trial )
	    if ( xmaxt .lt. xmaxi ) xmaxt = xmaxt + trial
	    if ( xmint + ndiv * trial * 1.0001 .ge. xmaxt ) then
		if ( ( xmaxi - xmini ) / ( xmaxt - xmint ) .gt. ratio ) then
		    ratio = ( xmaxi - xmini ) / ( xmaxt - xmint )
		    xmin = xmint
		    delx = trial
		    xmax = xmaxt
		end if
	    end if
   40	continue
	return
	end

This begins by finding the max and min of the data.  Delta is set to
the minimum allowed increment per axis division.  Ratio is the fraction
of the axis that would be occupied by the graph; it is initialized to
zero.  Then each increment of the form 1 * 10^j, 2 * 10^j, and 4 * 10^j
is examined.  First the j in 10^j is adjusted up or down so that 
table(i) * 10^(j-1) <= delta <= table(i) * 10^j (actually, we don't care
about j, only trial := table(i) * 10^j).  Using trial as an increment
per axis division, we find xmint = the largest multiple of trial that
is .le. xmini and xmaxt = the smallest multiple of trial that is .ge. xmaxt.
If ndiv * trial spans the interval from xmint to xmaxt then we check to see
how much of the axis would be occupied by the graph.  We accept trial if
it uses more axis that the previous.  When we have inspected all 4 forms
for trial, we return with the best one.

----------------------------------------------------------------------

Dave Dodson		                             dodson@convex.COM
Convex Computer Corporation      Richardson, Texas      (214) 497-4234

whit@milton.u.washington.edu (John Whitmore) (11/09/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:
>
>what I want to do is to have a subroutine which, given a range of real
>numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
>outside that range. obviously the first application is axis values for
>x-y plots.

	I had the same problem.  I don't like my solution, but...
I treated the axis endpoint as a range (obviously it has to include
the lowest data value, but no more than 20% of the range of the
plot should be empty.)  Then I took the axis endpoint and separated
it into mantissa_MSD, mantissa_LSD, and exponent.  'MSD' is for
'most significant digits', and 'LSD' for 'least significant digits.'
With the entirety of the value of the number in 'mantissa_LSD', 
and 'mantissa_MSD'= 0, one tests to see if the 'mantissa_MSD' and
exponent are sufficient (i.e. if that number is in the determined
range).  
	Every time the test fails, one takes one more digit from
'mantissa_LSD' and appends it to 'mantissa_MSD'.  Eventually one
comes up with some sort of success; the successful number will
be 'mantissa_MSD' or mantissa_MSD +/-1 times 'exponent'; one
discards the remaining digits.
	Additionally, one can build in a preference for even
numbers by simply adding to a table of the variations added to
'mantissa_MSD' the values 0.2, 0.4, 0.6, 0.8.

	Code follows (skip if you're not interested).

        SUBROUTINE BESCL9(X1IN,X2IN,RESULT,NTRIAL)
C
C       THIS SUBROUTINE FINDS THE MOST ROUNDED NUMBER IT CAN, IN
C       THE RANGE X TO XLIMIT, RETURNING THE ROUNDEST NUMBER IN
C       RESULT.  THE 'ROUNDNESS' CRITERION IS: FEWEST NONZERO DIGITS,
C       WITH PREFERENCE FOR LAST NONZERO DIGIT '5', OR ANY EVEN DIGIT.
C
C       X1IN, X2IN ARE REAL VARIABLES, INPUT, WHICH DELIMIT
C               THE RANGE IN WHICH A NUMBER IS TO BE FOUND
C       RESULT  REAL, OUTPUT, IS THE ROUND NUMBER RESULT
C       NTRIAL  INTEGER, OUTPUT, INDICATES THE NUMBER OF FRUITLESS ATTEMPTS
C               ATTEMPTS MADE TO MAKE THE NUMBER; MAY BE USEFUL IN
C               QUANTIFYING THE DEGREE OF ROUNDING ATTAINED
C
C               PREFERRED DIGITS ARE TABULATED IN ARRAY "BESTUN"
        PARAMETER(NB=7)
        REAL BESTUN(NB)
        DATA BESTUN/  0., 1., .5, .2, .4, .6, .8 /
C
        DATA TEN,ONE,TENTH/10., 1., .1/
C
        NTRIAL=0
        IF(X1IN .EQ. X2IN) THEN
C
C                       NO ROUNDING POSSIBLE OR REQUIRED
                RESULT=X1IN
                GO TO 5
        ELSE IF( X1IN .LT. 0. .EQV. X2IN .GT. 0.) THEN
C
C                       ZERO IS A ROUND NUMBER
                RESULT=0.
                GO TO 5
        ELSE IF(X1IN .GT. X2IN .EQV. X1IN .GT. 0.) THEN
C
C                       SEARCH 'UPWARD' FROM X1IN
                X=X2IN
                XLIMIT= X1IN
        ELSE
C
C                       SEARCH 'UPWARD' FROM X2IN
                X=X1IN
                XLIMIT=X2IN
        ENDIF
C
C       THE TRIPLE (P1,P2,P3) HOLDS THE TRIAL X ELEMENTS AS:
C       P1= MANTISSA'S LEAST SIGNIFICANT DIGITS
C       P2= POWER OF TEN, FACTORED FROM THE MANTISSA, WITH SIGN OF X
C       P3= MOST SIGNIFICANT DIGITS OF MANTISSA, WHICH ARE
C               HOPELESS TO TRY TO DO WITHOUT
C
        P1=ABS(X)
        P2=SIGN(ONE,X)
        P3=0.
C
C               MULTIPLY OR DIVIDE P2 BY POWER OF TEN TO PUT
C                P1 INTO THE INTERVAL (.1,1.]
C
 1      CONTINUE
        IF(P1 .LT. TENTH) THEN
                P1=P1*TEN
                P2=P2*TENTH
                GO TO 1
        ENDIF
 11     CONTINUE
        IF(P1 .GT. ONE) THEN
                P1=P1*TENTH
                P2=P2*TEN
                GO TO 11
        ENDIF
 2      CONTINUE
C
C               NOW TEST FOR ADEQUATE NEARNESS TO BEST NUMBERS
C               IF FIND LIKELY VALUE IN TABLE, WILL USE TABLE VALUE
C               AS A ROUNDING APPROXIMATION TO P1
        XLESS= MIN(X1IN,X2IN)
        XGRTR= MAX(X1IN,X2IN)
        DO 3 N=1,NB
                NTRIAL = NTRIAL + 1
                RESULT= (BESTUN(N)+P3) *P2
 3              IF( XLESS .LE. RESULT .AND. XGRTR .GE. RESULT) GO TO 5
C
C               THE SEARCH WAS FRUITLESS; PUSH DIGITS TO P3 AND TRY AGAIN
C
        A=   AINT(P1*TEN)
        P1=  P1*TEN-A
        P2=  P2*TENTH
        P3=  P3*TEN+A
        GO TO 2
C
 5      RETURN
        END


	John Whitmore
	Whit@milton.u.washington.edu

jmb@ukc.ac.uk (Mike Bremner) (11/10/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:

>what I want to do is to have a subroutine which, given a range of real
>numbers, for example  -123.445 -> 1033.02 , and choose whole numbers 
>outside that range. obviously the first application is axis values for
>x-y plots.

Applied Statistics algorithm AS 168 - Vol. 30 (1981), p.339 - deals with
scale selection for scatter plots, etc.

jlg@lanl.gov (Jim Giles) (11/10/90)

In article <U714092.90Nov1101629@eagle.larc.nasa.gov> u714092@eagle.larc.nasa.gov (prichard devon ) writes:
> [...]
> what I want to do is to have a subroutine which, given a range of real
> numbers, for example  -123.445 -> 1033.02 , and choose whole numbers
> outside that range. obviously the first application is axis values for
> x-y plots.

How about:

      Subroutine ticks(lbound,hbound,bottom,top,tick)
C       takes low and high bound of interval and returns
C       a 'tick' value for the interval between tick marks
C       on a graph axis.  The value of tick is always
C       either a power of 10, twice a power of ten or
C       five times a power of ten.  A new interval is
C       also returned which begins below the original
C       and ends above the original and is an integral
C       number of 'tick's long.  The number of ticks in
C       the new interval is (about) 5 to 13.
      real lbound,hbound,bottom,top,tick
      real lowtic, basetk
      integer tenp
      lowtic = (hbound-lbound)/12.5
      tenp=int(alog10(lowtick))
      if (tenp.gt.alog10(lowtick)) tenp=tenp-1
      basetk = 10**tenp
      if (tick.lt.lowtick) tick = 2 * basetk
      if (tick.lt.lowtick) tick = 5 * basetk
      if (tick.lt.lowtick) tick = 10 * basetk
      bottom = nint(lbound/tick)*tick
      if (bottom.gt.lbound) bottom = bottom - tick
      top = nint(hbound/tick)*tick
      if (top.lt.hbound) top = top + tick
      return

Given your numbers: -123.445 -> 1033.02, this routine returns
TICK equal to 100.0, BOTTOM equal to -200.0, and TOP equal to 1100.0.

The number of ticks in the result interval is always less than the
divisor for LOWTICK plus one (in this case, the divisor is 12.5,
so the maximum number of ticks is 13).  The minimum number of ticks
in the final interval is the LOWTICK divisor itself divided by 2.5.
So, if you want finer resolution, pick a number larger than 12.5.
If you want coarse resolution, pick a smaller number to divide by.
You can even make the divisor one of the subroutine arguments if
you want.  Given your numbers:-123.445 -> 1033.02, and a divisor
of 25, the answer is - TICK = 50.0, BOTTOM = -150.0, and TOP = 1050.0
giving 24 intervals.

Not elegant, but not bad off the top of my head.

J. Giles