--- a/scripts/specfun/betainc.m	Fri Feb 23 20:30:28 2018 +0100
+++ b/scripts/specfun/betainc.m	Sun Feb 25 00:02:44 2018 +0100
@@ -70,18 +70,18 @@
 
 function [y] = betainc (x, a, b, tail = "lower")
 
-  if (nargin > 4 || nargin < 3)
+  if ((nargin > 4) || (nargin < 3))
     print_usage ();
   endif
 
-  if (! isscalar (x) || ! isscalar (a) || ! isscalar (b))
+  if ((! isscalar (x)) || (! isscalar (a)) || (! isscalar (b)))
     [err, x, a, b] = common_size (x, a, b);
     if (err > 0)
       error ("betainc: x, a and b must be of common size or scalars");
     endif
   endif
 
-  if (iscomplex (x) || iscomplex (a) || iscomplex (b))
+  if ((iscomplex (x)) || (iscomplex (a)) || (iscomplex (b)))
     error ("betainc: inputs must be real or integer");
   endif
 
@@ -93,7 +93,7 @@
     error ("betainc: b must be strictly positive");
   endif
 
-  if (any (x > 1 | x < 0))
+  if (any ((x > 1) | (x < 0)))
     error ("betainc: x must be between 0 and 1");
   endif
 
@@ -117,7 +117,7 @@
 
   # If any of the argument is single, also the output should be
 
-  if (strcmpi (class (y), "single") || strcmpi (class (a), "single") || strcmpi (class (b), "single"))
+  if ((strcmpi (class (y), "single")) || (strcmpi (class (a), "single")) || (strcmpi (class (b), "single")))
     a = single (a);
     b = single (b);
     x = single (x);
@@ -125,21 +125,21 @@
   endif
 
   # In the following, we use the fact that the continued fraction we will
-  # use is more efficient when x <= a / (a+b).
+  # use is more efficient when x <= a / (a + b).
   # Moreover, to compute the upper version, which is defined as
   # I_x(a,b,"upper") = 1 - I_x(a,b) we used the property
   # I_x(a,b) + I_(1-x) (b,a) = 1.
 
   if (strcmpi (tail, "upper"))
-    flag = (x < a./(a+b));
-    x(!flag) = 1 - x(!flag);
-    [a(!flag), b(!flag)] = deal (b(!flag), a(!flag));
+    fflag = (x < (a ./ (a + b)));
+    x(! fflag) = 1 - x(! fflag);
+    [a(! fflag), b(! fflag)] = deal (b(! fflag), a(! fflag));
   elseif (strcmpi (tail, "lower"))
-    flag = (x > a./(a+b));
-    x (flag) = 1 - x(flag);
-    [a(flag), b(flag)] = deal (b(flag), a(flag));
+    fflag = (x > a./(a+b));
+    x (fflag) = 1 - x(fflag);
+    [a(fflag), b(fflag)] = deal (b(fflag), a(fflag));
   else
-    error ("betainc: invalid value for flag")
+    error ("betainc: invalid value for fflag")
   endif
 
   f = zeros (size (x), class (x));
@@ -155,7 +155,7 @@
   y = a .* log (x) + b .* log1p (-x) + gammaln (a + b) - ...
     gammaln (a) - gammaln (b) + log (f);
   y = real (exp (y));
-  y(flag) = 1 - y(flag);
+  y(fflag) = 1 - y(fflag);
   y = reshape (y, sz);
 
 endfunction

--- a/scripts/specfun/cosint.m	Fri Feb 23 20:30:28 2018 +0100
+++ b/scripts/specfun/cosint.m	Sun Feb 25 00:02:44 2018 +0100
@@ -114,7 +114,7 @@
   y(i_miss & flag_large) = yy;
 
   ## For values small in modulus, use the series expansion (also near (-oo, 0])
-  i_miss = ((i_miss) & (!flag_large));
+  i_miss = ((i_miss) & (! flag_large));
   if (iscomplex (x))
     ## indexing can lose imag part: if it was -0, we could end up on the
     ## wrong right side of the branch cut along the negative real axis.

--- a/scripts/specfun/expint.m	Fri Feb 23 20:30:28 2018 +0100
+++ b/scripts/specfun/expint.m	Sun Feb 25 00:02:44 2018 +0100
@@ -1,5 +1,5 @@
 ## Copyright (C) 2006-2017 Sylvain Pelissier
-## Copyright (C) 2017 Michele Ginesi
+## Copyright (C) 2018 Michele Ginesi
 ##
 ## This file is part of Octave.
 ##
@@ -107,7 +107,7 @@
   x = x(:);  # convert to column vector
 
   ## Initialize the result
-  if (isreal (x) && x >= 0)
+  if ((isreal (x)) && (x >= 0))
     E1 = zeros (size (x), class (x));
   else
     E1 = complex (zeros (size (x), class (x)));
@@ -115,11 +115,13 @@
   tol = eps (class (x));
 
   ## Divide the input into 3 regions and apply a different algorithm for each.
-  s_idx = (((real (x) + 19.5).^2 ./ (20.5^2) + imag (x).^2 ./ (10^2)) <= 1) ...
-          | (real (x) < 0 & abs (imag (x)) <= 1e-8);
-  cf_idx = ((((real (x) + 1).^2 ./ (38^2) + imag (x).^2 ./ (40^2)) <= 1) ...
+  s_idx = (((real (x) + 19.5) .^ 2 ./ (20.5 ^ 2) + ...
+          imag (x) .^ 2 ./ (10 ^ 2)) <= 1) ...
+          | (real (x) < 0 & abs (imag (x)) <= 1e-08);
+  cf_idx = ((((real (x) + 1) .^ 2 ./ (38 ^ 2) + ...
+            imag (x) .^ 2 ./ (40 ^ 2)) <= 1) ...
             & (! s_idx)) & (real (x) <= 35);
-  a_idx = ! s_idx & ! cf_idx;
+  a_idx = ((! s_idx) & (! cf_idx));
   x_s  = x(s_idx);
   x_cf = x(cf_idx);
   x_a  = x(a_idx);
@@ -133,22 +135,22 @@
   ## formula 5.1.11, p 229.
   ## FIXME: Why so long?  IEEE double doesn't have this much precision.
   gm = 0.577215664901532860606512090082402431042159335;
-  e1_s = -gm - log (x_s);
-  res = -x_s;
+  e1_s = - gm - log (x_s);
+  res = - x_s;
   ssum = res;
   k = 1;
   fflag = true (size (res));
-  while (k < 1e3 && any (fflag))
+  while ((k < 1e03) && (any (fflag)))
     res_tmp = res(fflag);
     x_s_tmp = x_s(fflag);
     ssum_tmp = ssum(fflag);
-    res_tmp .*= (k * -x_s_tmp/((k + 1)^2));
+    res_tmp .*= (k * (- x_s_tmp) / ((k + 1) ^ 2));
     ssum_tmp += res_tmp;
     k += 1;
     res(fflag) = res_tmp;
     ssum(fflag) = ssum_tmp;
     x_s(fflag) = x_s_tmp;
-    fflag = abs (res) > tol*abs (ssum);
+    fflag = (abs (res) > (tol * abs (ssum)));
   endwhile
   e1_s -= ssum;
 
@@ -160,7 +162,8 @@
   e1_cf = exp(- x_cf);
   e1_cf .*= __expint_lentz__ (x_cf, strcmpi (class (x_cf), "single"));
 
-  # Remove spurious imaginary part if needed
+  # Remove spurious imaginary part if needed (__expint_lentz__ works
+  # automathically with complex values)
 
   if (isreal (x_cf) && x_cf >= 0)
     e1_cf = real (e1_cf);
@@ -177,7 +180,7 @@
     res_tmp = res(fflag);
     oldres_tmp = res_tmp;
     x_a_tmp = x_a(fflag);
-    res_tmp ./= (-x_a_tmp / (k+1));
+    res_tmp ./= (- x_a_tmp / (k + 1));
     ssum(fflag) += res_tmp;
     k += 1;
     res(fflag) = res_tmp;

--- a/scripts/specfun/gammainc.m	Fri Feb 23 20:30:28 2018 +0100
+++ b/scripts/specfun/gammainc.m	Sun Feb 25 00:02:44 2018 +0100
@@ -1,7 +1,7 @@
 ## Copyright (C) 2016 Marco Caliari
 ## Copyright (C) 2016 Nir Krakauer
-## Copyright (C) 2017 Michele Ginesi
 ## Copyright (C) 2018 Stefan Schlögl
+## Copyright (C) 2018 Michele Ginesi
 ##
 ## This file is part of Octave.
 ##
@@ -94,7 +94,7 @@
 
 function y = gammainc (x, a, tail = "lower")
 
-  if (nargin >= 4 || nargin <= 1)
+  if ((nargin >= 4) || (nargin <= 1))
     print_usage ();
   endif
 
@@ -105,7 +105,7 @@
     endif
   endif
 
-  if (any (a < 0) || any (imag (a) != 0))
+  if ((any (a < 0)) || (any (imag (a) != 0)))
     error ("gammainc: a must be real and non negative");
   endif
 
@@ -118,7 +118,7 @@
     x = double (x);
   endif
 
-  if (strcmpi (class (a), "single") || strcmpi (class (x), "single"))
+  if ((strcmpi (class (a), "single")) || (strcmpi (class (x), "single")))
     x = single (x);
     a = single (a);
   endif
@@ -194,7 +194,7 @@
   flag_a_small = ((abs (a) < 2) & (abs(a) > 0) & (! i_done) & (x < 0));
   a(flag_a_small) += 2;
 
-  flag_s = (((x + 0.25 < a | x < 0 | a < 5) & (x > -20)) | (abs (x) < 1));
+  flag_s = ((((x + 0.25 < a) | (x < 0) | (a < 5)) & (x > -20)) | (abs (x) < 1));
 
   ## Case 8: x, a relatively small.
   ii = ((! i_done) & flag_s);
@@ -234,7 +234,7 @@
 
 ## x == 0, a == 0.
 function y = gammainc_00 (tail)
-  if (strcmpi (tail, "upper") || strcmpi (tail, "scaledupper"))
+  if ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledupper")))
     y = 0;
   else
     y = 1;
@@ -243,7 +243,7 @@
 
 ## x == 0.
 function y = gammainc_x0 (tail)
-  if (strcmpi (tail, "upper") || strcmpi (tail, "scaledlower"))
+  if ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledlower")))
     y = 1;
   elseif (strcmpi (tail, "lower"))
     y = 0;
@@ -256,7 +256,7 @@
 function y = gammainc_x_inf (tail)
   if (strcmpi (tail, "lower"))
     y = 1;
-  elseif (strcmpi (tail, "upper") || strcmpi (tail, "scaledupper"))
+  elseif ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledupper")))
     y = 0;
   else
     y = Inf;
@@ -267,7 +267,7 @@
 function y = gammainc_a_inf (tail)
   if (strcmpi (tail, "lower"))
     y = 0;
-  elseif (strcmpi (tail, "upper") || strcmpi (tail, "scaledlower"))
+  elseif ((strcmpi (tail, "upper")) || (strcmpi (tail, "scaledlower")))
     y = 1;
   else
     y = Inf;
@@ -329,7 +329,7 @@
 ## Numerical Recipes in Fortran 77 (6.2.5)
 ## series
 function y = gammainc_s (x, a, tail)
-  if (strcmpi (tail, "scaledlower") || strcmpi (tail, "scaledupper"))
+  if ((strcmpi (tail, "scaledlower")) || (strcmpi (tail, "scaledupper")))
     y = ones (size (x), class (x));
     term = x ./ (a + 1);
   else
@@ -339,7 +339,7 @@
     term = y .* x ./ (a + 1);
   endif
   n = 1;
-  while (any (abs (term(:)) > abs (y(:)) * eps))
+  while ((any (abs (term(:))) > (abs (y(:)) * eps)))
     ## y can be zero from the beginning (gammainc (1,1000))
     jj = abs (term) > abs (y) * eps;
     n += 1;

--- a/scripts/specfun/sinint.m	Fri Feb 23 20:30:28 2018 +0100
+++ b/scripts/specfun/sinint.m	Sun Feb 25 00:02:44 2018 +0100
@@ -84,7 +84,7 @@
 
   ## For values small in modulus we use the series expansion
 
-  i_miss = ((i_miss) & (!flag_large));
+  i_miss = ((i_miss) & (! flag_large));
   xx = x(i_miss);
   ssum = xx; # First term of the series expansion
   yy = ssum;