A Theoretical (vs. Proficient/Practical) Deterministic Primality Test

[a1call]

Originally Posted by danaj

Using different software:

n=100000001, candidate = 10000000199999989

0.00000081 seconds for proof (0.81 microseconds).

0.16 seconds trial division (about 4x faster than Batalov's similar function in Pari/GP, but the point being much faster than the simple method originally shown). There are presumably even faster solutions.

Calculating the 43 million digit primorial for n takes 3.3 seconds. Calculating the 756 million digit factorial takes about 31 seconds. It gives no benefit over the primorial for this purpose (it's the same primes but with lots of duplication -- it just makes everything slower).

The GCD with the primorial is only 0.00146s.

Summarizing, trial division using a gcd is faster than stepwise divisibility if you ignore the time to generate the primorial. If you include the time for the primorial (or factorial) generation, it's much slower. Both are massively slower than a proper primality test made for numbers this size.

Yes, the gcd is fast. But calculating the primorial is quite slow. I use this to speed up small trial division in my code -- calculate the big primorials once (three of them, and not close to this large), then I can do a single gcd to effectively do trial division by a number of small primes. I found this worked well for me with GMP, but I only use it out to 40k (the first 4203 primes), and it would depend on the underlying software. I also expect to amortize the primorial calculation time over every use (at that size it isn't a big loss if only used once, but the win comes with many uses).




With pc=100000000000000019999999999999971 example.

0.014 milliseconds BPSW (probable prime)

0.470 milliseconds proof (ECPP, but ends up being just a couple simple n+1 and n-1 tests)

Trial division whether by steps, by remainder tree, by blocked remainder tree, or by gcds, is going to be very slow in comparison.


For this thread, it's too bad Pari/GP doesn't give an error when the gcd arguments are integer + real. That would make it clear that it isn't doing at all what a1call expects. But I suppose there is some use case where these argument types makes sense.

Thank you very much for the free expert analysis danaj, it is greatly appreciated.

*I consider the time to calculate the factorial/primorial irrelevant at this smallish n point, because:

** It is not my immediate bottleneck
** The Prime-Candidate has a range of n<p<n^2 per any given n!/P_n#, hence each calculation can be used for more than 1 Prime-Candidate

As per factorials vs. primorials, primorials are much smaller and of course they are faster as long as it is practical to calculate them. However they have the disadvantage of being much harder to recognize patterns in their form than factorials, which I am (most probably naively) hoping to utilize.

As per the "other" software, I installed it on ubuntu according to the command here:
http://www.mersenneforum.org/showthread.php?t=21219

but can't find any reference to how to initiate the interface. Please provide this information as well.

I have vague memories of Pari-gp style floating point errors from the Cobol days, but haven't encountered them since, but then again haven't been doing much low level programming since either (which Pari-gp is not).

Thank you very much for your input.

[CRGreathouse]

Originally Posted by a1call

oracle(n)=print("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);


print(oracle(n));

I would either replace "print(oracle(n));" with "oracle(n);" or replace "print("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);" with "Str("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);". Either of these would get rid of the extraneous printed 0.

Fun fact: It's not really a 0 but the special value gnil which is what void functions return in GP.

[a1call]

Originally Posted by CRGreathouse

I would either replace "print(oracle(n));" with "oracle(n);" or replace "print("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);" with "Str("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);". Either of these would get rid of the extraneous printed 0.

Fun fact: It's not really a 0 but the special value gnil which is what void functions return in GP.

Thank you CRGreathouse,

/* Factorial - BAT-110-C - by Rashid Naimi */

n=10000000199999979;


factorialDigits(n)=if(n>3,ceil(lngamma(n+1)/log(10)),1);\\Credits CRGreathouse 
theNumberOfFollowing0(x)=
    {
        theCount=0;
        i=1;
        while(5^i<=x,theCount+=x\(5^i);i++;);
        return(theCount);
    }
oracle(n)=print("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);


Str(oracle(n));

[science_man_88]

Originally Posted by a1call

Thank you CRGreathouse,

/* Factorial - BAT-110-C - by Rashid Naimi */

n=10000000199999979;


factorialDigits(n)=if(n>3,ceil(lngamma(n+1)/log(10)),1);\\Credits CRGreathouse 
theNumberOfFollowing0(x)=
    {
        theCount=0;
        i=1;
        while(5^i<=x,theCount+=x\(5^i);i++;);
        return(theCount);
    }
oracle(n)=print("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",factorialDigits(n)," decimal digits.";);


Str(oracle(n));

are you using 5^i because 10^i doesn't work ? I get 2*(x\10^i) !=5^i

in every 10 you have z5*even = 10* something and y0 =10*y

that's two per 10 then for every 100 you have 200*500 = 100,000 adding one so in each 10 of 10^i you have an extra pairing added that adds a 0 on the end.

[science_man_88]

Originally Posted by science_man_88

are you using 5^i because 10^i doesn't work ? I get 2*(x\10^i) !=5^i

in every 10 you have z5*even = 10* something and y0 =10*y

that's two per 10 then for every 100 you have 200*500 = 100,000 adding one so in each 10 of 10^i you have an extra pairing added that adds a 0 on the end.

doh 2^i* 5^i =10^i

[a1call]

Originally Posted by science_man_88

are you using 5^i because 10^i doesn't work ? I get 2*(x\10^i) !=5^i

in every 10 you have z5*even = 10* something and y0 =10*y

that's two per 10 then for every 100 you have 200*500 = 100,000 adding one so in each 10 of 10^i you have an extra pairing added that adds a 0 on the end.

I just looked at the pattern, it seemed to add 1 per 5 and 1 per 25, and 1 per 125...
ran a for loop comparing the code to valuation values and breaking in case of mismatch. Worked fine till 10k, figured it's correct. I was expecting a \10 too, but didn't seem to be the case.

[science_man_88]

Originally Posted by a1call

I just looked at the pattern, it seemed to add 1 per 5 and 1 per 25, and 1 per 125...
ran a for loop comparing the code to valuation values and breaking in case of mismatch. Worked fine till 10k, figured it's correct. I was expecting a \10 too, but didn't seem to be the case.

well it's not completely:

per ten the last digits are:
0-> will create a multiple of ten when the last digit of a number
1
2
3
4
5 -> multiplied by even is a multiple of 10
6
7
8
9

but then within the 10's digit that pattern repeats there are two tens digits that when multiplied add another 0 to things etc.

[a1call]

Originally Posted by science_man_88

well it's not completely:

per ten the last digits are:
0-> will create a multiple of ten when the last digit of a number
1
2
3
4
5 -> multiplied by even is a multiple of 10
6
7
8
9

but then within the 10's digit that pattern repeats there are two tens digits that when multiplied add another 0 to things etc.

Well I thought about it and this is how I see it.
It takes a 2 and a 5. Way too many 2s in any factorial, far more than 5s. Every 5 multiplied will guarantee a 0 addition, every 25 counts you multiply two 5s rather than one, so it adds two 0s and so on.
Doesn't matter how many 2's you have(you always have to many and more than you need). It takes a 5 to complete one of the 2s.

[science_man_88]

Originally Posted by a1call

Well I thought about it and this is how I see it.
It takes a 2 and a 5. Way too many 2s in any factorial, far more than 5s. Every 5 multiplied will guarantee a 0 addition, every 25 counts you multiply two 5s rather than one, so it adds two 0s and so on.
Doesn't matter how many 2's you have(you always have to many and more than you need). It takes a 5 to complete one of the 2s.

in theory using what I talked about by PM ( which I somehow missed in this thread) you'll have a possibility of adding the n-th triangular number of 0's up to the nth power of 5.

[a1call]

So, here is an attempt at getting msd Most-Significant-Digits of a factorial returned.

/* Factorial - BAT-110-E - by Rashid Naimi */

n=100001;
msd=19;\\Number of Most-Significant-Digits to return
\p 100000;
\\
\\
factorialDigits(n)=if(n>3,ceil(lngamma(n+1)/log(10)),1);\\Credits CRGreathouse 
\\
theNumberOfFollowing0(x)=
    {
        theCount=0;
        i=1;
        while(5^i<=x,theCount+=x\(5^i);i++;);
        return(theCount);
    }
\\
mostSignificantDigits(n,msd)=
    {
        b=floor(lngamma(n+1)/log(10))-msd+1;
        a=(10^(lngamma(n+1)/log(10)));
        a=a/(10^b);
        a=round(a);
        return(a);
    }
\\
oracle(n,msd)=
    {
        a=factorialDigits(n);
        print("** ",n,"! has ",theNumberOfFollowing0(n)," following 0s and ",a," decimal digits.");
        if(msd>a,msd=a);
        mostSignificantDigits(n,msd);
    }
\\
\\
theFactorial=oracle(n,msd)
#digits(theFactorial)

realprecision = 100009 significant digits (100000 digits displayed)
** 100001! has 24999 following 0s and 456579 decimal digits.
2824257650254427478
19

There doesn't seem to be a way to input a variable to the Metacommand "\p", is there?

Thank you in advance.

[CRGreathouse]

You should really write definitions for what your functions do. So

factorialDigits(n)=if(n>3,ceil(lngamma(n+1)/log(10)),1);
addhelp(factorialDigits, "factorialDigits(n): Given an integer n >= 0, return the number of decimal digits of n!. Same as #digits(n!,10) but faster and with less memory used.");

This lets you type (in this case) ?factorialDigits to get the nice help entry. But more importantly, it makes the definition of the function clear to you and others reading your code. Ideally you would explain at least (1) what the inputs are, (2) what the inputs mean, (3) what the output is, and possibly (4) what algorithms are used.

For example, it's not at all clear to me what the second argument of mostSignificantDigits is supposed to mean.

Originally Posted by a1call

mostSignificantDigits(n,msd)=
    {
        b=floor(lngamma(n+1)/log(10))-msd+1;
        a=(10^(lngamma(n+1)/log(10)));
        a=a/(10^b);
        a=round(a);
        return(a);
    }

This could be

mostSignificantDigits(n,msd)=
{
  my(a,b,t); \\ always define your variables using my() or bad things can happen
  t = lngamma(n+1)/log(10); \\ Why compute this twice?
  b = floor(t) - msd + 1;
  a = 10^t;
  round(a / 10^b); \\ no need for return() at the end of a function
}
addhelp(mostSignificantDigits, "mostSignificantDigits(n,msd): Do something or other with n and msd, then return something else.");

Come to think of it, you could simply do

mostSignificantDigits(n,msd)=
{
  my(b,t);
  t = lngamma(n+1)/log(10);
  b = floor(t) - msd + 1;
  round(10^(t-b));
}

and... hmm... t - floor(t) is just frac(t), so

mostSignificantDigits(n,msd)=
{
  my(t = lngamma(n+1)/log(10));
  round(10^(frac(t) - msd + 1));
}

I still have no idea what it's doing, though.

Originally Posted by a1call

There doesn't seem to be a way to input a variable to the Metacommand "\p", is there?

Use

default(realprecision, x)