# HG changeset patch # User Michele Ginesi # Date 1519513364 -3600 # Node ID 40ab8be7d7ec78d9984fa74ec6274dbea48b2388 # Parent aece120e1ec108683ea56db226befd9c0f14a10c Fixed style in specfun scripts -- changed scripts/specfun/betainc.m changed scripts/specfun/cosint.m changed scripts/specfun/expint.m changed scripts/specfun/gammainc.m changed scripts/specfun/sinint.m diff -r aece120e1ec1 -r 40ab8be7d7ec scripts/specfun/betainc.m --- 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 diff -r aece120e1ec1 -r 40ab8be7d7ec scripts/specfun/cosint.m --- 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. diff -r aece120e1ec1 -r 40ab8be7d7ec scripts/specfun/expint.m --- 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; diff -r aece120e1ec1 -r 40ab8be7d7ec scripts/specfun/gammainc.m --- 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; diff -r aece120e1ec1 -r 40ab8be7d7ec scripts/specfun/sinint.m --- 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;