Bug report: TGraph::Eval()

From: Alexandr Malusek <Alexandr.Malusek_at_imv.liu.se>
Date: Tue, 26 Apr 2005 15:58:33 +0200


Hi,

There seems to be a bug in the TGraph::Eval() routine:

Test code:

{

  const Int_t dim = 5;
  Double_t x[dim] = {1.0, 2.0, 3.0, 4.0, 5.0};
  Double_t y[dim] = {1.0, 2.0, 3.0, 4.0, 5.0};
  TGraph *g = new TGraph(dim, x, y);
  cout << "y(4.99) = " << g->Eval(4.99) << endl;
  cout << "y(5.00) = " << g->Eval(5.00) << endl;
  cout << "y(5.01) = " << g->Eval(5.01) << endl;
}

Output:

y(4.99) = 4.99          # correct
y(5.00) = 1.01246e+07   # wrong
y(5.01) = 1.01449e+07   # wrong

I think the code in root/graf/src/TGraph.cxx:

Double_t TGraph::Eval(Double_t x, TSpline *spline, Option_t *option) const

  ...

//linear interpolation
//find point in graph immediatly below x
//In case x is < fX[0] or > fX[fNpoints-1] return the extrapolated point

      Int_t up, low = TMath::BinarySearch(fNpoints,fX,x);
      if (low == fNpoints-1) {up=low; low = up+1;}
      if (low == -1) {low=0; up=1;}
      if (fX[low] == fX[low+1]) return fY[low];
      Double_t yn = x*(fY[low]-fY[low+1]) +fX[low]*fY[low+1] - fX[low+1]*fY[low];
      return yn/(fX[low]-fX[low+1]);

should be

//linear interpolation
//find point in graph immediately below x
//In case x is < fX[0] or > fX[fNpoints-1] return the extrapolated point

      Int_t low = TMath::BinarySearch(fNpoints,fX,x);
      Int_t up = low+1;
      if (low == fNpoints-1) {up=low; low = up-1;}
      if (low == -1) {low=0; up=1;}
      if (fX[low] == fX[up]) return fY[low];
      Double_t yn = x*(fY[low]-fY[up]) +fX[low]*fY[up] - fX[up]*fY[low];
      return yn/(fX[low]-fX[up]);


The problem is with

if (low == fNpoints-1) {up=low; low = up+1;}

the rest are cosmetic changes to utilize the variable "up". The modified code gives correct results:

y(4.99) = 4.990000e+00
y(5.00) = 5.000000e+00
y(5.01) = 5.010000e+00

The problem was observed with root-4.03/02 and root-4.03/04 on Linux, other versions may also be affected.

Regards,
Alexandr Received on Tue Apr 26 2005 - 15:58:41 MEST

This archive was generated by hypermail 2.2.0 : Tue Jan 02 2007 - 14:45:07 MET