Merge pull request #1243 from cbm755/bernoulli

Bernoulli B_1 moving to 1/2

Author: cbm755

inst/@double/bernoulli.m

@@ -1,4 +1,4 @@
-%% Copyright (C) 2018-2019, 2022 Colin B. Macdonald
+%% Copyright (C) 2018-2019, 2022-2023 Colin B. Macdonald
 %%
 %% This file is part of OctSymPy.
 %%
@@ -32,6 +32,19 @@
 %% @end group
 %% @end example
 %%
+%% Note there are two conventions in the literature about the sign of B_1,
+%% but for certain the absolute value is one half:
+%% @example
+%% @group
+%% @c Would render with an extra zero on earlier Octave
+%% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
+%% abs (bernoulli (1))
+%%   @result{} 0.5000
+%% @end group
+%% @end example
+%% As of 2023, this numerical evaluation function is in a state of flux
+%% about which one it takes, @pxref{@@sym/bernoulli}.
+%%
 %% Polynomial example:
 %% @example
 %% @group
@@ -85,7 +98,14 @@
 
 %!assert (bernoulli (0), 1)
 %!assert (bernoulli (3), 0)
-%!assert (bernoulli (1), -0.5, -eps)
+
+%!test
+%! % two different definitions in literature
+%! assert (abs (bernoulli (1)), 0.5, -eps)
+
+%!xtest
+%! % we want to use B_1 = 1/2, possible with a version-specific filter
+%! assert (bernoulli (1), 0.5, -eps)
 
 %!test
 %! n = sym(88);
@@ -95,7 +115,7 @@
 %! assert (A, B, -eps);
 
 %!test
-%! m = [0 1; 2 4];
+%! m = [0 2; 3 4];
 %! n = sym(m);
 %! A = bernoulli (m);
 %! B = double (bernoulli (n));

inst/@sym/bernoulli.m

@@ -1,4 +1,4 @@
-%% Copyright (C) 2014-2016, 2018-2019, 2022 Colin B. Macdonald
+%% Copyright (C) 2014-2016, 2018-2019, 2022-2023 Colin B. Macdonald
 %%
 %% This file is part of OctSymPy.
 %%
@@ -22,12 +22,33 @@
 %% @deftypemethodx @@sym {@var{p} =} bernoulli (@var{n}, @var{x})
 %% Return symbolic Bernoulli numbers or Bernoulli polynomials.
 %%
-%% Examples:
+%% With a sufficiently recent SymPy version, the first seven
+%% Bernoulli numbers are:
 %% @example
 %% @group
-%% bernoulli(sym(6))
+%% @c doctest: +XFAIL_IF(pycall_sympy__ ('return Version(spver) < Version("1.12.dev")'))
+%% bernoulli (sym(0:6))
+%%   @result{} (sym) [1  1/2  1/6  0  -1/30  0  1/42]  (1×7 matrix)
+%% @end group
+%% @end example
+%%
+%% Note there are two different definitions in use which differ
+%% in the sign of the value of B_1.  As of 2023 and a sufficiently
+%% recent SymPy library, we use the definition with positive one half:
+%% @example
+%% @group
+%% @c doctest: +XFAIL_IF(pycall_sympy__ ('return Version(spver) < Version("1.12.dev")'))
+%% bernoulli (sym(1))
+%%   @result{} (sym) 1/2
+%% @end group
+%% @end example
+%%
+%% Other examples:
+%% @example
+%% @group
+%% bernoulli (sym(6))
 %%   @result{} (sym) 1/42
-%% bernoulli(sym(7))
+%% bernoulli (sym(7))
 %%   @result{} (sym) 0
 %% @end group
 %% @end example
@@ -36,7 +57,7 @@
 %% @example
 %% @group
 %% syms x
-%% bernoulli(2, x)
+%% bernoulli (2, x)
 %%   @result{} (sym)
 %%        2       1
 %%       x  - x + ─
@@ -69,9 +90,19 @@
 %! assert (isequal (bernoulli(3,x), x^3 - 3*x^2/2 + x/2))
 
 %!test
-%! m = sym([0 1; 8 888889]);
+%! % two different definitions in literature
+%! assert (isequal (abs (bernoulli (sym(1))), sym(1)/2))
+
+%!test
+%! % we use B_1 = 1/2
+%! if (pycall_sympy__ ('return Version(spver) >= Version("1.12.dev")'))
+%!   assert (isequal (bernoulli (sym(1)), sym(1)/2))
+%! end
+
+%!test
+%! m = sym([0 2; 8 888889]);
 %! A = bernoulli (m);
-%! B = [1 -sym(1)/2; -sym(1)/30 0];
+%! B = [1 sym(1)/6; -sym(1)/30 0];
 %! assert (isequal (A, B))
 
 %!test

inst/heaviside.m

@@ -1,5 +1,5 @@
 %% Copyright (C) 2006 Sylvain Pelissier <[email protected]>
-%% Copyright (C) 2015-2016 Colin B. Macdonald <[email protected]>
+%% Copyright (C) 2015-2016, 2019, 2022 Colin B. Macdonald <[email protected]>
 %%
 %% 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
@@ -35,6 +35,7 @@
 %% function returns 0.5 by default:
 %% @example
 %% @group
+%% @c Would render with an extra zero on earlier Octave
 %% @c doctest: +SKIP_IF(compare_versions (OCTAVE_VERSION(), '6.0.0', '<'))
 %% heaviside(0)
 %%   @result{} 0.5000