@vv

Interesting. So in fact this is a bug in simplify, not in sum:

simplify(polylog(1/2, -1)+1);
Error, (in polylog) numeric exception: division by zero

@Robert Israel

At least for trigs of Pi*rational there is a method that usually works:

ee:=-1-2*cos((8/83)*Pi)-2*cos((10/83)*Pi)+2*cos((9/83)*Pi)+2*cos((5/83)*Pi)+2*cos((7/83)*Pi)-
    2*cos((2/83)*Pi)+2*cos((1/83)*Pi)-2*cos((38/83)*Pi)+2*cos((41/83)*Pi)-2*cos((40/83)*Pi)-
    2*cos((6/83)*Pi)-2*cos((4/83)*Pi)-2*cos((36/83)*Pi)-2*cos((32/83)*Pi)-2*cos((34/83)*Pi)-
    2*cos((30/83)*Pi)-2*cos((26/83)*Pi)-2*cos((28/83)*Pi)-2*cos((24/83)*Pi)+2*cos((31/83)*Pi)+
    2*cos((33/83)*Pi)+2*cos((29/83)*Pi)+2*cos((37/83)*Pi)+2*cos((35/83)*Pi)+2*cos((39/83)*Pi)-
    2*cos((22/83)*Pi)-2*cos((18/83)*Pi)-2*cos((20/83)*Pi)-2*cos((14/83)*Pi)-2*cos((16/83)*Pi)-
    2*cos((12/83)*Pi)+2*cos((27/83)*Pi)+2*cos((25/83)*Pi)+2*cos((23/83)*Pi)+2*cos((13/83)*Pi)+
    2*cos((15/83)*Pi)+2*cos((11/83)*Pi)+2*cos((3/83)*Pi)+2*cos((19/83)*Pi)+2*cos((21/83)*Pi)+
    2*cos((17/83)*Pi):

(evala@convert)(ee, RootOf);
                               0

with(LinearAlgebra):
M := Matrix(5, 5,  (i, j) -> (10*i+j)*sin((1/180)*Pi*(10*i+j))):

timed := proc(e::uneval) local f, t, ans := NULL;
  if procname::indexed and [op](procname) = [real] then f := time[real] else f := time end if;
  t := f(proc() try assign('ans', eval(e)) catch: WARNING(lastexception[2 .. -1]) end try end proc());
  [t, ans] end proc:

Normalizer := evala:

M2 := (evala@convert)~(M, RootOf):
M3 := subs((evala@Algfield)(M2)[1], M2): # M2 will do as well, but will be slower

timed(MatrixInverse(M3));
Warning, singular matrix
                            [2.562]

forget(evala);
timed(Determinant(M3));
                           [1.625, 0]

forget(evala);
timed([LUDecomposition](M3)[3][5]);
        [2.531, Vector[row](5, [0, 0, 0, 0, 0])]

Normalizer:=normal is in fact something like Normalizer:=automatic, invoking LinearAlgebra:-DeduceNormalizer, which has a custom setting for algebraics. But Determinant apparently doesn't call DeduceNormalizer, so Normalizer := evala is needed.

@vv

For the same reason why it's needed for evalf/Int. Like the example from point #3 here: it doesn't work with higher settings of digits, because the tolerance setting also increases and is not attainable. But setting the tolerance (epsilon) to a more reasonable value works like a charm.

For evalf/Sum we cannot do a direct test, but, as some indirect evidence, at Digits:=19, evalf/Sum of 1/(k^2+ln(k)) takes about 40 seconds; for Digits:=20, it takes 110 seconds. The difference seems to be due mostly to the increase in the requested tolerance. But the point is that often one does not need the extra accuracy, only the extra digits to work with numerically unstable summands.

If Maple had a summation method based on the Euler–Maclaurin formula, the evalf/Int example would be directly relevant. Perhaps the issue is that Maple seems to have exactly one summation method. So it fails even for rather simple summands, like 1/(k^2*ln(k)) or sin(Pi*k/3)/k.

@Axel Vogt

The question is about computing a sum numerically when no closed form is available.

@vv

That's debatable. Take the sum of (1-10^(-20))^k/k and consider how many terms you need to add to get one correct digit. While for an appropriate summation method the example is easy (just evaluating the integral already gets you in the right ballpark).

Your example is just another illustration that evalf/Sum needs more user-controllable options like epsilon. Specifically, for your example, an option to control how many terms to take at the beginning.

Also:

foldl(seq, i, i = 1 .. 2);
                              i, i

Unlike in the example with apply, here i doesn't get properly localized.

ee := foldl(%seq, 0, i = 1 .. 2, j = 1 .. 2)

value(ee) # fails

eval(ee, %seq = seq) # fails

eval(subs(%seq = seq, ee)) # works

This one is less of an issue (it's safer to use [%seq]), but it looks like an accident of the implementation that the last one works and the other two don't.

@pchin

That works, thanks. ?plot/tickmarks could have been clearer, I guess, because it says "Restricted forms of the tickmarks and gridlines option may also be given directly to a plotting command" and then "The tickmark specification for all axes may be given through a single tickmarks=[s1, s2]", giving an impression that tickmarks=[xspec, yspec] should work as a suboption to axis as well, since it works as a top-level option.

@acer

Might be system-dependent? Because on Windows I'm seeing this: https://i.imgur.com/4b6pcoP.jpg. The details are definitely missing.

@acer

Impressive. One possible catch is that 'background' isn't pixel-perfect. In Embed(abs(fft2unmasked)^(1/4)) there's a prominent vertical midline, and nine dots visible (the ninth one is in the center), five of them really as tiny dots rather than crosses. In Preview2D, there's a much dimmer vertical midline, and only the four crosses visible. So I'd say this should be used in conjunction with Embed.

@Mariusz Iwaniuk

Well, there's nothing wrong with the antiderivative, it's even continuous. The problem seems to occur when int computes the limit of the antiderivative at infinity.

@vv

Indeed. So for S1 it's the other way around, :-asympt fails, but MultiSeries:-asympt knows how to deal with those elliptics:

MultiSeries:-asympt(subs(S(a, t) = int(s(tau), tau = a .. t), E(a, t)[2]), t);
    (1/108-1/108*I)*((9*I)*sqrt(9)*sqrt(6)*a^5+9*sqrt(9)*sqrt(6)*a^5+(8*I)*sqrt(3*I)*sqrt(-3*I)*
    EllipticK(I)*sqrt(9*a^4+1)+4*sqrt(3*I)*sqrt(-3*I)*EllipticK(I*sqrt(-2))*sqrt(9*a^4+1)+I*sqrt(9)*
    sqrt(6)*a+4*sqrt(9)*sqrt(1-(3*I)*a^2)*sqrt(1+(3*I)*a^2)*EllipticF((1/2+1/2*I)*a*sqrt(6), I)+
    sqrt(9)*sqrt(6)*a)*sqrt(6)/sqrt(9*a^4+1)+O(1/t)

evalf(subs(a = 1, %));
    2.261499331+0.*I+O(1/t)

But this is the position of the asymptote at -infinity. So as the asymptotic expansion of E(a, t) at +infinity, this is wrong.

@vv

I'm probably missing something here. This works:

:-asympt~(subs(S(a, t) = S__h, E(a, t)), t)

but not MultiSeries:-asympt. What are you applying MultiSeries:-asympt to?

@Carl Love

Then, if an expression with the head `seq` can be constructed using eval and suchlike, why not using composition?

apply works:

(f@curry(apply, seq))(i, i = 1 .. 2);
                            f(1, 2)

And composition works for time:

(f@time[real])(Threads:-Sleep(1));
                             f(0.)

It's like you say, the argument is evaluated before the special evaluation rules kick in.

(To digress, it also means that (`[]`@time)(e) is not the same as [time](e).)

In all fairness, there is this page: ?spec_eval_rules, and it does list evalf. But it only says about evalf that "it ensures that its second argument (if provided) is evaluated before its first argument".

Adding two other examples of the code crashing or not working correctly when several commands are put in the same cell. Both examples work fine if each input line is in a separate 2D Input cell.

Download 2d2.mw