mersenneforum.org

mersenneforum.org (https://www.mersenneforum.org/index.php)
-   PARI/GP (https://www.mersenneforum.org/forumdisplay.php?f=155)
-   -   PARI's commands (https://www.mersenneforum.org/showthread.php?t=13636)

science_man_88 2014-10-27 22:07

[QUOTE=EdH;386269]I'm not following a lot of what you're posting, but do have a query:

How does 0 mod 12 equate to (0 mod 4 || 0 mod 3)? (Or, am I too far removed from your code?)

Won't 0 mod 12 miss numbers like 9950971671152, which is 0 mod 4, 2 mod 3 and 8 mod 12?

Thanks for all the effort you' re putting in, though...[/QUOTE]

I broke the code down into two codes one handles those that are 0 mod 4 the other 0 mod 6 (Okay technically 6 mod 12 at the time, the combo of 0 mod 2 and 0 mod 3)

but right now the 0 mod 4 code has 3 times the numbers to test until I get the 0 mod 12 numbers over to the 0 mod 6 code and remove it from the need to be checked by the 0 mod 4 code. I also pointed out that to be abundant a number may have a divisor that is either perfect of abundant so I've made a script to test that before going on to test sigma(b)-b ( still ironing things out) only the ones that fail that script need to be tested to show they are abundant before proceeding with d=2 , I'm also thinking of further breaking down the findcyc4 script to avoid an if elseif like construct to find e mod 12 ( because once I move the other check over to the findcyc6 script inputs that are 0 mod 12 don't need to be considered on the findcyc4 script).[CODE]isaorp(x)=my(t=0);forstep(y=#x-1,2,-1,if(!isprimepower(x[y]) && vecsum(divisors(x[y]))>=x[y],return(1),t+=1));if(t==#x-2,return(0))[/CODE]

science_man_88 2014-10-28 21:15

okay so with some adjustment ( saw a code CRG posted in a sequence in the OEIS that helped as a new start):

[CODE]isaorp(n) = fordiv(n,d,if(5<d && d<sqrtint(n) && !isprimepower(d) && sigma(d)>=2*d,return(1)))[/CODE]

which along with other code alterations brought my attempts at both altering the code from scratch and altering my own versions of the code down from about 3 minutes down to under a minute again. don't know if I'd call isaorp an efficient sigma replacement but I did the math and if you assumed no overlap it can eliminate the need to check after the first sigma in potentially over 20% of cases I still end up having to do sigma at the moment because the way of doing it twice in one go slows down the code because instead of just calling sigma twice in a row it uses it three times the way I originally thought about skipping it. still haven't broken down findcyc4 to the two cases there in theory you could use isaorp before at least the first few sigma because it can can still eliminate potentially 20% of results needing sigma on the larger number itself.

danaj 2014-10-29 01:15

Note that fordiv does a full factorization, then creates the full divisor list, then loops calling the code block. So breaking early doesn't save much time. Measuring, calling fordiv(n,d,break()) takes almost 2x more time than sigma(n). I don't see how this would save time unless you did the checks before calling fordiv.

Your function says the number 100001024 is not abundant/perfect. Isn't it an abundant number?

science_man_88 2014-10-29 01:23

[QUOTE=danaj;386365]Note that fordiv does a full factorization, then creates the full divisor list, then loops calling the code block. So breaking early doesn't save [B]much[/B] time. Measuring, calling fordiv(n,d,break()) takes almost 2x more time than sigma(n). I don't see how this would save time unless you did the checks before calling fordiv.

Your function says the number 100001024 is not abundant/perfect. Isn't it an abundant number?[/QUOTE]

that's not what it's checking it's checking up to sqrtint(n) for abundant or perfect divisors of the number so the result of sigma of the full number doesn't need to be checked because it has properties that make it abundant. edit: sumdiv doesn't let you break out of it apparently so I used fordiv.

CRGreathouse 2014-10-29 01:50

[QUOTE=science_man_88;386169]finding out if the sum was prime seemed to take less time than the integer division it once held if I remember correct[/QUOTE]

No. Checking for primality in gp involves at least two GCD calculations, each of which take a number of integer divisions.

CRGreathouse 2014-10-29 01:56

[QUOTE=danaj;386365]Note that fordiv does a full factorization, then creates the full divisor list, then loops calling the code block. So breaking early doesn't save much time.[/QUOTE]

Right. If there was some good reason to do it, you should first factor and then use the factored result:

[code]f=factor(n);
fordiv(f, d,
if (blah(n) && foo(n) && bar(n) && sigma(d) >= 2*n, return(1))
);
sigma(f) ...[/code]

At least you don't have to factor it twice. But it would be better to code your own routine that uses early-abort factorization when the unfactored part is small enough, and the factor bound large enough, to prove that the number must be deficient, or else that the factored part alone suffices to show that the number is abundant. This should be orders of magnitude faster on average.

danaj 2014-10-29 02:53

[QUOTE=science_man_88;386366]that's not what it's checking it's checking up to sqrtint(n) for abundant or perfect divisors of the number so the result of sigma of the full number doesn't need to be checked because it has properties that make it abundant. edit: sumdiv doesn't let you break out of it apparently so I used fordiv.[/QUOTE]

But the number *is* abundant, and the function says it is not. We should continue looking at this number, but isaorp says we don't need to. Either the function is wrong or the name is wrong (because it is returning false when given input that "is abundant or perfect"). Assuming we are asking if n is abundant or perfect, then calling sigma(n) is *faster* than all the checks and early returns inside fordiv, because the setup for fordiv takes more time than all the work inside sigma. [By the way, it turns out that the ntheory library does the same thing for its divisor_sum and fordivisors loop, so no differences there. Not to mention the annoying restriction of no good way to early exit fordivisors. I have solutions for that but they're very inelegant in implementation.]

Based on EdH's program output, if the number is abundant, then we need the result of sigma(n). It seems what we really want is a fast "is this deficient" function, since that is the case where we get to exit early. Since the program limits the search to 4, I suppose a fast "is this perfect" could be used for the last call.

Charles's suggestion of inserting logic inside the factor routine would work since it could either return sigma(n) or truly exit early for deficient numbers (or both abundant and deficient for the last call).

Would it be possible to go from the other direction and see if we could skip the last sigma by constructing which abundant numbers have a sigma equal to 2x one of the perfect numbers? That may be completely infeasible.

As EdH mentioned, perhaps a new thread should be made for this, unless it is of enough interest regarding the Pari commands.

CRGreathouse 2014-10-29 05:29

I'm fine with a new thread.

I think the 'right' function would be one that takes n and k and returns -1, 0, or 1 based on whether sigma_{-1} (n) is less than, equal to, or greater than k, with early exit when possible.

science_man_88 2014-10-29 12:01

[QUOTE=danaj;386372]But the number *is* abundant, and the function says it is not. We should continue looking at this number, but isaorp says we don't need to. Either the function is wrong or the name is wrong (because it is returning false when given input that "is abundant or perfect"). Assuming we are asking if n is abundant or perfect, then calling sigma(n) is *faster* than all the checks and early returns inside fordiv, because the setup for fordiv takes more time than all the work inside sigma. [By the way, it turns out that the ntheory library does the same thing for its divisor_sum and fordivisors loop, so no differences there. Not to mention the annoying restriction of no good way to early exit fordivisors. I have solutions for that but they're very inelegant in implementation.]

Based on EdH's program output, if the number is abundant, then we need the result of sigma(n). It seems what we really want is a fast "is this deficient" function, since that is the case where we get to exit early. Since the program limits the search to 4, I suppose a fast "is this perfect" could be used for the last call.

Charles's suggestion of inserting logic inside the factor routine would work since it could either return sigma(n) or truly exit early for deficient numbers (or both abundant and deficient for the last call).

Would it be possible to go from the other direction and see if we could skip the last sigma by constructing which abundant numbers have a sigma equal to 2x one of the perfect numbers? That may be completely infeasible.

As EdH mentioned, perhaps a new thread should be made for this, unless it is of enough interest regarding the Pari commands.[/QUOTE]

maybe I should of named it isaorpdiv there are abundant numbers that have all divisors but themselves deficient it was a code to calculate them that I started with as they would be one of the lower divisors that could help exit the loop. to go from the other direction you would need to have all intermediate results in theory it's possible using a function like ali or aligen that I made earlier ( admittedly I may have to look them up.) but that doesn't really avoid a sigma and it involves partitions.

danaj 2014-10-31 07:07

I noticed this tonight. Works in 2.6.0, broken in 2.6.1 - 2.8 and tonight's git master. Looks like the commit that introduced it was 3bb8081166c9ac6915977323a4751669ae6c9ed0 in July 2013.

[code](23:05) gp > ispower(-167^10)
*** at top-level: ispower(-167^10)
*** ^----------------
*** ispower: domain error in mplog: argument <= 0
*** Break loop: type 'break' to go back to GP prompt
break> [/code]which ought to return 5.

LaurV 2014-10-31 12:46

I think is not broken, just fixed.
The powering has a higher priority than the unary minus.
It totally makes sense when you write "a^n-b^n", and it is arguable when you write just "-b^n", but to make it right, "(-b)^n" is more clear.


All times are UTC. The time now is 06:55.

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.