Mercurial > hg > octave-nkf
annotate scripts/polynomial/residue.m @ 16195:b52d2f9294b6
use class for reading lexer input
* input.h, input.cc (octave_base_input_reader, octave_terminal_reader,
octave_file_reader, octave_eval_string_reader): New classes.
(octave_gets): Now a member function of octave_base_reader.
(get_user_input): Delete function. Move actions to
octave_terminal_reader::get_input and
octave_eval_string_reader::get_input.
function. Call from octave_file_reader::get_input. Don't check
whether reading an eval string.
(current_eval_string, get_input_from_eval_string): Delete global
variables.
(get_debug_input): Check CURR_LEXER->input_from_eval_string instead of
global get_input_from_eval_string variable. Don't protect
get_input_from_eval_string.
* lex.h (octave_lexer::input_reader): New data member.
(octave_lexer::octave_lexer (void)): Initialize it.
(octave_lexer::octave_lexer (const std::string&),
octave_lexer::octave_lexer (FILE *)): New constructors.
* lex.h, lex.cc (octave_lexer::input_buffer::fill): New function.
(octave_lexer::input_buffer::read): Delete.
(octave_lexer::read): Rename from octave_read. Call reader to get
input, and then hand it to the input_buffer instead of asking the
input buffer to read input. Change all callers.
(octave_lexer::input_source, octave_lexer::input_from_eval_string):
New functions. Call octave_lexer::input_from_eval_string instead of
using get_input_from_eval_string.
* oct-parse.in.yy (octave_parser::octave_parser (FILE *),
octave_parser::octave_parser (const std::string&)): New constructors.
(parse_fcn_file): Pass FILE pointer to octave_parser constructor.
(eval_string): Pass string to octave_parser constructor instead of
setting global current_eval_string variable.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 05 Mar 2013 10:19:51 -0500 |
| parents | 5d3a684236b0 |
| children | 088d014a7fe2 |
| rev | line source |
|---|---|
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
14104
diff
changeset
|
1 ## Copyright (C) 1994-2012 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}. |
| 3368 | 26 ## @tex |
| 27 ## $$ | |
| 6978 | 28 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
| 3368 | 29 ## + \sum_{i=1}^N k_i s^{N-i}. |
| 30 ## $$ | |
| 31 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
32 ## @ifnottex |
| 3426 | 33 ## |
| 3368 | 34 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
35 ## @group |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
36 ## B(s) M r(m) N |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
37 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
38 ## 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
|
39 ## @end group |
| 3368 | 40 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
41 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
42 ## @end ifnottex |
| 3368 | 43 ## @noindent |
| 3499 | 44 ## where @math{M} is the number of poles (the length of the @var{r}, |
| 6978 | 45 ## @var{p}, and @var{e}), the @var{k} vector is a polynomial of order @math{N-1} |
| 46 ## representing the direct contribution, and the @var{e} vector specifies | |
| 8494 | 47 ## the multiplicity of the m-th residue's pole. |
| 3426 | 48 ## |
| 6964 | 49 ## For example, |
| 3426 | 50 ## |
| 3368 | 51 ## @example |
| 52 ## @group | |
| 6964 | 53 ## b = [1, 1, 1]; |
| 54 ## 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
|
55 ## [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
|
56 ## @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
|
57 ## @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
|
58 ## @result{} k = [](0x0) |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
59 ## @result{} e = [1; 2; 1] |
| 3368 | 60 ## @end group |
| 61 ## @end example | |
| 3426 | 62 ## |
| 3368 | 63 ## @noindent |
| 6978 | 64 ## which represents the following partial fraction expansion |
| 3368 | 65 ## @tex |
| 66 ## $$ | |
| 67 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} | |
| 68 ## $$ | |
| 69 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
70 ## @ifnottex |
| 3426 | 71 ## |
| 3368 | 72 ## @example |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
73 ## @group |
| 3368 | 74 ## s^2 + s + 1 -2 7 3 |
| 2311 | 75 ## ------------------- = ----- + ------- + ----- |
| 76 ## 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
|
77 ## @end group |
| 3368 | 78 ## @end example |
| 6964 | 79 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
80 ## @end ifnottex |
| 6978 | 81 ## |
|
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
|
82 ## The second calling form performs the inverse operation and computes |
|
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
|
83 ## the reconstituted quotient of polynomials, @var{b}(s)/@var{a}(s), |
|
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
|
84 ## from the partial fraction expansion; represented by the residues, |
|
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
|
85 ## poles, and a direct polynomial specified by @var{r}, @var{p} and |
|
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
|
86 ## @var{k}, and the pole multiplicity @var{e}. |
| 7011 | 87 ## |
| 88 ## 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
|
89 ## determined by the function @code{mpoles}. |
| 6978 | 90 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
91 ## For example: |
| 6964 | 92 ## |
| 93 ## @example | |
| 94 ## @group | |
| 7011 | 95 ## r = [-2; 7; 3]; |
| 96 ## p = [2; 2; 1]; | |
| 97 ## k = [1, 0]; | |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
98 ## [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
|
99 ## @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
|
100 ## @result{} a = [1, -5, 8, -4] |
| 7011 | 101 ## |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
102 ## where mpoles is used to determine e = [1; 2; 1] |
| 7011 | 103 ## @end group |
| 104 ## @end example | |
| 105 ## | |
| 106 ## Alternatively the multiplicity may be defined explicitly, for example, | |
| 107 ## | |
| 108 ## @example | |
| 109 ## @group | |
| 110 ## r = [7; 3; -2]; | |
| 111 ## p = [2; 1; 2]; | |
| 112 ## k = [1, 0]; | |
| 113 ## 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
|
114 ## [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
|
115 ## @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
|
116 ## @result{} a = [1, -5, 8, -4] |
| 6964 | 117 ## @end group |
| 118 ## @end example | |
| 119 ## | |
| 120 ## @noindent | |
| 6978 | 121 ## which represents the following partial fraction expansion |
| 6964 | 122 ## @tex |
| 123 ## $$ | |
| 6978 | 124 ## {-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 | 125 ## $$ |
| 126 ## @end tex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
127 ## @ifnottex |
| 6964 | 128 ## |
| 129 ## @example | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
130 ## @group |
|
14327
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
131 ## -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
|
132 ## ----- + ------- + ----- + s = -------------------------- |
|
4d917a6a858b
doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
133 ## (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
|
134 ## @end group |
| 6964 | 135 ## @end example |
|
10821
693e22af08ae
Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents:
10224
diff
changeset
|
136 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
137 ## @end ifnottex |
|
14104
614505385171
doc: Overhaul docstrings for polynomial functions.
Rik <octave@nomad.inbox5.com>
parents:
13963
diff
changeset
|
138 ## @seealso{mpoles, poly, roots, conv, deconv} |
| 3368 | 139 ## @end deftypefn |
| 1025 | 140 |
| 3202 | 141 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
| 6964 | 142 ## Author: Ben Abbott <bpabbott@mac.com> |
| 2312 | 143 ## Created: June 1994 |
| 144 ## Adapted-By: jwe | |
| 559 | 145 |
| 6978 | 146 function [r, p, k, e] = residue (b, a, varargin) |
| 559 | 147 |
| 7011 | 148 if (nargin < 2 || nargin > 4) |
| 6046 | 149 print_usage (); |
| 559 | 150 endif |
| 151 | |
| 6964 | 152 toler = .001; |
| 153 | |
| 7011 | 154 if (nargin >= 3) |
| 155 if (nargin >= 4) | |
| 156 e = varargin{2}; | |
| 157 else | |
| 158 e = []; | |
| 159 endif | |
| 6964 | 160 ## The inputs are the residue, pole, and direct part. Solve for the |
| 161 ## corresponding numerator and denominator polynomials | |
| 7011 | 162 [r, p] = rresidue (b, a, varargin{1}, toler, e); |
| 6964 | 163 return |
| 7011 | 164 endif |
| 559 | 165 |
| 2303 | 166 ## Make sure both polynomials are in reduced form. |
| 1025 | 167 |
| 168 a = polyreduce (a); | |
| 169 b = polyreduce (b); | |
| 559 | 170 |
| 1025 | 171 b = b / a(1); |
| 172 a = a / a(1); | |
| 559 | 173 |
| 1025 | 174 la = length (a); |
| 175 lb = length (b); | |
| 559 | 176 |
| 2303 | 177 ## Handle special cases here. |
| 1025 | 178 |
| 179 if (la == 0 || lb == 0) | |
| 559 | 180 k = r = p = e = []; |
| 181 return; | |
| 182 elseif (la == 1) | |
| 1025 | 183 k = b / a; |
| 559 | 184 r = p = e = []; |
| 185 return; | |
| 186 endif | |
| 187 | |
| 2303 | 188 ## Find the poles. |
| 1025 | 189 |
| 190 p = roots (a); | |
| 191 lp = length (p); | |
| 559 | 192 |
| 6964 | 193 ## Sort poles so that multiplicity loop will work. |
| 194 | |
| 195 [e, indx] = mpoles (p, toler, 1); | |
| 196 p = p (indx); | |
| 559 | 197 |
| 7398 | 198 ## 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
|
199 ## pole to the average of the group. This reduces the error in |
| 7398 | 200 ## the resulting poles. |
| 201 | |
| 202 p_group = cumsum (e == 1); | |
| 203 for ng = 1:p_group(end) | |
| 204 m = find (p_group == ng); | |
| 205 p(m) = mean (p(m)); | |
| 206 endfor | |
| 207 | |
| 2303 | 208 ## Find the direct term if there is one. |
| 1025 | 209 |
| 210 if (lb >= la) | |
| 6964 | 211 ## Also return the reduced numerator. |
| 1025 | 212 [k, b] = deconv (b, a); |
| 213 lb = length (b); | |
| 559 | 214 else |
| 215 k = []; | |
| 216 endif | |
| 217 | |
| 7398 | 218 ## Determine if the poles are (effectively) zero. |
| 219 | |
| 220 small = max (abs (p)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
221 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
222 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
|
223 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
224 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
|
225 endif |
| 7398 | 226 p(abs (p) < small) = 0; |
| 227 | |
| 228 ## Determine if the poles are (effectively) real, or imaginary. | |
| 229 | |
| 230 index = (abs (imag (p)) < small); | |
| 231 p(index) = real (p(index)); | |
| 232 index = (abs (real (p)) < small); | |
| 233 p(index) = 1i * imag (p(index)); | |
| 234 | |
| 235 ## The remainder determines the residues. The case of one pole | |
| 236 ## is trivial. | |
| 237 | |
| 1025 | 238 if (lp == 1) |
| 239 r = polyval (b, p); | |
| 559 | 240 return; |
| 241 endif | |
| 242 | |
| 6964 | 243 ## Determine the order of the denominator and remaining numerator. |
| 244 ## With the direct term removed the potential order of the numerator | |
| 245 ## is one less than the order of the denominator. | |
| 1025 | 246 |
| 6964 | 247 aorder = numel (a) - 1; |
| 248 border = aorder - 1; | |
| 1025 | 249 |
| 6964 | 250 ## Construct a system of equations relating the individual |
| 251 ## contributions from each residue to the complete numerator. | |
| 559 | 252 |
| 6964 | 253 A = zeros (border+1, border+1); |
| 254 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
|
255 for ip = 1:numel (p) |
| 6964 | 256 ri = zeros (size (p)); |
| 257 ri(ip) = 1; | |
| 258 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; | |
| 259 endfor | |
| 559 | 260 |
| 2303 | 261 ## Solve for the residues. |
| 1025 | 262 |
| 263 r = A \ B; | |
| 559 | 264 |
| 265 endfunction | |
| 6964 | 266 |
| 7011 | 267 function [pnum, pden, e] = rresidue (r, p, k, toler, e) |
| 6964 | 268 |
| 269 ## Reconstitute the numerator and denominator polynomials from the | |
| 270 ## residues, poles, and direct term. | |
| 271 | |
| 7011 | 272 if (nargin < 2 || nargin > 5) |
| 6964 | 273 print_usage (); |
| 274 endif | |
| 275 | |
| 7011 | 276 if (nargin < 5) |
| 277 e = []; | |
| 278 endif | |
| 279 | |
| 6964 | 280 if (nargin < 4) |
| 281 toler = []; | |
| 282 endif | |
| 283 | |
| 284 if (nargin < 3) | |
| 285 k = []; | |
| 286 endif | |
|
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
287 |
|
14552
86854d032a37
maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents:
14363
diff
changeset
|
288 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
|
289 indx = 1:numel (p); |
| 7011 | 290 else |
| 291 [e, indx] = mpoles (p, toler, 0); | |
| 292 p = p (indx); | |
| 293 r = r (indx); | |
| 294 endif | |
| 6964 | 295 |
|
14868
5d3a684236b0
maint: Use Octave coding conventions for cuddling parentheses in scripts directory
Rik <octave@nomad.inbox5.com>
parents:
14552
diff
changeset
|
296 indx = 1:numel (p); |
| 6964 | 297 |
| 298 for n = indx | |
| 299 pn = [1, -p(n)]; | |
|
14552
86854d032a37
maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents:
14363
diff
changeset
|
300 if (n == 1) |
| 6964 | 301 pden = pn; |
| 302 else | |
| 303 pden = conv (pden, pn); | |
| 304 endif | |
| 305 endfor | |
| 306 | |
| 307 ## D is the order of the denominator | |
| 308 ## K is the order of the direct polynomial | |
| 309 ## N is the order of the resulting numerator | |
| 310 ## pnum(1:(N+1)) is the numerator's polynomial | |
| 311 ## pden(1:(D+1)) is the denominator's polynomial | |
| 312 ## pm is the multible pole for the nth residue | |
| 313 ## pn is the numerator contribution for the nth residue | |
| 314 | |
| 315 D = numel (pden) - 1; | |
| 316 K = numel (k) - 1; | |
| 317 N = K + D; | |
| 318 pnum = zeros (1, N+1); | |
| 7011 | 319 for n = indx(abs (r) > 0) |
| 6964 | 320 p1 = [1, -p(n)]; |
| 7011 | 321 for m = 1:e(n) |
| 322 if (m == 1) | |
| 6964 | 323 pm = p1; |
| 324 else | |
| 325 pm = conv (pm, p1); | |
| 326 endif | |
| 327 endfor | |
| 328 pn = deconv (pden, pm); | |
| 329 pn = r(n) * pn; | |
| 7183 | 330 pnum = pnum + prepad (pn, N+1, 0, 2); |
| 6964 | 331 endfor |
| 332 | |
| 333 ## Add the direct term. | |
| 334 | |
| 335 if (numel (k)) | |
| 336 pnum = pnum + conv (pden, k); | |
| 337 endif | |
| 338 | |
| 339 ## 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
|
340 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
|
341 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
|
342 else |
|
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; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
344 endif |
| 6964 | 345 |
| 7011 | 346 pnum(abs (pnum) < small) = 0; |
| 347 pden(abs (pden) < small) = 0; | |
| 6964 | 348 |
| 349 pnum = polyreduce (pnum); | |
| 350 pden = polyreduce (pden); | |
| 351 | |
| 352 endfunction | |
| 6968 | 353 |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
354 |
| 6968 | 355 %!test |
| 356 %! b = [1, 1, 1]; | |
| 357 %! a = [1, -5, 8, -4]; | |
| 358 %! [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
|
359 %! 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
|
360 %! 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
|
361 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
362 %! assert (e, [1; 2; 1]); |
| 6994 | 363 %! k = [1 0]; |
| 7011 | 364 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); |
| 365 %! a = a; | |
| 366 %! [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
|
367 %! 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
|
368 %! assert (ar, a, 1e-12); |
| 7011 | 369 %! [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
|
370 %! 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
|
371 %! assert (ar, a, 1e-12); |
| 6994 | 372 |
| 373 %!test | |
| 374 %! b = [1, 0, 1]; | |
| 375 %! a = [1, 0, 18, 0, 81]; | |
| 7398 | 376 %! [r, p, k, e] = residue (b, a); |
| 6994 | 377 %! r1 = [-5i; 12; +5i; 12]/54; |
| 378 %! 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
|
379 %! 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
|
380 %! 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
|
381 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
382 %! assert (e, [1; 2; 1; 2]); |
| 6994 | 383 %! [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
|
384 %! 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
|
385 %! assert (ar, a, 1e-12); |
| 7011 | 386 |
| 387 %!test | |
| 388 %! r = [7; 3; -2]; | |
| 389 %! p = [2; 1; 2]; | |
| 390 %! k = [1 0]; | |
| 391 %! e = [2; 1; 1]; | |
| 392 %! [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
|
393 %! 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
|
394 %! assert (a, [1, -5, 8, -4], 1e-12); |
| 7011 | 395 %! [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
|
396 %! [jnk, n] = mpoles (p); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
397 %! 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
|
398 %! 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
|
399 %! 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
|
400 %! assert (er, e(n), 1e-12); |
| 7011 | 401 |
| 7188 | 402 %!test |
| 403 %! b = [1]; | |
| 404 %! a = [1, 10, 25]; | |
| 7398 | 405 %! [r, p, k, e] = residue (b, a); |
| 7188 | 406 %! r1 = [0; 1]; |
| 407 %! p1 = [-5; -5]; | |
|
14363
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
408 %! 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
|
409 %! 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
|
410 %! assert (isempty (k)); |
|
f3d52523cde1
Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents:
14327
diff
changeset
|
411 %! assert (e, [1; 2]); |
| 7188 | 412 %! [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
|
413 %! 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
|
414 %! assert (ar, a, 1e-12); |
|
13128
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
415 |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
416 ## 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
|
417 %!xtest |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
418 %! z1 = 7.0372976777e6; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
419 %! p1 = -3.1415926536e9; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
420 %! p2 = -4.9964813512e8; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
421 %! 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
|
422 %! 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
|
423 %! 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
|
424 %! r4 = z1/p2/p1; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
425 %! r = [r1; r2; r3; r4]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
426 %! p = [p1; p2; 0; 0]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
427 %! k = []; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
428 %! e = [1; 1; 1; 2]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
429 %! b = [1, z1]; |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
430 %! 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
|
431 %! [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
|
432 %! assert (br, b, 1e-8); |
|
d049192e5d15
Add test f for bug #34266
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13127
diff
changeset
|
433 %! 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
|
434 |