A PARI script to search for sociable numbers
 

I've been searching for sociable numbers in the 10^15 area for a while now, using the following PARI script:
[code]
pS=readvec(pariStatus)
a=pS[1]
for(b=a,a+10^9,i=b;for(n=1,200,m=n;i=sigma(i)-i;if(i<b||i>b*20,break()));if(m>199,print(b);write("pariCyclesFound",b));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b)))
[/code]The script uses the file pariStatus to keep track of how far along it is, such as:
[code]
1000010274050000
[/code]In this example, the script will start searching at 1000010274050000, find a cycle and keep searching:
[code]
...
1000010274060000 5.6430000000000000000000000000000000000 seconds/10000
1000010274063430
1000010274070000 5.5030000000000000000000000000000000000 seconds/10000
...
[/code]In truth. it doesn't guarantee it's a cycle, but if the sequence goes for 200 iterations without increasing to b*20, it is at least interesting. In the above example, that is a known amicable number. I do a lot more printing than I probably need to, but it shows activity. I'd prefer overwriting the pariStatus file, but couldn't figure out how to do that with PARI, so you can see the extern("rm pariStatus") command to get rid of the previous file prior to writing the current value. This is a time waster and is definitely hazardous if a power loss is timed just right, but I couldn't quite figure out a better way.

For a while I had nine machines set up to search 11*10^14, 12*10^14, ... 19*10^14, respectively. Each was set to start on power-up, so if the AC was lost, they would begin again, within a few numbers (seconds) of where they were interrupted. Occasionally, I would record their progress in case an interruption occurred at a time that caused pariStatus to be corrupted or lost.

I'm sure this could be improved upon and it may even have some flaw(s) I've not considered. So, I'm open to all critique...

[QUOTE=EdH;275857]I've been searching for sociable numbers in the 10^15 area for a while now, using the following PARI script:
[code]
pS=readvec(pariStatus)
a=pS[1]
for(b=a,a+10^9,i=b;for(n=1,200,m=n;i=sigma(i)-i;if(i<b||i>b*20,break()));if(m>199,print(b);write("pariCyclesFound",b));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b)))
[/code]The script uses the file pariStatus to keep track of how far along it is, such as:
[code]
1000010274050000
[/code]In this example, the script will start searching at 1000010274050000, find a cycle and keep searching:
[code]
...
1000010274060000 5.6430000000000000000000000000000000000 seconds/10000
1000010274063430
1000010274070000 5.5030000000000000000000000000000000000 seconds/10000
...
[/code]In truth. it doesn't guarantee it's a cycle, but if the sequence goes for 200 iterations without increasing to b*20, it is at least interesting. In the above example, that is a known amicable number. I do a lot more printing than I probably need to, but it shows activity. I'd prefer overwriting the pariStatus file, but couldn't figure out how to do that with PARI, so you can see the extern("rm pariStatus") command to get rid of the previous file prior to writing the current value. This is a time waster and is definitely hazardous if a power loss is timed just right, but I couldn't quite figure out a better way.

For a while I had nine machines set up to search 11*10^14, 12*10^14, ... 19*10^14, respectively. Each was set to start on power-up, so if the AC was lost, they would begin again, within a few numbers (seconds) of where they were interrupted. Occasionally, I would record their progress in case an interruption occurred at a time that caused pariStatus to be corrupted or lost.

I'm sure this could be improved upon and it may even have some flaw(s) I've not considered. So, I'm open to all critique...[/QUOTE]

add one more layer to this:

[QUOTE]for(x=2,1000000,if((sigma(x)-x)!=x && sigma(sigma(x)-x)-(sigma(x)-x)==x ,print(x","sigma(x)-x)))[/QUOTE] and you have sociable numbers. this checks for amicable pairs ( just adding another sigma and a copy of what is there should work for 3).

[QUOTE=science_man_88;275858]add one more layer to this:

and you have sociable numbers. this checks for amicable pairs ( just adding another sigma and a copy of what is there should work for 3).[/QUOTE]
I'm not aware of any order 3 sociable numbers, but there are order 4, 5, 6, 8, 9 and 28, so I'd like my search to include at least that many levels.

[QUOTE=EdH;275868]I'm not aware of any order 3 sociable numbers, but there are order 4, 5, 6, 8, 9 and 28, so I'd like my search to include at least that many levels.[/QUOTE]

I was thinking of making a loop to change the level of the sigma, but it might have to jump a lot to work.

[QUOTE=EdH;275868]I'm not aware of any order 3 sociable numbers, but there are order 4, 5, 6, 8, 9 and 28, so I'd like my search to include at least that many levels.[/QUOTE]

sorry for posting again so soon people. the following is what I come to.
[CODE]for(x=3,1000000,b=x;for(y=1,100,if(sigma(b)-b!=0,b=sigma(b)-b);if(b==x,print(x","y);break())))[/CODE]

[QUOTE=science_man_88;275873]sorry for posting again so soon people. the following is what I come to.
[CODE]for(x=3,1000000,b=x;for(y=1,100,if(sigma(b)-b!=0,b=sigma(b)-b);if(b==x,print(x","y);break())))[/CODE][/QUOTE]
Shouldn't be a problem with valid users posting "too soon!"

I like your approach, but I need to get a better picture of the whole. I think there should be a break if sigma(b)-b falls below x. I'm going to need to add my tracking and check timing.

Off to study, but maybe not too much tonight.

Unfortunately, a preliminary trial with edits here didn't find the test cycle from the previous example...:sad:, but it ran faster:smile:...

Thanks...

[QUOTE=EdH;275919]Shouldn't be a problem with valid users posting "too soon!"

I like your approach, but I need to get a better picture of the whole. I think there should be a break if sigma(b)-b falls below x. I'm going to need to add my tracking and check timing.

Off to study, but maybe not too much tonight.

Unfortunately, a preliminary trial with edits here didn't find the test cycle from the previous example...:sad:, but it ran faster:smile:...

Thanks...[/QUOTE]

how long does it take to cycle because I did checks up to 1000000 in length and it can't find it also I checked by hand that's not the full cycle if it is one:

[CODE](09:00)>sigma(1000010274070000)-1000010274070000
%154 = 1421114600505088
(09:00)>sigma(%)-%
%155 = 1415563371597376
(09:01)>sigma(%)-%
%156 = 1515576001115456
(09:01)>sigma(%)-%
%157 = 2075868183012544
(09:01)>sigma(%)-%
%158 = 2096174602984256
(09:01)>sigma(%)-%
%159 = 2063421874812754
(09:01)>sigma(%)-%
%160 = 1043754255742634
(09:01)>sigma(%)-%
%161 = 545264138006614
(09:01)>sigma(%)-%
%162 = 336474062378666
(09:01)>sigma(%)-%
%163 = 195012689531734
(09:01)>sigma(%)-%
%164 = 97506344765870
(09:01)>sigma(%)-%
%165 = 78039653764690
(09:01)>sigma(%)-%
%166 = 67385070876590
(09:01)>sigma(%)-%
%167 = 53908056701290
(09:01)>sigma(%)-%
%168 = 56557293070742
(09:01)>sigma(%)-%
%169 = 31204023763258[/CODE]

in fact:

[url]http://factordb.com/sequences.php?se=1&eff=2&aq=1000010274060000&action=all&fr=0&to=100[/url]

proves it's not a cycle.

[CODE](10:15)>for(x=1000010274060000,1000010274060000,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x || isprime(b),print(x","y);break())))
1000010274060000,228
(10:15)>##
*** last result computed in 1,235 ms.
(10:15)>for(x=1000010274060000,1000010274060000,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x,print(x","y);break())))
(10:22)>##
*** last result computed in 12,859 ms.[/CODE]

also finds an end with alteration.

I think you have confused my checkpoint, 1000010274060000, with the amicable number, [URL="http://factordb.com/sequences.php?se=1&eff=2&aq=1000010274063430&action=all&fr=0&to=100"]1000010274063430[/URL].:smile:

[code][85486] 45 Needham 2004
1000010274063430=2*7^2*5*31*53507*1230371
1161476072988602=2*7^2*17*251*571*727*6691[/code]are known amicable ([URL]http://amicable.adsl.dk/aliquot/c2/c2_16.txt[/URL]):

[code]
? sigma(1000010274063430)-1000010274063430
%1 = 1161476072988602
? sigma(1161476072988602)-1161476072988602
%2 = 1000010274063430
?
[/code]Note: In the following examples, the lines with time are checkpoints and the line without a time is the sociable result.

The new code now works, but when I add in all the print stuff, it bogs down. Although a slight improvement, this new method does not result in a significant savings:
[code]
for(b=1000010274050000,1000010274080000,i=b;for(n=1,100,i=sigma(i)-i;if(i<b||i>b*20,break());if(i==b,print(b);write("pariCyclesFound",b);break()));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b)))
[/code]gives:
[code]
...
1000010274060000 5.6750000000000000000000000000000000000 seconds/10000
1000010274063430
1000010274070000 5.4970000000000000000000000000000000000 seconds/10000
...
[/code]compared to:
[code]
for(b=1000010274050000,1000010274080000,i=b;for(n=1,200,m=n;i=sigma(i)-i;if(i<b||i>b*20,break()));if(m>199,print(b);write("pariCyclesFound",b));if(b%10000==0,print(b, " ",gettime()/1000+0.," seconds/10000");extern("rm pariStatus");write("pariStatus",b)))[/code]which gives:
[code]
...
1000010274060000 5.6780000000000000000000000000000000000 seconds/10000
1000010274063430
1000010274070000 5.5030000000000000000000000000000000000 seconds/10000
...
[/code]Actually, the loss is probably not directly the print statements, but the checks I do for whether to print. I also have received some info via PM that may help with my pariStatus reading/writing. I have to study it a bit. But, I have some other issues to attend to for now.

Thanks everyone for the assistance...

[QUOTE=EdH;275997]I think you have confused my checkpoint, 1000010274060000, with the amicable number, [URL="http://factordb.com/sequences.php?se=1&eff=2&aq=1000010274063430&action=all&fr=0&to=100"]1000010274063430[/URL].:smile:
[/QUOTE]

yeah probably but you should still find it with mine:

[CODE](10:29)>for(x=1000010274063430,1000010274063430,b=x;for(y=1,10000000,if((sigma(b)-b) != 0,b=sigma(b)-b);if(b==x,print(x","y);break())))
1000010274063430,2[/CODE]

[QUOTE=EdH;275857]II'd prefer overwriting the pariStatus file, but couldn't figure out how to do that with PARI, so you can see the extern("rm pariStatus") command to get rid of the previous file prior to writing the current value. This is a time waster and is definitely hazardous if a power loss is timed just right, but I couldn't quite figure out a better way.[/QUOTE]


GP doesn't have much in the way of I/O. I suggest
[code]extern("rm oldStatus; mv pariStatus oldStatus");
write("pariStatus", ...);
extern("rm oldStatus");[/code]