/*-----------------------------------------------------------------------------

Copyright (C) 2016

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

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

Function      rb2 - returns the rational base 2 number closest to x

Syntax        #include "sclfuncs.h"
              REAL rb2(REAL x)

Prototype in  sclfuncs.h

Description   x is any REAL on entry, positive or negative. The return
              value has the same sign as x. 

Remark        Invalid input is handled as follows:
              If x is larger than 2*REAL_MAX, then 2*REAL_MAX is returned.  
	      If x is smaller than -2*REAL_MAX, then -2*REAL_MAX is returned.
	      If scl::IsFinite(x) returns false, then x is returned.
	      
Return value  The nearest rational base 2 number. If x is halfway between 
              the two nearest rational base 2 numbers, the one farther from 
	      zero is returned..
	      
Functions     Library: (none)
called        libscl:  IsFinite

------------------------------------------------------------------------------*/

#include "sclfuncs.h"

namespace scl {

  REAL rb2(REAL x)
  {
    REAL rv = x;

    REAL sign = x < 0 ? -1.0 : 1.0;
    x *= sign;

    if (x == 0.0) return rv;

    const REAL upper_bound = REAL_MAX/2.0;

    if (x >= upper_bound) return sign*upper_bound;
    if (!IsFinite(x)) return rv;

    REAL xsave = x;

    rv = 1.0;

    if (x < 1.0) {
      REAL left = 1.0;
      while (x < 1.0) {
        x *= 2.0;
        left /= 2.0;
      }
      REAL right = 2.0*left;
      rv = (xsave - left) < (right - xsave) ? left : right;
    }
    else if (x > 1.0) {
      REAL right = 1.0;
      while (x > 1.0) {
        x /= 2.0;
        right *= 2.0;
      }
      REAL left = right/2.0;
      rv = (xsave - left) < (right - xsave) ? left : right;
    }

    rv *= sign;

    return rv;
  }

}
