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++;
  );
}