Dann is a software engineer at Microsoft. He can be contacted at a-cnadc@microsoft.com.
As a former calculus teacher, I know how hard it is to explain basic numerical-integration techniques. That's why I was especially pleased when Dann Corbit sent me this clear exposition of both the elementary methods and some new approaches that have just been developed. If you've ever wondered how numerical integration works, or even if you think you already understand it, you should also appreciate Dann's contribution.
-Tim Kientzle
Integration is an important step in many engineering and scientific calculations. Rather than interface to numerical Fortran libraries, I wrote a set of C++ classes to support numerical integration. The base class, Integral, implements the actual integration algorithm. It calls what is essentially a virtual method to evaluate the function to be integrated. You use it by subclassing and overriding this one virtual method.
To test the integration routines, I implemented classes of elliptic integrals; see Figure 1. Class EllipticIntegral inherits from Integral. Then the specific types of elliptic integrals inherit from class EllipticIntegral.
Elliptic integrals are a good test because there is more than one way to evaluate them, making it possible to check your results. Because they are used in a wide variety of problems (such as computing orbits of satellites), there are also books and other programs that provide values for comparison. I've checked my results against the Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, edited by Milton Abramowitz and Irene A. Stegun (U.S. Government Printing Office, 1964. LCCN 64060036), and against Mathematica.
Although my elliptic integrals of the first and second kind agree with Mathematica to 15 digits or so, my elliptic integrals of the third kind only agreed within 6 digits. Since the Handbook only lists six digits, I can't tell who is the more correct. Perhaps an intrepid reader who has access to high-precision results for elliptic integrals of the third kind will let me know if I have a flaw in my program.
The first method of integration I used was based on a hybrid of the Newton-Cotes formula and the "divide-and-conquer" strategy. The simplest example of this method is the trapezoidal rule. You divide the area to be integrated into a number of trapezoids; see Figure 2. You then add up the areas of the trapezoids. If your answer is not good enough, you can double the number of trapezoids and try again.
As you can see from Table 1, it takes a lot of trapezoids to get a reasonably accurate answer to the integration. Even with a million of them, you don't get double-precision accuracy. The problem is that straight line segments don't follow a curve very well. The shaded area in Figure 2(b) shows the error caused by using trapezoids to calculate the area.
Another way to think about trapezoidal-rule integration is that you're approximating your function with short line segments, then integrating the "piece-wise linear" function. By using more sophisticated approximations, you can get better results. For example, if you could fit polynomials through the points, then you would get a much better estimate, since polynomials follow curves much more closely than line segments. For example, the curve y=1/x looks much more like a bit of a parabola than a straight line segment.
While the trapezoidal rule approximates the function with little line segments, Simpson's rule (see Figure 6) approximates the function with little parabolas. This method involves calculating the function over uniform intervals and is an example of the closed Newton-Cotes formulas. It is accurate for smooth functions and for polynomials of low degree. Surprisingly, it provides exact answers for polynomials up to cubics. Simpson's rule assumes uniformly wide intervals and that the number of intervals, n, is even.
As Table 2 shows, it takes only 2048 little intervals to get a pretty good estimate using Simpson's rule, compared to over a million for the trapezoid rule.
You can use higher-order polynomials than parabolas to estimate the integral. Divide the interval [a,b] into n equal parts, let yn = 
(xn), and let h be the width of each part, (b-a)/n. Then find polynomials that fit through those points and integrate them to approximate the area under the curve. You can think of higher-order Newton-Cotes rules as integrals of quadratic equations, cubic equations, and so on, that approximate your original function. William Press' Numerical Recipes in C (Cambridge University Press, 1992) shows how to create your own Newton-Cotes formulas.
A metaphor for understanding Newton-Cotes (and other numerical-integration methods) is to imagine sets of boxes on a chain. Each box has exactly two properties: its position on the chain, and its multiplier or weight. For the Newton-Cotes rules, you assume that all of the boxes are evenly spaced along the chain, including one at each end. Now, suppose you're integrating a polynomial of degree n+1. The rules for integrating polynomials are very well known, and you can analytically determine the exact answer. This condition gives you a set of n+1 equations in n+1 unknowns. Solving this matrix gives you a Newton-Cotes rule. The next section shows how to derive Simpson's rule from the general-integration approach described previously. You can also derive methods that use higher-order polynomials. I used an eleven-point Newton-Cotes rule to make the integral class in Listing One. I recursively divided the interval into subintervals until each partition was small enough. Then I applied the Newton-Cotes rule. (Source code and executables that implement Newton-Cotes and other techniques discussed in this article are available electronically.)
Rather than reevaluating the integral at each level of recursion, I subdivide until the step size is small enough, then evaluate the integral once. I arrived at "small enough" heuristically, rather than using some sophisticated analysis. My real goal was to have a function that will accurately integrate "well-behaved" functions and be simple to implement.
The error for Newton-Cotes is related to step size and the nth derivative. For the eleven-point rule I used, the maximum possible error is proportional to h13 (h is the step size) times the 12th derivative of the function. A step size of .01 raised to the 13th power is 10-26. Even if the 12th derivative of (x) reaches 10 million, you still get 19 digits of accuracy, sufficient for long double results. Of course, that's the theoretical accuracy, which doesn't account for accumulated rounding errors due to limited precision. So I tried functions that are readily available as library functions-such as log and bessel functions-and compared the results with direct integration. I was satisfied that, for most purposes, I will get a very accurate answer by this method.
A much larger step size won't do. For instance, 0.05 to the 13th power is 10-17 (just double accuracy if the 12th derivative is 1). So, practically speaking, I will always need to recurse to about the same level. Certain ultra simple functions would not require this, but I won't integrate them numerically anyway.
What sort of integrations are good candidates for a Newton-Cotes method?
Newton-Cotes rules are fairly easy to derive from scratch. The basic idea is that you start with the values of (x) at n evenly spaced points x0, x1, ..., xn-1. The easiest thing to compute with these values is the weighted sum in Figure 3.
The problem is to find the weights wj. Newton-Cotes attempts to make the rule exact for low-degree polynomials. For example, it should be exact if you integrate 1, which gives you the equation in Figure 4(a). Similarly, if you want it exact for linear functions, you need to add the condition in Figure 4(b). Remember that the xj are predetermined, so you can only choose the n weights.
For concreteness, here's how to derive Simpson's rule: Simpson's rule is exact for quadratics, so you need three equations for 1, x, and x2, which allows you to use three weights. A few simplifications are useful. First, since you're integrating from a to b, you can rename x0 to a, x1 to (a+b)/2, and x2 to b. Also, these formulae are supposed to work for any values of a and b, so you can assume that a=0 and b=1. You then end up with Figures 5(a) and 5(b), corresponding to Figures 4(a) and 4(b). Figure 5(c) is the third equation you need.
Solving the three equations in Figure 5, you get w0 = 1/6, w1 = 2/3, w2 = 1/6, yielding the familiar forms of Simpson's rule given in Figure 6. By following the same outline, you can derive Newton-Cotes rules of any order.
Now, remember the simplifying assumption that the boxes be evenly spaced for the Newton-Cotes methods? It turns out that this is not the best possible place to put the boxes. If you instead allow the position on the chain to vary along with the multiplier, you get a grand total of 2n+2 equations in 2n+2 unknowns. This allows the solution to be exact for polynomials of degree <= 2n+1, with the same number of function evaluations (boxes). Alternatively, you can get the same accuracy with only half as many function evaluations. This is especially useful if the function is difficult to evaluate (for example, you might have to run a lab experiment for each value). Integration rules that use this derivation are called Gaussian rules. Listing Two is an implementation of a 22-point ''Gaussian rule," which is exact for 44th degree polynomials.
Each of the aforementioned techniques has weaknesses. Newton-Cotes rules require twice as many function evaluations as Gaussian rules for the same degree of accuracy, which makes them a poor choice if the function is difficult to evaluate. On the other hand, with Newton-Cotes rules, you can subdivide your intervals for better accuracy and reuse previous function values. With Gaussian rules, you can't reuse points because they are unevenly spaced. Higher-order Gaussian rules require you to throw away most of your previous effort.
There is a family of weighted chain rules that nicely solves these problems. These rules are called "Recursive Monotone Stable" (RMS) rules. They are defined recursively, starting with the Trapezoidal rule. Monotone means that each successive rule has smaller spaces between boxes. You also want none of the weights to be negative. A series of articles in ACM Transactions on Mathematical Software (June, 1991) describes such a sequence of rules.
RMS rules have remarkable properties. For instance, to evaluate at higher precision, all of the previous calculations can be reused. You only have to evaluate the function at some intermediate points to get higher precision. You can start with only a few function calls, adding more evaluations as needed. If the values fail to converge, you can call higher-order rules as needed, reusing our previous calculations. RMS rules compete successfully with other well-established rules for both precision and speed. Going back to our boxes-on-chains metaphor, with an RMS rule, you have lots of boxes on your chain, but most of them are turned off. You make one pass with only a few boxes turned on, then turn on some more and look at the answer. You can keep switching on more boxes until the answer is good enough. The RMS rule is shown in Listing Three.
Having an integral class is a handy abstraction for solving complex problems. Suppose, for instance, that you wanted to find the effort to push water aside from a sailboat hull (that is, the time rate of change of the momentum of the water). You could create a matrix of integrals to solve it. I plan to use this object to implement many of the formulae from the Abramowitz and Stegun Handbook of Mathematical Functions by direct integration, rather than series summations.
Curiously, a Fortran program from NETLIB was used to collect the constants and to study the algorithm used. Thus, my efforts at integration have come full circle. I was originally going to avoid Fortran by writing my own integration routine. If the best available algorithm is implemented in Fortran, why not learn from the numerical-analysis masters?
Baker, Louis, C Mathematical Function Handbook. McGraw-Hill, 1992. ISBN 0-07-911158-0.
Weisstein, Eric, Eric's Treasure Trove of Mathematics, http://www.gps.caltech.edu/ ~eww/math/math0.html
Since I wrote the article about my Integration class, I've made a few changes that may be of interest. I decided to evaluate improper integrals (those for which either or both of the limits of integration are infinite). There is a simple change of variables that can be used to integrate over an infinite range. From Eric's Treasure Trove of Mathematics web page, "Integrals (Infinite Limits)," we find the following transformation:
This formula can also be used to reduce a huge interval of integration to a small one. I use it for that purpose with my RMS rule. However, since the RMS rule includes the endpoints of the interval, I cannot use this transformation over infinite intervals with the RMS rule. I created a 200-point Gaussian rule that can be used for this purpose. The integral class now allows integration over enormous intervals, tiny intervals, and infinite intervals with equal ease. There are some caveats for integration over infinite intervals. If the function integrated decreases towards infinity slower than 1/x2, the answer provided could be far from the true value. If the function is infinite in only one direction, it must not pass through the origin (note the restriction ab > 0). There are still many improvements that could be made to these routines, for example, automatic asymptote detection.
--D.C.
Trapezoids Area Estimate Error 1 0.75 0.056852819 2 0.708333333 0.015186153 4 0.69702381 0.003876629 8 0.69412185037185 0.00097467 16 0.693391202207526 0.000244022 . . . . . . . . . 65536 0.693147180574497 1.5E-11 131072 0.693147180563583 3.6E-12 262144 0.693147180560854 9.1E-13 524288 0.693147180560172 2.3E-13 1048576 0.693147180560002 5.7E-14 Correct Answer 0.693147180559945
Intervals Integral by Simpson's rule
2 .6944444444444
4 .693253968254
8 .6931545306545
16 .6931476528194
32 .6931472102898
64 .6931471824215
128 6931471806763
256 6931471805672
512 .6931471805604
1024 .6931471805600
2048 .6931471805599
Figure 1: Elliptic integrals. (a) first kind; (b) second kind; (c) third kind.
Figure 3: Approximating an integral with a weighted sum of function values.
Figure 4: (a) Making the rule exact for 1; (b) making the rule exact for x.
Figure 5: (a) Making Simpson's rule exact for 1; (b) making Simpson's rule exact for x; (c) making Simpson's rule exact for x2.
Figure 6: Different forms of Simpson's rule.
Listing One
void EXPORT Integral::Evaluate() {
// We use this constant to divide the interval into two equal parts:
const long double ldHalf = 0.5e0L;
// Note that the maximum possible error is:
// -(1,346,350/326,918,592)*h^13*epsilon, where epsilon is some value
// of the 12th derivative of the supplied function in our sub-interval.
// See "Handbook of Mathematical Functions" by Abramowitz and Stegun,
// formula 25.4.20 on page 887.
// Empirically derived for long double stepsize:
const long double ldStepsizeForBestAccuracy = 0.01e0L;
// The stepsize used for the integral estimation
// The number .1 is due to the 11 terms in the Newton-Cotes formula
// creating 10 sub-intervals.
const long double ldH = ( ldX1 - ldX0 )*0.1e0L;
// If stepsize is too large, divide the integral into two equal sections.
if ( ldH > ldStepsizeForBestAccuracy ) {
long double ldMidpoint = ldX0 + ( ldX1 - ldX0 )*ldHalf;
Integral i0( ldX0, ldMidpoint, ldFprime, pldConstArray );
Integral i1( ldMidpoint, ldX1, ldFprime, pldConstArray );
ldIntegralResult=i0.ldGetIntegralResult()+i1.ldGetIntegralResult();
}
else { // Newton-Cotes constant for calculating the integral:
const long double A = 5.0e0L/299376.0e0L;
// Perform the summation.
ldIntegralResult = (
16067e0L * ( ldFprime( ldX0, pldConstArray )
+ldFprime( ldX0 + ldH*10, pldConstArray ) )
+ 106300e0L * ( ldFprime( ldX0 + ldH, pldConstArray )
+ldFprime( ldX0 + ldH*9, pldConstArray ) )
- 48525e0L * ( ldFprime( ldX0 + ldH*2, pldConstArray )
+ldFprime( ldX0 + ldH*8, pldConstArray ) )
+ 272400e0L * ( ldFprime( ldX0 + ldH*3, pldConstArray )
+ldFprime( ldX0 + ldH*7, pldConstArray ) )
- 260550e0L * ( ldFprime( ldX0 + ldH*4, pldConstArray )
+ldFprime( ldX0 + ldH*6, pldConstArray ) )
+ 427368e0L * ldFprime( ldX0 + ldH*5, pldConstArray )
) * ldH * A;
}
}
Listing Two
void EXPORT Integral::Evaluate() {
// We use this constant to divide the interval into two equal parts:
const long double ldHalf = 0.5e0L;
// Empirically derived for long double stepsize: The number .21 is due to
// the 22 weights in our Gaussian formula creating 21 subdivisions. This
// will give a stepsize of .01, when .21 is divided 21 times.
// Theoretically, the error should be less than that of Newton-Cotes.
const long double ldSmallEnoughInterval = 0.21e0L;
// If interval is too large, divide the integral into two equal sections.
if ( ( ldX1 - ldX0 ) > ldSmallEnoughInterval ) {
long double ldMidpoint = ldX0 + ( ldX1 - ldX0 )*ldHalf;
Integral i0( ldX0, ldMidpoint, ldFprime, pldConstArray );
Integral I1( ldMidpoint, ldX1, ldFprime, pldConstArray );
ldIntegralResult = i0.ldGetIntegralResult() + I1.ldGetIntegralResult();
}
else
ldIntegralResult = ldGaussian( ldFprime, ldX0, ldX1, pldConstArray );
}
// Internal Method, Gaussian rule for integral evaluation:
static long double ldGaussian(_ldIntegralFunction_t ldF,
long double ldA, long double ldB, long double *pldConstArr) {
const long double ldGLWeights[] = {
2.78342835580868333037e-002L, 6.27901847324523123016e-002L,
9.31451054638671257124e-002L, 1.16596882295995240001e-001L,
1.31402272255123331128e-001L, 1.36462543388950315360e-001L,
1.31402272255123331128e-001L, 1.16596882295995240001e-001L,
9.31451054638671257124e-002L, 6.27901847324523123016e-002L,
2.78342835580868333037e-002L,
};
const long double ldGLAbscissae[] = {
1.08856709269715036124e-002L, 5.64687001159523504590e-002L,
1.34923997212975337962e-001L, 2.40451935396594092037e-001L,
3.65228422023827513853e-001L, 5.00000000000000000000e-001L,
6.34771577976172486147e-001L, 7.59548064603405907963e-001L,
8.65076002787024662065e-001L, 9.43531299884047649541e-001L,
9.89114329073028496360e-001L,
};
long j; /* Array index */
long double ldXOffset; /* X scale factor for the interval */
long double ldSum; /* Partial summation term */
long double ldCenter = ( ldB + ldA ) * 0.5; /* Average of b and a */
long double ldH = ( ldB - ldA ) * 0.5; /* Half the width of the interval */
long double ldIntegralResult = 0.0e0L;
for ( j = 0; j < 11; j++ ) {
/* Calculate f(x) at two symmetrical offsets from the center: */
ldXOffset = ldH * ldGLAbscissae[j];
ldSum = ( *ldF )( ldCenter - ldXOffset, pldConstArr ) +
( *ldF )( ldCenter + ldXOffset, pldConstArr );
ldIntegralResult += ldGLWeights[j] * ldSum;
}
return ( ldIntegralResult * ldH );
}
Listing Three
#define HALF_KNOT_COUNT 21
#define ABS( a ) ( a ) > 0 ? ( a ) : -( a )
void EXPORT Integral::Evaluate() {
// Empirically derived for long double stepsize: const long double ldSmallEnoughInterval = 0.25e0L;
// If interval is too large, divide the integral into two equal sections.
if ( ( ldX1 - ldX0 ) > ldSmallEnoughInterval ) {
long double ldMidpoint = ldX0 + ( ldX1 - ldX0 ) * 0.5e0L;
Integral i0( ldX0, ldMidpoint, ldFprime, pldConstArray );
Integral I1( ldMidpoint, ldX1, ldFprime, pldConstArray );
ldIntegralResult=i0.ldGetIntegralResult()+I1.ldGetIntegralResult();
}
else {
const long double ldTol = LDBL_EPSILON*2.0e0L;
int I1 = -1;
int I2 = -1;
long double pldVectorOne[HALF_KNOT_COUNT] = { 0 };
long double pldVectorTwo[HALF_KNOT_COUNT] = { 0 };
int iRuleChoice = 0;
long double ldResult0;
long double ldResult1;
ldResult0 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice++, pldConstArray,
pldVectorOne, pldVectorTwo, &I1, &I2 );
ldResult1 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice++, pldConstArray,
pldVectorOne, pldVectorTwo, &I1, &I2 );
if ( ABS( ( ldResult0 - ldResult1 )/ldResult1 ) > ldTol ) {
ldResult0 = ldResult1;
ldResult1 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice++, pldConstArray,
pldVectorOne, pldVectorTwo, &I1, &I2 );
if ( ABS( ( ldResult0 - ldResult1 )/ldResult1 ) > ldTol )
ldResult1 = qxrul( ldFprime, ldX0, ldX1, iRuleChoice,
pldConstArray, pldVectorOne, pldVectorTwo, &I1, &I2 );
}
ldIntegralResult = ldResult1;
}
}
/* These routines are derived from function QXRUL from:
** ALGORITHM 691, COLLECTED ALGORITHMS FROM ACM.
** THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
** VOL. 17, NO. 2, PP. 218-232. JUNE, 1991.
** This was a component of the SLATEC FORTRAN library created by:
** LAWRENCE LIVERMORE NATIONAL LABORATORY
** To compute I = integral of f() over ( ldA, ldB ), with error estimate
** Parameters on entry:
** f() - long double function to be integrated
** ldA - long double lower limit of integration
** ldB - long double upper limit of integration
** iRuleChoice - int
** Selector for choice of local integration rule. Rule used with:
** 13 points if iRuleChoice = 0, 19 points if iRuleChoice = 1,
** 27 points if iRuleChoice = 2, 42 points if iRuleChoice = 3
** pldVectorOne - long double
** vector containing pI1 saved functional values of positive abscissas.
** pldVectorTwo - long double
** vector containing pI2 saved functional values of negative abscissas.
** pI1, pI2 - int number of elements in pldVectorOne, pldVectorTwo respectively.
** On return ** ldY - long double approximation to the integral I.
** pldVectorOne - long double
** vector containing pI1 saved functional values of positive abscissas.
** pldVectorTwo - long double
** vector containing pI2 saved functional values of negative abscissas.
** pI1, pI2 - int number of elements in pldVectorOne,pldVectorTwo respectively.
*/
static long double qxrul(_ldIntegralFunction_t ldF, long double ldA,
long double ldB, int iRuleChoice, long double *pldConstArray,
long double *pldVectorOne, long double *pldVectorTwo,
int *pI1, int *pI2) {
/* Starting positions into the weights array: */
int iStart[] = { 0, 7, 17, 31 };
/* Lengths of the weight sequence to process: */
int iLen[] = { 7, 10, 14, HALF_KNOT_COUNT };
/* Symmetrical offsets +/- for unit interval */
const long double ldAbscissae[HALF_KNOT_COUNT] = {
0.0000000e0L, 0.2500e0L, 0.5000e0L, 0.750000e0L, 0.87500e0L,
0.9375000e0L, 1.0000e0L, 0.3750e0L, 0.625000e0L, 0.96875e0L,
0.1250000e0L, 0.6875e0L, 0.8125e0L, 0.984375e0L, 0.18750e0L,
0.3125000e0L, 0.4375e0L, 0.5625e0L, 0.843750e0L, 0.90625e0L, 0.9921875e0L
};
const long double ldWeights[52] = { // Shortened to 25 digits for printing
0.1303262173284849021810473e0L, 0.2390632866847646220320330e0L,
0.2630626354774670227333506e0L, 0.2186819313830574175167853e0L,
0.0275789764664283686585960e0L, 0.1055750100538458443365035e0L,
0.0157119426059518225416843e0L, 0.1298751627936015783241174e0L,
0.2249996826462523640447835e0L, 0.1680415725925575286319047e0L,
0.1415567675701225879892812e0L, 0.1006482260551160175038684e0L,
0.0251060486072428247905834e0L, 0.0094029643600097471100311e0L,
0.0554269923329587516840678e0L, 0.0998673524740336752572038e0L,
0.0450752305681049246641588e0L, 0.0630094224964777393174617e0L,
0.1261383225537664703013000e0L, 0.1273864433581028272878710e0L,
0.0857650041431182051421409e0L, 0.0710288484231025339744731e0L,
0.0502638357285794240375983e0L, 0.0046836700106090938104326e0L,
0.1235837891364555000245005e0L, 0.1148933497158144016800200e0L,
0.0125257577422612263339148e0L, 0.1239572396231834242194190e0L,
0.0250130641375031057952595e0L, 0.0491595791814613009425885e0L,
0.0225916737495647471330203e0L, 0.0636276297878272455926934e0L,
0.0995006582734679464319326e0L, 0.0704822000271856536609874e0L,
0.0651229733939833564587270e0L, 0.0399822915031365972479053e0L,
0.0345651225708028750983205e0L, 0.0022121679758841144327603e0L,
0.0814032642594593804596783e0L, 0.0658321344760055290627354e0L,
0.0259291372645079254606423e0L, 0.1187141856692283347609436e0L,
0.0599994760538597198558967e0L, 0.0550093798019804173691026e0L,
0.0052644224217646559697603e0L, 0.0153312687405658695933837e0L,
0.0352715936975012310045570e0L, 0.0500055643165395512421280e0L,
0.0574416483117972010634072e0L, 0.0159882379728381343830125e0L,
0.0263566041022088499347248e0L, 0.0119600393794554109167011e0L
};
long double ldY = 0.0e0L;
long double ldOffset; long double ldMidpoint;
long double ldHalfInterval;
int iKnots;
int i;
iKnots = iLen[iRuleChoice];
ldHalfInterval = ( ldB - ldA ) * 0.5;
ldMidpoint = ldA + ldHalfInterval;
for ( i = 0; i < iKnots; i++ ) {
if ( i >= *pI1 || i >= *pI2 ) {
ldOffset = ldHalfInterval * ldAbscissae[i];
if ( i >= *pI1 )
pldVectorOne[i] = ( *ldF )( ldMidpoint + ldOffset, pldConstArray );
if ( i >= *pI2 )
pldVectorTwo[i] = ( *ldF )( ldMidpoint - ldOffset, pldConstArray );
}
ldY +=(pldVectorOne[i]+pldVectorTwo[i])*ldWeights[iStart[iRuleChoice]+i];
}
ldY *= ldHalfInterval;
if ( *pI1 <= iKnots ) *pI1 = iKnots;
if ( *pI2 <= iKnots ) *pI2 = iKnots;
return ldY;
}