1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
|
void Polynom_MakeLagrangeSerie(Float_t* aDst, UInt_t deg, const Float_t* aX, const Float_t* aY)
{
// Lagrange serie construction... thats "quite" easy... we just "add" deg+1 polynoms (of degree deg) together...
// For each polynom 'i', the polynom to add is the products of all (x - aX[j]) with (j != i).
// each polynom is then multiplied by aY[i] and devided by its value when x = aX[i], final result being aY[i] at aX[i], and 0 at aX[j],j!=i
// cleaning destination polynom...
RW_MEMCLEAR(aDst,sizeof(Float_t)*deg);
// allocating temporary polynom...
Float_t* aTmp = RW_ALLOC_ARRAY(Float_t,deg+1);
for (UInt_t i = 0; (i <= deg); ++i) {
Float_t val = aY[i];
if (val) {
Float_t tmpEval = 1.f;
RW_MEMCLEAR(aTmp,sizeof(Float_t)*deg);
aTmp[0] = 1.f; // ALWAYS !!
// ok... computing the polynom to add with 'i'
UInt_t n = 0;
for (UInt_t j = 0; (j <= deg); ++j) {
if (i != j) {
tmpEval *= (aX[i] - aX[j]);
++n;
for (UInt_t k = n; (k > 0); --k)
aTmp[k] = aTmp[k-1] - aX[j] * aTmp[k];
aTmp[0] *= -aX[j];
}
}
// Mult of the polynom !
Polynom_MAdd(aDst,aDst,deg,aTmp,deg,(val / tmpEval));
}
}
// WE ARE DONE !!!
RW_FREE(aTmp);
}
void Polynom_MAdd(Float_t* pOut, const Float_t* pSrc1, UInt_t deg1, const Float_t* pSrc2, UInt_t deg2, Float_t f)
{
UInt_t deg = RW_MIN2(deg1,deg2);
UInt_t i = 0;
for (;(i<=deg);++i) *pOut++ = *pSrc1++ + (*pSrc2++ * f);
for (;(i<=deg1);++i) *pOut++ = *pSrc1++;
for (;(i<=deg2);++i) *pOut++ = (*pSrc2++ * f);
} |
Partager