This site is supported by donations to The OEIS Foundation.

M. F. Hasler/A039986

From OeisWiki
Jump to: navigation, search

Code for this sequence:


A039986_row(n)={ /* compute terms of length n */
 if(n>1 /* criteria for the last digit apply only for length > 1 */
 , local(D=eval(Vec("0245681379")) /* all possible digits, with possible ending digits coming last: equivalent to vecsort([0..9],(a,b)->isprime(10+a)-isprime(10+b)) */
   , u=vectorv(n, i, 10^(n-i)) /* powers of 10: dot product d.u is slightly faster & shorter than fromdigits(d) */
   , nextperm()= /* change the "global" variable p to the "next" permutation of [1..n] */
       for(i=2, n, /* i = the position of the index we will decrease. (We permute "left first": 123..., 213..., 132... etc) */
         (t=p[i]-1)||next; /* can we decrease this index without getting 0? else goto next position */
         while(setsearch(Set(p[i+1..n]), t) /* is this decreased index available to the left? */
                 || d[t]==d[p[i]] /* ... and does this change the digit? ... */
            , t--||next(2) /* ... else try next smaller; goto next position if we run into 0 */
            );
         i<n /* rightmost position? */
         || bittest(650, d[t]) /* then digit must be among 1379 */
         || return /* else (a digit other than 1379 reached the right) we're done: all subsequent perms will be composite */;
         /* now all criteria are satisfied if we replace p[i] with the smaller t & reorder the other indices to the left of p[i] */
         p=concat([setminus(Set(p[1..i]), [t]) /* reorder the indices (other than t) to the left */
                  , t, p[i+1..n]]); /* then append t (= new p[i]), and the other indices to the right (unchanged) */
         return(1)
       ) /* end for(i); also the end of nextperm(): returns 0 if all permutations are done */
   , L=List() /*list of results*/, f /* "found" */, p /* the permutation */, d /* the chosen digits in increasing order */
   ) /* end of "local" variables and procedures; now follows the main loop */;
   forvec( i=vector(n, j, [7^(j==n), 10]) /* selection of the digits: indices 1..9 correspond to digits '0','2',...,'7','9' ; only increasing vectors i are generated, and last = n-th index must be >= 7, i.e., digit 1,3,7 or 9 */
   , vecsum(d=vecextract(D, i)/* list of digits corresponding to the list of indices i */
            )%3 || next /* sum of digits must not be divisible by 3 */; 
     f=0 /* nothing found so far */;
     p=[1..n] /* initial permutation = id */; 
     until(!nextperm() /* after doing the body of the loop, go to next relevant permutation; end if all done. */
     , isprime(vecextract(d, p) /* take selected digits in order of permutation p */ *u /*convert to number*/) /* a prime? */ 
       && (f && next(2) /* if we had already found another prime: no good; goto next selection of digits */
           || f=p /*else set "found" equal to the current permutation */)
     ) /* end of until() */ ; 
     f /*if we have found exactly one prime */ && d[f[1]] /* first digit nonzero */ 
     && listput(L, vecextract(d, f)*u) /* add the corresponding number to the list */
   , 1 /* only consider increasing indices */
   ); 
   Set(L) /* return the solutions in increasing order */
 , /* else (n=1): return the 1-digit primes (2,3,5,7) */
   primes(4) /* equivalent to [2,3,5,7] but OEIS inserts a space after each "," so the explicit list takes more space */
 )/*end if n>1*/
} \\ Returns all terms with n digits. - M. F. Hasler, Jul 01 2018