make dpchid into PP kernel (still in pchip.c as chia calls

lib/PDL/Primitive.pd

@@ -4891,19 +4891,62 @@ which should never happen.
 'Integrate (x,f(x)) over arbitrary limits.'
 );
 
-defslatec('pchip_chid',
-	  {D => 'dpchid'},
-	  'IntVal	n=$SIZE(n);
-           Mat          x       (n);
-           Mat          f       (n);
-           Mat          d       (n);
-           Logical [io] skip    ();
-           IntScalar    ia      ();
-           IntScalar    ib      ();
-           FuncRet [o]  ans     ();
-           Integer [o]  ierr    ();
-          ',
-'Evaluate the definite integral of a a piecewise
+pp_def('pchip_chid',
+  Pars => 'x(n);f(n);d(n);int [io]skip();longlong ia();longlong ib();[o]ans();longlong [o]ierr()',
+  RedoDimsCode => <<'EOF',
+if ($SIZE(n) < 2)
+  $CROAK("NUMBER OF DATA POINTS LESS THAN TWO");
+EOF
+  Code => <<'EOF',
+/*  Programming notes: */
+/*  1. This routine uses a special formula that is valid only for */
+/*   integrals whose limits coincide with data values.  This is */
+/*   mathematically equivalent to, but much more efficient than, */
+/*   calls to DCHFIE. */
+/* Local variables */
+PDL_Indx i, iup, low;
+/*  VALIDITY-CHECK ARGUMENTS. */
+if (!$skip()) {
+  loop (n=1) %{
+    if ($x() > $x(n=>n-1)) continue;
+    $ierr() = -1;
+    $CROAK("X-ARRAY NOT STRICTLY INCREASING");
+  %}
+}
+/*  FUNCTION DEFINITION IS OK, GO ON. */
+$skip() = 1;
+if ($ia() < 0 || $ia() > $SIZE(n)-1 || $ib() < 0 || $ib() > $SIZE(n)-1) {
+  $ierr() = -4;
+  $CROAK("IA OR IB OUT OF RANGE");
+}
+$ierr() = 0;
+/*  COMPUTE INTEGRAL VALUE. */
+$GENERIC() sum, value = 0.;
+if ($ia() != $ib()) {
+  low = PDLMIN($ia(),$ib());
+  iup = PDLMAX($ia(),$ib());
+  sum = 0.;
+  loop (n=low:iup) %{
+    $GENERIC() h = $x(n=>n+1) - $x();
+    sum += h * ($f() + $f(n=>n+1) + ($d() - $d(n=>n+1)) * (h / 6.));
+  %}
+  value = 0.5 * sum;
+  if ($ia() > $ib())
+    value = -value;
+}
+$ans() = value;
+EOF
+  GenericTypes => $F,
+  Doc => <<'EOF',
+=for ref
+
+Evaluate the definite integral of a piecewise cubic
+Hermite function over an interval whose endpoints are data
+points.
+
+=for usage
+
+Evaluate the definite integral of a a piecewise
 cubic Hermite function between C<x($ia)> and
 C<x($ib)>.
 
@@ -4957,10 +5000,7 @@ Error status returned by C<$ierr>:
 (VALUE will be zero in any of these cases.)
 NOTE: The above errors are checked in the order listed, and following
 arguments have B<NOT> been validated.
-',
-'Evaluate the definite integral of a piecewise cubic
-Hermite function over an interval whose endpoints are data
-points.'
+EOF
 );
 
 defslatec('pchip_chbs',