Source Code for Module hedge.interpolation

 1  """1D Newton interpolation helper module.""" 
 2   
 3  from __future__ import division 
 4   
 5  __copyright__ = "Copyright (C) 2007 Andreas Kloeckner" 
 6   
 7  __license__ = """ 
 8  This program is free software: you can redistribute it and/or modify 
 9  it under the terms of the GNU General Public License as published by 
10  the Free Software Foundation, either version 3 of the License, or 
11  (at your option) any later version. 
12   
13  This program is distributed in the hope that it will be useful, 
14  but WITHOUT ANY WARRANTY; without even the implied warranty of 
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
16  GNU General Public License for more details. 
17   
18  You should have received a copy of the GNU General Public License 
19  along with this program.  If not, see U{http://www.gnu.org/licenses/}. 
20  """ 
21   
22   
23   
24   
25  import pymbolic 
26   
27   
28   
29   
30 -def newton_interpolation_coefficients(x, y): 
31      assert len(x) == len(y) 
32      n = len(y) 
33      divided_differences = [y] 
34      last = y 
35   
36      for step in range(1, n): 
37          next = [(last[i+1]-last[i])/(x[i+step]-x[i]) 
38                  for i in range(n-step)] 
39          divided_differences.append(next) 
40          last = next 
41   
42      return [dd_col[-1] for dd_col in divided_differences] 
43   
44   
45   
46   
47 -def newton_interpolation_polynomial(x, y): 
48      coeff = newton_interpolation_coefficients(x, y) 
49   
50      var_x = pymbolic.var("x") 
51      linear_factors = [ 
52              pymbolic.Polynomial(var_x, ((0, pt), (1, 1))) 
53              for pt in x] 
54      pyramid_linear_factors = [1] 
55      for l in linear_factors: 
56          pyramid_linear_factors.append( 
57                  pyramid_linear_factors[-1]*l) 
58   
59      return pymbolic.linear_combination(coeff, pyramid_linear_factors) 
60   
61   
62   
63   
64 -def newton_interpolation_function(x, y): 
65      return pymbolic.compile(newton_interpolation_polynomial(x, y), ["x"]) 
66