I have extended the nub of this thread a little further, but it still seems to be computationally useless.

For example b = 3 and p = 61 so that mp = 2^61 - 1, I ran this:

Code:

? p=61;mp=2^p-1;D=mp-1;V=factor(D,1000000);forbigdiv(D,d -> r=lift(Mod(3,mp)^(D/d));if(2^logint(r,2)==r,print(">>>"d)))
>>>1
>>>61
>>>3
>>>183
>>>9
>>>549

Where "forbigdiv" is given

here
This means 3^((mp-1)/549) == 2^n for some n OR 3^4200078341008550 == 2^n for some n, thus implying M61 is prime.

I will experiment with other bases b, for example p ...