Smarandache prime(s)

[firejuggler]

890? what? what did I say? ( feel free to ignore ths post)

[danaj]

From what I can tell, you're finding patterns in the indices that will always result in an input divisible by a small prime. My point is that this is adding complexity to reduce time in an area that was already a minuscule portion of the total time. Reducing 2 hours down to 1 hour 59 minutes 59.999 seconds. Unless your isprime equivalent has no checks at all for small divisors built in, in which case you should add something that looks for more than just these cases. The point about being sure was the "only prime indexes need to be checked" not the divisibility by 3 test.

A filter on the index that was not related to n having a small factor would be of interest.

Clearly we were not aware of the Mathworld references. In 2012 the OEIS entry indicated it had been searched to 5000. CRG4 expanded that to 20k in April 2014, and I went to 64k a week later. I wish we'd known Weisstein had already looked to 64,728. Sloane mentioned this on 29 Sep this year, and 4 days later Alekseyev's result to 77k showed up.

I don't recall the process I used before, but it was likely something like do-simple-divisiblity-tests, output candidates to a file, run OpenPFGW on candidates, run anything left (doubtful) through BPSW. It may have just been running BPSW in GMP, but that's a lot slower than OpenPFGW for 10k+ digit inputs. Quick enough to double check the first 20k entries in a few hours though, and typically I run these to exercise and improve my software more than to get the actual result.

[danaj]

Originally Posted by LaurV

edit: pari goes quite fast to few thousands and stays there. Interesting that the "isprime" version is faster than other "probable prime" (modular) tests, because probably the builtin "isprime" filters the numbers which are obviously not primes (like even, etc). It can be seen how the "counter" progresses "in steps", jumping over big chunks of integers.

Almost as if it already doing some simple divisibility tests and they are really fast. Going to 3k I didn't see any speed difference between isprime and ispseudoprime using Pari 2.8.0. The former would be a bit problematic if you did find one however (when given a 500k digit input, will Pari's APRCL complain or go into effectively infinite computation?).

Pari/GP's ispseudoprime is 2-4x slower than my GMP is_prob_prime. OpenPFGW starts being faster than the GMP test somewhere in the 3k-10k range. It's more drudge work and requires a second step, but at 140k digits it is 5x faster than GMP M-R and 10x faster than Pari's M-R (ispseudoprime(n,1)) so seems well worth it. This task is easy to parallelize so parfor etc. don't help the performance argument (though they're convenient).

If someone has other results, I'd like to hear them. The last time I did a complete comparison of the speeds was in 2013, Pari 2.6.2.

[LaurV]

How about sieving?
That should be extremely efficient because it will not need multiple-precision arithmetic, all calculus is done modulo p (the prime which sieves the bunch). In C, someone can implement this to be 100 times more efficient (pari is not the king of bit-operations, and here 90% of time is screen printing).

Save this as "smake.gp", to run use "\r smake.gp", to resume see the end of code lines.
It will create a file with factors (2,3, and 5 make no sense to be saved), and a "sieve" vector from which the remaining candidates can be extracted and tested with pfgw, etc. It can be improved, this is just I could produce in the last 2 hours, in between tasks at job. It starts with the first 1 million candidates (n=1M) and it gets faster as p increases because it prints less into the file and the screen.

trialdiv(n=121,max_p=11)=
{
    forprime(p=2,max_p,
        if(n%p==0,
            return(p))
    );
    return(max_p+1);
}

smake_file_tf(max_n=100000, max_p=40000, min_factor=100, file="smake.txt", file_factors="smake_factors.txt")=
{   
    n=1;
    cnt=1;
    smarandache=1;
    until(n>max_n,
        n++;
        smarandache=eval(concat(Str(smarandache),Str(n)));
        if((a=trialdiv(smarandache,max_p))>max_p,
            write(file, smarandache);
            cnt++,
        /*else*/
            if(a>=min_factor,
                write(file_factors,n,",",a)
            )
        );
        printf("... %d, %d  (div=%d) ...             %c", n, cnt, a, 13)
    );
}

/*to create the global file discussed below (pari is silly when working with files)*/
smake_create_sieve_file(file="smake_sieve.txt")=
{
    v=vectorsmall(125000,x,255);
    v[1]=254;
    write(file,v);
}

/*this routine uses a global file, in case we want to resume sieving, this file must exist*/
smake_sieve(min_p=2, max_p=10^10, min_factor=7, file="smake_sieve.txt", file_factors="smake_sieve_factors.txt")=
{
    v=read(file);
    
    cnt=0;
    for(i=1,125000,
        for(j=0,7,
            if(bitand(v[i],1<<j),
                cnt++
            )
        )
    );
    
    p=nextprime(min_p);
    until(p>max_p,
        if(p>=min_factor,  /*write the factors if we need them*/
            write(file_factors,"Divisible by "p":")
        );

        /*compute a table with 10^n%p to speed up the things*/
        t=vector(6);
        t[1]=10%p;
        for(i=2,6,
            t[i]=(10*t[i-1])%p
        );
        
        /*start sieving*/
        n=1;
        s=1;
        i=1;
        j=1;
        x=1;
        while(n<1000000,
            /*go to the next candidate*/
            n++;
            if(n==10||n==100||n==1000||n==10000||n==100000,
                x++   /*faster than log*/
            );
            
            /*calculate s%p, where s=smarandache(n), all calculus is simple precision, small integers!!*/
            s=(s*t[x]+n)%p;
            
            /*calculate the index into the sieve (more complicate than C, because pari can't manipulate bits*/
            j<<=1;
            if(j>128,
                j=1;
                i++
            );
            
            /*check if we have a hit of the sieve*/
            if(s==0 && bitand(v[i],j),      /*stupid and slow bit testing in pari*/
                v[i]-=j;                    /*stupid bit reset in pari*/
                cnt--;                      /*one less candidate to play with*/
                if(p>=min_factor,
                    write(file_factors,n);  /*write the factors if we need them*/
                )
            );
            
            /*checkpoints (to see we are doing something :P)*/
            if(n%100000==0,
                printf("... Sieving with prime p=%d at n=%d, Remaining candidates=%d...      %c",p,n,cnt,13)
            )
        );
        
        /*print the checkpoints on screen, keep all of them visible*/
        print("Sieved with prime p="p", Remaining candidates="cnt"                               ");
    
        /*go to next prime*/
        p=nextprime(p+1)
    );
    write(concat(file,"_"),v);
}

smake_create_sieve_file();
smake_sieve(,10000);

/*
to resume: 
smake_sieve(10000,100000);
smake_sieve(100000,200000);
etc.
*/

[science_man_88]

Originally Posted by LaurV

Sum of digits of any 3 consecutive numbers is divisible by 3, so only 2/3 of the 1,3,5,7 (mod 10) have to be considered
This is what sm88 tries to tell you hehe, and you don't get it...

technically I'm saying 2/3rd of all indexes can be skipped because we can prove they are divisible by 3 without actually trying to divide them by it and then of the ones left over we can easily eliminate the ones that divide by 5 so we go from 100% needing to be checked to being able to code out the multiples of 3 and 5 from the start which only leaves 4 in 30 though this equates to the 1/3 elimination you talk about.

[science_man_88]

Originally Posted by danaj

The point about being sure was the "only prime indexes need to be checked" not the divisibility by 3 test.

I had to try to figure something out and going to prime indexes only seemed to be the likeliest thing that might work since it seems to work so often with these special forms of primes.

[LaurV]

Which is not true, or at least I don't understand which criteria would you use to eliminate, for example 427, or 1177. Those are indeed not prime, and their corresponding s's are not prime, but you don't know till you test. That is because s(427) is divisible by 3899827, and s(1177) by 392263, but what "rule" could you use to get these factors without testing effectively, I don't know. And I don't believe it is possible.

By analogy with mersenne numbers, 2^(n*7)-1 is not prime when n>1 because it can be written as 2^7-1 times something, and when we put all mersenne numbers in a row, we 'sieve out' all the multiplies of 7th, same as we do the sieve of Eratosthenes for integers. At the end, we are left with prime indexes (prime exponents) only. Same for 2^(n*11)-1, they are all divisible by 2^11-1. Same for 2^(n*14-1), they are all divisible by 2^14-1. Other series of numbers have more complicated relations that still allow to "sieve" them out.

But here you can not put the s(n) in a row and sieve, because except for 2 and 3, and a finite set of bigger primes that you can count on your fingers (which follow some modularity rule to 10, same as 3 does), for the other primes there is no rule. 2 will divide every even one, and 3 will divide every last two from each 3, and 5 is an exception too, because fits into "end in 0 or 5", i.e. it will divide every fifth.

But that's it. If you ran my pari script you may already know that, for example, the first indexes n for which s(n) has 7 as the smallest factor, are: 43, 109, 151, 193, 277, 319, 361 etc. Many of them are prime. And many composite indexes are missed. What's your rule for sieving them out?

And how do you get that s(427) and s(1177) are composite, without testing them? (i.e. trial division, to get the factors shown above).

Unless I am missing something, the talk about "only testing prime indexes" is false.

[ewmayer]

I've been working a bit on this sequence in complete oblivion of this thread, until sm88 made me aware of it after I posted my results in the Probability & Probabilistic Number Theory subforum. Long story short, I estimate #primes of just ~0.6 among the first million terms:

Expected number of primes in OEIS A007908

[wombatman]

I've tested terms up to n=80001, starting from n=77000, so far. No PRPs found yet. Continuing on.

[ATH]

I'm also working on this, so maybe we should start coordinating ranges? Are you sieving and then running PRP test with PFGW on the survivors?

I have sieved n=1 to 200,000 up to p=68,000,000 so far. Out of the 79,999 numbers from 3 to 199,999 which are 1,3,7,9 mod 10, I found factors for 76,924:
smarandachefactors.txt
53,332 of these factors are 3 which is exactly 2/3 for some reason, I included all those for completeness.


Here are the remaining 79999-76924 = 3075 n's up to 200,000: smarandacheremaining.txt

[ewmayer]

Originally Posted by ATH

53,332 of these factors are 3 which is exactly 2/3 for some reason, I included all those for completeness.

The divisibility-by-3 pattern is noted by Roderick MacPhee deep in the comments to the OEIS entry (thanks to sm88 for the point-out), and was one the first things I noted in my estimation of how-many-primes-do-we-expect-in-the-sequence.

BTW, if my estimate accurately captures the large-scale odds for primes in the sequence, checking probable primality (*please* only for the 4-of-every-30 terms which can possibly be prime, all else is a sheer waste of cycles) for n = 10^5 through 10^6 only boosts the odds of turning up a prime by ~10% over those for n < 10^5, from around 0.55 to around 0.60. So "you're gonna need a bigger boat", unless you get lucky.