File Coverage

lib/PDL/Primitive-pp-pchip_bvalu.c

Criterion	Covered	Total	%
statement	51	75	68.0
branch	47	246	19.1
condition			n/a
subroutine			n/a
pod			n/a
total	98	321	30.5

#line 453 "lib/PDL/PP.pm"
/*
 * THIS FILE WAS GENERATED BY PDL::PP from lib/PDL/Primitive.pd! Do not modify!
 */

#define PDL_FREE_CODE(trans, destroy, comp_free_code, ntpriv_free_code) \
    if (destroy) { \
	comp_free_code \
    } \
    if ((trans)->dims_redone) { \
	ntpriv_free_code \
    }

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "pdl.h"
#include "pdlcore.h"
#define PDL PDL_Primitive
extern Core* PDL; /* Structure hold core C functions */
#line 23 "lib/PDL/Primitive-pp-pchip_bvalu.c"
extern int pdl_srand_threads;
extern uint64_t *pdl_rand_state;
void pdl_srand(uint64_t **s, uint64_t seed, int n);
double pdl_drand(uint64_t *s);
#define PDL_MAYBE_SRAND \
  if (pdl_srand_threads < 0) \
    pdl_srand(&pdl_rand_state, PDL->pdl_seed(), PDL->online_cpus());
#define PDL_RAND_SET_OFFSET(v, thr, pdl) \
  if (v < 0) { \
    if (thr.mag_nthr >= 0) { \
      int thr_no = PDL->magic_get_thread(pdl); \
      if (thr_no < 0) return PDL->make_error_simple(PDL_EFATAL, "Invalid pdl_magic_get_thread!"); \
      v = thr_no == 0 ? thr_no : thr_no % PDL->online_cpus(); \
    } else { \
      v = 0; \
    } \
  }

#line 1857 "lib/PDL/PP.pm"
pdl_error pdl_pchip_bvalu_redodims(pdl_trans *__privtrans) {
  pdl_error PDL_err = {0, NULL, 0};
#line 45 "lib/PDL/Primitive-pp-pchip_bvalu.c"
  __privtrans->ind_sizes[0] = 3*(__privtrans->ind_sizes[2]-__privtrans->ind_sizes[1]);
  #ifndef PDL_DECLARE_PARAMS_pchip_bvalu_0
  #define PDL_DECLARE_PARAMS_pchip_bvalu_0(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ideriv,PDL_PPSYM_PARAM_ideriv,PDL_TYPE_PARAM_inbv,PDL_PPSYM_PARAM_inbv) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, t, (__privtrans->pdls[0]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, a, (__privtrans->pdls[1]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ideriv, ideriv, (__privtrans->pdls[2]), 0, PDL_PPSYM_PARAM_ideriv) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[3]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, ans, (__privtrans->pdls[4]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_inbv, inbv, (__privtrans->pdls[5]), 0, PDL_PPSYM_PARAM_inbv) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, work, (__privtrans->pdls[6]), 0, PDL_PPSYM_OP)
  #endif
  #define PDL_IF_BAD(t,f) f
  switch (__privtrans->__datatype) { /* Start generic switch */
    case PDL_F: {
      PDL_DECLARE_PARAMS_pchip_bvalu_0(PDL_Float,F,PDL_Indx,N,PDL_Indx,N)
        {PDL_Indx k = __privtrans->ind_sizes[2] - __privtrans->ind_sizes[1];
        if (k < 1)        return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "K DOES NOT SATISFY K.GE.1");
        if (__privtrans->ind_sizes[1] < k) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "N DOES NOT SATISFY N.GE.K");
        }
    } break;
    case PDL_D: {
      PDL_DECLARE_PARAMS_pchip_bvalu_0(PDL_Double,D,PDL_Indx,N,PDL_Indx,N)
        {PDL_Indx k = __privtrans->ind_sizes[2] - __privtrans->ind_sizes[1];
        if (k < 1)        return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "K DOES NOT SATISFY K.GE.1");
        if (__privtrans->ind_sizes[1] < k) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "N DOES NOT SATISFY N.GE.K");
        }
    } break;
    case PDL_LD: {
      PDL_DECLARE_PARAMS_pchip_bvalu_0(PDL_LDouble,E,PDL_Indx,N,PDL_Indx,N)
        {PDL_Indx k = __privtrans->ind_sizes[2] - __privtrans->ind_sizes[1];
        if (k < 1)        return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "K DOES NOT SATISFY K.GE.1");
        if (__privtrans->ind_sizes[1] < k) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "N DOES NOT SATISFY N.GE.K");
        }
    } break;
    default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_bvalu: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
  }
  #undef PDL_IF_BAD

  PDL_RETERROR(PDL_err, PDL->redodims_default(__privtrans));
  return PDL_err;
}


#line 1857 "lib/PDL/PP.pm"
pdl_error pdl_pchip_bvalu_readdata(pdl_trans *__privtrans) {
  pdl_error PDL_err = {0, NULL, 0};
#line 92 "lib/PDL/Primitive-pp-pchip_bvalu.c"
  register PDL_Indx __k3_size = __privtrans->ind_sizes[0];
  if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "broadcast.incs NULL");
  /* broadcastloop declarations */
  int __brcloopval;
  register PDL_Indx __tind0,__tind1; /* counters along dim */
  register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
  /* dims here are how many steps along those dims */
  register PDL_Indx __tinc0_t = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
  register PDL_Indx __tinc0_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
  register PDL_Indx __tinc0_ideriv = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
  register PDL_Indx __tinc0_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
  register PDL_Indx __tinc0_ans = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
  register PDL_Indx __tinc0_inbv = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
  register PDL_Indx __tinc0_work = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,0);
  register PDL_Indx __tinc1_t = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
  register PDL_Indx __tinc1_a = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
  register PDL_Indx __tinc1_ideriv = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
  register PDL_Indx __tinc1_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
  register PDL_Indx __tinc1_ans = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
  register PDL_Indx __tinc1_inbv = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
  register PDL_Indx __tinc1_work = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,1);
  #define PDL_BROADCASTLOOP_START_pchip_bvalu_readdata PDL_BROADCASTLOOP_START( \
    readdata, \
    __privtrans->broadcast, \
    __privtrans->vtable, \
    t_datap += __offsp[0]; \
    a_datap += __offsp[1]; \
    ideriv_datap += __offsp[2]; \
    x_datap += __offsp[3]; \
    ans_datap += __offsp[4]; \
    inbv_datap += __offsp[5]; \
    work_datap += __offsp[6]; \
    , \
    (  ,t_datap += __tinc1_t - __tinc0_t * __tdims0 \
    ,a_datap += __tinc1_a - __tinc0_a * __tdims0 \
    ,ideriv_datap += __tinc1_ideriv - __tinc0_ideriv * __tdims0 \
    ,x_datap += __tinc1_x - __tinc0_x * __tdims0 \
    ,ans_datap += __tinc1_ans - __tinc0_ans * __tdims0 \
    ,inbv_datap += __tinc1_inbv - __tinc0_inbv * __tdims0 \
    ,work_datap += __tinc1_work - __tinc0_work * __tdims0 \
    ), \
    (  ,t_datap += __tinc0_t \
    ,a_datap += __tinc0_a \
    ,ideriv_datap += __tinc0_ideriv \
    ,x_datap += __tinc0_x \
    ,ans_datap += __tinc0_ans \
    ,inbv_datap += __tinc0_inbv \
    ,work_datap += __tinc0_work \
    ) \
  )
  #define PDL_BROADCASTLOOP_END_pchip_bvalu_readdata PDL_BROADCASTLOOP_END( \
    __privtrans->broadcast, \
    t_datap -= __tinc1_t * __tdims1 + __offsp[0]; \
    a_datap -= __tinc1_a * __tdims1 + __offsp[1]; \
    ideriv_datap -= __tinc1_ideriv * __tdims1 + __offsp[2]; \
    x_datap -= __tinc1_x * __tdims1 + __offsp[3]; \
    ans_datap -= __tinc1_ans * __tdims1 + __offsp[4]; \
    inbv_datap -= __tinc1_inbv * __tdims1 + __offsp[5]; \
    work_datap -= __tinc1_work * __tdims1 + __offsp[6]; \
  )
  register PDL_Indx __inc_a_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_a_n;
  register PDL_Indx __inc_t_nplusk = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_t_nplusk;
  register PDL_Indx __inc_work_k3 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,6,0)]; (void)__inc_work_k3;
  #ifndef PDL_DECLARE_PARAMS_pchip_bvalu_1
  #define PDL_DECLARE_PARAMS_pchip_bvalu_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ideriv,PDL_PPSYM_PARAM_ideriv,PDL_TYPE_PARAM_inbv,PDL_PPSYM_PARAM_inbv) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, t, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, a, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ideriv, ideriv, (__privtrans->pdls[2]), 1, PDL_PPSYM_PARAM_ideriv) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, ans, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_inbv, inbv, (__privtrans->pdls[5]), 1, PDL_PPSYM_PARAM_inbv) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, work, (__privtrans->pdls[6]), 1, PDL_PPSYM_OP)
  #endif
  #define PDL_IF_BAD(t,f) f
  switch (__privtrans->__datatype) { /* Start generic switch */
    case PDL_F: {
      PDL_DECLARE_PARAMS_pchip_bvalu_1(PDL_Float,F,PDL_Indx,N,PDL_Indx,N)
      PDL_BROADCASTLOOP_START_pchip_bvalu_readdata    {
          #line 6145 "lib/PDL/Primitive.pd"
          PDL_Indx k = __privtrans->ind_sizes[2] - __privtrans->ind_sizes[1];
          if ((ideriv_datap)[0] < 0 || (ideriv_datap)[0] >= k)
            return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K");
          PDL_Indx i;
          int mflag;
          /* *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1) */
          /*   (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)). */
          do { /* inlined dintrv */
            PDL_Indx ihi = (inbv_datap)[0] + 1, lxt = __privtrans->ind_sizes[1];
            if (ihi >= lxt) {
              if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(lxt))]) {
                mflag = 1; i = lxt; break;
              }
              if (lxt <= 0) {
                mflag = -1; i = 0; break;
              }
              ihi = (inbv_datap)[0] = lxt;
            }
            char skipflag = 0;
            if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(ihi))]) {
              PDL_Indx inbv = (inbv_datap)[0];
              if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(inbv))]) {
                mflag = 0; i = (inbv_datap)[0]; break;
              }
          /* *** NOW X .LT. XT(IHI) . FIND LOWER BOUND */
              PDL_Indx istep = 1;
              while (1) {
                ihi = (inbv_datap)[0];
                (inbv_datap)[0] = ihi - istep;
                if ((inbv_datap)[0] <= 0) {
                  break;
                }
                PDL_Indx inbv = (inbv_datap)[0];
                if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(inbv))]) {
                  skipflag = 1;
                  break;
                }
                istep <<= 1;
              }
              if (!skipflag) {
                (inbv_datap)[0] = 0;
                if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(0))]) {
                  mflag = -1; i = 0; break;
                }
              }
              skipflag = 1;
          /* *** NOW X .GE. XT(ILO) . FIND UPPER BOUND */
            }
            if (!skipflag) {
              PDL_Indx istep = 1;
              while (1) {
                (inbv_datap)[0] = ihi;
                ihi = (inbv_datap)[0] + istep;
                if (ihi >= lxt) break;
                if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(ihi))]) {
                  skipflag = 1;
                  break;
                }
                istep <<= 1;
              }
              if (!skipflag) {
                if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(lxt))]) {
                  mflag = 1; i = lxt; break;
                }
                ihi = lxt;
              }
            }
          /* *** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL */
            while (1) {
              PDL_Indx middle = ((inbv_datap)[0] + ihi) / 2;
              if (middle == (inbv_datap)[0]) {
                mflag = 0; i = (inbv_datap)[0]; break;
              }
          /*   NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 */
              if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(middle))])
                ihi = middle;
              else
                (inbv_datap)[0] = middle;
            }
          } while (0); /* end dintrv inlined */
          if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(k-1))]) {
            return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "X IS N0T GREATER THAN OR EQUAL TO T(K)");
          }
          if (mflag != 0) {
            if ((x_datap)[0] > (t_datap)[0+(__inc_t_nplusk*(i))]) {
              return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "X IS NOT LESS THAN OR EQUAL TO T(N+1)");
            }
            while (1) {
              if (i == k-1) {
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "A LEFT LIMITING VALUE CANNOT BE OBTAINED AT T(K)");
              }
              --i;
              if ((x_datap)[0] != (t_datap)[0+(__inc_t_nplusk*(i))]) {
                break;
              }
            }
          /* *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES */
          /*   WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K */
          }
          PDL_Indx imk = i+1 - k, j;
          {/* Open k3=:k */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(k, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6246 "lib/PDL/Primitive.pd"
              (work_datap)[0+(__inc_work_k3*(k3))] = (a_datap)[0+(__inc_a_n*(imk+k3))];
          }} /* Close k3=:k */
          #line 6248 "lib/PDL/Primitive.pd"
          if ((ideriv_datap)[0] != 0) {
            for (j = 0; j < (ideriv_datap)[0]; ++j) {
              PDL_Indx kmj = k - j - 1;
              PDL_Float fkmj = kmj;
              {/* Open k3=0:kmj */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(kmj, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6253 "lib/PDL/Primitive.pd"
                  PDL_Indx ihi = i+1 + k3;
                  (work_datap)[0+(__inc_work_k3*(k3))] = ((work_datap)[0+(__inc_work_k3*(k3+1))] - (work_datap)[0+(__inc_work_k3*(k3))]) / ((t_datap)[0+(__inc_t_nplusk*(ihi))] - (t_datap)[0+(__inc_t_nplusk*(ihi-kmj))]) * fkmj;
                }} /* Close k3=0:kmj */
          #line 6256 "lib/PDL/Primitive.pd"
            }
          /* *** COMPUTE VALUE AT *X* IN (T(I),T(I+1)) OF IDERIV-TH DERIVATIVE, */
          /*   GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV). */
          }
          PDL_Indx km1 = k - 1;
          if ((ideriv_datap)[0] != km1) {
            PDL_Indx j, j1 = k, kpk = k + k, j2 = kpk, kmider = k - (ideriv_datap)[0];
            for (j = 0; j < kmider; ++j) {
              PDL_Indx ipj = i + j + 1;
              (work_datap)[0+(__inc_work_k3*(j1))] = (t_datap)[0+(__inc_t_nplusk*(ipj))] - (x_datap)[0];
              (work_datap)[0+(__inc_work_k3*(j2))] = (x_datap)[0] - (t_datap)[0+(__inc_t_nplusk*(i-j))];
              ++j1;
              ++j2;
            }
            for (j = (ideriv_datap)[0]; j < km1; ++j) {
              PDL_Indx kmj = k - j - 1, ilo = kpk + kmj - 1;
              {/* Open k3=0:kmj */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(kmj, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6273 "lib/PDL/Primitive.pd"
                  (work_datap)[0+(__inc_work_k3*(k3))] = ((work_datap)[0+(__inc_work_k3*(k3+1))] * (work_datap)[0+(__inc_work_k3*(ilo))] + (work_datap)[0+(__inc_work_k3*(k3))] *
                      (work_datap)[0+(__inc_work_k3*(k+k3))]) / ((work_datap)[0+(__inc_work_k3*(ilo))] + (work_datap)[0+(__inc_work_k3*(k+k3))]);
                  --ilo;
                }} /* Close k3=0:kmj */
          #line 6277 "lib/PDL/Primitive.pd"
            }
          }
          (ans_datap)[0] = (work_datap)[0+(__inc_work_k3*(0))];
          #line 313 "lib/PDL/Primitive-pp-pchip_bvalu.c"
          }PDL_BROADCASTLOOP_END_pchip_bvalu_readdata
    } break;
    case PDL_D: {
      PDL_DECLARE_PARAMS_pchip_bvalu_1(PDL_Double,D,PDL_Indx,N,PDL_Indx,N)
      PDL_BROADCASTLOOP_START_pchip_bvalu_readdata    {
          #line 6145 "lib/PDL/Primitive.pd"
          PDL_Indx k = __privtrans->ind_sizes[2] - __privtrans->ind_sizes[1];
          if ((ideriv_datap)[0] < 0 || (ideriv_datap)[0] >= k)
            return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K");
          PDL_Indx i;
          int mflag;
          /* *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1) */
          /*   (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)). */
          do { /* inlined dintrv */
            PDL_Indx ihi = (inbv_datap)[0] + 1, lxt = __privtrans->ind_sizes[1];
            if (ihi >= lxt) {
              if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(lxt))]) {
                mflag = 1; i = lxt; break;
              }
              if (lxt <= 0) {
                mflag = -1; i = 0; break;
              }
              ihi = (inbv_datap)[0] = lxt;
            }
            char skipflag = 0;
            if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(ihi))]) {
              PDL_Indx inbv = (inbv_datap)[0];
              if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(inbv))]) {
                mflag = 0; i = (inbv_datap)[0]; break;
              }
          /* *** NOW X .LT. XT(IHI) . FIND LOWER BOUND */
              PDL_Indx istep = 1;
              while (1) {
                ihi = (inbv_datap)[0];
                (inbv_datap)[0] = ihi - istep;
                if ((inbv_datap)[0] <= 0) {
                  break;
                }
                PDL_Indx inbv = (inbv_datap)[0];
                if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(inbv))]) {
                  skipflag = 1;
                  break;
                }
                istep <<= 1;
              }
              if (!skipflag) {
                (inbv_datap)[0] = 0;
                if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(0))]) {
                  mflag = -1; i = 0; break;
                }
              }
              skipflag = 1;
          /* *** NOW X .GE. XT(ILO) . FIND UPPER BOUND */
            }
            if (!skipflag) {
              PDL_Indx istep = 1;
              while (1) {
                (inbv_datap)[0] = ihi;
                ihi = (inbv_datap)[0] + istep;
                if (ihi >= lxt) break;
                if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(ihi))]) {
                  skipflag = 1;
                  break;
                }
                istep <<= 1;
              }
              if (!skipflag) {
                if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(lxt))]) {
                  mflag = 1; i = lxt; break;
                }
                ihi = lxt;
              }
            }
          /* *** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL */
            while (1) {
              PDL_Indx middle = ((inbv_datap)[0] + ihi) / 2;
              if (middle == (inbv_datap)[0]) {
                mflag = 0; i = (inbv_datap)[0]; break;
              }
          /*   NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 */
              if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(middle))])
                ihi = middle;
              else
                (inbv_datap)[0] = middle;
            }
          } while (0); /* end dintrv inlined */
          if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(k-1))]) {
            return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "X IS N0T GREATER THAN OR EQUAL TO T(K)");
          }
          if (mflag != 0) {
            if ((x_datap)[0] > (t_datap)[0+(__inc_t_nplusk*(i))]) {
              return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "X IS NOT LESS THAN OR EQUAL TO T(N+1)");
            }
            while (1) {
              if (i == k-1) {
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "A LEFT LIMITING VALUE CANNOT BE OBTAINED AT T(K)");
              }
              --i;
              if ((x_datap)[0] != (t_datap)[0+(__inc_t_nplusk*(i))]) {
                break;
              }
            }
          /* *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES */
          /*   WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K */
          }
          PDL_Indx imk = i+1 - k, j;
          {/* Open k3=:k */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(k, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6246 "lib/PDL/Primitive.pd"
              (work_datap)[0+(__inc_work_k3*(k3))] = (a_datap)[0+(__inc_a_n*(imk+k3))];
          }} /* Close k3=:k */
          #line 6248 "lib/PDL/Primitive.pd"
          if ((ideriv_datap)[0] != 0) {
            for (j = 0; j < (ideriv_datap)[0]; ++j) {
              PDL_Indx kmj = k - j - 1;
              PDL_Double fkmj = kmj;
              {/* Open k3=0:kmj */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(kmj, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6253 "lib/PDL/Primitive.pd"
                  PDL_Indx ihi = i+1 + k3;
                  (work_datap)[0+(__inc_work_k3*(k3))] = ((work_datap)[0+(__inc_work_k3*(k3+1))] - (work_datap)[0+(__inc_work_k3*(k3))]) / ((t_datap)[0+(__inc_t_nplusk*(ihi))] - (t_datap)[0+(__inc_t_nplusk*(ihi-kmj))]) * fkmj;
                }} /* Close k3=0:kmj */
          #line 6256 "lib/PDL/Primitive.pd"
            }
          /* *** COMPUTE VALUE AT *X* IN (T(I),T(I+1)) OF IDERIV-TH DERIVATIVE, */
          /*   GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV). */
          }
          PDL_Indx km1 = k - 1;
          if ((ideriv_datap)[0] != km1) {
            PDL_Indx j, j1 = k, kpk = k + k, j2 = kpk, kmider = k - (ideriv_datap)[0];
            for (j = 0; j < kmider; ++j) {
              PDL_Indx ipj = i + j + 1;
              (work_datap)[0+(__inc_work_k3*(j1))] = (t_datap)[0+(__inc_t_nplusk*(ipj))] - (x_datap)[0];
              (work_datap)[0+(__inc_work_k3*(j2))] = (x_datap)[0] - (t_datap)[0+(__inc_t_nplusk*(i-j))];
              ++j1;
              ++j2;
            }
            for (j = (ideriv_datap)[0]; j < km1; ++j) {
              PDL_Indx kmj = k - j - 1, ilo = kpk + kmj - 1;
              {/* Open k3=0:kmj */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(kmj, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6273 "lib/PDL/Primitive.pd"
                  (work_datap)[0+(__inc_work_k3*(k3))] = ((work_datap)[0+(__inc_work_k3*(k3+1))] * (work_datap)[0+(__inc_work_k3*(ilo))] + (work_datap)[0+(__inc_work_k3*(k3))] *
                      (work_datap)[0+(__inc_work_k3*(k+k3))]) / ((work_datap)[0+(__inc_work_k3*(ilo))] + (work_datap)[0+(__inc_work_k3*(k+k3))]);
                  --ilo;
                }} /* Close k3=0:kmj */
          #line 6277 "lib/PDL/Primitive.pd"
            }
          }
          (ans_datap)[0] = (work_datap)[0+(__inc_work_k3*(0))];
          #line 461 "lib/PDL/Primitive-pp-pchip_bvalu.c"
          }PDL_BROADCASTLOOP_END_pchip_bvalu_readdata
    } break;
    case PDL_LD: {
      PDL_DECLARE_PARAMS_pchip_bvalu_1(PDL_LDouble,E,PDL_Indx,N,PDL_Indx,N)
      PDL_BROADCASTLOOP_START_pchip_bvalu_readdata    {
          #line 6145 "lib/PDL/Primitive.pd"
          PDL_Indx k = __privtrans->ind_sizes[2] - __privtrans->ind_sizes[1];
          if ((ideriv_datap)[0] < 0 || (ideriv_datap)[0] >= k)
            return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K");
          PDL_Indx i;
          int mflag;
          /* *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1) */
          /*   (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)). */
          do { /* inlined dintrv */
            PDL_Indx ihi = (inbv_datap)[0] + 1, lxt = __privtrans->ind_sizes[1];
            if (ihi >= lxt) {
              if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(lxt))]) {
                mflag = 1; i = lxt; break;
              }
              if (lxt <= 0) {
                mflag = -1; i = 0; break;
              }
              ihi = (inbv_datap)[0] = lxt;
            }
            char skipflag = 0;
            if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(ihi))]) {
              PDL_Indx inbv = (inbv_datap)[0];
              if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(inbv))]) {
                mflag = 0; i = (inbv_datap)[0]; break;
              }
          /* *** NOW X .LT. XT(IHI) . FIND LOWER BOUND */
              PDL_Indx istep = 1;
              while (1) {
                ihi = (inbv_datap)[0];
                (inbv_datap)[0] = ihi - istep;
                if ((inbv_datap)[0] <= 0) {
                  break;
                }
                PDL_Indx inbv = (inbv_datap)[0];
                if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(inbv))]) {
                  skipflag = 1;
                  break;
                }
                istep <<= 1;
              }
              if (!skipflag) {
                (inbv_datap)[0] = 0;
                if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(0))]) {
                  mflag = -1; i = 0; break;
                }
              }
              skipflag = 1;
          /* *** NOW X .GE. XT(ILO) . FIND UPPER BOUND */
            }
            if (!skipflag) {
              PDL_Indx istep = 1;
              while (1) {
                (inbv_datap)[0] = ihi;
                ihi = (inbv_datap)[0] + istep;
                if (ihi >= lxt) break;
                if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(ihi))]) {
                  skipflag = 1;
                  break;
                }
                istep <<= 1;
              }
              if (!skipflag) {
                if ((x_datap)[0] >= (t_datap)[0+(__inc_t_nplusk*(lxt))]) {
                  mflag = 1; i = lxt; break;
                }
                ihi = lxt;
              }
            }
          /* *** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL */
            while (1) {
              PDL_Indx middle = ((inbv_datap)[0] + ihi) / 2;
              if (middle == (inbv_datap)[0]) {
                mflag = 0; i = (inbv_datap)[0]; break;
              }
          /*   NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 */
              if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(middle))])
                ihi = middle;
              else
                (inbv_datap)[0] = middle;
            }
          } while (0); /* end dintrv inlined */
          if ((x_datap)[0] < (t_datap)[0+(__inc_t_nplusk*(k-1))]) {
            return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "X IS N0T GREATER THAN OR EQUAL TO T(K)");
          }
          if (mflag != 0) {
            if ((x_datap)[0] > (t_datap)[0+(__inc_t_nplusk*(i))]) {
              return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "X IS NOT LESS THAN OR EQUAL TO T(N+1)");
            }
            while (1) {
              if (i == k-1) {
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_bvalu:" "A LEFT LIMITING VALUE CANNOT BE OBTAINED AT T(K)");
              }
              --i;
              if ((x_datap)[0] != (t_datap)[0+(__inc_t_nplusk*(i))]) {
                break;
              }
            }
          /* *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES */
          /*   WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K */
          }
          PDL_Indx imk = i+1 - k, j;
          {/* Open k3=:k */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(k, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6246 "lib/PDL/Primitive.pd"
              (work_datap)[0+(__inc_work_k3*(k3))] = (a_datap)[0+(__inc_a_n*(imk+k3))];
          }} /* Close k3=:k */
          #line 6248 "lib/PDL/Primitive.pd"
          if ((ideriv_datap)[0] != 0) {
            for (j = 0; j < (ideriv_datap)[0]; ++j) {
              PDL_Indx kmj = k - j - 1;
              PDL_LDouble fkmj = kmj;
              {/* Open k3=0:kmj */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(kmj, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6253 "lib/PDL/Primitive.pd"
                  PDL_Indx ihi = i+1 + k3;
                  (work_datap)[0+(__inc_work_k3*(k3))] = ((work_datap)[0+(__inc_work_k3*(k3+1))] - (work_datap)[0+(__inc_work_k3*(k3))]) / ((t_datap)[0+(__inc_t_nplusk*(ihi))] - (t_datap)[0+(__inc_t_nplusk*(ihi-kmj))]) * fkmj;
                }} /* Close k3=0:kmj */
          #line 6256 "lib/PDL/Primitive.pd"
            }
          /* *** COMPUTE VALUE AT *X* IN (T(I),T(I+1)) OF IDERIV-TH DERIVATIVE, */
          /*   GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV). */
          }
          PDL_Indx km1 = k - 1;
          if ((ideriv_datap)[0] != km1) {
            PDL_Indx j, j1 = k, kpk = k + k, j2 = kpk, kmider = k - (ideriv_datap)[0];
            for (j = 0; j < kmider; ++j) {
              PDL_Indx ipj = i + j + 1;
              (work_datap)[0+(__inc_work_k3*(j1))] = (t_datap)[0+(__inc_t_nplusk*(ipj))] - (x_datap)[0];
              (work_datap)[0+(__inc_work_k3*(j2))] = (x_datap)[0] - (t_datap)[0+(__inc_t_nplusk*(i-j))];
              ++j1;
              ++j2;
            }
            for (j = (ideriv_datap)[0]; j < km1; ++j) {
              PDL_Indx kmj = k - j - 1, ilo = kpk + kmj - 1;
              {/* Open k3=0:kmj */ PDL_EXPAND2(register PDL_Indx k3=0, __k3_stop=PDLMIN(kmj, (__k3_size))); for(; k3<__k3_stop; k3+=1) {
            #line 6273 "lib/PDL/Primitive.pd"
                  (work_datap)[0+(__inc_work_k3*(k3))] = ((work_datap)[0+(__inc_work_k3*(k3+1))] * (work_datap)[0+(__inc_work_k3*(ilo))] + (work_datap)[0+(__inc_work_k3*(k3))] *
                      (work_datap)[0+(__inc_work_k3*(k+k3))]) / ((work_datap)[0+(__inc_work_k3*(ilo))] + (work_datap)[0+(__inc_work_k3*(k+k3))]);
                  --ilo;
                }} /* Close k3=0:kmj */
          #line 6277 "lib/PDL/Primitive.pd"
            }
          }
          (ans_datap)[0] = (work_datap)[0+(__inc_work_k3*(0))];
          #line 609 "lib/PDL/Primitive-pp-pchip_bvalu.c"
          }PDL_BROADCASTLOOP_END_pchip_bvalu_readdata
    } break;
    default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_bvalu: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
  }
  #undef PDL_IF_BAD
  return PDL_err;
}

static pdl_datatypes pdl_pchip_bvalu_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
static PDL_Indx pdl_pchip_bvalu_vtable_realdims[] = { 1, 1, 0, 0, 0, 0, 1 };
static char *pdl_pchip_bvalu_vtable_parnames[] = { "t","a","ideriv","x","ans","inbv","work" };
static short pdl_pchip_bvalu_vtable_parflags[] = {
  0,
  0,
  PDL_PARAM_ISTYPED,
  0,
  PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
  PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE,
  PDL_PARAM_ISCREAT|PDL_PARAM_ISTEMP|PDL_PARAM_ISWRITE
};
static pdl_datatypes pdl_pchip_bvalu_vtable_partypes[] = { -1, -1, PDL_IND, -1, -1, PDL_IND, -1 };
static PDL_Indx pdl_pchip_bvalu_vtable_realdims_starts[] = { 0, 1, 2, 2, 2, 2, 2 };
static PDL_Indx pdl_pchip_bvalu_vtable_realdims_ind_ids[] = { 2, 1, 0 };
static char *pdl_pchip_bvalu_vtable_indnames[] = { "k3","n","nplusk" };
pdl_transvtable pdl_pchip_bvalu_vtable = {
  PDL_TRANS_DO_BROADCAST, 0, pdl_pchip_bvalu_vtable_gentypes, 4, 7, NULL /*CORE21*/,
  pdl_pchip_bvalu_vtable_realdims, pdl_pchip_bvalu_vtable_parnames,
  pdl_pchip_bvalu_vtable_parflags, pdl_pchip_bvalu_vtable_partypes,
  pdl_pchip_bvalu_vtable_realdims_starts, pdl_pchip_bvalu_vtable_realdims_ind_ids, 3,
  3, pdl_pchip_bvalu_vtable_indnames,
  pdl_pchip_bvalu_redodims, pdl_pchip_bvalu_readdata, NULL,
  NULL,
  0,"PDL::Primitive::pchip_bvalu"
};


pdl_error pdl_run_pchip_bvalu(pdl  *t,pdl  *a,pdl  *ideriv,pdl  *x,pdl  *ans,pdl  *inbv) {
  pdl_error PDL_err = {0, NULL, 0};
  if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
  pdl_trans *__privtrans = PDL->create_trans(&pdl_pchip_bvalu_vtable);
  if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
  __privtrans->pdls[0] = t;
  __privtrans->pdls[1] = a;
  __privtrans->pdls[2] = ideriv;
  __privtrans->pdls[3] = x;
  __privtrans->pdls[4] = ans;
  __privtrans->pdls[5] = inbv;
  PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
  PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
  return PDL_err;
}