mersenneforum.org Some random calculations on (MM19+2)/3
 Register FAQ Search Today's Posts Mark Forums Read

 2016-09-27, 20:34 #1 Raman Noodles     "Mr. Tuch" Dec 2007 Chennai, India 3·419 Posts Some random calculations on (MM19+2)/3 (2524287+1)/3 is a composite number with no known prime factors. I had some doubt over it when factordb.com web site page marked it as 'Status Unknown', but I myself wrote out a Modular Exponentiation function, not a function, but rather a one single line script in PARI/GP tiny enough to check it out! Ten years back, in the year 2006, when I first implemented Modular Exponentiation, I wrote a big Java function iteratively. It is also possible to do it recursively, but modPow function is readily available already and implemented in built in Java. For learning purposes only. Code: static java.math.BigInteger modPow(BigInteger base,BigInteger power,BigInteger modulo) { String binary=new String(); BigInteger copy=power; while(copy.compareTo(BigInteger.valueOf(0))>0) { if((copy.remainder(BigInteger.valueOf(2))).compareTo(BigInteger.valueOf(0))==0) { binary=binary+"0"; } if((copy.remainder(BigInteger.valueOf(2))).compareTo(BigInteger.valueOf(0))==1) { binary=binary+"1"; } copy=copy.divide(BigInteger.valueOf(2)); } BigInteger multiplier=base; int position=0; BigInteger result=BigInteger.valueOf(1); while(position0,x=x*j);y=y*2;j=j^2);x For recursively doing Modular Exponentiation in PARI/GP, it throws away a nesting error of stack overflow exception. PARI/GP does not have an in built Modular Exponentiation function already available or readily implemented, but indeed we can use up with the Mod(x, y) operator. Code: pow(a,b)=print(log(b)/log(10));if(b==1,a,if(b%2==0,pow(a^2,b/2),a*pow(a^2,(b-1)/2))) N=(2^8191+1)/3 lift(Mod(3,N)^(N-1)) lift(pow(Mod(3,N),N-1)) *** deep recursion: print(log(b)/log(10));if(b ^-------------------- N=(2^524287+1)/3 lift(Mod(3,N)^(N-1)) lift(pow(Mod(3,N),N-1)) *** if: the PARI stack overflows ! current stack size: 4000000 (3.815 Mbytes) [hint] you can increase GP stack with allocatemem() For a PRP test of (2524287+1)/3, my script executed in 4 hours and 15 minutes. But Mod(3, N)^(N-1) where N = (2524287+1)/3 executed in 1 hour and 45 minutes which amounts to be about 142.86% as efficient than when compared to my own implementation. What causes this speed up even when PARI/GP programming code is only used up? Code: N=(2^524287+1)/3 lift(Mod(3,N)^(N-1)) lift(power(Mod(3,N),N-1))
2016-09-28, 03:44   #2
LaurV
Romulan Interpreter

"name field"
Jun 2011
Thailand

17·19·31 Posts

Quote:
 Originally Posted by Raman For a PRP test of (2524287+1)/3, my script executed in 4 hours and 15 minutes. But Mod(3, N)^(N-1) where N = (2524287+1)/3 executed in 1 hour and 45 minutes which amounts to be about 142.86% as efficient than when compared to my own implementation. What causes this speed up even when PARI/GP programming code is only used up?
Well, at last, an interesting post from you..
(actually second today, after the children's puzzle, which I actually enjoyed solving it, even at my old age).

One reason why your function is slower is the printing itself. Printing does strange things with the command prompt window, like scrolling it, etc, which take CPU time, because they are executed in the same context thread like of the cmd prompt. Every time you print you lose precious milliseconds, which I would hazard a guess, it takes more than 50% of the time. .

You can test this very easy, by editing the window properties and make the buffer 3000 lines tall, instead of the default 300 or so (right click on the window's title, chose properties, layout, screen buffer size height). Your function will be even slower.

A solution would be to display the progress only every 100, or 1000 iterations or so (why do you think P95 does not display the progress every iteration?).

Something like:
Code:
gp > power(a,b)=x=1;y=1;cnt=0;j=a;z=floor(log(b)/log(2))+1;for(i=1,z,if(cnt++>=100, cnt=0; printf("...Iteration: %d / %d (%7.3f%% )...%c",i,z,100.0*i/z,13));if(bitand(b,y)>0,x=x*j);y=y*2;j=j^2);x
This avoids scrolling the screen, by printing all on the same line, easier to read too, as the useful text does not scroll away, and also, to increment the counter is faster than printing, and faster than modular i%100 (in case you don't want to use additional counter variable), and it does printing every 100 iterations.

Other improvements you could do: pari/gp "bitand" (and all bitwise functions) are terrible slow. Also, using modular variable is slower than using integers. If you compute every iteration, you should make the calculus in integers, and call the function with integer parameters. Like using the % operator after every squaring/rotation. Oh, rotation, by the way, is faster than multiplication, so you can also gain some time by writing (integer!) operations like y<<=1 instead of y=y*2. In fact, y*=2 is faster than y=y*2 (because it evaluate y only once), and y=y<<1 is faster then both (because it only shifts the integer to the left, there is no calculus done), but y<<=1 is much faster than both. Etc. Then, if you work in integers, you also avoid the lifting operation every time.

Edit: even if the "floor(log blah blah)" is executed only once, this is a tricky operation, whose result depends on the default precision. You would be better off (and faster, and more accurate) with "z=#binary(b)".

Last fiddled with by LaurV on 2016-09-28 at 04:01

 2016-09-28, 05:57 #3 LaurV Romulan Interpreter     "name field" Jun 2011 Thailand 234358 Posts Ha, replying to my own post... man, see what you are doing to me... I got this disease from you... I put this together during my lunch break, which is a bit faster that the "default" one and also shows the progress, it does not uses booleans, but only shifts and subtractions. Now the only thing you need to make is to make the powering step faster (that is still dealing with double number of digits, and pari does not implement FFT). (of course, that was a joke; compiled stuff will always be much faster that this!) Code: /* it raises 'a^b (mod m)' showing the progress every 'step' iterations */ modpow(a=2,b=10260,m=10261,step=10)= { if(b<0, return(0)); if(b==0, return(1)); if(b==1, return(a%m)); my(size=#binary(b)); my(mask=1<<(size-1)); my(cnt=0); my(pow=a); my(iter=0); b-=mask; while(mask>1, /* pow*=pow; */ /* pow%=m; */ pow=(pow*pow)%m; mask>>=1; if(mask<=b, b-=mask; /* pow*=a; */ /* pow%=m */ pow=(pow*a)%m ); if(cnt++>=step, cnt=0; printf("...Iteration: %d / %d (%7.3f%% )...%c",iter+=step,size,100.0*iter/size,13) ) ); return(pow); } For the same definition of N that you have, here is the timing: Code: gp > Mod(3,N)^1000000000; time = 564 ms. gp > Mod(3,N)^1000000000; time = 498 ms. gp > Mod(3,N)^1000000000; time = 514 ms. gp > Mod(3,N)^1000000000; time = 503 ms. gp > modpow(3,1000000000,N,100); time = 523 ms. gp > modpow(3,1000000000,N,100); time = 510 ms. gp > modpow(3,1000000000,N,100); time = 491 ms. gp > modpow(3,1000000000,N,100); time = 471 ms. gp > modpow(3,1000000000,N,100); time = 502 ms. (the laptop is doing other things in the same time, which can not be stopped, therefore the multiple runs, with different times). Edit: comparison: (this is the "fast" one which does not print every iteration) Code: gp > lift(power(Mod(3,N),1000000000)); time = 971 ms. gp > lift(power(Mod(3,N),1000000000)); time = 923 ms. gp > lift(power(Mod(3,N),1000000000)); time = 906 ms. gp > lift(power(Mod(3,N),1000000000)); time = 1006 ms. gp > lift(power(Mod(3,N),1000000000)); time = 965 ms. gp > lift(power(Mod(3,N),1000000000)); time = 955 ms. gp > Last fiddled with by LaurV on 2016-09-28 at 06:08
2016-09-28, 07:03   #4
CRGreathouse

Aug 2006

10111010110112 Posts

Quote:
 Originally Posted by Raman Code: power(a,b)=x=1;y=1;j=a;z=floor(log(b)/log(2))+1;for(i=1,z,print("Iteration: "i" / "z);if(bitand(b,y)>0,x=x*j);y=y*2;j=j^2);x
Some quick notes:
• It is better to use logint(b,2) rather than floor(log(b)/log(2)), as it avoids rounding issues. For example, on my 64-bit machine with default precision (38 decimal places), power(Mod(2,7), 2^33) gives the wrong result.
• Similarly, prefer x*=j and y*=2 (or even y<<=1) to x=x*j and y=y*2. These avoid aliasing issues* and are easier to read.
• For readability, functions like this should be written over several lines.
• You should always declare scratch variables with one of my(), local(), or inline(). If you don't know which you need, my() is the one you want.

Together with LaurV's suggestion on printing (which is, in most languages, horribly slow) you get

Code:
power(a,b)=
{
for(i=1,z,
if(cbnt++>99,
print("Iteration: "i" / "z);
cnt=0
);
if(bitand(b,y)>0, x*=j);
y <<= 1;
j=j^2
);
x;
}
I find this much easier to read (and you probably do, too). It also lets you take stock of what's going on. For example, you're walking through the bits of b by making a large variable, y, which has only one set bit. It would be easier if you could access the words individually so you didn't have to keep such a large number around (while needlessly computing ANDs across so many bits). And once you do this, it might be more natural to compute one word of the exponent before checking to see if you should print, saving 63 additions and 63 comparisons (with 64-bit words).

The next natural step would be to take this cleaned-up code and pass it to gp2c, where it would generate C source. This would be a faster version of the code, which you could then edit to make it faster still.

Under the hood, gp does modular exponentiation in the function Fp_pow, which is called by powgi which is called from gpow which is essentially the GP ^ operator.

* This makes them somewhat faster. For example, t=7;for(i=1,1e7,t*=1) takes 437 ms vs. 1296 ms for t=7;for(i=1,1e7,t=t*1) on my machine.

2016-09-28, 07:37   #5
axn

Jun 2003

5,387 Posts

Quote:
 Originally Posted by CRGreathouse * This makes them somewhat faster. For example, t=7;for(i=1,1e7,t*=1) takes 437 ms vs. 1296 ms for t=7;for(i=1,1e7,t=t*1) on my machine.
Code:
GP/PARI CALCULATOR Version 2.7.3 (development 16369-68ebfde)
i686 running mingw (ix86/GMP-5.1.3 kernel) 32-bit version
compiled: Jan  4 2015, gcc version 4.9.1 (GCC)
(readline v6.2 enabled, extended help enabled)
Code:
? #
timer = 1 (on)
?  t=7;for(i=1,1e7,t*=1)
time = 1,451 ms.
?  t=7;for(i=1,1e7,t=t*1)
time = 1,357 ms.
?  t=7;for(i=1,1e7,t*=1)
time = 1,467 ms.
?  t=7;for(i=1,1e7,t=t*1)
time = 1,342 ms.

EDIT:- Nope
Code:
          GP/PARI CALCULATOR Version 2.8.0 (alpha)
amd64 running mingw (x86-64/GMP-6.0.0 kernel) 64-bit version
compiled: Aug  1 2016, gcc version 5.4.0 20160609 (GCC)
(readline v6.2 enabled, extended help enabled)
Code:
? #
timer = 1 (on)
? t=7;for(i=1,1e7,t*=1)
time = 1,357 ms.
? t=7;for(i=1,1e7,t=t*1)
time = 1,279 ms.
? t=7;for(i=1,1e7,t=t*1)
time = 1,326 ms.
? t=7;for(i=1,1e7,t*=1)
time = 1,405 ms.

Last fiddled with by axn on 2016-09-28 at 07:42

 2016-09-28, 22:11 #6 danaj   "Dana Jacobsen" Feb 2011 Bangkok, TH 32×101 Posts My brief experiment, using the n=(2^524287+1)/3 example: 31m Perl/ntheory powmod(3,n-1,n) (GMP's mpz_powm) 33m Pari/GP Mod(3,n)^(n-1) 66m lift(power(Mod(3,n),n-1) with CRG's power and cnt++ > 999. For the first, we could also just call is_pseudoprime($n,3). The first two are close enough to be equal given the simple test. I'm not sure what the timings are for others with this function, but it doesn't look good here. 2016-09-29, 04:39 #7 R. Gerbicz "Robert Gerbicz" Oct 2005 Hungary 62516 Posts Quote:  Originally Posted by danaj My brief experiment, using the n=(2^524287+1)/3 example: 31m Perl/ntheory powmod(3,n-1,n) (GMP's mpz_powm) 33m Pari/GP Mod(3,n)^(n-1) 66m lift(power(Mod(3,n),n-1) with CRG's power and cnt++ > 999. And those are still nowhere from a decent prp implementation for these type of numbers. Modular reduction by 3*N=2^e+1 is easy using bit arithmetic, and the bit pattern of (2^e+1)/3 is trivial. But just to make the code cleaner we calculate 3^(3*N-1) mod N. Here is my code in Pari-Gp: (using only two divisions, but you can easily eliminate those also) Code: H(e)={if(e%2==0,return(0)); n=2^e+1;m=n-2;r=3%n;for(i=1,e,r=sqr(r);hi=shift(r,-e);lo=bitand(r,m); r=lo-hi;if(r<0,r+=n));return((r-9)%(n/3)==0)} your one liner was: Code: F(e)=n=(2^e+1)/3;return(Mod(3,n)^(n-1)==Mod(1,n)) and the tests: Code: ? H(524287) %12 = 0 ? ## *** last result computed in 8min, 25,077 ms. ? F(524287) %13 = 0 ? ## *** last result computed in 31min, 27,651 ms. 2016-09-29, 10:10 #8 Raman Noodles "Mr. Tuch" Dec 2007 Chennai, India 125710 Posts My submitted post has been moved here from its original place of Cunningham tables 2+ table discussion thread? It is always a good thing to think about ideas for some time, ponder around some thoughts, before submitting it for thread posting purposes and give away some time over there by providing room to wards for an other people's inputs! My new blurb under my mersenne forum user name has been given as Noodles in white colour font text? Is it really appropriate to be called as blurb? Or mersenne forum position, as some one else would call it in to an other bulletin boards. Post counts in mersenne forum show as majority of times as prime factorizations and minority of times as showing number representation in different bases. The latter only seems to use appearing it in to base 2, 8, 10, 16 in mersenne forum. These are really only the commonly used bases even in real life. Not the other bases at all? Should also consider with an other novel representations such as sum of two squares representations and an other normal representations such as with normal number representations with out any changes in to number text representations at all. A variety of representation, each time showing as different things would be extremely nicer! Coming back to the subject, It would be good for me to have a PARI/GP script for calculating ECM curve's group order. Here is the MAGMA script. Code: // courtesy of Allan Steel // first save this in a file (say go.magma) // then start magma and type: // load "go.magma"; // then type: // FindGroupOrder(p, sigma); // where p is the found factor (must be prime) // and sigma is the curve parameter (Suyama's parametrization) FindGroupOrder := function (p, sigma) K := GF(p); v := K ! (4*sigma); u := K ! (sigma^2-5); x := u^3; b := 4*x*v; a := (v-u)^3*(3*u+v); A := a/b-2; x := x/v^3; b := x^3 + A*x^2 + x; E := EllipticCurve([0,b*A,0,b^2,0]); return FactoredOrder(E); end function; To convert it in to PARI/GP script from MAGMA script, I need to know a few things. What is ! operator in to MAGMA language? Are there any equivalents of EllipticCurve([0,b*A,0,b^2,0]) and GF(p) in to PARI/GP language? Quote:  Originally Posted by GP2 If you google the number 493613348917766417426006591803225943649556875046256846631045789294232301761, you find this old thread: "Faking factors with Complex Multiplication", and within that thread post #4 discusses that number. It does seem it was found artificially, if I understand that post correctly. It is very trivial to construct extremely smooth prime factors which are very larger for p-1 and p+1. But for ECM curves, constructing extremely smooth prime factors with a chosen σ parameter of the value will be a real challenge! By the way observe and remember that stage 1 of ECM curves or p-1 or p+1 will only always take up with LCM of all numbers below B1 bound, till then. Quote:  Originally Posted by Raman How about posting all of my opinions over here, rather than in their appropriate threads? It is possible for us to predict solar eclipses and lunar eclipses, Mercury transits and Venus transits, but can be ever able to predict earthquakes? Quote:  Originally Posted by LaurV Well, at last, an interesting post from you.. (actually second today, after the children's puzzle, which I actually enjoyed solving it, even at my old age). How about two quotes here? And many quotes from different places? Every quote tag can choose to have a pointer tag or not, it is always up to the thread poster, but only never for the private messages quote tags and code tags! To post threads replies and things immediately after wards of the thinking processes and thought processes! Quote:  Originally Posted by LaurV This is by far the best post of Raman in the last year... Last year? During the last year, I had only been active in this mersenne forum for the past one month. To wards this direction, after about Monday 17 March 2014, I had stayed away from this mersenne forum for nearly about 2½ years (all most about a 'once in a blue moon'!) till about Friday 26 August 2016, except for two posts up on Monday 23 May 2016. (It would also be a some good thing if mersenne forum thread header for every post would too display day of week along with date and time in to the title bar and IP address does get tagged with every thread post invisible in to the normal users of mersenne forum along with in to the ISP provider address). Any way that I will be celebrating with my ten years of spent time in to the mersenne forum, by using some time sooner right now, by using including my older mersenne forum ID, involving my extended time periods of inactivity! Quote:  Originally Posted by Raman How can I insert line break in to the mersenne forum thread post? Quote tags and besides of code tags can be used as point separators and besides of line breaks. People living in the border of the South Indian state of the Tamiɻ Nadu only know two languages namely English and Tamiɻ while on the other hand, people living in the border of the South Indian state of the Andhra Pradesh know four languages namely English, Tamiɻ, Telugu and Hindi. No Jawahar Navodaya Vidyalaya schools are opened up at all ultimately with in to the state of Tamiɻ Nadu! If Tamiɻ people respect Hindi people, then Hindi people will respect Tamiɻ people. Just as always simple and easy rule of life! New languages can only be picked up as a child, needs a fresh mind and family and social environment and interactions, unlikely as an adult. But, that in to human body, may be that languages is in built as a child, by using parent's genes or some thing like that type of thing. By the way that why is 18 years of oldness is considered in to be age of majority with majority of countries and nations? But where as up on the other hand, learning many an other programming languages from one single programming language is always easier, and also interesting too, most of them have been got with similar and user friendly syntax in to the scripting language, except only when programming code syntax differs significantly and need to think about novel different and strange ways to write out programming code with in to the exotic scripting language. Consider with in to the population growth rate during the last past few years ultimately and untimely and the probability of an other people being born with in to the last past few years ultimately and untimely! Population decline and later marriages and weddings on wards, too much of the obsession with in to the higher studies, etc., and so on, things as these, also exists out right at the country and nation of the Japan and Brazil at all. Right? Just as always, things as these. Quote:  Originally Posted by Raman How can I insert line break in to the mersenne forum thread post? In the factordb.com web site page, we can be able to report prime factors of composite numbers. But that is it possible for us to mark up now known and analyzed numbers like those such as, as (2524287+1)/3, etc., and so on, as composite, which are right now marked up like those such as, as 'Status Unknown'? It also accepts too primality certificates for some type of tests? Does it? For what types of tests? But the largest known prime numbers like those such as, as the Mersenne prime numbers exponents are marked as 'Status Unknown'. It does not accept the Lucas Lehmer Test residue as a primality certificate. Does it? It also automatically too runs and executes away workers out for the smallest unfactored composite numbers that are added up by using an other people regularly and routinely from time to time and from now and then, right now. Does it? Also smallest probable prime numbers (smallest PRP numbers) and smallest 'Status Unknown' numbers too. Does it? Three 32 digit prime factors for computation from Aliquot Sequence of 4788 Iteration Number 10662! Higher powers of 2 and an other numbers would not last off standing out for long time periods with in to the computation from Aliquot Sequences Iteration Numbers! Mersenne Forum number of users registered: 14968 totally over all! Mersenne Forum number of attachments IDs: 14969 totally over all! Last fiddled with by Raman on 2016-09-29 at 10:41 2016-09-29, 12:08 #9 LaurV Romulan Interpreter "name field" Jun 2011 Thailand 17×19×31 Posts Quote:  Originally Posted by danaj My brief experiment, using the n=(2^524287+1)/3 example: 31m Perl/ntheory powmod(3,n-1,n) (GMP's mpz_powm) 33m Pari/GP Mod(3,n)^(n-1) 66m lift(power(Mod(3,n),n-1) with CRG's power and cnt++ > 999. For the first, we could also just call is_pseudoprime($n,3). The first two are close enough to be equal given the simple test. I'm not sure what the timings are for others with this function, but it doesn't look good here.
If I don't dare much, would you like to run my code also? I am curious how it compares with the Mod, on the same computer where you ran the other stuff.
You would need to run modpow(3,n-1,n,100) or so.

Last fiddled with by LaurV on 2016-09-29 at 12:10

2016-09-29, 17:58   #10
Batalov

"Serge"
Mar 2008
Phi(4,2^7658614+1)/2

9,901 Posts

Quote:
 Originally Posted by R. Gerbicz And those are still nowhere from a decent prp implementation for these type of numbers.
Exactly right!

The 'golden standard' running time is ~53 seconds (with mprime, of course, and on 1 thread to make playing field even; it is even faster with mprime's || code).
Splitting hairs between 31 or 33 minutes is pretty much the same on this benchmark, but 8 minutes is a good demonstration by Robert. The exponentiations need to run mod(2^p+1) with one final mod((2^p+1)/3); if there is a trick in Pari to recognize special mod(2^p+1) and pass it to GMP, I haven't found it.

But for the naked libgmp test, there should be a way to implement this, no?

2016-09-29, 18:35   #11
danaj

"Dana Jacobsen"
Feb 2011
Bangkok, TH

38D16 Posts

I debated mentioning using gwnum, which would be much faster even just doing a straightforward modular power similar to the others. Robert of course points out that applying some thought can reduce the brute force needs.

Quote:
 Originally Posted by Batalov ]But for the naked libgmp test, there should be a way to implement this, no?
Do you mean applying Robert's algorithm in GMP? I could try that on the same machine (but might be a few days before I have time, someone else could do before).

 Similar Threads Thread Thread Starter Forum Replies Last Post jasong Lounge 46 2017-05-09 12:32 Sam Kennedy Programming 2 2012-12-27 08:52 Greenk12 Factoring 1 2008-11-15 13:56 Googol Information & Answers 2 2007-07-27 03:09 alienz Lounge 19 2003-07-06 03:11

All times are UTC. The time now is 10:12.

Sat Aug 13 10:12:11 UTC 2022 up 37 days, 4:59, 2 users, load averages: 0.95, 1.02, 1.04