Mercurial > octave
view scripts/specfun/cosint.m @ 25054:6652d3823428 stable
maint: Update copyright dates in all source files.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 30 Mar 2018 09:19:05 -0400 |
parents | f80e68529bbf |
children | 6dd232798997 |
line wrap: on
line source
## Copyright (C) 2017-2018 Michele Ginesi ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## Author: Michele Ginesi <michele.ginesi@gmail.com> ## -*- texinfo -*- ## @deftypefn {} {} cosint (@var{x}) ## Compute the cosine integral function: ## @tex ## $$ ## {\rm Ci} (x) = - \int_x^\infty {{\cos (t)} \over t} dt ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## +oo ## / ## Ci (x) = - | (cos (t)) / t dt ## / ## x ## @end group ## @end example ## ## @end ifnottex ## An equivalent definition is ## @tex ## $$ ## {\rm Ci} (x) = \gamma + \log (x) + \int_0^x {{\cos (t) - 1} \over t} dt ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## x ## / ## | cos (t) - 1 ## Ci (x) = gamma + log (x) + | ------------- dt ## | t ## / ## 0 ## @end group ## @end example ## ## @end ifnottex ## Reference: ## ## @nospell{M. Abramowitz and I.A. Stegun}, ## @cite{Handbook of Mathematical Functions} ## 1964. ## ## @seealso{sinint, expint, cos} ## ## @end deftypefn function y = cosint (x) if (nargin != 1) print_usage (); endif if (! isnumeric (x)) error ("cosint: X must be numeric"); endif ## Convert to floating point if necessary if (isinteger (x)) x = double (x); endif ## Convert to column vector orig_sz = size (x); if (iscomplex (x)) ## Work around reshape which narrows to real (bug #52953) x = complex (real (x)(:), imag (x)(:)); else x = x(:); end ## Initialize the result y = zeros (size (x), class (x)); tol = eps (class (x)); todo = true (size (x)); ## Special values y(x == Inf) = 0; y((x == -Inf) & !signbit (imag (x))) = 1i * pi; y((x == -Inf) & signbit (imag (x))) = -1i * pi; todo(isinf (x)) = false; ## For values large in modulus, but not in the range (-oo,0), we use the ## relation with expint. flag_large = (abs (x) > 2); xx = x(flag_large); ## Abramowitz, relation 5.2.20 ii_sw = (real (xx) <= 0 & imag (xx) < 0); xx(ii_sw) = conj (xx(ii_sw)); ii_nw = (real (xx) < 0); xx(ii_nw) *= -1; yy = -0.5 * (expint (1i * xx) + expint (-1i * xx)); yy(ii_nw) += 1i * pi; yy(ii_sw) = conj (yy(ii_sw)); y(todo & flag_large) = yy; todo(flag_large) = false; ## For values small in modulus, use the series expansion (also near (-oo, 0]) 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. xx = complex (real (x)(todo), imag (x)(todo)); else xx = x(todo); end ssum = - xx .^ 2 / 4; # First term of the series expansion ## FIXME: This is way more precision than a double value can hold. gma = 0.57721566490153286060651209008; # Euler gamma constant yy = gma + log (complex (xx)) + ssum; # log(complex(...) handles signed zero flag_sum = true (nnz (todo), 1); it = 0; maxit = 300; while (any (flag_sum) && (++it < maxit)); ssum .*= - xx .^ 2 * (2 * it) / ((2 * it + 2) ^ 2 * (2 * it + 1)); yy(flag_sum) += ssum (flag_sum); flag_sum = (abs (ssum) >= tol); endwhile y(todo) = yy; ## Clean up values which are purely real flag_neg_zero_imag = (real (x) < 0) & (imag (x) == 0) & signbit (imag (x)); y(flag_neg_zero_imag) = complex (real (y(flag_neg_zero_imag)), -pi); ## Restore original shape y = reshape (y, orig_sz); endfunction %!assert (cosint (1.1), 0.38487337742465081550, 2 * eps); %!test %! x = [2, 3, pi; exp(1), 5, 6]; %! A = cosint (x); %! B = [0.422980828774864996, 0.119629786008000328, 0.0736679120464254860; ... %! 0.213958001340379779, -0.190029749656643879, -0.0680572438932471262]; %! assert (A, B, -5e-15); %!assert (cosint (0), - Inf) %!assert (cosint (-0), -inf + 1i*pi) %!assert (cosint (complex (-0, 0)), -inf + 1i*pi) %!assert (cosint (complex (-0, -0)), -inf - 1i*pi) %!assert (cosint (inf), 0) %!assert (cosint (-inf), 1i * pi) %!assert (cosint (complex (-inf, -0)), -1i * pi) %!assert (isnan (cosint (nan))) %!assert (class (cosint (single (1))), "single") ## tests against maple %!assert (cosint (1), 0.337403922900968135, -2*eps) %!assert (cosint (-1), 0.337403922900968135 + 3.14159265358979324*I, -2*eps) %!assert (cosint (pi), 0.0736679120464254860, -2e-15) %!assert (cosint (-pi), 0.0736679120464254860 + 3.14159265358979324*I, -2*eps) %!assert (cosint (300), -0.00333219991859211178, -2*eps) %!assert (cosint (1e4), -0.0000305519167244852127, -2*eps) %!assert (cosint (20i), 1.28078263320282944e7 + 1.57079632679489662*I, -2*eps) %!test %! x = (0:4).'; %! y_ex = [-Inf %! 0.337403922900968135 %! 0.422980828774864996 %! 0.119629786008000328 %! -0.140981697886930412]; %! assert (cosint (x), y_ex, -3e-15); %!test %! x = -(1:4).'; %! y_ex = [0.337403922900968135 + pi*1i %! 0.422980828774864996 + pi*1i %! 0.119629786008000328 + pi*1i %! -0.140981697886930412 + pi*1i]; %! assert (cosint (x), y_ex, -4*eps); %!test %! x = complex (-(1:4).', 0); %! y_ex = [0.337403922900968135 + pi*1i %! 0.422980828774864996 + pi*1i %! 0.119629786008000328 + pi*1i %! -0.140981697886930412 + pi*1i]; %! assert (cosint (x), y_ex, -4*eps); %!test %! x = complex (-(1:4).', -0); %! y_ex = [0.337403922900968135 - pi*1i %! 0.422980828774864996 - pi*1i %! 0.119629786008000328 - pi*1i %! -0.140981697886930412 - pi*1i]; %! assert (cosint (x), y_ex, -4*eps); %!test %! x = 1i * (0:4).'; %! y_ex = [-Inf %! 0.837866940980208241 + 1.57079632679489662*I %! 2.45266692264691452 + 1.57079632679489662*I %! 4.96039209476560976 + 1.57079632679489662*I %! 9.81354755882318556 + 1.57079632679489662*I]; %! assert (cosint (x), y_ex, -4*eps); %!test %! x = -1i * (1:4).'; %! y_ex = [0.837866940980208241 - 1.57079632679489662*I %! 2.45266692264691452 - 1.57079632679489662*I %! 4.96039209476560976 - 1.57079632679489662*I %! 9.81354755882318556 - 1.57079632679489662*I]; %! assert (cosint (x), y_ex, -4*eps); %!test %! x = [1+2i; -2+5i; 2-5i; 100; 10i; -1e-4 + 1e-6*1i; -20-1i]; %! A = [ 2.03029639329172164 - 0.151907155175856884*I %! 1.61538963829107749 + 19.7257540553382650*I %! 1.61538963829107749 + 16.5841614017484717*I %! -0.00514882514261049214 %! 1246.11448604245441 + 1.57079632679489662*I %! -8.63307471207423322 + 3.13159298695312800*I %! 0.0698222284673061493 - 3.11847446254772946*I ]; %! B = cosint (x); %! assert (A, B, -3*eps); %! B = cosint (single (x)); %! assert (A, B, -3*eps ("single")); ## Fails along negative real axis %!test %! x = [-25; -100; -1000]; %! yex = [-0.0068485971797025909189 + pi*1i %! -0.0051488251426104921444 + pi*1i %! 0.000826315511090682282 + pi*1i]; %! y = cosint (x); %! assert (y, yex, -5*eps); ## FIXME: Need a test for bug #52953 %#!test <*52953> ## Test input validation %!error cosint () %!error cosint (1,2) %!error <X must be numeric> cosint ("1")