Agree, with 2 additions.

First, the "recent" pages are not much useful not only because of spam, but also because they don't reflect comments to the answers.

Second, this site hasn't been much usable starting from its "upgrade" (which actually was a large "downgrade") to its 2nd version (about a year ago).

The only (but big) improvement that I see in this 2nd version is the ability to edit our posts even after somebody replied to them. The problem is that you won't know that somebody replied until you check the thread, and if your post in that thread was not an answer, but a comment, it's not that easy to find that thread if it's more than one or two days old.

It may be a little bit easier for people receiving email notifications. But who uses email nowadays? I don't. It's not as it used to be 15 years ago.

This site used to be unique by being by far the best compared to other competing systems; and now it is unique by having the most problems, also by far.

Alec

However, convert/exp produces the answer in a different form.

For the first term, that can be corrected by using the following conversion to myexp instead of conversion to exp,

`convert/myexp`:=proc(e)
local a,b,c,f;
f:=proc(d)
if typematch(d,`&*`(a::nonreal,specfunc(c::algebraic,exp))) then b:=1 
else typematch(d,`&*`(a::nonreal,b::algebraic,specfunc(c::algebraic,exp)))
fi;
abs(a)*b*exp(evalc(Re(c)+I*(Im(c)+argument(a))))
end;
evalindets(combine(convert(e,exp)),
{`&*`(nonreal,algebraic,specfunc(algebraic,exp)),
`&*`(nonreal,specfunc(algebraic,exp))},f)
end:

convert(MultiSeries:-asympt(HankelH1(t,s),s,1),compose,polynom,myexp) assuming t>=0;

             1/2 / 1  \1/2
            2    |----|    exp((-1/2 Pi t + s - 1/4 Pi) I)
                 \Pi s/

For the second and other terms, it would be convenient (after using something like collect(...,exp,factor)) to match the following type, `&*`(nonreal,seq(algebraic),specfunc(algebraic,exp)), but the seq(algebraic) type can be used only in procedures parameters.

Alec

Edit: 2nd attempt,

`convert/myexp`:=proc(e) local k;
evalindets(maptype(`+`,factor,combine(convert(e,exp))),
And(`*`,satisfies(x->op(1,x)::nonreal and membertype(exp(algebraic),[op(2..-1,x)],'k'))),
d->applyop(abs,1,applyop(c->evalc(Re(c)+I*(Im(c)+argument(op(1,d)))),[k+1,1],d)))
end:

It works not only for the first term, but for the second and other terms, too, as

convert(MultiSeries:-asympt(HankelH1(t,s),s,2),compose,polynom,myexp) assuming t>=0;

       1/2 / 1  \1/2
  1/8 2    |----|    exp((-1/2 Pi t + s + 1/4 Pi) I) (2 t - 1) (2 t + 1)/s
           \Pi s/

                       1/2 / 1  \1/2
                    + 2    |----|    exp((-1/2 Pi t + s - 1/4 Pi) I)
                           \Pi s/

but the order of terms still might be improved.

Alec

However, convert/exp produces the answer in a different form.

For the first term, that can be corrected by using the following conversion to myexp instead of conversion to exp,

`convert/myexp`:=proc(e)
local a,b,c,f;
f:=proc(d)
if typematch(d,`&*`(a::nonreal,specfunc(c::algebraic,exp))) then b:=1 
else typematch(d,`&*`(a::nonreal,b::algebraic,specfunc(c::algebraic,exp)))
fi;
abs(a)*b*exp(evalc(Re(c)+I*(Im(c)+argument(a))))
end;
evalindets(combine(convert(e,exp)),
{`&*`(nonreal,algebraic,specfunc(algebraic,exp)),
`&*`(nonreal,specfunc(algebraic,exp))},f)
end:

convert(MultiSeries:-asympt(HankelH1(t,s),s,1),compose,polynom,myexp) assuming t>=0;

             1/2 / 1  \1/2
            2    |----|    exp((-1/2 Pi t + s - 1/4 Pi) I)
                 \Pi s/

For the second and other terms, it would be convenient (after using something like collect(...,exp,factor)) to match the following type, `&*`(nonreal,seq(algebraic),specfunc(algebraic,exp)), but the seq(algebraic) type can be used only in procedures parameters.

Alec

Edit: 2nd attempt,

`convert/myexp`:=proc(e) local k;
evalindets(maptype(`+`,factor,combine(convert(e,exp))),
And(`*`,satisfies(x->op(1,x)::nonreal and membertype(exp(algebraic),[op(2..-1,x)],'k'))),
d->applyop(abs,1,applyop(c->evalc(Re(c)+I*(Im(c)+argument(op(1,d)))),[k+1,1],d)))
end:

It works not only for the first term, but for the second and other terms, too, as

convert(MultiSeries:-asympt(HankelH1(t,s),s,2),compose,polynom,myexp) assuming t>=0;

       1/2 / 1  \1/2
  1/8 2    |----|    exp((-1/2 Pi t + s + 1/4 Pi) I) (2 t - 1) (2 t + 1)/s
           \Pi s/

                       1/2 / 1  \1/2
                    + 2    |----|    exp((-1/2 Pi t + s - 1/4 Pi) I)
                           \Pi s/

but the order of terms still might be improved.

Alec

@Axel Vogt 

The answer with 23 digits is especially interesting,

Digits:=23:
frem(10^23,evalf(2*Pi));

                                 -10.

That answer can be also achieved without changing the Digits setting as

restart;
frem(10^10,evalf(2*Pi));

                                 -10.

Also,

frem(10^10-1,evalf(2*Pi));

                                 -11.

but

frem(10^10+1,evalf(2*Pi));

                                 -10.

2*Pi is not an exception for that. For example, here is what happens with e,

for i from 11 to 15 do
    frem(10^i,evalf(exp(1)))
od;
                                  10.
                                 100.
                                1000.
                               10000.
                              100000.

(Corollary: e ≥ 200,000).

Alec

@Axel Vogt 

The answer with 23 digits is especially interesting,

Digits:=23:
frem(10^23,evalf(2*Pi));

                                 -10.

That answer can be also achieved without changing the Digits setting as

restart;
frem(10^10,evalf(2*Pi));

                                 -10.

Also,

frem(10^10-1,evalf(2*Pi));

                                 -11.

but

frem(10^10+1,evalf(2*Pi));

                                 -10.

2*Pi is not an exception for that. For example, here is what happens with e,

for i from 11 to 15 do
    frem(10^i,evalf(exp(1)))
od;
                                  10.
                                 100.
                                1000.
                               10000.
                              100000.

(Corollary: e ≥ 200,000).

Alec

This could be also written as

t:=freeze(ln(x)):
thaw(MultiSeries:-asympt(LambertW(exp(t)),t,3));

                                                       2
                                ln(ln(x))     ln(ln(x))
            ln(x) - ln(ln(x)) + --------- + O(----------)
                                  ln(x)              2
                                                ln(x)

It's not shorter or better - just looks different.

The difference with Axel Vogt's reply is in the O-term. MultiSeries:-asympt has the correct form of it, and the top level asympt doesn't.

Alec

This could be also written as

t:=freeze(ln(x)):
thaw(MultiSeries:-asympt(LambertW(exp(t)),t,3));

                                                       2
                                ln(ln(x))     ln(ln(x))
            ln(x) - ln(ln(x)) + --------- + O(----------)
                                  ln(x)              2
                                                ln(x)

It's not shorter or better - just looks different.

The difference with Axel Vogt's reply is in the O-term. MultiSeries:-asympt has the correct form of it, and the top level asympt doesn't.

Alec

Not only for real x. Also, evalc would give the same result in the original post even without assumptions (evalc automatically assumes that all variables in it are real).

To do this properly, one has to look at the branch cuts of sqrt(z^2+1), sqrt(z+I), and sqrt(z-I).

z:=x+y*I:
solve({Re(z^2+1)<=0,Im(z^2+1)=0}) assuming real;

                  {x = 0, y <= -1}, {x = 0, 1 <= y}

solve({Re(z+I)<=0,Im(z+I)=0}) assuming real;

                           {y = -1, x <= 0}

solve({Re(z-I)<=0,Im(z-I)=0}) assuming real;

                           {y = 1, x <= 0}

They divide complex plane in 3 domains, in one of which (containing the real axis) the identity holds, and in 2 others - doesn't. That can be visualized with densityplot,

plots:-densityplot(abs(sqrt(z^2+1)-sqrt(z+I)*sqrt(z-I)),
x=-4..4,y=-4..4,colorstyle=HUE,style=patchnogrid,axes=box);

_______________
Alec Mihailovs, PhD

It can be done without guessing noticing that in this example

                                   x
                                  /
                                 |
          phi[k + 1] - phi[k] =  |   phi[k] - phi[k - 1] dx
                                 |
                                /
                                  0

so every next difference is the integral of the previous one, and φ_k is the sum of the differences, so the limit of φ_k is

sum(diff(x^2,x$(-n)),n=1..infinity);

                                             2
                       2 exp(x) - 2 - 2 x - x

(Knowing the Taylor series for exp(x) and the integral of x^m, this sum can be calculated mentally, without Maple.)

_______________
Alec Mihailovs, PhD

@Preben Alsholm 

The correct one is

int(exp(-x)/(1-exp(-x))-exp(-x)/x,x=0..infinity);

                                gamma

Whoever wrote the table, skipped the second term in the formula. The more general bug is

int(exp(-t*x)/(1-exp(-x)),x=0..infinity) assuming t>0;

                               -Psi(t)

Same as above it would be correct for Re(t)>0 with the same second term,

int(exp(-t*x)/(1-exp(-x))-exp(-x)/x,x=0..infinity) assuming t>0;

                    infinity
                   /
                  |           exp(-t x)    exp(-x)
                  |          ----------- - ------- dx
                  |          1 - exp(-x)      x
                 /
                   0

which Maple returns unevaluated instead.

Alec

PS For some strange reason, my comments are getting posted as answers. Sometimes they get changed later to comments , and sometimes - not. -Alec

Yes, it's that weird case in which approaching [-1, -1/2, -2] by [a, b, c] along different curves (in the (a,b,c)-space) gives different limits for the hypergeom.

In particular,

evalf(Limit(hypergeom([-k,1/2-k],[-2*k],0.99),k=1));

                             1.663750000

evalf(Limit(hypergeom([-k,1/2-k],[3*k-5],0.99),k=1));

                             0.1450000000

while

eval(3/4+z^2/4,z=0.1);

                             0.7525000000

The singularity at z=0 for non-integer k is not actually a pole in its usual sense - because the function is even. There is a pole on every branch though (there are 2 branches), with residues of the opposite sign, and a branch cut coming through them.

The series given by Maple for P is correct only in the right half-plane (with Re(z)>0 or Re(z)=0 and Im(z)>0 - in other words, with csgn(z)=1) and should have the opposite sign at z^-1 and other odd powers of z in the left half plane (with Re(z)<0 or Re(z)=0 and Im(z)<0 - i.e. with csgn(z)=-1).

The residue defined through the contour integral (with a half of the contour on one branch and another half - on another branch, so it is not actually a contour) is 0 - the integral on the right hand side of it cancels with the integral on the left hand side, and the value given by Maple,

residue(P,z=0);

                              (-2 k - 1)
                             2

is just plain wrong - either residues of both poles should be listed - this one and one with the opposite sign, or their sum which is 0. Otherwise the total sum of residues will be wrong.

By the way, the series for sqrt(z^2) and 1/sqrt(z^2) etc. are also wrong,

series(sqrt(z^2),z);

                                  z

Simplifying a function helps,

series(simplify(sqrt(z^2)),z);

                              csgn(z) z

series(simplify(P),z,1);

                         (-k)          -1      0
                    1/2 4     csgn(z) z   + O(z )

simplify(P);

                   (-k)                        (2 k + 1)
                  4     csgn(z) (1 + csgn(z) z)
              1/2 --------------------------------------
                                    z

(which is again wrong for positive integer k).

But no help for the residue, it's actually getting even worse, including z in the answer,

residue(simplify(P),z=0);

                                 csgn(z)
                             1/2 -------
                                    k
                                   4

Also, that gives a proof that 0=1 in Maple,

f:=(sqrt(z^2)-z)/2:
eval(series(f,z)=series(simplify(f),z),z=-1);

                                0 = 1

The following is interesting and puzzling a little bit,

eval(simplify(P),k=1) assuming z>0:   #limit on one branch
eval(simplify(P),k=1) assuming z<0:   #limit on another branch
expand(%+%%);                         #the sum of limits on both branches

                                      2
                                     z
                              3/4 + ----
                                     4

Alec

@Alejandro Jakubi 

In case if somebody is going to report that (I am not going to)...

In addition to this particular bug in series/RootOf, there is a bigger problem in it (which may exist even in earlier Maple versions).

In this example, for x and y with large absolute values, the curve is close to 16*x^4-y^4 = 0, which has 4 solutions, y=2*x, y=-2*x, y=2*I*x, and y=-2*I*x.

So the asympt and corresponding series/RootOf should produce 4 different solutions for 4 different roots. Hovewer, it doesn't - and adding index=1, index=2 etc. in the RootOf doesn't help - it produces the same answer.

Alec

@Alejandro Jakubi 

In case if somebody is going to report that (I am not going to)...

In addition to this particular bug in series/RootOf, there is a bigger problem in it (which may exist even in earlier Maple versions).

In this example, for x and y with large absolute values, the curve is close to 16*x^4-y^4 = 0, which has 4 solutions, y=2*x, y=-2*x, y=2*I*x, and y=-2*I*x.

So the asympt and corresponding series/RootOf should produce 4 different solutions for 4 different roots. Hovewer, it doesn't - and adding index=1, index=2 etc. in the RootOf doesn't help - it produces the same answer.

Alec

Maple traditionally ignores special cases.

One of the well-known examples is

int(x^k,x);
                                (k + 1)
                               x
                               --------
                                k + 1

which is wrong for k=-1.

In this example with hypergeom, not only conversion to StandardFunctions is wrong. Conversion to elementary is wrong even with assumptions,

convert(P,elementary) assuming k::posint;

                   (-2 k - 1)        2 1/2 (2 k + 1)
                  2           (1 + (z )   )
                  ----------------------------------
                                 2 1/2
                               (z )

eval(%,k=1);

                                   2 1/2 3
                            (1 + (z )   )
                            --------------
                                  2 1/2
                              8 (z )

while the correct answer is

convert(eval(P,k=1),elementary);

                                      2
                                     z
                              3/4 + ----
                                     4

Without assumptions, series are also wrong for positive integers k,

series(P,z,1);

             (-2 k - 1)  -1    (-2 k - 1)
            2           z   + 2           (2 k + 1) + O(z)

series(P,z,3) assuming k::posint;

                       (-2 k)                2
                      2       (2 k + 1) + O(z )

Alec

I've just voted up all you posted thus far, so hopefully you are getting there (where you'll be able to vote up :)

Alec