@janhardo 

@Carl Love The replacement is necessary for r>0 too.

@Carl Love 
Thanks. I wonder why the formula was not included somewhere (e.g. in `simplify/Psi` or some convert); maybe it was forgotten.

Your extension is wrong for Psi(r) when r is a negative rational; probably because the original Psi has 2 arguments (or recursive overload problem); it would be interesting to investigate.

Edit. The mistake: (n-1)*(q+1) must be replaced with (n-1)*q+p

Just curious: what's wrong with a .mla update? It used to work well in the past. Is it not the case in "modern Maple"?

@tomleslie Maybe because you are using a 65 bit OS  :-)

@nguyenhuyenag 
The question is what information about f is contained in discriminant(f,[k])?
When such a computation is useful for a polynomial in several variables?

Also, if you try  g:=eval(f, a=d+e) in your last example, it does not terminate.

What are you trying to obtain?
E.g. in your last example, what does discriminant(f, [k]) represent for f?

@afernande1008 

It's a bug in Physics:-diff. Use diff instead.

@Axel Vogt I don't think so. The function is

f := f0 + r:  # where
f0 := x -> add(- i * sin((i + 1) * x + i), i=1..5):
r := rand(-1. .. 1.):
plot(f, -10 .. 10, numpoints=5000);

and it is practically nowhere continuous. DirectSearch does not try to "clean-up" the noise r().

The program is written in Algol. The syntax is close to Maple, so, if you produce from the pdf a correct text file (with indentation) it should be easy to obtain a Maple version.

@Thomas Richard 

For k=64/81 the polynomial ir reducible:

((3*w - 2)*(162*w^4 - 297*w^3 + 72*w^2 + 48*w + 32))/81

For solvable irreducibles, the values of k=m/n, |m|, |n| <= 100  are: -16/9, 1/81, 25/9, 80/81.

@Thomas Richard 

The value k=1/81 is interesting because the polynomial is irreducible and solve with option explicit works. 

@Carl Love 

Actually Im( ln(1+ z) ) does not have a limit at z=-1. The limit of ln(1+z) at z=-1 is "complex infinity" (denoted by infinity+I*infinity) in Maple.

Note that assuming does not make much difference between `<` and `<=`.

E.g.

limit(z^n, n=infinity) assuming abs(z)<=1; #   0

@David Sycamore 

Probably you have not noticed the new name for the procedure.

Download SP.mw

SP:= proc(N::posint, p1::prime:= 2)
local S:= Array([p1], datatype= integer[8]), Used:= table([p1= ()]), s, k, p:= p1, i,j,ok;
    for k from 2 to N do
        s:= `+`(convert(p, base, 10)[]);
        p:= s; ok:=false;
        for i do  for j by 2 to 10^i-1 do
            if p::prime and Used[p]<>() then ok:=true; break fi;
            p:= Scale10(s, i) + j;
        od; if ok then break fi;
        od;
        S(k):= p; Used[p]:= ()
    od;
    [seq(S)]
end proc: