* Bug 16620 fixed: polyint() introduced
[scilab.git] / scilab / modules / polynomials / macros / polyint.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2020 - Samuel GOUGEON - Le Mans Universit√©
3 //
4 // This file is hereby licensed under the terms of the GNU GPL v2.0,
5 // pursuant to article 5.3.4 of the CeCILL v.2.1.
6 // This file was originally licensed under the terms of the CeCILL v2.1,
7 // and continues to be available under such terms.
8 // For more information, see the COPYING file which you should have received
9 // along with this program.
10
11 function M = polyint(M, cstes)
12     // M : array (from scalar to hypermatrix) of polynomials with real or
13     //     complex coefficients
14     // Cstes: Integration constants (real or complex): scalar, array of size
15     //     size(M), or list of mixed scalar or sized(M) arrays or
16     //     undefined elements (valued as 0).
17
18     // EXAMPLE
19     // =======
20     if argn(2)==0 then
21         m = mode(); mode(6)
22         M = [ -1+%z, 1+2*%z+%z^2 ; 2-%z^3, 6]
23         mprintf("\n%s", "polyint(M)")
24         I = polyint(M)
25         mprintf("\n%s", "polyint(M, 1)")
26         I = polyint(M, 1)
27         mprintf("\n%s", "polyint(M, [1 2 ; 3 4])")
28         I = polyint(M, [1 2 ; 3 4])
29         mode(m)
30         M = I
31         return
32     end
33
34     s = size(M)
35     n = length(M)
36
37     // CHECKING ARGUMENTS
38     // ==================
39     rhs = argn(2)
40     if type(M)<>2 then
41         msg = _("%s: Argument #%d: Polynomial expected.\n")
42         error(msprintf(msg, "polyint",1))
43     end
44     if ~isdef("cstes","l")
45         cstes = 0
46     end
47     if and(type(cstes)<>[1 15])
48         msg = _("%s: Argument #%d: numbers or list of numbers expected.\n")
49         error(msprintf(msg, "polyint",2))
50     end
51     if type(cstes) <> 15 then
52         cstes = list(cstes)
53     end
54     for i = 1:length(cstes)
55         k = cstes(i)
56         if isdef("k","l")
57             if type(k) <> 1
58                 msg = _("%s: Argument #%d(%d): Decimal or complex number expected.\n")
59                 error(msprintf(msg, "polyint",2,i))
60             end
61             if length(k)<>1 & or(size(k)<>s)
62                 msg = _("%s: Argument #%d(%d): Array of size [%s] expected.\n")
63                 ts = strcat(msprintf("%d\n",s(:))," ")
64                 error(msprintf(msg, "polyint", 2, i, ts))
65             end
66         end
67     end
68
69     // PROCESSING
70     // ==========
71     vn = poly(0, varn(M))
72     M = matrix(M, -1, 1)
73     for i = cstes
74         if ~isdef("i","l")
75             i = zeros(n,1)
76         elseif length(i)==1
77             i = i * ones(n,1)
78         else
79             i = matrix(i,-1,1)
80         end
81         M = M * vn
82         d = max(max(degree(M),0))
83         p = (0:d) .*. ones(M)
84         c = coeff(M)
85         k = find(p > 0)
86         c(k) = c(k) ./ p(k)
87         c(:,1) = c(:,1) + i
88         M = inv_coeff(c, d, varn(M))
89     end
90
91     M = matrix(M, s)
92 endfunction