Mercurial > octave

changeset 24913:9da779d2f029 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
cosint improve signed zero imag input near branch cut Stop special casing the origin, instead have log deal with signed zeros (real and complex). Ensure we approach the branch cut (negative real axis) correctly based on signed 0 for imaginary part. cosint.m: handle complex signed zero, add and modify BIST.
author Colin Macdonald <cbm@m.fsf.org>
date Tue, 23 Jan 2018 23:14:50 -0800
parents 08c871c4147b
children d7293106945c
files scripts/specfun/cosint.m
diffstat 1 files changed, 38 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/specfun/cosint.m	Tue Jan 23 22:43:03 2018 -0800
+++ b/scripts/specfun/cosint.m	Tue Jan 23 23:14:50 2018 -0800
@@ -77,18 +77,25 @@
   endif
 
   sz = size (x);
-  x = x(:);
+  #x = reshape (x, numel (x), 1);
+  if (iscomplex (x))
+    ## workaround reshape narrowing to real (#52953)
+    x = complex (real (x)(:), imag (x)(:));
+  else
+    x = x(:);
+  end
+
   y = zeros (size (x), class (x));
   tol = eps (class (x));
 
   i_miss = true (length (x), 1);
 
   ## special values
-  y(x == 0) = - Inf;
   y(x == Inf) = 0;
-  y(x == - Inf) = 1i * pi;
+  y((x == -Inf) & !signbit (imag(x))) = 1i * pi;
+  y((x == -Inf) &  signbit (imag(x))) = -1i * pi;
 
-  i_miss = ((i_miss) & (x != 0) & (x != Inf) & (x != - Inf));
+  i_miss = ((i_miss) & (x != Inf) & (x != -Inf));
 
   ## For values large in modulus and not in (-oo,0), we use the relation
   ## with expint
@@ -108,10 +115,16 @@
 
   ## For values small in modulus, use the series expansion (also near (-oo, 0])
   i_miss = ((i_miss) & (!flag_large));
-  xx = x(i_miss);
+  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)(i_miss), imag (x)(i_miss));
+  else
+    xx = x(i_miss);
+  end
   ssum = - xx .^ 2 / 4; # First term of the series expansion
   gma = 0.57721566490153286060651209008; # Euler gamma constant
-  yy = gma + log (xx) + ssum;
+  yy = gma + log (complex (xx)) + ssum; # log(complex(...) handles signed zero
   flag_sum = true (nnz (i_miss), 1);
   it = 1;
   maxit = 300;
@@ -139,8 +152,12 @@
 %! 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")
@@ -164,12 +181,19 @@
 %! assert (cosint(x), y_ex, -3e-15);
 
 %!test
-%! x = -(0:4) ';
-%! y_ex = [-Inf
-%!         0.337403922900968135 + 3.14159265358979324*I
-%!         0.422980828774864996 + 3.14159265358979324*I
-%!         0.119629786008000328 + 3.14159265358979324*I
-%!         -0.140981697886930412 + 3.14159265358979324*I];
+%! 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
@@ -182,9 +206,8 @@
 %! assert (cosint (x), y_ex, -4*eps);
 
 %!test
-%! x = - 1i * (0:4).';
-%! y_ex = [- Inf
-%!         0.837866940980208241 - 1.57079632679489662*I
+%! 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];