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