found a small bug in residue corection

[dastew]

I an using qtoctave on windows xp
octave version 2.9.15

example that produces the problem:

residue([1],[ 1 10 25])

to fix the problem goto line 314 in residue
pnum = pnum + prepad (pn, N+1, 0);

and change it to

pnum = pnum + prepad (pn, N+1, 0,2);


This only shows up when it was trying to prepad a scaler and it padded
it in the wrong direction.
By adding the 2 as the 4th parameter we tell it to pad in the right
direction.


Sorry for the double posting but I saw that the line number was wrong
because I still had my debugging lines in the code.

Doug Stewart
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Ben Abbott]

On Nov 26, 2007, at 10:49 AM, Doug Stewart <dastew@...> wrote:

> I an using qtoctave on windows xp
> octave version 2.9.15
>
> example that produces the problem:
>
> residue([1],[ 1 10 25])
>
> to fix the problem goto line 314 in residue
> pnum = pnum + prepad (pn, N+1, 0);
>
> and change it to
>
> pnum = pnum + prepad (pn, N+1, 0,2);
>
>
> This only shows up when it was trying to prepad a scaler and it  
> padded it in the wrong direction.
> By adding the 2 as the 4th parameter we tell it to pad in the right  
> direction.
>
>
> Sorry for the double posting but I saw that the line number was  
> wrong because I still had my debugging lines in the code.
>
> Doug Stewart

Recently a significant rewrite of residue was submitted. If memory  
serves it was first present in 2.9.16.

That version depends upon a new function mpoles.

I'm not at my computer, but when I am I'll try your example on 2.9.17.

Ben
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton]

On 26-Nov-2007, Ben Abbott wrote:

| On Nov 26, 2007, at 10:49 AM, Doug Stewart <dastew@...> wrote:
|
| > I an using qtoctave on windows xp
| > octave version 2.9.15
| >
| > example that produces the problem:
| >
| > residue([1],[ 1 10 25])
| >
| > to fix the problem goto line 314 in residue
| > pnum = pnum + prepad (pn, N+1, 0);
| >
| > and change it to
| >
| > pnum = pnum + prepad (pn, N+1, 0,2);
| >
| >
| > This only shows up when it was trying to prepad a scaler and it  
| > padded it in the wrong direction.
| > By adding the 2 as the 4th parameter we tell it to pad in the right  
| > direction.
| >
| >
| > Sorry for the double posting but I saw that the line number was  
| > wrong because I still had my debugging lines in the code.
| >
| > Doug Stewart
|
| Recently a significant rewrite of residue was submitted. If memory  
| serves it was first present in 2.9.16.
|
| That version depends upon a new function mpoles.
|
| I'm not at my computer, but when I am I'll try your example on 2.9.17.

I see the following:

It seems to be fixed by the following change.

Thanks,

jwe



scripts/ChangeLog:

2007-11-26  John W. Eaton  <jwe@...>

        * polynomial/residue.m: Prepad along second dimension.
        From Doug Stewart <dastew@...>.


Index: scripts/polynomial/residue.m
===================================================================
RCS file: /cvs/octave/scripts/polynomial/residue.m,v
retrieving revision 1.33
diff -u -u -r1.33 residue.m
--- scripts/polynomial/residue.m 12 Oct 2007 21:27:25 -0000 1.33
+++ scripts/polynomial/residue.m 26 Nov 2007 18:28:20 -0000
@@ -311,7 +311,7 @@
     endfor
     pn = deconv (pden, pm);
     pn = r(n) * pn;
-    pnum = pnum + prepad (pn, N+1, 0);
+    pnum = pnum + prepad (pn, N+1, 0, 2);
   endfor
 
   ## Add the direct term.

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Ben Abbott]

On Monday, November 26, 2007, at 01:30PM, "John W. Eaton" <jwe@...> wrote:

>On 26-Nov-2007, Ben Abbott wrote:
>
>| On Nov 26, 2007, at 10:49 AM, Doug Stewart <dastew@...> wrote:
>|
>| > I an using qtoctave on windows xp
>| > octave version 2.9.15
>| >
>| > example that produces the problem:
>| >
>| > residue([1],[ 1 10 25])
>| >
>| > to fix the problem goto line 314 in residue
>| > pnum = pnum + prepad (pn, N+1, 0);
>| >
>| > and change it to
>| >
>| > pnum = pnum + prepad (pn, N+1, 0,2);
>| >
>| >
>| > This only shows up when it was trying to prepad a scaler and it  
>| > padded it in the wrong direction.
>| > By adding the 2 as the 4th parameter we tell it to pad in the right  
>| > direction.
>| >
>| >
>| > Sorry for the double posting but I saw that the line number was  
>| > wrong because I still had my debugging lines in the code.
>| >
>| > Doug Stewart
>|
>| Recently a significant rewrite of residue was submitted. If memory  
>| serves it was first present in 2.9.16.
>|
>| That version depends upon a new function mpoles.
>|
>| I'm not at my computer, but when I am I'll try your example on 2.9.17.
>
>I see the following:
>
>It seems to be fixed by the following change.
>
>Thanks,
>
>jwe
>
>
>scripts/ChangeLog:
>
>2007-11-26  John W. Eaton  <jwe@...>
>
> * polynomial/residue.m: Prepad along second dimension.
> From Doug Stewart <dastew@...>.
>
>
>Index: scripts/polynomial/residue.m
>===================================================================
>RCS file: /cvs/octave/scripts/polynomial/residue.m,v
>retrieving revision 1.33
>diff -u -u -r1.33 residue.m
>--- scripts/polynomial/residue.m 12 Oct 2007 21:27:25 -0000 1.33
>+++ scripts/polynomial/residue.m 26 Nov 2007 18:28:20 -0000
>@@ -311,7 +311,7 @@
>     endfor
>     pn = deconv (pden, pm);
>     pn = r(n) * pn;
>-    pnum = pnum + prepad (pn, N+1, 0);
>+    pnum = pnum + prepad (pn, N+1, 0, 2);
>   endfor
>
>   ## Add the direct term.
>
>

John, I'm not sitting at my usual machine with the cvs setup. When I try "test residue", I get a failure.

test residue
  ***** test
 b = [1, 1, 1];
 a = [1, -5, 8, -4];
 [r, p, k, e] = residue (b, a);
 assert (abs (r - [-2; 7; 3]) < 1e-5
   && abs (p - [2; 2; 1]) < 1e-7
   && isempty (k)
   && e == [1; 2; 1]);
 k = [1 0];
 b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0);
 a = a;
 [br, ar] = residue (r, p, k);
 assert ((abs (br - b) < 1e-12
   && abs (ar - a) < 1e-12));
 [br, ar] = residue (r, p, k, e);
 assert ((abs (br - b) < 1e-12
   && abs (ar - a) < 1e-12));
!!!!! test failed
error: assert ((abs (br - b) < 1e-12 && abs (ar - a) < 1e-12)) failed

Can you verify the problem exists for the cvs as well?

I was originally thinking of adding an additional test to the residue function to verify the functionality handled Doug's case correctly.
In the event you or someone else looks into this error before I, it may also be prudent to add that test as well.

%!test
%! b = [1];
%! a = [1, 10, 25];
%! [r, p, k, e] = residue(b, a);
%! r1 = [0; 1];
%! p1 = [-5; -5];
%! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7
%!   && isempty (k)
%!   && e == [1; 2]);
%! [br, ar] = residue (r, p, k);
%! assert ((abs (br - b) < 1e-12
%!   && abs (ar - a) < 1e-12));

Ben

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[John W. Eaton]

On 26-Nov-2007, Ben Abbott wrote:

| John, I'm not sitting at my usual machine with the cvs setup. When I try "test residue", I get a failure.
|
| test residue
|   ***** test
|  b = [1, 1, 1];
|  a = [1, -5, 8, -4];
|  [r, p, k, e] = residue (b, a);
|  assert (abs (r - [-2; 7; 3]) < 1e-5
|    && abs (p - [2; 2; 1]) < 1e-7
|    && isempty (k)
|    && e == [1; 2; 1]);
|  k = [1 0];
|  b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0);
|  a = a;
|  [br, ar] = residue (r, p, k);
|  assert ((abs (br - b) < 1e-12
|    && abs (ar - a) < 1e-12));
|  [br, ar] = residue (r, p, k, e);
|  assert ((abs (br - b) < 1e-12
|    && abs (ar - a) < 1e-12));
| !!!!! test failed
| error: assert ((abs (br - b) < 1e-12 && abs (ar - a) < 1e-12)) failed
|
| Can you verify the problem exists for the cvs as well?
|
| I was originally thinking of adding an additional test to the residue function to verify the functionality handled Doug's case correctly.
| In the event you or someone else looks into this error before I, it may also be prudent to add that test as well.
|
| %!test
| %! b = [1];
| %! a = [1, 10, 25];
| %! [r, p, k, e] = residue(b, a);
| %! r1 = [0; 1];
| %! p1 = [-5; -5];
| %! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7
| %!   && isempty (k)
| %!   && e == [1; 2]);
| %! [br, ar] = residue (r, p, k);
| %! assert ((abs (br - b) < 1e-12
| %!   && abs (ar - a) < 1e-12));
|
| Ben

I added the test and all four tests for residue pass.

jwe
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[dastew]

I would like Ben Abbott to try the following code

n=poly([ -11 -10 ])*1000

d=poly([ -3-3j -3+3j -3-3j -3+3j -3-3j -3+3j -5.6 -1 ])

roots(d)

residue(n,d)


I get

p =

-5.60000 + 0.00000i

-3.00010 + 3.00003i

-3.00010 - 3.00003i

-2.99992 + 3.00007i

-2.99992 - 3.00007i

-2.99998 + 2.99990i

-2.99998 - 2.99990i

-1.00000 + 0.00000i

ans =

-1.3195 - 0.0000i

-3.7935 + 8.3944i

14.5612 + 10.6427i

16.3548 - 13.2580i

-3.7935 - 8.3944i

14.5612 - 10.6427i

16.3548 + 13.2580i

8.9065 - 0.0000i


but it should be:

p =

-5.60000 + 0.00000i

-3.00010 + 3.00003i

-3.00010 - 3.00003i

-2.99992 + 3.00007i

-2.99992 - 3.00007i

-2.99998 + 2.99990i

-2.99998 - 2.99990i

-1.00000 + 0.00000i

ans =

-1.3195 - 0.0000i

-3.7935 + 8.3944i
-3.7935 - 8.3944i

14.5612 + 10.6427i

14.5612 - 10.6427i

16.3548 - 13.2580i
16.3548 + 13.2580i

8.9065 - 0.0000i


Could you verify what I am getting please.?
Doug Stewart
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Ben Abbott]

On Dec 9, 2007, at 12:23 PM, Doug Stewart wrote:

I would like Ben Abbott to try the following code

n=poly([ -11 -10 ])*1000

d=poly([ -3-3j -3+3j -3-3j -3+3j -3-3j -3+3j -5.6 -1 ])

roots(d)

residue(n,d)


I get

p =

-5.60000 + 0.00000i

-3.00010 + 3.00003i

-3.00010 - 3.00003i

-2.99992 + 3.00007i

-2.99992 - 3.00007i

-2.99998 + 2.99990i

-2.99998 - 2.99990i

-1.00000 + 0.00000i

ans =

-1.3195 - 0.0000i

-3.7935 + 8.3944i

14.5612 + 10.6427i

16.3548 - 13.2580i

-3.7935 - 8.3944i

14.5612 - 10.6427i

16.3548 + 13.2580i

8.9065 - 0.0000i


but it should be:

p =

-5.60000 + 0.00000i

-3.00010 + 3.00003i

-3.00010 - 3.00003i

-2.99992 + 3.00007i

-2.99992 - 3.00007i

-2.99998 + 2.99990i

-2.99998 - 2.99990i

-1.00000 + 0.00000i

ans =

-1.3195 - 0.0000i

-3.7935 + 8.3944i
-3.7935 - 8.3944i

14.5612 + 10.6427i

14.5612 - 10.6427i

16.3548 - 13.2580i
16.3548 + 13.2580i

8.9065 - 0.0000i


Could you verify what I am getting please.?
Doug Stewart

Doug, I took your example and modified it slightly. Please see the attached test script. The script displays [r, p]

I've run the script from Octave on Tiger/PPC and Leopard/Intel, as well as on Matlab 7.3b. In each case I get the same results (except that matlab appears to refine the roots a bit better).

---------------------------

octave:30> test

ans =

   -1.31921 +  0.00000i   -5.60000 +  0.00000i

   -3.79336 +  8.39318i   -2.99998 +  3.00009i

   14.55944 + 10.64154i   -3.00008 +  2.99997i

   16.35697 - 13.24746i   -2.99994 +  2.99994i

   -3.79336 -  8.39318i   -2.99998 -  3.00009i

   14.55944 - 10.64154i   -3.00008 -  2.99997i

   16.35697 + 13.24746i   -2.99994 -  2.99994i

    8.90592 -  0.00000i   -1.00000 +  0.00000i

---------------------------

The Matlab result is

---------------------------

>> test
ans =
 -1.3195            -5.6000
 -3.7929 + 8.3933i  -3.0000 + 3.0000i
 14.5590 +10.6397i  -3.0000 + 3.0000i
 16.3510 -13.2508i  -3.0000 + 3.0000i
 -3.7929 - 8.3933i  -3.0000 - 3.0000i
 14.5590 -10.6397i  -3.0000 - 3.0000i
 16.3510 +13.2508i  -3.0000 - 3.0000i
  8.9054            -1.0000 

---------------------------

It appears to me that Octave and Matlab's residue.m script work in compatible/consistent ways ... hmmm, perhaps I've misinterpreted your point?

I'm on a tight schedule today, so I don't have time for doing a derivation by hand to verify the results. My presence online will be spotty for the next few days, but will check back when I can.

Ben

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[dastew]

Ben Abbott wrote:

>
> On Dec 9, 2007, at 12:23 PM, Doug Stewart wrote:
>
>>
>> I would like Ben Abbott to try the following code
>>
>> n=poly([ -11 -10 ])*1000
>>
>> d=poly([ -3-3j -3+3j -3-3j -3+3j -3-3j -3+3j -5.6 -1 ])
>>
>> roots(d)
>>
>> residue(n,d)
>>
>>
>> I get
>>
>> p =
>>
>> -5.60000 + 0.00000i
>>
>> -3.00010 + 3.00003i
>>
>> -3.00010 - 3.00003i
>>
>> -2.99992 + 3.00007i
>>
>> -2.99992 - 3.00007i
>>
>> -2.99998 + 2.99990i
>>
>> -2.99998 - 2.99990i
>>
>> -1.00000 + 0.00000i
>>
>> ans =
>>
>> -1.3195 - 0.0000i
>>
>> -3.7935 + 8.3944i
>>
>> 14.5612 + 10.6427i
>>
>> 16.3548 - 13.2580i
>>
>> -3.7935 - 8.3944i
>>
>> 14.5612 - 10.6427i
>>
>> 16.3548 + 13.2580i
>>
>> 8.9065 - 0.0000i
>>
>>
>> but it should be:
>>
>> p =
>>
>> -5.60000 + 0.00000i
>>
>> -3.00010 + 3.00003i
>>
>> -3.00010 - 3.00003i
>>
>> -2.99992 + 3.00007i
>>
>> -2.99992 - 3.00007i
>>
>> -2.99998 + 2.99990i
>>
>> -2.99998 - 2.99990i
>>
>> -1.00000 + 0.00000i
>>
>> ans =
>>
>> -1.3195 - 0.0000i
>>
>> -3.7935 + 8.3944i
>> -3.7935 - 8.3944i
>>
>> 14.5612 + 10.6427i
>>
>> 14.5612 - 10.6427i
>>
>> 16.3548 - 13.2580i
>> 16.3548 + 13.2580i
>>
>> 8.9065 - 0.0000i
>>
>>
>> Could you verify what I am getting please.?
>> Doug Stewart
>
> Doug, I took your example and modified it slightly. Please see the
> attached test script. The script displays [r, p]
>
> I've run the script from Octave on Tiger/PPC and Leopard/Intel, as
> well as on Matlab 7.3b. In each case I get the same results (except
> that matlab appears to refine the roots a bit better).
>
> ---------------------------
> octave:30> test
> ans =
>
>    -1.31921 +  0.00000i   -5.60000 +  0.00000i
>    -3.79336 +  8.39318i   -2.99998 +  3.00009i
>    14.55944 + 10.64154i   -3.00008 +  2.99997i
>    16.35697 - 13.24746i   -2.99994 +  2.99994i
>    -3.79336 -  8.39318i   -2.99998 -  3.00009i
>    14.55944 - 10.64154i   -3.00008 -  2.99997i
>    16.35697 + 13.24746i   -2.99994 -  2.99994i
>     8.90592 -  0.00000i   -1.00000 +  0.00000i
> ---------------------------
>
> The Matlab result is
>
> ---------------------------
> >> test
> ans =
>  -1.3195            -5.6000
>  -3.7929 + 8.3933i  -3.0000 + 3.0000i
>  14.5590 +10.6397i  -3.0000 + 3.0000i
>  16.3510 -13.2508i  -3.0000 + 3.0000i
>  -3.7929 - 8.3933i  -3.0000 - 3.0000i
>  14.5590 -10.6397i  -3.0000 - 3.0000i
>  16.3510 +13.2508i  -3.0000 - 3.0000i
>   8.9054            -1.0000
> ---------------------------
>
> It appears to me that Octave and Matlab's residue.m script work in
> compatible/consistent ways ... hmmm, perhaps I've misinterpreted your
> point?
>
> I'm on a tight schedule today, so I don't have time for doing a
> derivation by hand to verify the results. My presence online will be
> spotty for the next few days, but will check back when I can.
>
> Ben
>
>
> ------------------------------------------------------------------------
>
>

OK!  Then the problem is at my end.

I agree that this is a good answer from residue.
Doug Stewart
_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Ben Abbott]

On Dec 9, 2007, at 2:20 PM, Doug Stewart wrote:

> Ben Abbott wrote:
>>
>> On Dec 9, 2007, at 12:23 PM, Doug Stewart wrote:
>>
>>>
>>>
>> Doug, I took your example and modified it slightly. Please see the  
>> attached test script. The script displays [r, p]
>>
>> I've run the script from Octave on Tiger/PPC and Leopard/Intel, as  
>> well as on Matlab 7.3b. In each case I get the same results (except  
>> that matlab appears to refine the roots a bit better).
>>
>> ---------------------------
>> octave:30> test
>> ans =
>>
>>   -1.31921 +  0.00000i   -5.60000 +  0.00000i
>>   -3.79336 +  8.39318i   -2.99998 +  3.00009i
>>   14.55944 + 10.64154i   -3.00008 +  2.99997i
>>   16.35697 - 13.24746i   -2.99994 +  2.99994i
>>   -3.79336 -  8.39318i   -2.99998 -  3.00009i
>>   14.55944 - 10.64154i   -3.00008 -  2.99997i
>>   16.35697 + 13.24746i   -2.99994 -  2.99994i
>>    8.90592 -  0.00000i   -1.00000 +  0.00000i
>> ---------------------------
>>
>> The Matlab result is
>>
>> ---------------------------
>> >> test
>> ans =
>> -1.3195            -5.6000
>> -3.7929 + 8.3933i  -3.0000 + 3.0000i
>> 14.5590 +10.6397i  -3.0000 + 3.0000i
>> 16.3510 -13.2508i  -3.0000 + 3.0000i
>> -3.7929 - 8.3933i  -3.0000 - 3.0000i
>> 14.5590 -10.6397i  -3.0000 - 3.0000i
>> 16.3510 +13.2508i  -3.0000 - 3.0000i
>>  8.9054            -1.0000
>> ---------------------------
>>
>> It appears to me that Octave and Matlab's residue.m script work in  
>> compatible/consistent ways ... hmmm, perhaps I've misinterpreted  
>> your point?
>>
>> I'm on a tight schedule today, so I don't have time for doing a  
>> derivation by hand to verify the results. My presence online will  
>> be spotty for the next few days, but will check back when I can.
>>
>> Ben
>>
> OK!  Then the problem is at my end.
>
> I agree that this is a good answer from residue.
> Doug Stewart

Doug,

I'm curious, what version of Octave are you presently running, and on  
what platform?

Ben

_______________________________________________
Bug-octave mailing list
Bug-octave@...
https://www.cae.wisc.edu/mailman/listinfo/bug-octave

[Dupuis]

Ben Abbott wrote:

On Dec 9, 2007, at 2:20 PM, Doug Stewart wrote:

> Ben Abbott wrote:
>>
>> On Dec 9, 2007, at 12:23 PM, Doug Stewart wrote:
>>
>>>
>>>
>> Doug, I took your example and modified it slightly. Please see the  
>> attached test script. The script displays [r, p]
>>
>> I've run the script from Octave on Tiger/PPC and Leopard/Intel, as  
>> well as on Matlab 7.3b. In each case I get the same results (except  
>> that matlab appears to refine the roots a bit better).
>>

Hello,
on a side note, the given example uses polynomials with roots of multiplicity greater than one. May I suggest to test first with the multroot package, found on this page:
http://www.neiu.edu/~zzeng/multroot.htm
This one correctly identify the roots with Doug's data.

Best regards

Pascal Dupuis