*Bug #13511 fixed - Polynomials: lcm with large values 39/14839/2
Paul Bignier [Wed, 9 Jul 2014 15:04:22 +0000 (17:04 +0200)]
gcd and bezout fixed as well (lcm calls them both).

Change-Id: Id88f4bba882f1b10bc36b5ae0d11ba8d246ae26b

17 files changed:
scilab/CHANGES_5.5.X
scilab/modules/elementary_functions/help/en_US/discrete/gcd.xml
scilab/modules/elementary_functions/help/en_US/discrete/lcm.xml
scilab/modules/elementary_functions/help/fr_FR/discrete/gcd.xml
scilab/modules/elementary_functions/help/fr_FR/discrete/lcm.xml
scilab/modules/overloading/macros/%s_bezout.sci [new file with mode: 0644]
scilab/modules/overloading/macros/%s_gcd.sci [new file with mode: 0644]
scilab/modules/overloading/macros/%s_lcm.sci [new file with mode: 0644]
scilab/modules/polynomials/help/en_US/bezout.xml
scilab/modules/polynomials/help/fr_FR/bezout.xml
scilab/modules/polynomials/macros/gcd.sci
scilab/modules/polynomials/macros/lcm.sci
scilab/modules/polynomials/sci_gateway/fortran/sci_f_bezout.f
scilab/modules/polynomials/tests/nonreg_tests/bug_12679.dia.ref
scilab/modules/polynomials/tests/nonreg_tests/bug_12679.tst
scilab/modules/polynomials/tests/nonreg_tests/bug_13511.dia.ref [new file with mode: 0644]
scilab/modules/polynomials/tests/nonreg_tests/bug_13511.tst [new file with mode: 0644]

index 06c3667..6b67a1c 100644 (file)
@@ -189,6 +189,8 @@ Scilab Bug Fixes
 
 * Bug #13510 fixed - Datatip callback cleared 'd'.
 
+* Bug #13511 fixed - lcm returned incorrect results when used with doubles.
+
 * Bug #13512 fixed - dae crashed if the evaluation function had wrong prototype.
 
 * Bug #13515 fixed - There were wrong results for matrix/hypermatrix with bitset function.
@@ -197,7 +199,7 @@ Scilab Bug Fixes
 
 * Bug #13527 fixed - hilb did not check correctly the type of first input argument.
 
-* Bug #13554 fixed - rubberbox returns wrong values.
+* Bug #13554 fixed - rubberbox returned wrong values.
 
 * Bug #13585 fixed - suitesparse 4.3.1 was not supported.
 
index b87d940..f07a8ef 100644 (file)
@@ -67,9 +67,6 @@
             The greatest common divisor of an array <literal>p</literal> of real numbersof real numbers can be obtained by
             converting it to a polynomial before calling <literal>gcd</literal>, through <code>p = inv_coeff(p, 0)</code>.
         </para>
-        <para>
-            If <literal>p</literal> is given as an integer double (type 1), then it is treated as an <literal>int32</literal>.
-        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
index 4021a39..7e3fe4a 100644 (file)
@@ -60,9 +60,6 @@
             The least common multiple of an array <literal>p</literal> of real numbers can be obtained by
             converting it to a polynomial before calling <literal>lcm</literal>, through <code>p = inv_coeff(p, 0)</code>.
         </para>
-        <para>
-            If <literal>p</literal> is given as an integer double (type 1), then it is treated as an <literal>int32</literal>.
-        </para>
     </refsection>
     <refsection>
         <title>Examples</title>
@@ -74,11 +71,11 @@ p = [s s*(s+1)^2 s^2*(s+2)];
 p.*fact, pp
 
 // Integer case
-V = int32([2^2*3^5, 2^3*3^2,2^2*3^4*5]);
+V = int32([2^2*3^5 2^3*3^2 2^2*3^4*5]);
 lcm(V)
 
 // Double case
-V = [2^2*3^5, 2^3*3^2,2^2*3^4*5];
+V = [2^2*3^5  2^3*3^2 2^2*3^4*5];
 lcm(V)
  ]]></programlisting>
     </refsection>
index 16d5123..0a06e7e 100644 (file)
@@ -49,9 +49,6 @@
             Le PGCD d'une matrice <literal>p</literal> de réels peut s'obtenir en la convertissant en polynôme
             avant d'appeler <literal>gcd</literal>, grâce à la commande <literal>p = inv_coeff(p, 0)</literal>.
         </para>
-        <para>
-            Si <literal>p</literal> est donné comme un flottant entier (type 1), alors il est traité comme un <literal>int32</literal>.
-        </para>
     </refsection>
     <refsection>
         <title>Exemples</title>
index 6c3e019..97a7305 100644 (file)
@@ -52,9 +52,6 @@
             Le plus petit commun multiple d'une matrice <literal>p</literal> de réels peut s'obtenir en la convertissant
             en polynôme avant d'appeler <literal>lcm</literal>, grâce à la commande <literal>p = inv_coeff(p, 0)</literal>.
         </para>
-        <para>
-            Si <literal>p</literal> est donné comme un flottant entier (type 1), alors il est traité comme un <literal>int32</literal>.
-        </para>
     </refsection>
     <refsection>
         <title>Exemples</title>
@@ -66,11 +63,11 @@ p = [s s*(s+1)^2 s^2*(s+2)];
 p.*fact, pp
 
 // Cas des entiers
-V = int32([2^2*3^5, 2^3*3^2,2^2*3^4*5]);
+V = int32([2^2*3^5 2^3*3^2 2^2*3^4*5]);
 lcm(V)
 
 // Cas des doubles
-V = [2^2*3^5, 2^3*3^2,2^2*3^4*5];
+V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
 lcm(V)
  ]]></programlisting>
     </refsection>
diff --git a/scilab/modules/overloading/macros/%s_bezout.sci b/scilab/modules/overloading/macros/%s_bezout.sci
new file mode 100644 (file)
index 0000000..641c26b
--- /dev/null
@@ -0,0 +1,29 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [g,u] = %s_bezout(a,b)
+    //   g = bezout(a,b) is the greatest common divisor of a and b.
+    //       a and b must contain non-negative integer scalars.
+    //   [g,U] = bezout(a,b) also returns a (2x2) unimodular matrix U such that:
+    //   [a,b]*U = [g,0].
+    //   These are useful for solving Diophantine equations and computing
+    //   Hermite transformations.
+
+    u = [1 0 a];
+    v = [0 1 b];
+    while v(3)<>0
+        q = floor(u(3)/v(3));
+        t = u - v*q;
+        u = v;
+        v = t;
+    end
+    g = u(3);
+    u=[u(1) -v(1);u(2) -v(2)]
+
+endfunction
diff --git a/scilab/modules/overloading/macros/%s_gcd.sci b/scilab/modules/overloading/macros/%s_gcd.sci
new file mode 100644 (file)
index 0000000..2480524
--- /dev/null
@@ -0,0 +1,40 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [x,uu]=%s_gcd(p)
+    //Given a polynomial vector p, [pgcd,u]=gcd(p) computes the gcd
+    //of components and a unimodular matrix (with polynomial inverse) u,
+    //with minimal degree such that [p1 p2]*u=[0 ... 0 pgcd]
+    //!
+
+    [lhs,rhs]=argn(0)
+    [m,n]=size(p);
+    mn=m*n
+    p=matrix(p,1,mn)
+    x=p(1);
+    uu=1
+    for l=2:mn,
+        [x,u]=bezout(x,p(l)),
+        if lhs==2 then
+            uu=[uu(:,1:l-2) uu(:,l-1)*u(1,[2 1])];uu(l,l-1:l)=u(2,[2 1]);
+        end
+    end,
+    if lhs==1 then return,end
+    for l=mn:-1:2
+        pivot=uu(l,l-1);
+        for k=l:mn
+            q=floor(uu(l,k)/pivot)
+            r=uu(l,k)-q*pivot
+            if q<>0 then
+                uu(1:l-1,k)=uu(1:l-1,k)-q*uu(1:l-1,l-1)
+                uu(l,k)=r;
+            end
+        end
+    end
+endfunction
diff --git a/scilab/modules/overloading/macros/%s_lcm.sci b/scilab/modules/overloading/macros/%s_lcm.sci
new file mode 100644 (file)
index 0000000..1ea9d93
--- /dev/null
@@ -0,0 +1,24 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - Serge Steer
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
+
+function [q,fact]=%s_lcm(p)
+    //p=lcm(p) computes the lcm of polynomial vector p
+    //[pp,fact]=lcm(p) computes besides the vector fact of factors
+    //such that  p.*fact=pp*ones(p)
+    //!
+
+    k=find(p==0)
+    if k<>[] then q=p(k(1)),fact=0*ones(p),fact(k)=1,return,end
+
+    q=p(1);
+    for k=2:size(p,"*")
+        q=floor(q/%s_gcd([q,p(k)]))*p(k);
+    end
+    fact=floor(q./p)
+endfunction
index fdbe231..01fb490 100644 (file)
@@ -25,7 +25,7 @@
             <varlistentry>
                 <term>p1, p2</term>
                 <listitem>
-                    <para>two real polynomials or two integer scalars (type equal to 8)</para>
+                    <para>two real polynomials or two integer scalars (type equal to 1 or 8)</para>
                 </listitem>
             </varlistentry>
             <varlistentry>
     <refsection>
         <title>Examples</title>
         <programlisting role="example"><![CDATA[
-// Polynomial case
-x = poly(0, 'x');
+// Polynomials
+x = poly(0,'x');
 p1 = (x+1)*(x-3)^5;
 p2 = (x-2)*(x-3)^3;
-[thegcd,U] = bezout(p1, p2)
+[pgcd,U] = bezout(p1, p2)
 det(U)
 clean([p1 p2]*U)
-thelcm = p1*U(1,2)
+ppcm = p1*U(1,2)
 lcm([p1 p2])
 
-// Integer case
+// Integers
 i1 = int32(2*3^5);
 i2 = int32(2^3*3^2);
 [thegcd,U] = bezout(i1, i2)
@@ -83,6 +83,15 @@ V = int32([2^2*3^5 2^3*3^2 2^2*3^4*5]);
 [thegcd,U] = gcd(V)
 V*U
 lcm(V)
+
+// Doubles
+i1 = 2*3^5;
+i2 = 2^3*3^2;
+[thegcd,U] = bezout(i1, i2)
+V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
+[thegcd,U] = gcd(V)
+V*U
+lcm(V)
  ]]></programlisting>
     </refsection>
     <refsection role="see also">
index 0611f09..56e7ec3 100644 (file)
@@ -14,7 +14,7 @@
             <varlistentry>
                 <term>p1, p2  </term>
                 <listitem>
-                    <para>deux polynômes réels ou deux entiers (type égal à 8)
+                    <para>deux polynômes réels ou deux entiers (type égal à 1 ou 8)
                     </para>
                 </listitem>
             </varlistentry>
     </refsection>
     <refsection>
         <title>Exemples</title>
-        <programlisting role="example"><![CDATA[ 
-// cas des polynômes
-x=poly(0,'x');
-p1=(x+1)*(x-3)^5;p2=(x-2)*(x-3)^3;
-[pgcd,U]=bezout(p1,p2) 
+        <programlisting role="example"><![CDATA[
+// Polynômes
+x = poly(0,'x');
+p1 = (x+1)*(x-3)^5;
+p2 = (x-2)*(x-3)^3;
+[pgcd,U] = bezout(p1, p2)
 det(U)
-clean([p1,p2]*U)
-ppcm=p1*U(1,2)
-lcm([p1,p2])
+clean([p1 p2]*U)
+ppcm = p1*U(1,2)
+lcm([p1 p2])
 
-// cas des entiers
-i1=int32(2*3^5); i2=int32(2^3*3^2);
-[thegcd,U]=bezout(i1,i2) 
-V=int32([2^2*3^5, 2^3*3^2,2^2*3^4*5]);
-[thegcd,U]=gcd(V)
+// Entiers
+i1 = int32(2*3^5);
+i2 = int32(2^3*3^2);
+[thegcd,U] = bezout(i1, i2)
+V = int32([2^2*3^5 2^3*3^2 2^2*3^4*5]);
+[thegcd,U] = gcd(V)
+V*U
+lcm(V)
+
+// Doubles
+i1 = 2*3^5;
+i2 = 2^3*3^2;
+[thegcd,U] = bezout(i1, i2)
+V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
+[thegcd,U] = gcd(V)
 V*U
 lcm(V)
  ]]></programlisting>
index ee9f47a..a957735 100644 (file)
@@ -18,15 +18,17 @@ function [x, uu] = gcd(p)
         error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "gcd", 1));
     end
 
+    [lhs,rhs]=argn(0);
+
     if type(p)==1 then
         if floor(p)<>p then
             error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "gcd", 1));
         else
-            p = iconvert(p,4);
+            if lhs==2 then [x,uu]=%s_gcd(p), else x=%s_gcd(p), end
+            return
         end
     end
 
-    [lhs,rhs]=argn(0)
     if type(p)==8 then
         if lhs==2 then [x,uu]=%i_gcd(p), else x=%i_gcd(p), end
         return
index f079000..35af760 100644 (file)
@@ -22,7 +22,8 @@ function [p, fact] = lcm(p)
         if floor(p)<>p then
             error(msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "lcm", 1));
         else
-            p = iconvert(p,4);
+            if argn(1)==2 then [p, fact] = %s_lcm(p), else p = %s_lcm(p), end
+            return
         end
     end
 
index 21cc216..f5509c1 100644 (file)
@@ -34,7 +34,7 @@ c
       ilb=iadr(lstk(top))
       ilbr=ilb
       if(istk(ilb).lt.0) ilb=iadr(istk(ilb+1))
-      if(istk(ilb).gt.2) then
+      if(istk(ilb).eq.1.or.istk(ilb).gt.2) then
          fun=-1
          call funnam(ids(1,pt+1),'bezout',ilb)
          return
index 1a295ae..e7ab07e 100644 (file)
@@ -27,7 +27,7 @@ V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
 V_int = int32(V);
 [thegcd, U] = gcd(V);
 [thegcd, U_int] = gcd(V);
-assert_checkequal(V*U, int32([0 0 36]));
+assert_checkequal(V*U, [0 0 36]);
 assert_checkequal(V_int*U_int, int32([0 0 36]));
 assert_checkequal(gcd(uint8([15 20])), uint8(5));
 assert_checkequal(gcd([iconvert(15, 4) iconvert(20, 4)]), int32(5));
@@ -46,7 +46,7 @@ assert_checkequal(pp, [2*s^2 + 5*s^3 + 4*s^4 + s^5]);
 // Normal behavior, with integers
 V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
 V_int = int32(V);
-assert_checkequal(lcm(V), int32(9720));
+assert_checkequal(lcm(V), 9720);
 assert_checkequal(lcm(V_int), int32(9720));
 // Trying to use booleans, strings or decimals should yield an error
 refMsg4 = msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "lcm", 1);
index 88f98dd..fbefa58 100644 (file)
@@ -30,7 +30,7 @@ V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
 V_int = int32(V);
 [thegcd, U] = gcd(V);
 [thegcd, U_int] = gcd(V);
-assert_checkequal(V*U, int32([0 0 36]));
+assert_checkequal(V*U, [0 0 36]);
 assert_checkequal(V_int*U_int, int32([0 0 36]));
 assert_checkequal(gcd(uint8([15 20])), uint8(5));
 assert_checkequal(gcd([iconvert(15, 4) iconvert(20, 4)]), int32(5));
@@ -52,7 +52,7 @@ assert_checkequal(pp, [2*s^2 + 5*s^3 + 4*s^4 + s^5]);
 // Normal behavior, with integers
 V = [2^2*3^5 2^3*3^2 2^2*3^4*5];
 V_int = int32(V);
-assert_checkequal(lcm(V), int32(9720));
+assert_checkequal(lcm(V), 9720);
 assert_checkequal(lcm(V_int), int32(9720));
 // Trying to use booleans, strings or decimals should yield an error
 refMsg4 = msprintf(_("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n"), "lcm", 1);
diff --git a/scilab/modules/polynomials/tests/nonreg_tests/bug_13511.dia.ref b/scilab/modules/polynomials/tests/nonreg_tests/bug_13511.dia.ref
new file mode 100644 (file)
index 0000000..dee3450
--- /dev/null
@@ -0,0 +1,17 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2014 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 13511 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=13511
+//
+// <-- Short Description -->
+// lcm yield incorrect results when used with doubles.
+assert_checkequal(lcm([96 6250 10000 18700]), 56100000);
diff --git a/scilab/modules/polynomials/tests/nonreg_tests/bug_13511.tst b/scilab/modules/polynomials/tests/nonreg_tests/bug_13511.tst
new file mode 100644 (file)
index 0000000..94c81a6
--- /dev/null
@@ -0,0 +1,18 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2014 - Scilab Enterprises - Paul Bignier
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+//
+// <-- CLI SHELL MODE -->
+//
+// <-- Non-regression test for bug 13511 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=13511
+//
+// <-- Short Description -->
+// lcm yield incorrect results when used with doubles.
+
+assert_checkequal(lcm([96 6250 10000 18700]), 56100000);