[QUOTE=axn;238387][CODE]lucaslehmer(p)= s=Mod(4, 2^p-1);for(x=1,p-2,s=s^2-2);if(s==0,1,0);[/CODE]
Some small changes.[/QUOTE]

thanks axn you're pretty nice. I have a second implementation to catch all p that have sumdigits(last result)==9 now, as well as one for all #'s that have sumdigits(last result)==9.Neither sequence is in the OEIS but I haven't found much of a pattern yet either so they may not be of interest. I may check back soon if i can find a pattern.

That wasn't what I did before,as far as I can tell. but it is useful in my learning.Once again I spoke a bit too soon,

[CODE]lucaslehmer2(p)= s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2 && sumdigits(s)==9,print1(p","))[/CODE]

axn what are the possible improvements on this? This leads me to more curiosity than the first one I made. This one can help show all the exponents that have a sumdigits(last term)==9 not just the 2^p-1 is prime ones. If we can predict this one, maybe we can find a pattern in the p that make 2^p-1 prime (I'm guessing based on position in this supersuquence).

[QUOTE=science_man_88;238413][CODE]lucaslehmer2(p)= s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2 && sumdigits(s)==9,print1(p","))[/CODE]

axn what are the possible improvements on this? This leads me to more curiosity than the first one I made. This one can help show all the exponents that have a sumdigits(last term)==9 not just the 2^p-1 is prime ones. If we can predict this one, maybe we can find a pattern in the p that make 2^p-1 prime (I'm guessing based on position in this supersuquence).[/QUOTE]

Good luck with that..

This is too vague to work with.

Perhaps, if you explained this a bit more clearly, I could either help or shoot this "conjecture" to hell, where it likely belongs.

[QUOTE=3.14159;238416]Good luck with that..

This is too vague to work with.

Perhaps, if you explained this a bit more clearly, I could either help or shoot this "conjecture" to hell, where it likely belongs.[/QUOTE]

Put a for loop around it not much to explain!

[QUOTE=science_man_88;238417]Put a for loop around it not much to explain![/QUOTE]

What you posted is too vague. I either need math that explains what you posted, or a clear explanation of whatever it is that you posted, with minimal ambiguity.

[QUOTE=3.14159;238420]What you posted is too vague. I either need math that explains what you posted, or a clear explanation of whatever it is that you posted, with minimal ambiguity.[/QUOTE]

It's a lucas-lehmer test as denoted by the title. The main difference between this and the first one is the fact that the test checks to see if sumdigits(last result)==9.

[QUOTE=science_man_88;238421]It's a lucas-lehmer test as denoted by the title. The main difference between this and the first one is the fact that the test checks to see if sumdigits(last result)==9.[/QUOTE]

.. You could just download an implementation of it..

And, why would the latter part (Digits of.. what?, summing to 9), be necessary?

[QUOTE=3.14159;238423].. You could just download an implementation of it..[/QUOTE]

Why? By following the basic pseudocode on wikipedia I was able to create my own basic Pari script. I'll admit I haven't made the most efficient script but I'm happy with it, Also why do you feel you need to throw me down? To boost yourself up? If so you're sinking like a rock.

[QUOTE=3.14159;238423].. You could just download an implementation of it..

And, why would the latter part (Digits of.. what?, summing to 9), be necessary?[/QUOTE]

sumdigits(last result)==9, and the fact it's labelled lucaslehmer2(p) didn't make it obvious ? if not I feel very sad for you.

[QUOTE=science_man_88;238424]Why? By following the basic pseudocode on wikipedia I was able to create my own basic Pari script. I'll admit I haven't made the most efficient script but I'm happy with it, Also why do you feel you need to throw me down? To boost yourself up? If so you're sinking like a rock.[/QUOTE]

I dunno. A downloaded implementation is normally faster, for any serious searcher.

Throw you down? When? How? By pointing out that you're seeing "patterns", where there are none? There were no putdowns. It was an observation of mine, and nothing else.

A putdown would be something along the lines of, "You're wrong, because you're a damn moron who doesn't have any idea what they're talking about."

[QUOTE=3.14159;238416]

Perhaps, if you explained this a bit more clearly, I could either help or [B][I][U][SIZE="7"]shoot this "conjecture" to hell[/SIZE][/U][/I][/B], where it likely belongs.[/QUOTE]

You know nothing of a put down or even an attempt there, can't see one for the life of me.

[QUOTE=science_man_88;238432]You know nothing of a put down or even an attempt there, can't see one for the life of me.[/QUOTE]

Sorry, guy, no putdowns there. No matter how you may try perverting what I posted, there were no putdowns.

Now, let's prevent this thread from further derailing, and get back on track.

The original request still holds: Some math, or a clear explanation without ambiguity.

[QUOTE=3.14159;238433]Sorry, guy, no putdowns there. No matter how you may try perverting what I posted, there were no putdowns.

Now, let's prevent this thread from further derailing, and get back on track.

The original request still holds: Some math, or a clear explanation without ambiguity.[/QUOTE]

you lie through your teeth. on the other note I've already stated sumdigits(last result)==9, and the lucaslehmer title explain it all, also i never asked you, I asked the one who seemed not to put me down.

[QUOTE=science_man_88;238435]you lie through your teeth. on the other note I've already stated sumdigits(last result)==9, and the lucaslehmer title explain it all, also i never asked you I asked the one who seemed no to put me down.[/QUOTE]

[spoiler]Lied? When? This is too much fun!

Let's look at the definition of "lie";

[QUOTE=Dictionary][spoiler](Noun) [B]A false statement made with deliberate intent to deceive; an intentional untruth; a falsehood.[/B]
(Verb) [B]To speak falsely or utter untruth knowingly, as with intent to deceive.[/spoiler][/B]
[/QUOTE]

.. No lies told, no putdowns made. Hurl mud to your heart's desire, but it does no good.[/spoiler]

And, now: What does "sumdigits" do?

You claim to shoot me down is not a put down that's is false.If you can and are willing to shoot them down you think you're smarter. This is like calling them dumb, which is pretty much saying you moron.

"You're wrong, because you're a damn moron who doesn't have any idea what they're talking about."


Moron appears in both paths they are equal in this case.

[QUOTE=3.14159;238438]

And, now: What does "sumdigits" do?[/QUOTE]

if you don't remember look in this thread I made it here.

[QUOTE]you claim to shoot me down is not a put down that's is false.If you can and are willing to shoot them down you think you're smarter. This is like calling them dumb, which is pretty much saying you moron.
[/QUOTE]

I attacked your ideas. You are not your ideas. There was no reason to take this as a "personal attack". There were no putdowns. Like I said some time ago, if you insist on singing that there are, where there are none, you are free to do so.

[QUOTE]if you don't remember look in this thread I made it here.
[/QUOTE]

Where at?

[QUOTE=science_man_88;226878]I also did this with your knowledge:

[CODE]sumdigits(n) = if(n%9!=0,return(n%9),return(n%9+9))[/CODE]

but you might think it's kinda lame.[/QUOTE]

post 882 I believe.

[QUOTE=science_man_88;238443]post 882 I believe.[/QUOTE]

This can be simplified to x mod 9.

[QUOTE=3.14159;238444]This can be simplified to x mod 9.[/QUOTE]

No it can't as 0 becomes 9 so you must adapt it.

[QUOTE=science_man_88;238413][CODE]lucaslehmer2(p)= s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2 && sumdigits(s)==9,print1(p","))[/CODE]

axn what are the possible improvements on this? This leads me to more curiosity than the first one I made. This one can help show all the exponents that have a sumdigits(last term)==9 not just the 2^p-1 is prime ones. If we can predict this one, maybe we can find a pattern in the p that make 2^p-1 prime (I'm guessing based on position in this supersuquence).[/QUOTE]

You write
[CODE]if(x=p-2 && ...[/CODE]
which means "set x to the value p-2, then if that value is nonzero continue to the next condition". This probably isn't what you intended.

Also, what is sumdigits? [url=http://oeis.org/A010888]A010888[/url]? [url=http://oeis.org/A007953]A007953[/url]? Something else?

[QUOTE=CRGreathouse;238496]Also, what is sumdigits? [url=http://oeis.org/A010888]A010888[/url]? [url=http://oeis.org/A007953]A007953[/url]? Something else?[/QUOTE]

sumdigits is a function in the thread already pointed back to by Pi's questioning.

[QUOTE=science_man_88;238497]sumdigits is a function in the thread already pointed back to by Pi's questioning.[/QUOTE]

So, neither. It's like the digital root but with sumdigits(0) = 9.

Faster code:
[code]sumdigits(n)=n=n%9;if(n,n,9)[/code]

[QUOTE=CRGreathouse;238501]So, neither. It's like the digital root but with sumdigits(0) = 9.

Faster code:
[code]sumdigits(n)=n=n%9;if(n,n,9)[/code][/QUOTE]

all i have to do is take out the return() on each to cut the time by 66% or so. and yours didn't beat mine by much on my machine ( I'll double check with multiple runs)

[CODE](12:13) gp > sumdigits1(n)=if(n%9!=0,n%9,9)
%5 = (n)->if(n%9!=0,n%9,9)
(12:15) gp > sumdigits2(n)=if(n%9!=0,return(n%9),return(9))
%6 = (n)->if(n%9!=0,return(n%9),return(9))
(12:15) gp > sumdigits3(n)=n=n%9;if(n,n,9)
%7 = (n)->n=n%9;if(n,n,9)
(12:15) gp > for(x=1,10000000,sumdigits1(x))
(12:16) gp > ##
*** last result computed in 6,172 ms.
(12:16) gp > for(x=1,10000000,sumdigits2(x))
(12:16) gp > ##
*** last result computed in 16,563 ms.
(12:16) gp > for(x=1,10000000,sumdigits3(x))
(12:17) gp > ##
*** last result computed in 15,141 ms.[/CODE]

unfortunately my first sumdigits seems most effective when added into lucaslehmers()

Interesting. I wonder what Pari is doing under the hood here.

Mine is still faster for large numbers, but I can't imagine why yours is faster for small numbers.

The optimal Pari (as opposed to GP) code is surely
[code]long
digital_root(GEN n) {
if (typ(n) != t_INT)
pari_err(typeer, "digital_root");
register long ret = smodis(n, 9);
if (ret)
return ret;
return 9;
}[/code]
and I can't guess why your version would be closer to this than mine. Maybe I'll gp2c them later and find out.

[QUOTE=science_man_88;238510]unfortunately my first sumdigits seems most effective when added into lucaslehmers()[/QUOTE]

Do you need the Lucas-Lehmer residue, or just its digital root? You can compute the latter much faster.

[QUOTE=CRGreathouse;238530]Do you need the Lucas-Lehmer residue, or just its digital root? You can compute the latter much faster.[/QUOTE]

Digital root of the lucas lehmer residue ... really I didn't know that, Thanks. I'm just trying to limit it to a list of digital root 9 because we know the exponents must be in that group to work.

[QUOTE=CRGreathouse;238501]So, neither. It's like the digital root but with sumdigits(0) = 9.

Faster code:
[code]sumdigits(n)=n=n%9;if(n,n,9)[/code][/QUOTE]

what about
[code]sumdigits(n)=(n-1)%9+1[/code]?

[QUOTE=axn;238542]what about
[code]sumdigits(n)=(n-1)%9+1[/code]?[/QUOTE]

using flat out print(sumdigits(x)) it ties my fastest on my machine. I'll have to check it out on the lucas test.

[CODE](15:57) gp > for(x=1,10000000,sumdigits1(x))
(16:01) gp > ##
*** last result computed in 6,172 ms.
(16:01) gp > for(x=1,10000000,sumdigits2(x))
(16:01) gp > ##
*** last result computed in 16,891 ms.
(16:01) gp > for(x=1,10000000,sumdigits3(x))
(16:01) gp > ##
*** last result computed in 15,312 ms.
(16:01) gp > for(x=1,10000000,sumdigits4(x))
(16:02) gp > ##
*** last result computed in 5,266 ms.
(16:02) gp >[/CODE]

we have a new winner.

[CODE](16:04) gp > for(x=1,100000,lucaslehmer2(x))
3,5,7,13,17,19,23,31,33,51,61,71,89,101,107,127,139,191,271,273,305,331,347,351,367,397,405,407
*** _^s: user interrupt after 12,562 ms.
(16:04) gp > for(x=1,1000,lucaslehmer2(x))
3,5,7,13,17,19,23,31,33,51,61,71,89,101,107,127,139,191,271,273,305,331,347,351,367,397,405,407
(16:04) gp > ##
*** last result computed in 3,922 ms.
(16:04) gp > lucaslehmer2(p)= s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2 && sumdigits2(s)==9,p
%21 = (p)->s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2&&sumdigits2(s)==9,print1(p","))
(16:05) gp > for(x=1,1000,lucaslehmer2(x))
3,5,7,13,17,19,23,31,33,51,61,71,89,101,107,127,139,191,271,273,305,331,347,351,367,397,405,407
(16:05) gp > ##
*** last result computed in 3,953 ms.
(16:05) gp > lucaslehmer2(p)= s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2 && sumdigits3(s)==9,p
%22 = (p)->s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2&&sumdigits3(s)==9,print1(p","))
(16:05) gp > for(x=1,1000,lucaslehmer2(x))
3,5,7,13,17,19,23,31,33,51,61,71,89,101,107,127,139,191,271,273,305,331,347,351,367,397,405,407
(16:05) gp > ##
*** last result computed in 3,969 ms.
(16:05) gp > lucaslehmer2(p)= s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2 && sumdigits4(s)==9,p
%23 = (p)->s=4;for(x=1,p-2,s=(s^2-2)%(2^p-1));if(x=p-2&&sumdigits4(s)==9,print1(p","))
(16:05) gp > for(x=1,1000,lucaslehmer2(x))
3,5,7,13,17,19,23,31,33,51,61,71,89,101,107,127,139,191,271,273,305,331,347,351,367,397,405,407
(16:05) gp > ##
*** last result computed in 4,000 ms.
(16:05) gp >[/CODE]

this is weirder and weirder lol.

[QUOTE=science_man_88;238535]Digital root of the lucas lehmer residue ... really I didn't know that, Thanks. I'm just trying to limit it to a list of digital root 9 because we know the exponents must be in that group to work.[/QUOTE]

Actually, I'm not sure about my claim anymore. It would certainly be easy to compute [TEX](2+\sqrt3)^{2^{p-2}}[/TEX] mod 9 rather than mod [TEX]2^p-1,[/TEX] but that won't give you what you want since [TEX]9\not\mid2^p-1.[/TEX]

Of course, to me, that suggests that this is not really going to be meaningful... but that's your issue, not mine.

[CODE][COLOR="Red"]3,5,7,13,17,19,[/COLOR]23,[COLOR="red"]31,[/COLOR]33,51,[COLOR="red"]61,[/COLOR]71,[COLOR="red"]89,[/COLOR]101,[COLOR="red"]107,127,[/COLOR]139,191,271,273,305,331,347,351,367,397,405,407,427,435,457,467,489,[COLOR="red"]521,[/COLOR]525,539,543,549,559,565,577,583,589,597,601,[COLOR="red"]607,[/COLOR]611,613,617,619,641,643,661,693,717,729,787,793,809,817,819,837,871,879,891,899,983,987,991,[/CODE]

this is the sequence I got for lucaslehmer2()

[QUOTE=science_man_88;238535]Digital root of the lucas lehmer residue ... really I didn't know that, Thanks. I'm just trying to limit it to a list of digital root 9 because we know the exponents must be in that group to work.[/QUOTE]

Exponents? 0 mod 9? 9 does not divide any prime, so this condition is impossible to meet.

You can try 1 mod 9, 2 mod 9, 4 mod 9, 5 mod 9, 7 mod 9, and 8 mod 9.

[QUOTE=3.14159;238561]Exponents? 0 mod 9? 9 does not divide any prime, so this condition is impossible to meet.

You can try 1 mod 9, 2 mod 9, 4 mod 9, 5 mod 9, 7 mod 9, and 8 mod 9.[/QUOTE]

I'm not saying anything about a prime dividing by 9 lol. If 0=9 mod 9 then you know the mersenne prime exponent residues all fit in the group of exponents with residues that are 9 mod 9, they are a subgroup.

[QUOTE=science_man_88;238577]I'm not saying anything about a prime dividing by 9 lol. If 0=9 mod 9 then you know the mersenne prime exponent residues all fit in the group of exponents with residues that are 9 mod 9, they are a subgroup.[/QUOTE]

Okay.. what about it?

And;

Adding Fermat's factoring algorithm was rather easy;

[code]fmtfcs(a, x, m) = {
for(n=a,x,
if(issquare(n^2+m), print(n))
);
}[/code]

Okay; A few edits to that..

Revision:[code]fmtfcs(a,x,m)=for(n=a,x,if(issquare(n^2+m),print(sqrtint(n^2+m)+n, " * ", sqrtint(n^2+m)-n)));[/code]

Bad Pi! :wink: Functions should return, not print. What if you wanted to use that data from another program?

[QUOTE=CRGreathouse;238588]Bad Pi! :wink: Functions should return, not print. What if you wanted to use that data from another program?[/QUOTE]

Okay, fine. I've changed it to "return". Nope, that's unable to return the factors. Back to printing it goes.

I'm busy collecting relations for a c106. I'm only about ... 31% done. 159k relations are needed. So, I estimate that this will take two or three days. (I've only collected 49105 relations, so far, given eight hours.)

If adding Fermat's was easy.. How would one go about in adding an amateur QS?

(Yes, yes, I know that PARI already uses MPQS.)

My knowledge of it so far is;

1. Set bound.
2. Set factor base.
3. b-smooth congruence hunting.
4. Try to set ≥2 congruences into a congruence of squares.
5. gcd.
6. Factors.

I did a comparison, and it seems that trial division by far outspeeds the Fermat script.

I think I didn't have the most efficient script there is.. But I'll try improving on it.

I'm now 2/3 done with the collection of relations for the c106 I'm trying to split (Hopefully it does not crash.)|

Update: 72% done.

Record split; Part 3!

[code]starting SIQS on c106: 4475950135356613778937951617143741209311215873146646798431494279597781200242083780419779316018232365322301

==== sieving in progress ( 4 threads): 159064 relations needed ====
==== Press ctrl-c to abort and save state ====
159188 rels found: 35794 full + 123394 from 2418802 partial, ( 48.30 rels/sec)

SIQS elapsed time = 51034.7770 seconds.


***factors found***

PRP53 = 59032375916416525973964206982565204661985227334884669
PRP53 = 75821954747240330698482729033441429926359516470800129
[/code]

[QUOTE=3.14159;238632]I did a comparison, and it seems that trial division by far outspeeds the Fermat script. [/QUOTE]

That's actually usual for Fermat's method, so this doesn't indicate that your version is bad.

[QUOTE=CRGreathouse;238648]That's actually usual for Fermat's method, so this doesn't indicate that your version is bad.[/QUOTE]

It might be a more difficult task to code a factoring method which is faster than trial division.

I did some reading on some of the highly composite numbers; Is there any way in which one could compute a specific one, without checking every integer smaller than it?

Ex: I wanted to compute the 100th term in the sequence. Any way to do that without checking every integer before it?

[QUOTE=3.14159;238679]It might be a more difficult task to code a factoring method which is faster than trial division.[/QUOTE]

I recommend Pollard's rho, it's pretty easy.

[QUOTE=3.14159;238685]I did some reading on some of the highly composite numbers; Is there any way in which one could compute a specific one, without checking every integer smaller than it?

Ex: I wanted to compute the 100th term in the sequence. Any way to do that without checking every integer before it?[/QUOTE]

The fastest way would be to compute the first 100 members of the sequence, which does not require checking all numbers up to that point (or even checking the products of primorials up to that point). There's extensive discussion at a link from the Sloane sequence; it's from the person who calculated lots of terms there.

I am stumped. Must be something really trivial:
[CODE]gp > polcyclo(15,10)
*** variable name expected: polcyclo(15,10)
^---
gp > p=polcyclo(15)
%1 = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
gp > p(10)
*** unused characters: p(10)
^----
[/CODE]
How do I eval? I understand that I create a POL object. But how do I marry it to the argument? [FONT=Courier New]subst(polcyclo(15),x,10)[/FONT] ? Looks ugly.

[QUOTE=Batalov;238709]I am stumped. Must be something really trivial:
[CODE]gp > polcyclo(15,10)
*** variable name expected: polcyclo(15,10)
^---
gp > p=polcyclo(15)
%1 = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
gp > p(10)
*** unused characters: p(10)
^----
[/CODE]
How do I eval? I understand that I create a POL object. But how do I marry it to the argument? [FONT=Courier New]subst(polcyclo(15),x,10)[/FONT] ? Looks ugly.[/QUOTE]

Doesn't work?

I get a return of 90090991.

[QUOTE=CRGreathouse;238691]I recommend Pollard's rho, it's pretty easy.[/QUOTE]

Well, I'll have to learn how Pollard's Rho works first.

[QUOTE=Batalov;238709]
How do I eval?[/QUOTE]

eval() lol

[QUOTE=Batalov;238709]I am stumped. Must be something really trivial:
[CODE]gp > polcyclo(15,10)
*** variable name expected: polcyclo(15,10)
^---
gp > p=polcyclo(15)
%1 = x^8 - x^7 + x^5 - x^4 + x^3 - x + 1
gp > p(10)
*** unused characters: p(10)
^----
[/CODE]
How do I eval? I understand that I create a POL object. But how do I marry it to the argument? [FONT=Courier New]subst(polcyclo(15),x,10)[/FONT] ? Looks ugly.[/QUOTE]

substpol (or subst) are usual. If all you want is to evaluate the polynomial (not to take derivatives or extract coefficients, etc.) then it's best to create it as a closure instead:
[code]p=n->polcyclo(15,n);
p(15)[/code]

[QUOTE=CRGreathouse;238752]substpol (or subst) are usual. If all you want is to evaluate the polynomial (not to take derivatives or extract coefficients, etc.) then it's best to create it as a closure instead:
[code]p=n->polcyclo(15,n);
p(15)[/code][/QUOTE]
Thanks! That's what I was looking for.
[SIZE=1](and yes, subst() works, too, but I was looking for a hint to best practices)[/SIZE]

Wait; How do you code Pollard's Rho?

[QUOTE=3.14159;238811]Wait; How do you code Pollard's Rho?[/QUOTE]

[url]http://en.wikipedia.org/wiki/Pollard's_rho_algorithm[/url] may help.

If I knew a few more things maybe I could make it work,However I don't so I won't.

[QUOTE=Batalov;238805]Thanks! That's what I was looking for.
[SIZE=1](and yes, subst() works, too, but I was looking for a hint to best practices)[/SIZE][/QUOTE]

Yes, thats certainly substpol. Honestly, I haven't tested this vs. subst for performance or anything, but I imagine there are reasons to split the special case out.

Here's the script by D. Broadhurst in all its beauty minus some lines that I couldn't understand or that were not entirely by the book and cost nothing to be by the book. The script's result appears to tell us that this
27810-digit number is prime by Konyagin-Pomerance N-1 test with 30.2% factored N-1.

I was looking for a way not to put the 7334-digit prime there literally.

If this is correct by the book, I could proceed adding it to PFGW. Does it look proper to you, Charles? Thanks.

P.S. That previous number was actually 32.99% factored (but still short of B-L-S N-1 proof; N+1 side doesn't help enough). Here's the [I]30.2[/I]% factored N-1-based KP proof. Runs really fast, because the number is fairly small and within reach of Primo anyway.

[QUOTE=Batalov;238857]If this is correct by the book, I could proceed adding it to PFGW. Does it look proper to you, Charles? Thanks.[/QUOTE]

It looks fine. There are things I would do differently, to wit:
[code]test()={
print("D Broadhurst, KP N-1 test, 8 Jul 2001");
\\ checked for correctness as a Konyagin-Pomerance N-1 test

\\ Origin: http://tech.groups.yahoo.com/group/primeform/files/KP/KonPom.gp
\\ KP Combined test seemed to be a work in progress
\\ it would need more operations as were in the original script
\\ (See Theorem 4.2.9, p.186 in Crandall-Pomerance PN-ACP)

my(N=(8*10^27811-71)/9,p7334,p557,p202,lsm,F,R,c1,c4);

\\ Supply N-1 factors in the lsm[] array:

p7334 = subst(polcyclo(13905),x,10)/83431/333721;
p557 = subst(polcyclo(927),x,10)/46351/181693/1405333/43604227/75178674467317/664958944842271489;
p202 = subst(polcyclo(618),x,10)/619;

lsm=[2,2,2,2,3,3,3,5,7,11,13,19,31,37,41,211,241,271,619,1031,1237,2161,7211,9091,29611,46351,52579,181693,238681,333667,380071,\
497491,704521,1358983,1405333,2906161,3762091,7034077,\
43604227,44092859,102860539,569836171,984385009,2013681931,5905014721,326345481191,8595591045751,8985695684401,75178674467317,\
301917502615411,612053256358933,2702394989404991,41343860637475219,167798066344417387,664958944842271489,4185502830133110721,\
244261115576975427841,290934235666553153131,3462948230491376473027,182725114866521155647161,567664121498896654896380203841,\
896048585318577702680084550566846611,1471865453993855302660887614137521979,7073532436575439739956079653594194521,\
85517237338860102283865323577960208481,5294796903161592416528456780680376286484870226446771978908657527791,\
153211620887015423991278431667808361439217294295901387715486473457925534859044796980526236853,\
292815973319957863643663090772618219691707118842359217315010640397456279188940711108958040881306849,\
757,5563,6481,83431,333721,1577071,16357951,22053331,70541929,310362841,2016447481,\
14175966169,18245696041,258360989311,577603663291,31023833790241,440334654777631,\
483418418597220677238517353915231961831,8610583349234340055547908764091017276717091,\
p202,p557,p7334];

F=prod(k=1,length(lsm),lsm[k]);
R=(N-1)/F;

if(denominator(R)==1
&& gcd(F,R)==1
&& F^(1/3)/6>1
&& 1>3*N/F^(10/3),
,
error("Fail")
);

\\ Here, log(F)/log(N) =~ 0.3299 > 3/10, but less than 1/3
print("Factored part is: ", log(F)/log(N));

c1=divrem(R,F);
c4=c1[1];
c1=c1[2];

print();

\\ Unsure why the square test was not performed six times in the original script
\\ This is part 1 from Theorem 4.1.6 (or 4.2.9):
for(t=0,5,
if(issquare((c1+t*F)^2+4*t-4*c4),error("Fail @ "t))
);

my(b=contfrac(c1/F), v=0, vn=1, u=1, un=b[1], n=1);
while(F^(1/3)>vn,
bn=b[n++];
uo=u;
u=un;
vo=v;
v=vn;
un=bn*un+uo;
vn=bn*vn+vo
);

\\ Here, u/v is as per Part 2 of Theorem 4.1.6

my(d=round(c4*v/F), zs=[v,u*F-c1*v,c4*v-d*F+u,-d], q=0);
for(k=1,4,
print("zs["k"] = ", zs[k]);
q += zs[k]*x^(4-k)
);

print();

\\ required precision is dependent on the size of N
default(realprecision, 10000);
my(ps=polroots(q), r, q1, q2);

for(k=1,3,
r=floor(real(ps[k]));
print(r);
q1=substpol(q, x, r);
q2=substpol(q, x, r+1);
if(q1*q2>=0,error("Fail root ",k))
);

print("OK, no integral roots");
\\ if there's a root a, check: if aF+1 is a nontrivial factor of N, then N is composite

\\ Note to myself: for implementation in PFGW: one root is very close to 0 and can be found by Newton iter, with initial 0; then reduce to quadratic.
print(ps[2]);
};

test()[/code]
but these are largely stylistic differences.

[QUOTE=CRGreathouse;238648]That's actually usual for Fermat's method, so this doesn't indicate that your version is bad.[/QUOTE]

Meh. I think I've read too many kook sites claiming to have some drastic improvement to Fermat's algorithm in one way or another, and wished to entertain the thought of programming those in, then trying to see if the algorithm, with these "improvements, actually competes with trial factoring, for a fairly large semiprime number.

Here's [URL="http://www.naturalnumbers.org/fixfermatfact.html"]one[/URL] of them.

Some sort of iterated algorithm to try, for some fairly large semiprime n (Ex: 5403615307171909)

.. Oh, dear. This tutorial of theirs is a nightmare.

[QUOTE=3.14159;238905]Meh. I think I've read too many kook sites claiming to have some drastic improvement to Fermat's algorithm in one way or another, and wished to entertain the thought of programming those in, then trying to see if the algorithm, with these "improvements, actually competes with trial factoring, for a fairly large semiprime number.

Here's [URL="http://www.naturalnumbers.org/fixfermatfact.html"]one[/URL] of them.

Some sort of iterated algorithm to try, for some fairly large semiprime n (Ex: 5403615307171909)

.. Oh, dear. This tutorial of theirs is a nightmare.[/QUOTE]

the main part of the site I don't get without the code is the table they list off. I should know vb.net a bit though it's only been under 2 years.

Figured it out I think. Want my take on it Pi ?

[QUOTE=science_man_88;238917]Figured it out I think. Want my take on it Pi ?[/QUOTE]

Sure.

[QUOTE=3.14159;238919]Sure.[/QUOTE]

the important part is this:

[CODE]1. Determine the integer square root of N.

sqrN = Int(Sqr(n))

2. Determine the square of Step 1.

sqrN2 = sqrN * sqrN

(Note that sqrN2 is initially the greatest perfect square less than N.)

3. Determine the integer square root of Step 2 minus N.

PCS = Int(Sqr(Abs(sqrN2 - n)))

(Let's say PCS means 'perfected counting square'.)

4. Determine if the sum of Step 2 and Step 3 is a factor of N.

If mod(n, PCS + sqrN) = 0 then we have a factor, and we're done.

(Similarly, PCS minus sqrN would produce the same result.)

5. If not, add 1 to Step 1 (that is, sqrN+1) and repeat.[/CODE]


I think you can start with a for loop setting the minimum to floor(sqrt(n)) and looping until the number if needed (don't worry there will be a break command). sqrN2 = sqrN * sqrN then set this up as another variable just in the loop( it depends on the loop variable). PCS = Int(Sqr(Abs(sqrN2 - n))) also set in the loop(in th order shown, replacing Int() with floor()). Make and if statement(in the loop) to check the modulus. If it is false do nothing if it is 0 break(1) ( however best practice I think is to use the reverse jmp of asm so may have to reverse something).

[QUOTE=science_man_88;238921]the important part is this:

[CODE]1. Determine the integer square root of N.

sqrN = Int(Sqr(n))

2. Determine the square of Step 1.

sqrN2 = sqrN * sqrN

(Note that sqrN2 is initially the greatest perfect square less than N.)

3. Determine the integer square root of Step 2 minus N.

PCS = Int(Sqr(Abs(sqrN2 - n)))

(Let's say PCS means 'perfected counting square'.)

4. Determine if the sum of Step 2 and Step 3 is a factor of N.

If mod(n, PCS + sqrN) = 0 then we have a factor, and we're done.

(Similarly, PCS minus sqrN would produce the same result.)

5. If not, add 1 to Step 1 (that is, sqrN+1) and repeat.[/CODE]


I think you can start with a for loop setting the minimum to floor(sqrt(n)) and looping until the number if needed (don't worry there will be a break command). sqrN2 = sqrN * sqrN then set this up as another variable just in the loop( it depends on the loop variable). PCS = Int(Sqr(Abs(sqrN2 - n))) also set in the loop(in th order shown, replacing Int() with floor()). Make and if statement(in the loop) to check the modulus. If it is false do nothing if it is 0 break(1) ( however best practice I think is to use the reverse jmp of asm so may have to reverse something).[/QUOTE]

Without all their fluffy language, I understand it as:

1. sqrtint(n).
2. (sqrtint(n))^2.
3. n - (sqrtint(n))^2.
4. (sqrtint(n))^2 + (n -(sqrtint(n))^2) mod n.
5. If (sqrtint(n))^2 + (n -(sqrtint(n))^2) mod n = 0, factors found.
6. If (sqrtint(n))^2 + (n -(sqrtint(n))^2) mod n !=0, repeat steps 1-5 until factors are found.

I understand it as:

[CODE]test(n)=for(sqrN=floor(sqrt(n)),n,sqrN2=sqrn^2;PCS=floor(sqrt(abs(sqrN2-n)));if(Mod(sqrN2+PCS,n)!=0,,break()))[/CODE]

but Pari tells me:

[CODE] *** floor: incorrect type in gfloor.[/CODE]

[QUOTE=science_man_88;238924]I understand it as:

[CODE]test(n)=for(sqrN=floor(sqrt(n)),n,sqrN2=sqrn^2;PCS=floor(sqrt(abs(sqrN2-n)));if(Mod(sqrN2+PCS,n)!=0,,break()))[/CODE]

but Pari tells me:

[CODE] *** floor: incorrect type in gfloor.[/CODE][/QUOTE]

Just use sqrtint(n), and n. There's no need to toss in unnecessary variables.

Since they call it "Improved squares", I'll title it, "impsqs".

I'll have to follow the six steps.

[CODE]test(n)=for(sqrN=sqrtint(n),n,sqrN2=sqrN^2;PCS=sqrtint(sqrN2-n);if(Mod(sqrN2+PCS,n)!=0,,break()))[/CODE]

now I get:

[CODE]*** sqrtint: negative integer in sqrtint.[/CODE]

[QUOTE=science_man_88;238927][CODE]test(n)=for(sqrN=sqrtint(n),n,sqrN2=sqrN^2;PCS=sqrtint(sqrN2-n);if(Mod(sqrN2+PCS,n)!=0,,break()))[/CODE]

now I get:

[CODE]*** sqrtint: negative integer in sqrtint.[/CODE][/QUOTE]

Take a closer look at the script. There's a mistake there, and it might be right under your nose.

[QUOTE=3.14159;238928]Take a closer look at the script. There's a mistake there, and it might be right under your nose.[/QUOTE]

I don't see it. they have PCS= sqrtint(sqrN2-n).

oh reversed the modulo dah lol.

[QUOTE=science_man_88;238929]I don't see it. they have PCS= sqrtint(sqrN2-n).[/QUOTE]

Here's how to do it without using a script, in PARI/GP;

[CODE](14:03) gp > 3926696294764459
%203 = 3926696294764459
(14:03) gp > sqrtint(%)
%204 = 62663356
(14:03) gp > %204^2
%205 = 3926696185182736
(14:03) gp > %203-%205
%206 = 109581723
(14:04) gp > Mod(%203, %206 + %204)
%207 = Mod(81019925, 172245079)[/CODE]

Okay, I just realized that was incorrect.

[QUOTE=3.14159;238930]Here's how to do it without using a script, in PARI/GP;

[CODE](14:03) gp > 3926696294764459
%203 = 3926696294764459
(14:03) gp > sqrtint(%)
%204 = 62663356
(14:03) gp > %204^2
%205 = 3926696185182736
(14:03) gp > %203-%205
%206 = 109581723
(14:04) gp > Mod(%203, %206 + %204)
%207 = Mod(81019925, 172245079)[/CODE][/QUOTE]

Yeah they point out subtraction should work must of confused me lol.

[QUOTE]Int(Sqr(Abs(sqrN2 - n)))[/QUOTE]

[QUOTE=science_man_88;238931]Yeah they point out subtraction should work must of confused me lol.[/QUOTE]

The "Abs", means, "The absolute value of".

Therefore, there are no negative numbers involved.

Perhaps you can develop a script from that piece of information.

I've only got one factor properly, for the number they originally give the example with. The the other factor became 37 though.

got it working:
[CODE]test(n)=for(sqrN=sqrtint(n),n,sqrN2=sqrN^2;PCS=sqrtint(abs(sqrN2 - n));if(Mod(n,PCS+sqrN)!=0,,print(sqrN+PCS"*"sqrN-PCS);break()));[/CODE]

[QUOTE=science_man_88;238936]got it working:
[CODE]test(n)=for(sqrN=sqrtint(n),n,sqrN2=sqrN^2;PCS=sqrtint(abs(sqrN2 - n));if(Mod(n,PCS+sqrN)!=0,,print(sqrN+PCS"*"sqrN-PCS);break()));[/CODE][/QUOTE]

It works, but trial division outspeeds it by a factor of 40.

Nevermind.. It gains the edge when it comes to a c18.

[QUOTE=3.14159;238937]It works, but trial division outspeeds it by a factor of 40.

I knew these guys were lying kooks.[/QUOTE]

Yeah the largest Mersenne prime I felt like letting it finish was 2^19-1, anything over that I stopped at about 15 seconds in.

Also; For the rest of you out there;

I thought of this type of prime; k * (n! * p(n)#) + 1, or, if you'd like to generalize further, + b, where b is an integer.

There should be a k such that k * (n! * p(n)#) + 1 is composite with every n > 1.

I'll try looking for some that fit that description, to some extent.

[QUOTE=3.14159;238939]Also; For the rest of you out there;

I thought of this type of prime; k * (n! * p(n)#) + 1, or, if you'd like to generalize further, + b, where b is an integer.

There should be a k such that k * (n! * p(n)#) + 1 is composite with every n > 1.

I'll try looking for some that fit that description, to some extent.[/QUOTE]

p(n)# means what again ? I forgot.

[QUOTE=science_man_88;238940]p(n)# means what again ? I forgot.[/QUOTE]

p(n)#.

p(5)# = 2 * 3 * 5 * 7 * 11 = 2310.

Of course, the smallest prime of the type I mentioned is 3.

[QUOTE=3.14159;238941]p(n)#.

p(5)# = 2 * 3 * 5 * 7 * 11 = 2310.[/QUOTE]

so a primorial up to prime(n).

[QUOTE=science_man_88;238943]so a primorial up to prime(n).[/QUOTE]

Yes. Product of the first n primes.

And; Patch this for me.

[CODE]cs(a,x,m)=for(n=a,x,print(m*fp(n)+1));if(isprime(m*fp(n)+1),break)[/CODE]

Note that fp(n) = n! * p(n)#.

I got an error message; [CODE]*** gtos expected an integer, got 'n'.[/CODE]

[QUOTE=3.14159;238944]Yes. Product of the first n primes.

And; Patch this for me.

[CODE]cs(a,x,m)=for(n=a,x,print(m*fp(n)+1));if(isprime(m*fp(n)+1),break)[/CODE]

Note that fp(n) = n! * p(n)#.

I got an error message; [CODE]*** gtos expected an integer, got 'n'.[/CODE][/QUOTE]

only thing I see wrong is p(n)[COLOR="Red"]#[/COLOR] lol

never mind I fixed it embed the if and it worked for me lol.

[QUOTE=science_man_88;238959]only thing I see wrong is p(n)[COLOR="Red"]#[/COLOR] lol

never mind I fixed it embed the if and it worked for me lol.[/QUOTE]

.. Wait.. What revisions did you make to the code?

Also; Prime: 590 * (400! * p(400)#) + 1.

[QUOTE=3.14159;238966].. Wait.. What revisions did you make to the code?

Also; Prime: 590 * (400! * p(400)#) + 1.[/QUOTE]

[CODE]cs(a,x,m)=for(n=a,x,print(m*fp(n)+1);if(isprime(m*fp(n)+1),break(1)))[/CODE]

Still baffled by the error in your first code Pi ? If so I can tell you it's simple, in fact I over thought it.

[QUOTE= Pi][COLOR="lime"]cs(a,x,m)=for(n=a,x,print(m*fp(n)+1)[/COLOR]);[COLOR="lime"]if(isprime(m*fp(n)+1),break)[/COLOR][/QUOTE]
[QUOTE=science_man_88;238970][CODE][COLOR="Lime"]cs(a,x,m)=for(n=a,x,print(m*fp(n)+1)[/COLOR];[COLOR="lime"]if(isprime(m*fp(n)+1),break([COLOR="Black"]1[/COLOR]))[/COLOR])[/CODE][/QUOTE]

figure it out...

[QUOTE=science_man_88;238976]figure it out...[/QUOTE]

I already did; I was merely inattentive. Don't flatter yourself.

[QUOTE=3.14159;238978]I already did; I was merely inattentive. Don't flatter yourself.[/QUOTE]

Sorry I just see no reason for me to correct you if I don't,as you have more ability than me.

[QUOTE=science_man_88;238979]Sorry I just see no reason for me to correct you if I don't,as you have more ability than me.[/QUOTE]

I have an idea, in order to search for those k's with no primes. I'll do an amateur sieve, but it keeps [B]primes[/B] out instead.

Back to the original topic for a moment:

I've been doing some work over at Rosetta Code, adding examples of simple tasks in Pari/GP. Here's the link:
[url]http://rosettacode.org/wiki/Category:PARI/GP[/url]

You may find this useful if you're trying to do something simple in Pari that seems like it should be easy, but you just don't know the commands or can't think of a good way to do it. I found it very instructive to write these examples; you may be able to either learn from it or teach me (feel free to post here, PM me, or simply change one of the RC entries).

I have 132 tasks solved at the moment,* ranging from simple ones like language basics (loops, arithmetic) and primality by trial division to complicated ones like [url=http://rosettacode.org/wiki/Roots_of_a_function#PARI.2FGP]root finding[/url].

* Actually, I imagine that a few of those were written by others, but I certainly wrote the bulk. You can look at the edit history if you're interested in authorship (but why would you be?).

[QUOTE=CRGreathouse;238992]Back to the original topic for a moment:

I've been doing some work over at Rosetta Code, adding examples of simple tasks in Pari/GP. Here's the link:
[url]http://rosettacode.org/wiki/Category:PARI/GP[/url]

You may find this useful if you're trying to do something simple in Pari that seems like it should be easy, but you just don't know the commands or can't think of a good way to do it. I found it very instructive to write these examples; you may be able to either learn from it or teach me (feel free to post here, PM me, or simply change one of the RC entries).

I have 132 tasks solved at the moment,* ranging from simple ones like language basics (loops, arithmetic) and primality by trial division to complicated ones like [url=http://rosettacode.org/wiki/Roots_of_a_function#PARI.2FGP]root finding[/url].

* Actually, I imagine that a few of those were written by others, but I certainly wrote the bulk. You can look at the edit history if you're interested in authorship (but why would you be?).[/QUOTE]

thanks maybe I can use my book on system commands of ms-dos
that I own to help system() is a function just can't remember the coding for some of the things on [URL="http://rosettacode.org/wiki/Reports:Tasks_not_implemented_in_PARI/GP"]here[/URL] and some of these are talked of in number freak maybe that can help lol.

Also stack is on that list. I've actually wanted to for a long time try to make pari resemble asm lol.