Smarandache prime(s)

[Batalov]

The right thing to do is to throw the sieve on a GPU. (A wizard like Anand 'axn' could make it fly!)

$ gp -q
? F(n,p)=x=Mod(10,p);((999990001+(((9999001+(((120999998998*111^2*x^179-1099022)/121*x^2702-1000)/111^2-81)*1111^2)*x^36000-x^4)/1111^2-81)*11111^2)*x^(5*n-49995)-99999*n-x^5)/99999^2;
? v=readvec("../M/smarandacheremaining.txt");
? for(i=1,#v,n=v[i];if(n<90000,next);if(n>100000,break);forprime(p=20000000,10^8,if(F(n,p),,print(n" "p);break)))
91213 93402091
91789 89739989
...
#obviously this is not sieving but a proof-of-concept

After all, the n>10^5 tests take ~4hrs on a Phenom (but probably only a couple hours on a bunch of Haswell cores which anyone can quickly order for less than 2 cents per physical-core-hour on Amazon; 4 cents per test, cheap; one can try even 36 parallel tests per c4.8xlarge but it is likely that 18 tests bound to individual cores will be faster).

((((((((1490840987654358*10^179-1099022)/121*10^2699-1)/111^2*1111^2-89981)*10^35999-1)/1111^2*11111^2-899981)*10^449999-1)/11111^2*111111^2-8999981)*10^(6*100093-599989)-999999*100093-10^6)/999999^2 is composite: RES64: [9A5D7922B39D12C5] (14699.5044s+0.8407s)
((((((((1490840987654358*10^179-1099022)/121*10^2699-1)/111^2*1111^2-89981)*10^35999-1)/1111^2*11111^2-899981)*10^449999-1)/11111^2*111111^2-8999981)*10^(6*100213-599989)-999999*100213-10^6)/999999^2 is composite: RES64: [F386E35B6397B949] (14991.0264s+0.2702s)

[danaj]

Originally Posted by Antonio

A quick Perl script only took 50min. to reduce the range 100k-200k to a PFGW ABC file with 2662 indices to check, using mod 30 test plus a bit of extra trial division.

I edited gmp_main.c, bumped up the BGCD3 number to 20k, made is_prob_prime just do the pretest without extra trial division. Total 5 lines of change (bump BGCD3_PRIMES and BGCD3_NEXTPRIME, #if 0/#endif around the tree sieve calls, return res immediately after pretest).

perl -Mntheory=:all -E 'my $s=""; for my $n (1..200000) { $s .= $n; say $n if $n >= 100000 && is_prob_prime($s); }'

took under 9 minutes to reduce to 2096 candidates (trial division to ~200k using mpz_gcd). Bumping to 148k (trial division to 2M) pushes the time to 25 minutes but only 1779 candidates left.

Of course there are better ways to sieve this, but this is why I don't understand all the earlier talk about checking divisibility by 3 when it takes a fraction of a single PRP test to do all factors under 1M using naive methods (the code above does division by 2,3,5 directly on the string, so never even turns into a mpz_t).

[ewmayer]

Just curious, how deep are y'all sieving, and how long that take for a typical work interval? I banged together some simple dumb-sieving code (test each odd term of the sequence in turn) using my own own multiword-int library and running 1-threaded on my 2GHz Core2 I get from sieving using the first 2M odd primes:

10^5 - 10^5+1000: 17 candidates survived, ~200 sec to TF the range

Factoring in the speedup of the key MUL instruction and typical clock rates, that translates to ~60 sec on a **y Bridge or later CPU.

[danaj]

Originally Posted by ewmayer

Just curious, how deep are y'all sieving, and how long that take for a typical work interval? I banged together some simple dumb-sieving code (test each odd term of the sequence in turn) using my own own multiword-int library and running 1-threaded on my 2GHz Core2 I get from sieving using the first 2M odd primes:

10^5 - 10^5+1000: 17 candidates survived, ~200 sec to TF the range

Factoring in the speedup of the key MUL instruction and typical clock rates, that translates to ~60 sec on a **y Bridge or later CPU.

I answered above, but more details. 3770k running 12 other full-cpu tasks. For your range: 3.3 seconds to sieve to 40k (26 candidates). 8.3 seconds to sieve to 2M (17 candidates). 11.8s to sieve to 10M (16 candidates). Almost 4 minutes to sieve to 100M (15 candidates) so we'd need to rethink the method somewhere past 10M.

- take the text string, handle cases of length 1 with a switch statement. Return 0 if last character is even or 5. sum digits and return 0 if divisible by 3. I've tried doing the same for divisibility by 11 but it didn't seem worth the trouble.

- convert to mpz_t.

- Do simple tests if really small. Irrelevant for our case.

- mpz_gcd_ui with 3-53 and 59-101. (in theory we could skip 3 and 5 here, but this is reachable by other paths)

- mpz_gcd with a precalculated primorial of primes to 997. Again some wasted work from the primes above. Oh well.

- If between 300 and 700 bits, another mpz_gcd with product of primes to 10k removing up to 997. Not our case.

- If over 700 bits, an mpz_gcd with product of primes to 40k removing up to 997.

For this example I return the result here. Normally I would go on with possibly some simple tree sieving if over 1600 bits, followed by ES BPSW.

[Batalov]

Here is a primitive sieve.

/**************************************************/
/* A Smarandache prime simple sieve  (from stdin) */
/*                      (copyleft) 2015 S.Batalov */
/* Example use:                                   */
/*   gcc -O3 -o Smar_sieve Smar_sieve.c -lgmp     */
/*   ./Smar_sieve 1E6 1E8 < in.list > out.list    */
/**************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

typedef unsigned long long uint64_t;

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t c) {
  uint64_t d; /* to hold the result of a*b mod c */
  /* calculates a*b mod c, stores result in d */
  asm ("mov %1, %%rax;"        /* put a into rax */
       "mul %2;"               /* mul a*b -> rdx:rax */
       "div %3;"               /* (a*b)/c -> quot in rax remainder in rdx */
       "mov %%rdx, %0;"        /* store result in d */
       :"=r"(d)                /* output */
       :"r"(a), "r"(b), "r"(c) /* input */
       :"%rax", "%rdx"         /* clobbered registers */
      );
  return d;
}

uint64_t powmod(uint64_t i, uint64_t j, uint64_t k) { /* a general, well-known code */
  uint64_t r = 1;
  while(j > 1){
    if(j&1) r=mulmod(r,i,k);
    i=mulmod(i,i,k);
    j>>=1;
  }
  return mulmod(r,i,k);
}

int ncmp(const void* a, const void* b) {return *(int *)a - *(int *)b; }

#define submod(a,b) {uint64_t bb=b%f; if(a < bb) a += f-bb; else a -= bb;}
#define divmod(a,d) mpz_set_ui(R, d); mpz_invert(R, R, P); a = mulmod(a,mpz_get_ui(R),f)

main(int argc,char **argv) {
  uint64_t f,b,c,r,k,t,C,M,p,p1,p2,T[6];
  char buf[64], *s;
  int *v, nv=0,n,n0=999999,N=0,dn;
  mpz_t P, R;

  if(argc<3) { printf(" Use: %s p P < input_file > output_file\n"
               "\t 10^6 < p < P are numeric parameters\n"
               "\t a Smarandache simple sieve, v.1.1\n"
               "\t (copyleft) 2015 S.Batalov\n",argv[0]); exit(0); }
  p1 = atof(argv[1]);
  p2 = atof(argv[2]);
  if(p1<1000000) { printf(" can only sieve for primes p > 1E6\n", p1); exit(0); }
  mpz_init_set_ui(P, p1); mpz_init (R); 
  v = calloc(150000,sizeof(int));
  while((s=fgets(buf,64,stdin))) {
    n=strtol(s,NULL,10); if(!n) continue;
    if((n%6)!=1) { printf(" %d in input file is not =1 (mod 6)\n", n); exit(0); }
    if(n<100000 || n>999999) { printf(" %d in input file is not a 6-digit decimal\n", n); exit(0); }
    v[nv++] = n;
  }
  if(!nv) exit(-1);
  qsort(v,nv,sizeof(int),ncmp);
  printf("ABC Sm($a) // sieved from %lld to %lld\n",p1,p2);
  while(nv) {
    mpz_nextprime (P, P);
    if ((f = mpz_get_ui (P)) > p2) break;
    t = powmod(10,179,f);
    C = mulmod(t,15208068915062105958ULL%f,f);
    submod(C,11211123422ULL);
    t = powmod(10,2699,f);
    C = mulmod(C, t, f);
    submod(C,1109890222ULL);
    divmod(C, 12321ULL);
    t = powmod(10,35999,f);
    C = mulmod(C, t, f);
    C = mulmod(C, 123454321ULL, f);
    submod(C,1110988902222ULL);
    divmod(C, 1234321ULL);
    C = mulmod(C,12345654321ULL % f, f);
    t = powmod(10,449999,f);
    C = mulmod(C, t, f);
    submod(C,1111098889022222ULL);
    divmod(C, 123454321ULL);
    t = powmod(10,6*v[0]-599989,f);
    C = mulmod(C, t, f);
    M = (v[0]*999999ULL+1000000)%f;
    r = 1000000000000000000ULL%f;
    T[1] = mulmod(r, r, f);
    T[2] = mulmod(T[1], T[1], f);
    T[3] = mulmod(T[1], T[2], f);
    T[4] = mulmod(T[2], T[2], f);
    T[5] = mulmod(T[2], T[3], f);
    for(n=0; ; ) {
        int dn = v[n+1]-v[n];
        if (C == M) {
            fprintf(stderr,"%lld | Sm(%d)\n",f,v[n]);
            if(--nv==n) break;
            for(k=n;k<nv;k++) v[k]=v[k+1];
        } else if(++n>=nv) break;
        if(dn<=30) C = mulmod(C, T[dn/6], f);
        else {
          t = powmod(r, dn/3, f);
          C = mulmod(C, t, f);
        }
        M = (M + 999999ULL*dn) % f;
    }
  }
  for(n=0; n<nv; n++) printf("%d\n",v[n]);
  exit(0);
}

[ec2-user@ip-172-30-3-137 Sm]$ ./Smar_sieve 85592030 200000000 < smarandacherem100.txt |more
85592939 | Sm(199159)
86964313 | Sm(171073)
87543821 | Sm(144109)
87588569 | Sm(193891)
88385917 | Sm(140263)
89333941 | Sm(151069)
89472527 | Sm(131929)
89784251 | Sm(179041)
90260101 | Sm(132127)
94177561 | Sm(187441)
94314593 | Sm(127789)
94341713 | Sm(143227)
99484159 | Sm(198163)
99897437 | Sm(159007)
102931133 | Sm(162397)
106940389 | Sm(112051)
107352643 | Sm(186403)
109119211 | Sm(180811)
111669431 | Sm(196639)
112264637 | Sm(192307)
113959421 | Sm(153019)
113973677 | Sm(133267)
114024563 | Sm(162091)
114205153 | Sm(179239)
114597209 | Sm(122389)
114764879 | Sm(168709)
116039461 | Sm(152719)
117367157 | Sm(106543)
118750189 | Sm(129883)
119512619 | Sm(153487)
123252557 | Sm(173323)
126016399 | Sm(175561)
126510031 | Sm(175039)
126631811 | Sm(160057)
127430837 | Sm(165529)
129924467 | Sm(197413)
132909737 | Sm(158833)
135164857 | Sm(189811)
135234367 | Sm(180097)
135277343 | Sm(102469)
136059929 | Sm(164683)
137179997 | Sm(172543)
...

[wombatman]

Originally Posted by Batalov

I guess I'll start sieving 10⁵ < n < 10⁶, modestly as usual.
I've already given the obligatory hint in post #8: in each 10^k <= n < 10^k+1 region there is a simple polynomial formula. It is much easier to sieve that formula than pile up "numbers as strings".

Exercises:
A) calculate x+2x²+3x³+4x⁴+...+kx^k
B) calculate x^k-1+2x^k-2+3x^k-3+4x^k-4+...+(k-1)x+k (where x=10 and k=1..9). Hint: it is = (10[SUP]k+1[/SUP]-C)/81
C) calculate sum{j=10..k} j x^k-j (where x=100).
D) etc (where x=1000 and j=100..k), etc etc etc.
E) sum them up with proper weights

*sigh*

Alright, I know this is probably very basic, but I cannot for the life of me figure out why these are steps. I can see how you get most of the formula for , but I don't understand why it's divided by 81 = (10-1)^2. Can you explain a bit more and/or point me to something I can read to understand it better?

Thanks.

[Batalov]

I'll be glad to.

Here is a path that helps with a lot of sums with similar properties (e.g. coefficients at x^j depend on j linearly):
Let
A(k) = x+2x²+3x³+4x⁴+...+kx^k
mutliply by x =>
xA(k) = x²+2x³+3x⁴+...+(k-1)x^k+kx^k+1
Subtract one from the other, and we get
(x-1)A(k) = -(x+x²+x³+x⁴+...+x^k)+kx^k+1 = -(x^k+1-x)/(x-1)+kx^k+1 -- because we know easily how to sum up the "blue" part.
Then, A(k) = ((x-1)kx^k+1-x^k+1+x)/(x-1)^2, where x could be 10 (with k=1..9). Or for the two-digit chunk part, it is 100 and the denominator is 99^2; and for the two-digit chunk part, it is 100 and the denominator is 999^2; and so on.

This sum is related to sum that we need:
F1(k) = x^k-1+2x^k-2+3x^k-3+4x^k-4+...+(k-1)x+k, because
A(k)+F1(k) = k * (1+x+x²+x³+x⁴+...+x^k) = k(x^k+1-1)/(x-1).
So we can combine and then simplify the expression (or we could have started right away with F1(k) and compared it to xF1(k) and so on). That's the basic idea. I think I missed a few terms above, but I only wanted to show how it works. Try it on a sheet of paper and fix my little typos along the way (they are only in this message because it is a 'free finger painting'; the formulae and the siever above are of course debugged and are clean).

The other way to see this is simply with numbers (with x gone and being explicitly 10 and simply serve as a position in decimal):
Suppose I want to find the short form for:

 12345678 = Sm(8). I'll multiply it by 10 and then subtract Sm(8) from 10*Sm(8):

123456780 = 10*Sm(8)
-12345678 = Sm(8)
---------------
111111102 = 9*Sm(8) = (10^9-82)/9 ==> Sm(8) = (10^9-82)/81

Now, for two-digit chunks, we want to deal with two-digit chunks alone (i.e. starting with 10, let's call it Sm'(k) )
  10111213141516171819202122232425 = Sm'(25)

1011121314151617181920212223242500 = Sm'(25)*100 
- 10111213141516171819202122232425 = Sm'(25)
------------------------------------------------------------------------
1001010101010101010101010101010075 = Sm'(25)*99 which will look much better 
if I multiply it by 99...
..then with a few more operations, the denominator for that part will be 99^2 and the numerator will be some (C*10^N-D)

and so on, similarly for 3-digit chunks, 4-digit, etc

Phew! :-)

[R. Gerbicz]

Originally Posted by Batalov

       if(dn<=30) C = mulmod(C, T[dn/6], f);
        else {
          t = powmod(r, dn/3, f);
          C = mulmod(C, t, f);
        }

You will compute a lot of times the same r^(dn/3) mod p if dn>30.

There is a better way: let maxe=maximum(v[n+1]-v[n])/6 and sh is an integer for that 4^(sh-1)<=maxe<4^sh.
With (less than) 2^(sh+1) mulmods precompute:
low[i]=u^i mod p for i<2^sh and
high[i]=u^(i*2^sh) mdo p for i<2^sh. (where u=10^36).

In array E store in increasing order the possible values of (v[n+1]-v[n])/6 and in one run
compute T[E[i]]=mulmod(high[E[i]>>sh],low[E[i]&hash],p) (where hash=(2^sh)-1).

in this way you have precomputed every needed 1000000^(6*E[i])=u^E[i] mod p values with cnt+2^(sh+1) mulmods (if there are cnt different v[n+1]-v[n] values).

(And you need to update the maxe,sh,E array,cnt only if you've found at least one hit for p).

[Batalov]

This sounds good. I'll fit this in. Thanks!

What I did for starters was I checked the distribution of gaps and they were heavily skewed towards gaps <=30, so I wrote that simple prep. The ugly math is done only once, and then the sieve pass is run on the survivors (which is down from ~1450 to 1145 *).

I will sieve to 20G and then run 100000<n<200000 on the cloud. On Haswells, tests take 2.5hrs 1.25 hrs (actually two times faster; I ran them as a test of speed but forgot that all cores on that node are already busy with 'Arithmetic progressions, L=3' LLR processes which will soon finish. I only wanted to find one for my collection, and now I already have) and later will take much more.

Sm(100117) is composite: RES64: [D5D4C7CFBFB43BAE] (8958.0309s+0.0551s)
Sm(100327) is composite: RES64: [932DED2A9750E421] (8902.2297s+0.0556s)
...
Generic modular reduction using generic reduction FMA3 FFT length 160K, Pass1=640, Pass2=256 on A 1640037-bit number
Sm(100801) is composite: RES64: [D185691DF429494D] (9154.2415s+0.0555s)
Generic modular reduction using generic reduction FMA3 FFT length 160K, Pass1=640, Pass2=256 on A 1641591-bit number
Sm(100879) is composite: RES64: [9B63C0DCF33CC4D5] (9161.4165s+0.0572s)

______________
*not bad, LaurV, eh? Iirc, you expected "18000 (just an estimation) candidates left. You may sieve to 10^13 and have 17500. So?" Well, your estimation is way off. I've already sieved away 21% (...approaching 10G).

P.P.S. Actually looked back and noticed that I posted the code before I've made one final trivial acceleration:
t = powmod(r, dn/3, f); --> t = powmod(T[1], dn/6, f);

[LaurV]

Originally Posted by Batalov

*not bad, LaurV, eh?

grrr.. well... the last few posts from you and Robert are very instructive...

[Batalov]

Just in case, looked up the OEIS sequence comments today and all the candidates under 10^5 are reported done by maxal.

I am approaching 140000 (and in parallel am sieving up to 10^11 for 140000<n<200000); I will run the complete 100000<n<200000 range.