rootsof1(p,n)= { local(a); /* Primitive solutions to a^(p-1) % p^n == 1 for odd prime p */ a = vector((p-3)/2, i, i+1); for(j=2, n, a = vector(#a, i, lift(Mod(lift(Mod(a[i],p^j)^(p-1))\p^(j-1),p)/Mod(a[i],p)^(p-2))*p^(j-1)+a[i]); ); vecsort(concat([a, vector(#a, i, p^n-a[i]), [1, p^n-1]])) } A125609(n,p=3)= { local(v,i); v = rootsof1(p,n); i = 0; while (1, for (j=1, #v, if (isprime(i*p^n + v[j]), return(i*p^n + v[j])); ); i++; ); }