* Bug #12186 fixed - Performances of the function 'horner' improved. 33/10233/3
Alain Lamy [Wed, 16 Jan 2013 17:02:29 +0000 (18:02 +0100)]
Change-Id: I2b6c02e510abaf342905277b4d51f95f45c921b4

scilab/CHANGES_5.4.X
scilab/modules/polynomials/macros/horner.sci
scilab/modules/polynomials/tests/unit_tests/horner.dia.ref [new file with mode: 0644]
scilab/modules/polynomials/tests/unit_tests/horner.tst [new file with mode: 0644]

index ec61a18..8ec5059 100644 (file)
@@ -271,6 +271,8 @@ Bug fixes
 
 * Bug #12184 fixed - Performances of the function 'derivat' improved.
 
+* Bug #12186 fixed - Performances of the function 'horner' improved.
+
 * Bug #12190 fixed - Description of sprspn updated in help page.
 
 * Bug #12204 fixed - Fix a typo in the French localization
index 1d719f2..5f663b7 100644 (file)
@@ -8,7 +8,7 @@
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 
-function r=horner(p,x)
+function [r] = horner(p,x)
 // horner(P,x) evaluates the polynomial or rational matrix P = P(s)
 // when the variable s of the polynomial is replaced by x
 // x can be a scalar or polynomial or rational matrix.
@@ -18,42 +18,103 @@ function r=horner(p,x)
 // To evaluate a rational matrix at given frequencies use
 // preferably the freq primitive ;
 // See also: freq, repfreq.
+// Improvements:
+// Special cases aded to improve efficiency:
+// - p = row vector, x = column vector
+// - p = column vector, x = row vector
+// - x = scalar
 //!
 //
-  if argn(2)<>2 then 
+  if (argn(2) <> 2) then
     error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),'horner',2))
   end
-  if x==[]|p==[] then r=[],return,end
-  tp=type(p)
-  if tp<=2 then
-    [m,n]=size(p)
-    if m==-1 then indef=%t,m=1,n=1,p=p+0;else indef=%f,end
-    if m*n==1 then //special case for a single polynomial
-      C=coeff(p)
-      r=C($)*ones(x)
-      for kk=degree(p):-1:1,r=r.*x+C(kk);end;
+
+  if (x == [] | p == []) then
+    r = []
+    return
+  end
+
+  tp = type(p)
+
+  if (tp <= 2) then
+    // tp <= 2 <=> matrix of reals, complexes or polynomials
+    [m,n] = size(p)
+
+    if (m == -1) then
+      indef=%t, m=1, n=1, p=p+0
     else
-      r=[]
-      for l=1:m
-       rk=[]
-       for k=1:n
-         plk=p(l,k)
-         d=degree(plk)
-         rlk=coeff(plk,d)*ones(x); // for the case horner(1,x)
-         for kk=1:d,
-           rlk=rlk.*x+coeff(plk,d-kk);
-         end;
-         rk=[rk rlk]
-       end
-       r=[r;rk]
+      indef=%f
+    end
+
+    [mx,nx] = size(x)
+
+    if (m*n == 1) then
+      // special case: p = 1x1 polynomial, x = matrix
+      cp = coeff(p)
+      r = cp($) * ones(x)
+      for (k = degree(p) : -1 : 1)
+        r = r .* x + cp(k)
+      end
+
+    elseif (n*mx == 1)
+      // p = one column, x = one row
+      nd = max(degree(p));
+      r = zeros(p) * x;
+      for (k = nd : -1: 0)
+        c = coeff(p, k);
+        r = r .* (ones(p) * x) + c * ones(x);
+      end
+
+    elseif (m*nx == 1)
+      // p = one row, x = one column
+      nd = max(degree(p));
+      r = x * zeros(p);
+      for (k = nd : -1: 0)
+        c = coeff(p, k);
+        r = r .* (x * ones(p))+ ones(x) * c;
+      end
+
+    elseif (mx*nx == 1)
+      // p = matrix, x = scalar
+      nd = max(degree(p));
+      r = zeros(p);
+      for (k = nd : -1: 0)
+        c = coeff(p, k);
+        r = r * x + c;
+      end
+
+    else
+      // other cases
+      r = []
+      for (l = 1 : m)
+        rk = []
+             for (k = 1 : n)
+               plk = p(l,k)
+               d = degree(plk)
+               rlk = coeff(plk,d) * ones(x); // for the case horner(1,x)
+               for (kk = 1 : d)
+                 rlk = rlk .* x + coeff(plk,d-kk)
+               end
+               rk = [rk, rlk]
+             end
+             r = [r; rk]
       end
     end
-    if indef then r=r*eye(),end
-  elseif typeof(p)=='rational' then
-    r=horner(p(2),x)./horner(p(3),x),
-  elseif tp==129 then
-    r=horner(p(:),x);r=r(1):r(2):r(3)
+
+    if (indef) then
+      r = r * eye()
+    end
+
+  elseif (typeof(p) == "rational") then
+    r = horner(p(2),x) ./ horner(p(3),x)
+
+  elseif (tp == 129) then
+    // implicit polynomial for indexing
+    r = horner(p(:),x)
+    r = r(1) : r(2) : r(3)
+
   else
     error(msprintf(gettext("%s: Unexpected type for input argument #%d.\n"),'horner',1))
   end
+
 endfunction
diff --git a/scilab/modules/polynomials/tests/unit_tests/horner.dia.ref b/scilab/modules/polynomials/tests/unit_tests/horner.dia.ref
new file mode 100644 (file)
index 0000000..acef3c3
--- /dev/null
@@ -0,0 +1,72 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - CNES - Alain Lamy
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- CLI SHELL MODE -->
+// --------------------------------
+// p = polynomial
+s=poly(0,"s");
+p0 = [1+s+3*s^2];
+p = p0;
+r1 = horner(p,[1]);
+assert_checkequal(r1, 5);
+r1 = horner(p,[2]);
+assert_checkequal(r1, 15);
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5 15]);
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31 53]);
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15;31,53]);
+// --------------------------------
+// p = row vector
+p = [p0,p0^2];
+r1 = horner(p,[1]);
+assert_checkequal(r1, [5 25]);
+r1 = horner(p,[2]);
+assert_checkequal(r1, [15 225]);
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5,15,25,225]);
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31,53,961,2809]);
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15,25,225;31,53,961,2809]);
+// --------------------------------
+// p = column vector
+p = [p0;p0^2];
+r1 = horner(p,[1]);
+assert_checkequal(r1, [5; 25]);
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5,15;25,225]);
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31,53;961,2809]);
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15;31,53;25,225;961,2809]);
+// --------------------------------
+// p = matrix
+p = [p0,p0^2; p0, p0^3; p0, p0^4];
+r1 = horner(p,[1]);
+assert_checkequal(r1, [5,25;5,125;5,625]);
+r1 = horner(p,[2]);
+assert_checkequal(r1, [15,225;15,3375;15,50625]);
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5,15,25,225;5,15,125,3375;5,15,625,50625]);
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31,53,961,2809;31,53,29791,148877;31,53,923521,7890481]);
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15,25,225;31,53,961,2809;5,15,125,3375;31,53,29791,148877;5,15,625,50625;31,53,923521,7890481]);
+// --------------------------------
+// p = rational
+p = (1 + p0 ) / (2 + p0^2 + p0^3);
+r1 = horner(p,[1]);
+assert_checkalmostequal(r1, 0.0394737, 10^-6);
+r1 = horner(p,[2]);
+assert_checkalmostequal(r1, 0.0044420, 10^-5);
+r1 = horner(p,[1,2]);
+assert_checkalmostequal(r1, [0.0394737,0.0044420], 10^-5);
+r1 = horner(p,[3,4]);
+assert_checkalmostequal(r1, [0.0010405,0.0003560], 10^-4);
+r1 = horner(p,[1,2;3,4]);
+assert_checkalmostequal(r1, [0.0394737,0.0044420;0.0010405,0.0003560], 10^-4);
diff --git a/scilab/modules/polynomials/tests/unit_tests/horner.tst b/scilab/modules/polynomials/tests/unit_tests/horner.tst
new file mode 100644 (file)
index 0000000..0c8a06c
--- /dev/null
@@ -0,0 +1,105 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2013 - CNES - Alain Lamy
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- CLI SHELL MODE -->
+
+// --------------------------------
+// p = polynomial
+s=poly(0,"s");
+p0 = [1+s+3*s^2];
+p = p0;
+
+r1 = horner(p,[1]);
+assert_checkequal(r1, 5);
+
+r1 = horner(p,[2]);
+assert_checkequal(r1, 15);
+
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5 15]);
+
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31 53]);
+
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15;31,53]);
+
+
+// --------------------------------
+// p = row vector
+p = [p0,p0^2];
+
+r1 = horner(p,[1]);
+assert_checkequal(r1, [5 25]);
+
+r1 = horner(p,[2]);
+assert_checkequal(r1, [15 225]);
+
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5,15,25,225]);
+
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31,53,961,2809]);
+
+
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15,25,225;31,53,961,2809]);
+
+// --------------------------------
+// p = column vector
+p = [p0;p0^2];
+
+r1 = horner(p,[1]);
+assert_checkequal(r1, [5; 25]);
+
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5,15;25,225]);
+
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31,53;961,2809]);
+
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15;31,53;25,225;961,2809]);
+
+// --------------------------------
+// p = matrix
+p = [p0,p0^2; p0, p0^3; p0, p0^4];
+
+r1 = horner(p,[1]);
+assert_checkequal(r1, [5,25;5,125;5,625]);
+
+r1 = horner(p,[2]);
+assert_checkequal(r1, [15,225;15,3375;15,50625]);
+
+r1 = horner(p,[1,2]);
+assert_checkequal(r1, [5,15,25,225;5,15,125,3375;5,15,625,50625]);
+
+r1 = horner(p,[3,4]);
+assert_checkequal(r1, [31,53,961,2809;31,53,29791,148877;31,53,923521,7890481]);
+
+r1 = horner(p,[1,2;3,4]);
+assert_checkequal(r1, [5,15,25,225;31,53,961,2809;5,15,125,3375;31,53,29791,148877;5,15,625,50625;31,53,923521,7890481]);
+
+// --------------------------------
+// p = rational
+p = (1 + p0 ) / (2 + p0^2 + p0^3);
+
+r1 = horner(p,[1]);
+assert_checkalmostequal(r1, 0.0394737, 10^-6);
+
+r1 = horner(p,[2]);
+assert_checkalmostequal(r1, 0.0044420, 10^-5);
+
+r1 = horner(p,[1,2]);
+assert_checkalmostequal(r1, [0.0394737,0.0044420], 10^-5);
+
+r1 = horner(p,[3,4]);
+assert_checkalmostequal(r1, [0.0010405,0.0003560], 10^-4);
+
+r1 = horner(p,[1,2;3,4]);
+assert_checkalmostequal(r1, [0.0394737,0.0044420;0.0010405,0.0003560], 10^-4);
+