Mercurial > hg > octave-nkf
annotate scripts/polynomial/residue.m @ 20818:9d2023d1a63c
binoinv.m: Implement binary search algorithm for 28X performance increase (bug #34363).
* binoinv.m: Call new functions scalar_binoinv or vector_binoinv to calculate
binoinv. If there are still uncalculated values then call bin_search_binoinv
to perform binary search for remaining values. Add more BIST tests.
* binoinv.m (scalar_binoinv): New subfunction to calculate binoinv for scalar x.
Stops when x > 1000.
* binoinv.m (vector_binoinv): New subfunction to calculate binoinv for scalar x.
Stops when x > 1000.
| author | Lachlan Andrew <lachlanbis@gmail.com> |
|---|---|
| date | Sun, 11 Oct 2015 19:49:40 -0700 |
| parents | 83792dd9bcc1 |
| children |
| rev | line source |
|---|---|
|
19898
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17744
diff
changeset
|
1 ## Copyright (C) 1994-2015 John W. Eaton |
| 6964 | 2 ## Copyright (C) 2007 Ben Abbott |
| 2313 | 3 ## |
| 4 ## This file is part of Octave. | |
| 5 ## | |
| 6 ## Octave is free software; you can redistribute it and/or modify it | |
| 7 ## under the terms of the GNU General Public License as published by | |
| 7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
| 9 ## your option) any later version. | |
| 2313 | 10 ## |
| 11 ## Octave is distributed in the hope that it will be useful, but | |
| 12 ## WITHOUT ANY WARRANTY; without even the implied warranty of | |
| 13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 14 ## General Public License for more details. | |
| 15 ## | |
| 16 ## You should have received a copy of the GNU General Public License | |
| 7016 | 17 ## along with Octave; see the file COPYING. If not, see |
| 18 ## <http://www.gnu.org/licenses/>. | |
| 904 | 19 |
| 3368 | 20 ## -*- texinfo -*- |
|
13929
9cae456085c2
Grammarcheck of documentation before 3.6.0 release.
Rik <octave@nomad.inbox5.com>
parents:
13128
diff
changeset
|
21 ## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a}) |
|
13127
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
22 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
23 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
24 ## The first calling form computes the partial fraction expansion for the |
|
435d1b905e31
Update residue.m's docstring to more modern form, all calling forms on top
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
11587
diff
changeset
|
25 ## quotient of the polynomials, @var{b} and @var{a}. |
|
20375
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
26 ## |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
27 ## The quotient is defined as |
| 3368 | 28 ## @tex |
| 29 ## $$ | |
| 6978 | 30 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
| 3368 | 31 ## + \sum_{i=1}^N k_i s^{N-i}. |
| 32 ## $$ | |
| 33 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
34 ## @ifnottex |
| 3426 | 35 ## |
| 3368 | 36 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
37 ## @group |
|
20375
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
38 ## B(s) M r(m) N |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
39 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
40 ## A(s) m=1 (s-p(m))^e(m) i=1 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
41 ## @end group |
| 3368 | 42 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
43 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
44 ## @end ifnottex |
| 3368 | 45 ## @noindent |
|
20375
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
46 ## where @math{M} is the number of poles (the length of the @var{r}, @var{p}, |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
47 ## and @var{e}), the @var{k} vector is a polynomial of order @math{N-1} |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
48 ## representing the direct contribution, and the @var{e} vector specifies the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
49 ## multiplicity of the m-th residue's pole. |
| 3426 | 50 ## |
| 6964 | 51 ## For example, |
| 3426 | 52 ## |
| 3368 | 53 ## @example |
| 54 ## @group | |
| 6964 | 55 ## b = [1, 1, 1]; |
| 56 ## a = [1, -5, 8, -4]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
57 ## [r, p, k, e] = residue (b, a) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
58 ## @result{} r = [-2; 7; 3] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
59 ## @result{} p = [2; 2; 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
60 ## @result{} k = [](0x0) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
61 ## @result{} e = [1; 2; 1] |
| 3368 | 62 ## @end group |
| 63 ## @end example | |
| 3426 | 64 ## |
| 3368 | 65 ## @noindent |
| 6978 | 66 ## which represents the following partial fraction expansion |
| 3368 | 67 ## @tex |
| 68 ## $$ | |
| 69 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} | |
| 70 ## $$ | |
| 71 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
72 ## @ifnottex |
| 3426 | 73 ## |
| 3368 | 74 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
75 ## @group |
| 3368 | 76 ## s^2 + s + 1 -2 7 3 |
| 2311 | 77 ## ------------------- = ----- + ------- + ----- |
| 78 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
79 ## @end group |
| 3368 | 80 ## @end example |
| 6964 | 81 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
82 ## @end ifnottex |
| 6978 | 83 ## |
|
20375
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
84 ## The second calling form performs the inverse operation and computes the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
85 ## reconstituted quotient of polynomials, @var{b}(s)/@var{a}(s), from the |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
86 ## partial fraction expansion; represented by the residues, poles, and a direct |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
87 ## polynomial specified by @var{r}, @var{p} and @var{k}, and the pole |
|
f1d0f506ee78
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
20038
diff
changeset
|
88 ## multiplicity @var{e}. |
| 7011 | 89 ## |
| 90 ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is | |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
91 ## determined by the function @code{mpoles}. |
| 6978 | 92 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
93 ## For example: |
| 6964 | 94 ## |
| 95 ## @example | |
| 96 ## @group | |
| 7011 | 97 ## r = [-2; 7; 3]; |
| 98 ## p = [2; 2; 1]; | |
| 99 ## k = [1, 0]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
100 ## [b, a] = residue (r, p, k) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
101 ## @result{} b = [1, -5, 9, -3, 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
102 ## @result{} a = [1, -5, 8, -4] |
| 7011 | 103 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
104 ## where mpoles is used to determine e = [1; 2; 1] |
| 7011 | 105 ## @end group |
| 106 ## @end example | |
| 107 ## | |
| 108 ## Alternatively the multiplicity may be defined explicitly, for example, | |
| 109 ## | |
| 110 ## @example | |
| 111 ## @group | |
| 112 ## r = [7; 3; -2]; | |
| 113 ## p = [2; 1; 2]; | |
| 114 ## k = [1, 0]; | |
| 115 ## e = [2; 1; 1]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
116 ## [b, a] = residue (r, p, k, e) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
117 ## @result{} b = [1, -5, 9, -3, 1] |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
118 ## @result{} a = [1, -5, 8, -4] |
| 6964 | 119 ## @end group |
| 120 ## @end example | |
| 121 ## | |
| 122 ## @noindent | |
| 6978 | 123 ## which represents the following partial fraction expansion |
| 6964 | 124 ## @tex |
| 125 ## $$ | |
| 6978 | 126 ## {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s = {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} |
| 6964 | 127 ## $$ |
| 128 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
129 ## @ifnottex |
| 6964 | 130 ## |
| 131 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
132 ## @group |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
133 ## -2 7 3 s^4 - 5s^3 + 9s^2 - 3s + 1 |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
134 ## ----- + ------- + ----- + s = -------------------------- |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
135 ## (s-2) (s-2)^2 (s-1) s^3 - 5s^2 + 8s - 4 |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
136 ## @end group |
| 6964 | 137 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
138 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
139 ## @end ifnottex |
|
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
13963
diff
changeset
|
140 ## @seealso{mpoles, poly, roots, conv, deconv} |
| 3368 | 141 ## @end deftypefn |
| 1025 | 142 |
| 3202 | 143 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
| 6964 | 144 ## Author: Ben Abbott <bpabbott@mac.com> |
| 2312 | 145 ## Created: June 1994 |
| 146 ## Adapted-By: jwe | |
| 559 | 147 |
| 6978 | 148 function [r, p, k, e] = residue (b, a, varargin) |
| 559 | 149 |
| 7011 | 150 if (nargin < 2 || nargin > 4) |
| 6046 | 151 print_usage (); |
| 559 | 152 endif |
| 153 | |
| 6964 | 154 toler = .001; |
| 155 | |
| 7011 | 156 if (nargin >= 3) |
| 157 if (nargin >= 4) | |
| 158 e = varargin{2}; | |
| 159 else | |
| 160 e = []; | |
| 161 endif | |
| 6964 | 162 ## The inputs are the residue, pole, and direct part. Solve for the |
| 163 ## corresponding numerator and denominator polynomials | |
| 7011 | 164 [r, p] = rresidue (b, a, varargin{1}, toler, e); |
|
17312
088d014a7fe2
Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents:
14868
diff
changeset
|
165 return; |
| 7011 | 166 endif |
| 559 | 167 |
| 2303 | 168 ## Make sure both polynomials are in reduced form. |
| 1025 | 169 |
| 170 a = polyreduce (a); | |
| 171 b = polyreduce (b); | |
| 559 | 172 |
|
20441
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20375
diff
changeset
|
173 b /= a(1); |
|
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20375
diff
changeset
|
174 a /= a(1); |
| 559 | 175 |
| 1025 | 176 la = length (a); |
| 177 lb = length (b); | |
| 559 | 178 |
| 2303 | 179 ## Handle special cases here. |
| 1025 | 180 |
| 181 if (la == 0 || lb == 0) | |
| 559 | 182 k = r = p = e = []; |
| 183 return; | |
| 184 elseif (la == 1) | |
| 1025 | 185 k = b / a; |
| 559 | 186 r = p = e = []; |
| 187 return; | |
| 188 endif | |
| 189 | |
| 2303 | 190 ## Find the poles. |
| 1025 | 191 |
| 192 p = roots (a); | |
| 193 lp = length (p); | |
| 559 | 194 |
| 6964 | 195 ## Sort poles so that multiplicity loop will work. |
| 196 | |
| 197 [e, indx] = mpoles (p, toler, 1); | |
|
20038
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
198 p = p(indx); |
| 559 | 199 |
| 7398 | 200 ## For each group of pole multiplicity, set the value of each |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
201 ## pole to the average of the group. This reduces the error in |
| 7398 | 202 ## the resulting poles. |
| 203 | |
| 204 p_group = cumsum (e == 1); | |
| 205 for ng = 1:p_group(end) | |
| 206 m = find (p_group == ng); | |
| 207 p(m) = mean (p(m)); | |
| 208 endfor | |
| 209 | |
| 2303 | 210 ## Find the direct term if there is one. |
| 1025 | 211 |
| 212 if (lb >= la) | |
| 6964 | 213 ## Also return the reduced numerator. |
| 1025 | 214 [k, b] = deconv (b, a); |
| 215 lb = length (b); | |
| 559 | 216 else |
| 217 k = []; | |
| 218 endif | |
| 219 | |
| 7398 | 220 ## Determine if the poles are (effectively) zero. |
| 221 | |
| 222 small = max (abs (p)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
223 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
224 small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
225 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
226 small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
227 endif |
| 7398 | 228 p(abs (p) < small) = 0; |
| 229 | |
| 230 ## Determine if the poles are (effectively) real, or imaginary. | |
| 231 | |
| 232 index = (abs (imag (p)) < small); | |
| 233 p(index) = real (p(index)); | |
| 234 index = (abs (real (p)) < small); | |
| 235 p(index) = 1i * imag (p(index)); | |
| 236 | |
| 237 ## The remainder determines the residues. The case of one pole | |
| 238 ## is trivial. | |
| 239 | |
| 1025 | 240 if (lp == 1) |
| 241 r = polyval (b, p); | |
| 559 | 242 return; |
| 243 endif | |
| 244 | |
| 6964 | 245 ## Determine the order of the denominator and remaining numerator. |
| 246 ## With the direct term removed the potential order of the numerator | |
| 247 ## is one less than the order of the denominator. | |
| 1025 | 248 |
| 6964 | 249 aorder = numel (a) - 1; |
| 250 border = aorder - 1; | |
| 1025 | 251 |
| 6964 | 252 ## Construct a system of equations relating the individual |
| 253 ## contributions from each residue to the complete numerator. | |
| 559 | 254 |
| 6964 | 255 A = zeros (border+1, border+1); |
| 256 B = prepad (reshape (b, [numel(b), 1]), border+1, 0); | |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14552
diff
changeset
|
257 for ip = 1:numel (p) |
| 6964 | 258 ri = zeros (size (p)); |
| 259 ri(ip) = 1; | |
| 260 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; | |
| 261 endfor | |
| 559 | 262 |
| 2303 | 263 ## Solve for the residues. |
| 1025 | 264 |
| 265 r = A \ B; | |
| 559 | 266 |
| 267 endfunction | |
| 6964 | 268 |
| 7011 | 269 function [pnum, pden, e] = rresidue (r, p, k, toler, e) |
| 6964 | 270 |
| 271 ## Reconstitute the numerator and denominator polynomials from the | |
| 272 ## residues, poles, and direct term. | |
| 273 | |
| 7011 | 274 if (nargin < 2 || nargin > 5) |
| 6964 | 275 print_usage (); |
| 276 endif | |
| 277 | |
| 7011 | 278 if (nargin < 5) |
| 279 e = []; | |
| 280 endif | |
| 281 | |
| 6964 | 282 if (nargin < 4) |
| 283 toler = []; | |
| 284 endif | |
| 285 | |
| 286 if (nargin < 3) | |
| 287 k = []; | |
| 288 endif | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
289 |
|
14552
86854d032a37
maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents:
14363
diff
changeset
|
290 if (numel (e)) |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14552
diff
changeset
|
291 indx = 1:numel (p); |
| 7011 | 292 else |
| 293 [e, indx] = mpoles (p, toler, 0); | |
|
20038
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
294 p = p(indx); |
|
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19898
diff
changeset
|
295 r = r(indx); |
| 7011 | 296 endif |
| 6964 | 297 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14552
diff
changeset
|
298 indx = 1:numel (p); |
| 6964 | 299 |
| 300 for n = indx | |
| 301 pn = [1, -p(n)]; | |
|
14552
86854d032a37
maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents:
14363
diff
changeset
|
302 if (n == 1) |
| 6964 | 303 pden = pn; |
| 304 else | |
| 305 pden = conv (pden, pn); | |
| 306 endif | |
| 307 endfor | |
| 308 | |
| 309 ## D is the order of the denominator | |
| 310 ## K is the order of the direct polynomial | |
| 311 ## N is the order of the resulting numerator | |
| 312 ## pnum(1:(N+1)) is the numerator's polynomial | |
| 313 ## pden(1:(D+1)) is the denominator's polynomial | |
| 314 ## pm is the multible pole for the nth residue | |
| 315 ## pn is the numerator contribution for the nth residue | |
| 316 | |
| 317 D = numel (pden) - 1; | |
| 318 K = numel (k) - 1; | |
| 319 N = K + D; | |
| 320 pnum = zeros (1, N+1); | |
| 7011 | 321 for n = indx(abs (r) > 0) |
| 6964 | 322 p1 = [1, -p(n)]; |
| 7011 | 323 for m = 1:e(n) |
| 324 if (m == 1) | |
| 6964 | 325 pm = p1; |
| 326 else | |
| 327 pm = conv (pm, p1); | |
| 328 endif | |
| 329 endfor | |
| 330 pn = deconv (pden, pm); | |
| 331 pn = r(n) * pn; | |
|
20441
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20375
diff
changeset
|
332 pnum += prepad (pn, N+1, 0, 2); |
| 6964 | 333 endfor |
| 334 | |
| 335 ## Add the direct term. | |
| 336 | |
| 337 if (numel (k)) | |
|
20441
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20375
diff
changeset
|
338 pnum += conv (pden, k); |
| 6964 | 339 endif |
| 340 | |
| 341 ## Check for leading zeros and trim the polynomial coefficients. | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
342 if (isa (r, "single") || isa (p, "single") || isa (k, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
343 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single"); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
344 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
345 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
346 endif |
| 6964 | 347 |
| 7011 | 348 pnum(abs (pnum) < small) = 0; |
| 349 pden(abs (pden) < small) = 0; | |
| 6964 | 350 |
| 351 pnum = polyreduce (pnum); | |
| 352 pden = polyreduce (pden); | |
| 353 | |
| 354 endfunction | |
| 6968 | 355 |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
356 |
| 6968 | 357 %!test |
| 358 %! b = [1, 1, 1]; | |
| 359 %! a = [1, -5, 8, -4]; | |
| 360 %! [r, p, k, e] = residue (b, a); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
361 %! assert (r, [-2; 7; 3], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
362 %! assert (p, [2; 2; 1], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
363 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
364 %! assert (e, [1; 2; 1]); |
| 6994 | 365 %! k = [1 0]; |
| 7011 | 366 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); |
| 367 %! a = a; | |
| 368 %! [br, ar] = residue (r, p, k); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
369 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
370 %! assert (ar, a, 1e-12); |
| 7011 | 371 %! [br, ar] = residue (r, p, k, e); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
372 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
373 %! assert (ar, a, 1e-12); |
| 6994 | 374 |
| 375 %!test | |
| 376 %! b = [1, 0, 1]; | |
| 377 %! a = [1, 0, 18, 0, 81]; | |
| 7398 | 378 %! [r, p, k, e] = residue (b, a); |
| 6994 | 379 %! r1 = [-5i; 12; +5i; 12]/54; |
| 380 %! p1 = [+3i; +3i; -3i; -3i]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
381 %! assert (r, r1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
382 %! assert (p, p1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
383 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
384 %! assert (e, [1; 2; 1; 2]); |
| 6994 | 385 %! [br, ar] = residue (r, p, k); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
386 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
387 %! assert (ar, a, 1e-12); |
| 7011 | 388 |
| 389 %!test | |
| 390 %! r = [7; 3; -2]; | |
| 391 %! p = [2; 1; 2]; | |
| 392 %! k = [1 0]; | |
| 393 %! e = [2; 1; 1]; | |
| 394 %! [b, a] = residue (r, p, k, e); | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
395 %! assert (b, [1, -5, 9, -3, 1], 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
396 %! assert (a, [1, -5, 8, -4], 1e-12); |
| 7011 | 397 %! [rr, pr, kr, er] = residue (b, a); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
398 %! [jnk, n] = mpoles (p); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
399 %! assert (rr, r(n), 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
400 %! assert (pr, p(n), 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
401 %! assert (kr, k, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
402 %! assert (er, e(n), 1e-12); |
| 7011 | 403 |
| 7188 | 404 %!test |
| 405 %! b = [1]; | |
| 406 %! a = [1, 10, 25]; | |
| 7398 | 407 %! [r, p, k, e] = residue (b, a); |
| 7188 | 408 %! r1 = [0; 1]; |
| 409 %! p1 = [-5; -5]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
410 %! assert (r, r1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
411 %! assert (p, p1, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
412 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
413 %! assert (e, [1; 2]); |
| 7188 | 414 %! [br, ar] = residue (r, p, k); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
415 %! assert (br, b, 1e-12); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
416 %! assert (ar, a, 1e-12); |
|
13128
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
417 |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
418 ## The following test is due to Bernard Grung (bug #34266) |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
419 %!xtest |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
420 %! z1 = 7.0372976777e6; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
421 %! p1 = -3.1415926536e9; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
422 %! p2 = -4.9964813512e8; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
423 %! r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
424 %! r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
425 %! r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
426 %! r4 = z1/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
427 %! r = [r1; r2; r3; r4]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
428 %! p = [p1; p2; 0; 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
429 %! k = []; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
430 %! e = [1; 1; 1; 2]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
431 %! b = [1, z1]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
432 %! a = [1, -(p1 + p2), p1*p2, 0, 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
433 %! [br, ar] = residue (r, p, k, e); |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
434 %! assert (br, b, 1e-8); |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
435 %! assert (ar, a, 1e-8); |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
436 |