Find minimal 3x3 square of distinct primes, with all row/col/diag sums being mersenne prime exponents and total sum prime

srch.3x3.mpe.gp

readvec("mpe.gp");
srch(P)={
  my(M=oo);
  my(S=[]);
  my(used=Vec(0, 9*P[#P]));
  my(ismpe=Vec(0,mpe[#mpe]));
  foreach(mpe,p,ismpe[p]=1);

  foreach(P,a,used[a]=1;
   foreach(P,b,if(!used[b],used[b]=1;
    foreach(P,c,if(!used[c]&&ismpe[a+b+c],used[c]=1;
     foreach(P,d,if(!used[d],used[d]=1;
      foreach(P,e,if(!used[e]&&ismpe[d+b+e],used[e]=1;
       foreach(P,f,if(!used[f]&&ismpe[a+e+f],used[f]=1;
        foreach(P,g,if(!used[g]&&ismpe[f+b+g]&&ismpe[c+d+g],used[g]=1;
         foreach(P,h,if(!used[h]&&ismpe[a+h+g],used[h]=1;
          foreach(P,i,if(!used[i]&&ismpe[h+b+i]&&ismpe[f+i+c],
           if(isprime(a+b+c+d+e+f+g+h+i),
            if(a+b+c+d+e+f+g+h+i<M,M=a+b+c+d+e+f+g+h+i;S=[a,b,c,d,e,f,g,h,i];);
           );
          ));
         used[h]=0;));
        used[g]=0;));
       used[f]=0;));
      used[e]=0;));
     used[d]=0;));
    used[c]=0;));
   used[b]=0;));
  used[a]=0;);
  S
};

s(i)=concat(if(i<10,"_",""),Str(i));

{
 m=eval(getenv("m"));

 [a,b,c,d,e,f,g,h,i]=srch(primes([3,m]));

 print("{[",a,",",e,",",f,";\n",h,",",b,",",i,";\n",g,",",d,",",c,"];}");
 print("\n_________",g+b+f);
 print("________/");
 print(s(a)," ",s(e)," ",s(f)," __",a+e+f);
 print(s(h)," ",s(b)," ",s(i)," __",h+b+i);
 print(s(g)," ",s(d)," ",s(c)," __",g+d+c);
 print("________\\");
 print(a+h+g," ",e+b+d," ",f+i+c," ",a+b+c);
 print("\nsum: ",a+b+c+d+e+f+g+h+i);
}

Read mpe.gp from this gist:
https://gist.github.com/Hermann-SW/f44c3500b8db158eb87ae8c2fd010e21

hermann@7950x:~$ time m=239 gp -q < srch.3x3.mpe.gp 
{[5,37,47;
31,7,23;
53,17,19];}

_________107
________/
_5 37 47 __89
31 _7 23 __61
53 17 19 __89
________\
89 61 89 31

sum: 239

real	0m1.985s
user	0m1.425s
sys	0m0.560s
hermann@7950x:~$ 

Discussed here:
https://www.mersenneforum.org/node/1055343?p=1061301#post1061301