From f1e979f880e7d5892e1dc9ba1a7d1e863b40730d Mon Sep 17 00:00:00 2001 From: "Colin B. Macdonald" Date: Thu, 21 Apr 2016 17:25:06 -0700 Subject: [PATCH 1/5] Add sinint, copied from specfun pkg Original author seems to be Sylvain Pelissier, already on our contributors list. Others also made changes---as always, please let us know if you feel you should be on the contributors list for Symbolic! --- inst/sinint.m | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 inst/sinint.m diff --git a/inst/sinint.m b/inst/sinint.m new file mode 100644 index 000000000..f1f3e0a82 --- /dev/null +++ b/inst/sinint.m @@ -0,0 +1,45 @@ +## Copyright (C) 2006 Sylvain Pelissier +## +## This program 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. +## +## This program 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 +## this program; if not, see . + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{y} =} sinint (@var{x}) +## Compute the sine integral defined by: +## @verbatim +## x +## / +## sinint(x) = | sin(t)/t dt +## / +## 0 +## @end verbatim +## @seealso{cosint, expint} +## @end deftypefn + +function y = sinint (x) + if (nargin != 1) + print_usage; + endif + y = zeros(size(x)); + if prod(size(x)) < 101 + for k = 1:prod(size(x)) + y(k) = sum(besselj([0:100]+0.5,(x(k)/2)).^2); + endfor + y = y.*pi; + else + for k=0:100 + y += besselj(k+0.5,x/2).^2; + endfor + y = y.*pi; + endif +endfunction From 0c18a44ea9b3d85f5bee037d99825d335a8280dd Mon Sep 17 00:00:00 2001 From: "Colin B. Macdonald" Date: Thu, 21 Apr 2016 21:43:54 -0700 Subject: [PATCH 2/5] adopt matlab-compatible code --- inst/sinint.m | 68 +++++++++++++++++++++++++-------------------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/inst/sinint.m b/inst/sinint.m index f1f3e0a82..3757632b3 100644 --- a/inst/sinint.m +++ b/inst/sinint.m @@ -1,45 +1,45 @@ -## Copyright (C) 2006 Sylvain Pelissier -## -## This program 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. -## -## This program 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 -## this program; if not, see . +%% Copyright (C) 2006 Sylvain Pelissier +%% +%% This program 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. +%% +%% This program 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 +%% this program; if not, see . -## -*- texinfo -*- -## @deftypefn {Function File} {@var{y} =} sinint (@var{x}) -## Compute the sine integral defined by: -## @verbatim -## x -## / -## sinint(x) = | sin(t)/t dt -## / -## 0 -## @end verbatim -## @seealso{cosint, expint} -## @end deftypefn +%% -*- texinfo -*- +%% @deftypefn {Function File} {@var{y} =} sinint (@var{x}) +%% Compute the sine integral defined by: +%% @verbatim +%% x +%% / +%% sinint(x) = | sin(t)/t dt +%% / +%% 0 +%% @end verbatim +%% @seealso{cosint, expint} +%% @end deftypefn function y = sinint (x) - if (nargin != 1) - print_usage; - endif + if (nargin ~= 1) + print_usage (); + end y = zeros(size(x)); if prod(size(x)) < 101 for k = 1:prod(size(x)) y(k) = sum(besselj([0:100]+0.5,(x(k)/2)).^2); - endfor + end y = y.*pi; else for k=0:100 - y += besselj(k+0.5,x/2).^2; - endfor + y = y + besselj(k+0.5,x/2).^2; + end y = y.*pi; - endif -endfunction + end +end From 8fce3f18153dd4713d75dfc14336bee66bfe05af Mon Sep 17 00:00:00 2001 From: "Colin B. Macdonald" Date: Mon, 25 Apr 2016 13:51:21 -0700 Subject: [PATCH 3/5] update docs --- inst/sinint.m | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/inst/sinint.m b/inst/sinint.m index 3757632b3..5eb7c1341 100644 --- a/inst/sinint.m +++ b/inst/sinint.m @@ -14,8 +14,11 @@ %% this program; if not, see . %% -*- texinfo -*- -%% @deftypefn {Function File} {@var{y} =} sinint (@var{x}) -%% Compute the sine integral defined by: +%% @documentencoding UTF-8 +%% @defun sinint (@var{x}) +%% Compute the sine integral function. +%% +%% The sine integral function is defined for complex inputs @var{x} by: %% @verbatim %% x %% / @@ -23,8 +26,16 @@ %% / %% 0 %% @end verbatim -%% @seealso{cosint, expint} -%% @end deftypefn +%% +%% Example: +%% @example +%% @group +%% sinint (1) +%% @result{} 0.94608 +%% @end group +%% @end example +%% @seealso{cosint, expint, @@sym/sinint} +%% @end defun function y = sinint (x) if (nargin ~= 1) From 0bdc4232adae3ef6b5f7ce129404a2bbe3f4f861 Mon Sep 17 00:00:00 2001 From: "Colin B. Macdonald" Date: Mon, 25 Apr 2016 13:51:32 -0700 Subject: [PATCH 4/5] Add BIST, many of which fail --- inst/sinint.m | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/inst/sinint.m b/inst/sinint.m index 5eb7c1341..c9335de34 100644 --- a/inst/sinint.m +++ b/inst/sinint.m @@ -54,3 +54,51 @@ y = y.*pi; end end + + +%!assert (isequal (sinint (0), 0)) + +%!xtest +%! % FIXME: could relax to within eps +%! assert (isequal (sinint (inf), pi/2)) +%! assert (isequal (sinint (-inf), -pi/2)) + +% FIXME: these should be tighter +%!assert (sinint (1), 0.9460830703671830149414, 10*eps) +%!assert (sinint (-1), -0.9460830703671830149414, 10*eps) +%!assert (sinint (pi), 1.851937051982466170361, 20*eps) +%!assert (sinint (-pi), -1.851937051982466170361, 20*eps) + +%!xtest +%! % FIXME: when these pass, split these out into individual asserts +%! assert (sinint (300), 1.570881088213749519252, 2*eps) +%! assert (sinint (1e4), 1.570233121968771218148, 2*eps) +%! assert (sinint (20i), 1.280782633202829445942e8*1i, -5*eps) + +%!test +%! % random complex points in disc radius 4 +%! a = 0; b = 4; N = 30; +%! t = floor(9*rand(1,N))+1; +%! rs = (b-a) * sym(t)/9 + a; +%! r = (b-a) * t/9 + a; +%! ths = linspace(sym(0), 2*sym(pi), N); +%! th = linspace(0, 2*pi, N); +%! zs = rs.*exp(1i*ths); +%! z = r.*exp(1i*th); +%! A = sinint(z); +%! B = double(sinint(zs)); +%! assert (A, B, -20*eps) % should be tighter + +%!xtest +%! % random complex points in annulus [4, 25] +%! a = 5; b = 25; N = 30; +%! t = floor(9*rand(1,N))+1; +%! rs = (b-a) * sym(t)/9 + a; +%! r = (b-a) * t/9 + a; +%! ths = linspace(sym(0), 2*sym(pi), N); +%! th = linspace(0, 2*pi, N); +%! zs = rs.*exp(1i*ths); +%! z = r.*exp(1i*th); +%! A = sinint(z); +%! B = double(sinint(zs)); +%! assert (A, B, -20*eps) From 9ae596358c64ed3d2081c1693f516511bfe7235e Mon Sep 17 00:00:00 2001 From: "Colin B. Macdonald" Date: Mon, 25 Apr 2016 15:42:54 -0700 Subject: [PATCH 5/5] Add msg to docs about poor accuracy --- inst/sinint.m | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/inst/sinint.m b/inst/sinint.m index c9335de34..745ace45f 100644 --- a/inst/sinint.m +++ b/inst/sinint.m @@ -34,6 +34,16 @@ %% @result{} 0.94608 %% @end group %% @end example +%% +%% @strong{WARNING} This implementation is rather poor. +%% @itemize +%% @item It has very large errors for @code{abs(x) > 100}. +%% @item Even for small @var{x}, the relative error is often +%% as large as @code{500*eps} (accurate to only about +%% 13 decimals). +%% @end itemize +%% If you can, please help improve Octave by contributing to an +%% improved implementation! %% @seealso{cosint, expint, @@sym/sinint} %% @end defun @@ -89,6 +99,14 @@ %! B = double(sinint(zs)); %! assert (A, B, -20*eps) % should be tighter +%!test +%! % random real points in [0 100] +%! xs = linspace(0, sym(100), 20); +%! x = linspace(0, 100, 20); +%! A = sinint(x); +%! B = double(sinint(xs)); +%! assert (A, B, 500*eps) % should be tighter + %!xtest %! % random complex points in annulus [4, 25] %! a = 5; b = 25; N = 30;