* Bug 15058 fixed: gcd and lcm result could be puzzingly <0
[scilab.git] / scilab / modules / polynomials / macros / lcm.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) ????-2008 - INRIA
3 // Copyright (C) INRIA - Serge Steer
4 // Copyright (C) 2017 - Samuel GOUGEON : http://bugzilla.scilab.org/15058
5 //
6 // Copyright (C) 2012 - 2016 - Scilab Enterprises
7 //
8 // This file is hereby licensed under the terms of the GNU GPL v2.0,
9 // pursuant to article 5.3.4 of the CeCILL v.2.1.
10 // This file was originally licensed under the terms of the CeCILL v2.1,
11 // and continues to be available under such terms.
12 // For more information, see the COPYING file which you should have received
13 // along with this program.
14
15 function [p, fact] = lcm(p)
16     //p=lcm(p) computes the lcm of polynomial vector p
17     //[pp,fact]=lcm(p) computes besides the vector fact of factors
18     //such that  p.*fact=pp*ones(p)
19     //!
20
21     if type(p)<>1 & type(p)<>2 & type(p)<>8 then
22         msg = _("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n");
23         error(msprintf(msg, "lcm", 1));
24     end
25
26     if type(p)==1 then
27         if floor(p) <> p then
28             msg = _("%s: Wrong type for argument #%d: Integer array or Polynomial expected.\n");
29             error(msprintf(msg, "lcm", 1));
30         else
31             p = iconvert(p, 4);
32         end
33     end
34
35     // Integers:
36     if type(p)==8 then
37         if argn(1)==2 then
38             [p, fact] = %i_lcm(p);
39         else
40             p = %i_lcm(p);
41         end
42         if p~=[] & p<0 then
43             p = -p
44             fact = -fact;
45         end
46         return
47     end
48
49     // Polynomials:
50     [m, n] = size(p),
51     p = matrix(p, m*n, 1),
52     p0 = p(1);
53     fact = 1;
54     for l = 2:m*n,
55         [u, v] = simp(p0, p(l)),
56         p0 = p0*v,
57         fact = [v*fact, u],
58     end,
59     fact = matrix(fact, m, n),
60     p = p0;
61 endfunction
62
63 // ----------------------------------------------------------------------------
64
65 function [q, fact] = %i_lcm(p)
66     // p = lcm(p) computes the lcm of polynomial vector p
67     // [pp,fact]=lcm(p) computes besides the vector fact of factors
68     // such that  p.*fact=pp*ones(p)
69
70     k = find(p==0);
71     if k <> [] then
72         q = p(k(1));
73         fact = 0*ones(p);
74         fact(k) = 1;
75         return
76     end
77
78     q = p(1);
79     for k = 2:size(p,"*")
80         q = q / %i_gcd([q,p(k)]) * p(k);
81     end
82     fact = q ./ p;
83 endfunction