a365295(pmax, nmax=34) = {my(ndim=max(1110,floor(1.11*sqrt(pmax))), p=vectorsmall(ndim), n=1); p[1] = 1; for (k=2, pmax, if (ispower(k), p[n++]=k; if(n==ndim,break))); my(v=vectorsmall(pmax)); for (i=1, n-1, for (j=i+1, n, my(s=p[i]+p[j]); if (s<=#v, v[s]++))); for (r=1, nmax, for (k=1, pmax, if (v[k]==r, print([r,k]); break)))}; a365295(10^7); Results for a365295(2*10^9, 50) (32 GB RAM needed): [1, 5] [2, 17] [3, 129] [4, 468] [5, 1025] [6, 2628] [7, 12025] [8, 32045] [9, 27625] [10, 138125] [11, 430625] [12, 204425] [13, 160225] [14, 2010025] [15, 2348125] [16, 801125] [17, 1743625] [18, 2082925] [19, 4978025] [20, 4005625] [21, 12325625] [22, 30525625] [23, 73046025] [24, 5928325] [25, 13287625] [26, 46437625] [27, 45177925] [28, 35409725] [29, 120737825] [30, 52073125] [31, 66438125] [32, 29641625] [33, 32846125] [34, 956974625] [36, 77068225] [37, 232210225] [38, 162109025] [40, 148208125] [41, 577500625] [42, 1301828125] [45, 885243125] [48, 243061325] [49, 426999625]