## Alec Mihailovs

Dr. Aleksandrs Mihailovs

## 4470 Reputation

20 years, 25 days
Mihailovs, Inc.
Owner, President, and CEO
Tyngsboro, Massachusetts, United States

## Social Networks and Content at Maplesoft.com

I received my Ph.D. from the University of Pennsylvania in 1998 and I have been teaching since then at SUNY Oneonta for 1 year, at Shepherd University for 5 years, at Tennessee Tech for 2 years, at Lane College for 1 year, and this year I taught at the University of Massachusetts Lowell. My research interests include Representation Theory and Combinatorics.

## Agree...

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

## Different form of an answer...

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

## Different form of an answer...

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

## with 23 digits...

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

## with 23 digits...

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

## freeze and thaw...

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

## freeze and thaw...

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

## Gem 12...

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

## Gem 11...

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

## The correct one...

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
```

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...

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

## Another series/RootOf problem...

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

## Another series/RootOf problem...

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

## 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 )
```

```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