The R Project SVN R

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2002--2025 The R Core Team
 *
 *  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, a copy is available at
 *  https://www.R-project.org/Licenses/
 *
 * The interv() used to be Fortran in ../library/modreg/src/bvalue.f
 * and part of Hastie and Tibshirani's public domain GAMFIT package.
 *
 * Translated by f2c (version 20010821), cleaned up and extended by
 * Martin Maechler.
 */

#include <R_ext/Applic.h>
#include <R_ext/Boolean.h>
#include <R_ext/Utils.h>

/* Was in API header R_ext/Utils.h

   This is called from Fortran in stats/src/bvalue.f,
   3x stats/src/s*.f for smooth.spline() 
   and packages gam and mda.

   Was Rboolean *. but that is not possible from Fortran.  Since it is
   inteneded to be called from C it should not be in a C header.
*/
int F77_SUB(interv)(double *xt, int *n, double *x,
            int *rightmost_closed, int *all_inside,
            int *ilo, int *mflag)
{
    // This has an int return value which is ignored here
    return findInterval(xt, *n, *x,
            (Rboolean) *rightmost_closed,
            (Rboolean) *all_inside,
            *ilo, mflag);
}

/* In API header R_ext/Utils.h */
int findInterval2(double *xt, int n, double x,
          Rboolean rightmost_closed, Rboolean all_inside,
          Rboolean left_open, // <- new in findInterval2()
          int ilo, int *mflag)
{
    int istep, middle, ihi;

/* computes  `left' := max( i ; 1 <= i <= n   &&  xt[i] <= x )  .

 ******  i n p u t  ******

  xt    numeric vector of length  n , assumed to be nondecreasing
  n length(xt)
  x the point whose location with respect to the sequence  xt  is
    to be determined.
  rightmost_closed {logical} indicating if the rightmost xt[] interval
    should be closed, i.e. result := n-1 if x == x[n]
    (when left_open, the *leftmost* interval should be closed.)
  all_inside {logical} indicating if result should be coerced
        to lie inside {1, n-1}
  left_open  {logical} use intervals (s, t] instead of [s, t)
  ilo   typically the result of the last call to findInterval(.)
    `ilo' used to be a static variable (in Fortran) which is not
    desirable in R anymore (threads!).
    Instead, you *should* use a reasonable value, in the first call.

 ******  o u t p u t  ******

  left, mflag  both integers, whose value is

   0     -1      if            x <  xt[1]
   i      0      if  xt[i]  <= x <  xt[i+1]
   n      1      if  xt[n]  <= x

    in particular,  mflag = 0 is the 'usual' case.  mflag != 0
    indicates that  x  lies outside the halfopen interval
    xt[1] <= y < xt[n] . the asymmetric treatment of the
    interval is due to the decision to make all pp functions cont-
    inuous from the right.

   Note that if all_inside, left is 1 instead of 0 and n-1 instead of n;
   and if rightmost_closed and x == xt[n],  left is    n-1 instead of n.


 ******  m e t h o d  ******

  the program is designed to be efficient in the common situation that
  it is called repeatedly, with  x  taken from an increasing or decreasing
  sequence. this will happen, e.g., when a pp function is to be graphed.
  The first guess for  left  is therefore taken to be the value returned at
  the previous call and stored in the  l o c a l   variable  ilo .
  a first check ascertains that  ilo < n (this is necessary since the
  present call may have nothing to do with the previous call).

  then, if  xt[ilo] <= x < xt[ilo+1], we set  left = ilo
  and are done after just three comparisons.
  otherwise, we repeatedly double the difference  istep = ihi - ilo
  while also moving  ilo  and  ihi  in the direction of  x , until
              xt[ilo] <= x < xt[ihi] ,
  after which we use bisection to get, in addition, ilo+1 = ihi .
  left = ilo  is then returned.
*/

#define left_boundary  { *mflag = -1; \
    return((all_inside || (rightmost_closed && x == xt[1])) ? 1 : 0); }

#define right_boundary { *mflag = +1;                   \
    return((all_inside || (rightmost_closed && x == xt[n]))     \
        ? (n - 1) : n); }

#define X_grtr(XT_v) x > XT_v || (!left_open && x >= XT_v)
#define X_smlr(XT_v) x < XT_v ||  (left_open && x <= XT_v)

    if(n == 0) { *mflag = 0 ; return 0; }
    /* 1-indexing : */
    --xt;

    if(ilo <= 0) {
    if (X_smlr(xt[1]))      left_boundary;
    ilo = 1;
    }
    ihi = ilo + 1;
    if (ihi >= n) {
    if (X_grtr(xt[n]))      right_boundary;
    if (n <= 1) /* x < xt[1] */ left_boundary;
    ilo = n - 1;
    ihi = n;
    }

    if (X_smlr(xt[ihi])) {
    if (X_grtr(xt[ilo])) {
        /* `lucky': same interval as last time */
        *mflag = 0;    return ilo;
    }
    /* **** now x < xt[ilo] .   decrease  ilo  to capture  x */
    if(!left_open) for(istep = 1; ; istep *= 2) {
        ihi = ilo;
        ilo = ihi - istep;
        if (ilo <= 1)
        break;
        if (x >= xt[ilo])       goto L50;
    } else for(istep = 1; ; istep *= 2) {
        ihi = ilo;
        ilo = ihi - istep;
        if (ilo <= 1)
        break;
        if (x > xt[ilo])        goto L51;
    }
    ilo = 1;
    if (X_smlr(xt[1]))      left_boundary;
    }
    else {
    /* **** now x >= xt[ihi] .  increase  ihi  to capture  x */
    if(!left_open) for(istep = 1; ; istep *= 2) {
        ilo = ihi;
        ihi = ilo + istep;
        if (ihi >= n)
        break;
        if (x < xt[ihi])        goto L50;
    }
    else for(istep = 1; ; istep *= 2) {
        ilo = ihi;
        ihi = ilo + istep;
        if (ihi >= n)
        break;
        if (x <= xt[ihi])       goto L51;
    }
    if (X_grtr(xt[n]))      right_boundary;
    ihi = n;
    }
    
    if (left_open) goto L51; /* There _is_ a path to here, avoiding return and goto */

L50: // ! left_open
    /* **** now xt[ilo] <= x < xt[ihi] . narrow the interval. */
    for(;;) {
    middle = (ilo + ihi) / 2;
    if (middle == ilo) {
        *mflag = 0;    return ilo;
    }
    /* note. it is assumed that middle = ilo in case ihi = ilo+1 . */
    if (x >= xt[middle])
        ilo = middle;
    else
        ihi = middle;
    }

L51: // left_open
    /* **** now xt[ilo] < x <= xt[ihi] . narrow the interval. */
    for(;;) {
    middle = (ilo + ihi) / 2;
    if (middle == ilo) {
        *mflag = 0;    return ilo;
    }
    /* note. it is assumed that middle = ilo in case ihi = ilo+1 . */
    if (x > xt[middle])
        ilo = middle;
    else
        ihi = middle;
    }
} /* findInterval2 */

// has been in API -- keep for compatibility:
int findInterval(double *xt, int n, double x,
         Rboolean rightmost_closed,  Rboolean all_inside, int ilo,
         int *mflag)
{
    return findInterval2(xt, n, x, rightmost_closed, all_inside, FALSE, ilo, mflag);
}