20160927, 20:34  #1 
Noodles
"Mr. Tuch"
Dec 2007
Chennai, India
3·419 Posts 
Some random calculations on (MM19+2)/3
(2^{524287}+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(position<binary.length()) { if(binary.charAt(position)=='1') { result=result.multiply(multiplier); result=result.remainder(modulo); } multiplier=multiplier.multiply(multiplier); multiplier=multiplier.remainder(modulo); position=position+1; } return result; } Code:
static java.math.BigInteger modPow(BigInteger base,BigInteger power,BigInteger modulo) { if(power.remainder(BigInteger.valueOf(2)).compareTo(BigInteger.valueOf(0))==0) { return modPow((base.multiply(base)).remainder(modulo),power.divide(BigInteger.valueOf(2)),modulo); } else { return (base.multiply(modPow((base.multiply(base)).remainder(modulo),(power.subtract(BigInteger.valueOf(1))).divide(BigInteger.valueOf(2)),modulo))).remainder(modulo); } } certain operation and we can be able to power it. But, it shows no progress for bigger numbers. I included showing the progress in my written one line script in PARI/GP. 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 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,(b1)/2))) N=(2^8191+1)/3 lift(Mod(3,N)^(N1)) lift(pow(Mod(3,N),N1)) *** deep recursion: print(log(b)/log(10));if(b ^ N=(2^524287+1)/3 lift(Mod(3,N)^(N1)) lift(pow(Mod(3,N),N1)) *** if: the PARI stack overflows ! current stack size: 4000000 (3.815 Mbytes) [hint] you can increase GP stack with allocatemem() 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)^(N1)) lift(power(Mod(3,N),N1)) 
20160928, 03:44  #2  
Romulan Interpreter
"name field"
Jun 2011
Thailand
17·19·31 Posts 
Quote:
(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 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 20160928 at 04:01 

20160928, 05:57  #3 
Romulan Interpreter
"name field"
Jun 2011
Thailand
23435_{8} 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<<(size1)); 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); } 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. 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 20160928 at 06:08 
20160928, 07:03  #4  
Aug 2006
1011101011011_{2} Posts 
Quote:
Together with LaurV's suggestion on printing (which is, in most languages, horribly slow) you get Code:
power(a,b)= { my(x=1,y=1,j=a,z=logint(b,2)+1,cnt); 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; } The next natural step would be to take this cleanedup 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. 

20160928, 07:37  #5  
Jun 2003
5,387 Posts 
Quote:
Code:
GP/PARI CALCULATOR Version 2.7.3 (development 1636968ebfde) i686 running mingw (ix86/GMP5.1.3 kernel) 32bit version compiled: Jan 4 2015, gcc version 4.9.1 (GCC) threading engine: single (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 (x8664/GMP6.0.0 kernel) 64bit version compiled: Aug 1 2016, gcc version 5.4.0 20160609 (GCC) threading engine: single (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 20160928 at 07:42 

20160928, 22:11  #6 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
3^{2}×101 Posts 
My brief experiment, using the n=(2^524287+1)/3 example:

20160929, 04:39  #7  
"Robert Gerbicz"
Oct 2005
Hungary
625_{16} Posts 
Quote:
and the bit pattern of (2^e+1)/3 is trivial. But just to make the code cleaner we calculate 3^(3*N1) mod N. Here is my code in PariGp: (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=n2;r=3%n;for(i=1,e,r=sqr(r);hi=shift(r,e);lo=bitand(r,m); r=lohi;if(r<0,r+=n));return((r9)%(n/3)==0)} Code:
F(e)=n=(2^e+1)/3;return(Mod(3,n)^(n1)==Mod(1,n)) 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. 

20160929, 10:10  #8  
Noodles
"Mr. Tuch"
Dec 2007
Chennai, India
1257_{10} 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^25); x := u^3; b := 4*x*v; a := (vu)^3*(3*u+v); A := a/b2; 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; 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:
By the way observe and remember that stage 1 of ECM curves or p1 or p+1 will only always take up with LCM of all numbers below B1 bound, till then. Quote:
Quote:
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! 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 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. 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 (2^{524287}+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 20160929 at 10:41 

20160929, 12:08  #9  
Romulan Interpreter
"name field"
Jun 2011
Thailand
17×19×31 Posts 
Quote:
You would need to run modpow(3,n1,n,100) or so. Last fiddled with by LaurV on 20160929 at 12:10 

20160929, 17:58  #10  
"Serge"
Mar 2008
Phi(4,2^7658614+1)/2
9,901 Posts 
Quote:
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? 

20160929, 18:35  #11 
"Dana Jacobsen"
Feb 2011
Bangkok, TH
38D_{16} 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.
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). 
Thread Tools  
Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
random comments, random questions and thread titles made for Google  jasong  Lounge  46  20170509 12:32 
Could someone please check my calculations?  Sam Kennedy  Programming  2  20121227 08:52 
About random number (random seed) in Msieve  Greenk12  Factoring  1  20081115 13:56 
Is based 10 used for calculations?  Googol  Information & Answers  2  20070727 03:09 
HOW MANY CALCULATIONS AT ONCE??  alienz  Lounge  19  20030706 03:11 