forked from fricas/fricas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rinterp.spad
144 lines (119 loc) · 4.62 KB
/
rinterp.spad
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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
)if false
\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra rinterp.spad}
\author{Martin Rubey}
\maketitle
\begin{abstract}
Rational Interpolation
\end{abstract}
\eject
\section{Introduction}
This file contains a crude na\"ive implementation of rational interpolation,
where the coefficients of the rational function are in any given field.
\section{Questions and Outlook}
\begin{itemize}
\item Maybe this file should be joined with pinterp.spad, where polynomial
Lagrange interpolation is implemented. I have a second version that parallels
the structure of pinterp.spad closely.
\item There are probably better ways to implement rational interpolation. Maybe
{http://www.cs.ucsb.edu/~omer/personal/abstracts/rational.html} contains
something useful, but I don't know.
\item Comments welcome!
\end{itemize}
\section{RationalInterpolation}
)endif
)abbrev package RINTERP RationalInterpolation
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolation(xx, F) : Exports == Implementation where
xx : Symbol
F : Field
Exports == with
interpolate : (List F, List F, NonNegativeInteger,
NonNegativeInteger) -> Fraction Polynomial F
-- The implementation sets up a system of linear equations and solves it.
Implementation == add
interpolate(xlist, ylist, m, k) ==
)if false
First we check whether we have the right number of points and values. Clearly
the number of points and the number of values must be identical. Note that we
want to determine the numerator and denominator polynomials only up to a
factor. Thus, we want to determine $m+k+1$ coefficients, where $m$ is the degree
of the polynomial in the numerator and $k$ is the degree of the polynomial in
the denominator.
In fact, we could also leave -- for example -- $k$ unspecified and determine it
as $k=[[#xlist]]-m-1$: I don't know whether this would be better.
)endif
#xlist ~= #ylist =>
error "Different number of points and values."
#xlist ~= m+k+1 =>
error "wrong number of points"
)if false
The next step is to set up the matrix. Suppose that our numerator polynomial is
$p(x)=a_0+a_1x+\dots+a_mx^m$ and that our denominator polynomial is
$q(x)=b_0+b_1x+\dots+b_mx^m$. Then we have the following equations, writing $n$
for $m+k+1$:
\noindent
$$
\begin{array}{rl}
p(x_1)-y_1q(x_1)&=a_0+a_1x_1+\dots +a_mx_1^m-y_1(b_0+b_1x_1+\dots +b_kx_1^k)=0\\
p(x_2)-y_2q(x_2)&=a_0+a_1x_2+\dots +a_mx_2^m-y_2(b_0+b_1x_2+\dots +b_kx_2^k)=0\\
&\;\;\vdots\\
p(x_n)-y_nq(x_n)&=a_0+a_1x_n+\dots +a_mx_n^m-y_n(b_0+b_1x_n+\dots +b_kx_n^k)=0
\end{array}
$$
This can be written as
$$
\left[
\begin{array}{cccccccc}
1&x_1&\dots&x_1^m&-y_1&-y_1x_1&\dots&-y_1x_1^k\\
1&x_2&\dots&x_2^m&-y_2&-y_2x_2&\dots&-y_2x_2^k\\
&&&\vdots&&&&\\
1&x_n&\dots&x_n^m&-y_n&-y_nx_n&\dots&-y_nx_2^k
\end{array}
\right]
\left[
\begin{array}{c}
a_0\\a_1\\\vdots\\a_m\\b_0\\b_1\\\vdots\\b_k
\end{array}
\right]
=\mathbf 0
$$
We generate this matrix columnwise:
)endif
tempvec : List F := [1 for i in 1..(m+k+1)]
collist : List List F := cons(tempvec,
[(tempvec := [tempvec.i * xlist.i _
for i in 1..(m+k+1)]) _
for j in 1..max(m, k)])
collist := append([collist.j for j in 1..(m+1)], _
[[- collist.j.i * ylist.i for i in 1..(m+k+1)] _
for j in 1..(k+1)])
-- Now we can solve the system:
res : List Vector F := nullSpace((transpose matrix collist) _
::Matrix F)
)if false
Note that it may happen that the system has several solutions. In this case,
some of the data points may not be interpolated correctly. However, the
solution is often still useful, thus we do not signal an error.
)endif
if #res ~= 1 then output("Warning: unattainable points!"
)$OutputPackage
)if false
In this situation, all the solutions will be equivalent, thus we can always
simply take the first one:
)endif
reslist : List List Polynomial F := _
[[(res.1).(i+1)*(xx::Polynomial F)^i for i in 0..m], _
[(res.1).(i+m+2)*(xx::Polynomial F)^i for i in 0..k]]
-- Finally, we generate the rational function:
reduce((_+), reslist.1)/reduce((_+), reslist.2)
)if false
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}
)endif