- --- In primenumbers@yahoogroups.com, "David Broadhurst" <d.broadhurst@...> wrote:
>

I can shed some light on this behaviour.

> --- In primenumbers@yahoogroups.com,

> "richard_in_reading" <richard_in_reading@> wrote:

>

> > some primes pop in for an iteration and go again

> > and some come in and stick around

>

> How to predict which sort of behaviour?

>

> For example the divisor

> 93408839 | 137^(137^(137^(137^(137^137)))) + 184

> pops in and then out:

>

> {pmod(a,m)=local(k,q,t);k=#a;q=[m];t=a[k];

> for(j=2,k-1,q=concat(eulerphi(q[1]),q));

> for(j=1,k-1,t=Mod(a[k-j],q[j])^lift(t));t}

>

> for(k=5,7,print(pmod(vector(k,j,137),93408839)+184))

>

> Mod(40934825, 93408839)

> Mod(0, 93408839)

> Mod(10826080, 93408839)

>

For any prime q, define b(q,n) by the recurrence relation

b(q,n+1)=q^b(q,n), with initial condition b(q,0)=1.

[So b(q,1)=q, b(q,2)=q^q, and so on.]

Theorem: For k any integer, if a prime p divides (b(q,n)+k) and

(b(q(n+1)+k) then p divides (b(q,n+2)+k).

Proof: Assume the conditions of the theorem hold.

Then b(q,n+1) = b(q,n) mod p,

so q^b(q,n+1) = q^b(q,n) mod p,

i.e. b(q,n+2) = b(q,n+1) mod p,

so (b(q,n+2)+k) = (b(q,n+1)+k) mod p = 0.

Q.E.D.

In words: once any prime p has "popped in" to the divisor list for /two consecutive/ terms of the sequence, then it stays there for ever.

(So those 20 mins of computer time described in my post of yesterday were unnecessary, except as experimental confirmation.)

-Mike Oakes - Dear all,

I come back to this topic, looking for a working "pow_mod".

This :

{mypmod(b,n,p)=local(m=[p],f=0);

while(n=n-1,m=concat(eulerphi(m[1]),m));

for(p=n=1,#m,n=lift(Mod(b,m[p])^n);

if(f=f+(n*log(b)>=log(m[p])),n=n+m[p]));n%m[#m];}

yields:

mypmod(2,3,5)

= 2

but 2^2^2 % 5 = 16 % 5 = 1.

(and idem for 2^2^2^2 - of course.).

David gave some other code, in

https://groups.yahoo.com/neo/groups/primenumbers/conversations/topics/20526

:

"

Below is a Pari-GP procedure "pmod(a,m)" to compute

a[1]^(a[2]^(a[3]^ ... ^(a[k-1]^a[k]) ... ) modulo m

where the modulus "m" need not be prime.

(...)

{pmod(a,m)=local(k,q,t);k=#a;q=[m];t=a[k];

for(j=2,k-1,q=concat(eulerphi(q[1]),q));

for(j=1,k-1,t=Mod(a[k-j],q[j])^lift(t));t}

"

Although this may work for prime m

(at least, pmod([2,2,2],5) = Mod(1,5) as should),

it gives:

pmod([2,2,2],4)

= Mod(1, 4)

which is most certainly wrong.

So, to put it short, has anyone a working pmod() in his "library" ?

Thanks in advance!

Maximilian

--- In primenumbers@yahoogroups.com, "David Broadhurst"

<d.broadhurst@...> wrote:

> fail

(...)

I consider that his code should not be trusted.

David