Perfect number verb
Sanju Dutt
Trying to write a perfect number verb using control flow structures
What needs to be altered here?
Thanks,
Sanju
```
is_perfect_number=:3 : 0
sum=.0
i=.0
n=.y
for_i. n do.
if n | i -:0
sum=.sum+i
return. sum -:n
end.
)
```
Raul Miller
So I would be tempted to implement something like this:
divisors=: {{,*/@>{(^ i.@>:)&.>/ __ q: y}}
proper_divisors=: _1 }. divisors
is_perfect_number=: = +/@proper_divisors"0
(#~ is_perfect_number) 1+i.10000
6 28 496 8128
Note that this divisors verb will throw an error if you try to use it
on a number (like 0) which is not a positive integer. To work around
that, I would be tempted to change the definition of
is_perfect_number, like this:
is_perfect_number=: = +/@proper_divisors ::6:"0
With this change, a slightly simpler example would be:
I. is_perfect_number i.10000
6 28 496 8128
I hope this makes sense to you,
--
Raul
Henry Rich
Michael Day
I was reminded of perfect numbers coming up in a pub-quiz a few years ago,
and that I'd tried to impress fellow team-members with the following:
NB. Euler proved that even perfect numbers are of form (2^(p-1)) x (2^p
- 1)
NB. if both p and 2^p - 1 are prime, ie if 2^p - 1 is a "Mersenne Prime"
NB. eg, 2 is prime, and 3 = 2^2 - 1 is prime, (J-ers note:
Conventional maths notation)
NB. 2^(2-1) is 2, so 2 x 3 is perfect
NB. Similarly, 3 and 7 = 2^3 - 1 are prime, so 28 = (2^2) x (2^3 - 1)
= 4 x 7 is perfect...
NB. Odd perfect numbers have not been found!
NB. find all perfect numbers <: y
fast_perfect =: 3 : 0 NB. exploit Euler's lemma
N =. y
l =. '' NB. seed for list of perfect numbers
p =. 1
NB. loop over all primes until N < 2^(p-1) x (2^p - 1)
while. N >: cand =. <.@*/ n =. 0 _1 + 2^ _1 0 + p =. 4 p: p do.
if. 1 p: {: n do. NB. ok if 2^p - 1 is prime
l =. l, cand
end.
end.
l NB. letter l, not modulus !
)
This is pretty fast - much faster than my "slow-perfect" fn, which
I won't bother you with. No doubt a vector approach would be
even quicker, but this is what I had lying around, ready-made.
Not quite perfect, but good enough?
ts'fast_perfect 1000'
2.53e_5 3392
ts'fast_perfect 100000000'
2.57e_5 3456
fast_perfect 100000000
6 28 496 8128 33550336
Mike
> an email to forum+un...@jsoftware.com.
>
This email has been checked for viruses by Avast antivirus software.
www.avast.com
Sanju Dutt
Michael Day
incoming
mail yet!
> an email to forum+un...@jsoftware.com.
>
Michael Day
I don't think we actually said what needed altering in Sanju's suggested function.
Here goes: I'm using fixed width, and hope the use of NB. statements
renders the following suitable for copy & paste into an edit window.
NB. First of all, an annotated listing of your verb:
is_perfect_number=:3 : 0
sum=.0
i=.0 NB. not necessary as i is (re)defined in a for_i.
loop
n=.y
for_i. n do. NB. problem - this loop has only one element, ie
n itself
if n | i -:0 NB. problems - a) "if" is not a control word -
use if. exp do. action end.
b) in J you need the test
expression to be 0 = i | n
sum=.sum+i NB. problem - this action is always performed as
it's not within an if. ... end.
return. NB. problem - this return. is a valid control
statement, but terminates before the next action
sum -:n NB. problem - this action is never performed as
it follows return.
end.
)
NB. Second, fairly minimal changes to get it to work...
sum=. 0
n =. y
for_i. >: i. <: n do. NB. need proper divisors, ie those < n itself
if. 0 = i | n do. NB. I like indenting inner loops
sum =. sum + i
end.
end.
sum = n
)
NB. A few "improvements" - a matter of
opinion, perhaps
NB. A vectorised version of ipnsd ... a fairly typical J/APL way
of avoiding an explicit for loop
ipnsd_v =: {{
y = +/ vec * 0 = y |~ vec =. >: i. <: y
}}
NB. This uses the property that if i is a
divisor of n, then so is n%i
NB. A(nother) vectorised version of ipnsd looping only up to
sqrt of n ...
ipnsd_va =: {{
if. y > 1 do.
sum =. +/ vec * ok =. 0 = y |~ vec =. >: i. <. %: y
y = sum + +/ vec * ok =. 0 = y |~ vec =. y <.@% }.
ok#vec
else.
0
end.
}}
NB. Performance of the last is similar to my previous offering, slow_perfect, except
NB. that the latter returns all perfect numbers <: n , eg
slow_perfect 1000
6 28 496
I.ipnsd_va"0 i.1000
6 28 496
NB. Here's my "slow_perfect" again, for
completeness
NB. seek perfect numbers < or = N, by brute force
slow_perfect =: 3 : 0 NB. There are better ways to do this!
(see below)
N =. y
l =. '' NB. seed for list of perfect numbers
for_n. }.>:i.N do. NB. looping for each 2 <= n
<= N
f =. }.i.>:>.%:n NB. all numbers from 1 to square
root n
f =. f#~ 0 = f|n NB. select those which are factors of
n
f =. ~.f, n%}.f NB. if f divides n, then so does n%f
l =. l, n#~ n = +/f NB. append n to l if it's perfect, ie
if n = sum(f)
end.
l
)
NB. Time & space verb:
ts =: 6!:2 , 7!:2@]
NB. test performance testing
ts'ipnsd"0 ]i.1000'
0.0416124 20256
ts'ipnsd_v"0 ]i.1000' NB. better to vectorise
0.0053995 30368
ts'ipnsd_va"0 ]i.1000'
0.0044173 14880
ts'slow_perfect 1000'
0.0040656 11104
NB. BUT ... don't try doing the following calculation with any of the above versions!
,.fast_perfect 1e326
28
496
8128
33550336
137438691328
2305843008139952128
ts'fast_perfect 1e32' NB. quite
nippy!
0.000408 5040
ts'fast_perfect 1e128' NB. and so
on!!!
0.0057781 5632
NB. Here's fast_perfect, again for completeness:
fast_perfect =: 3 : 0 NB. exploit Euler's lemma
NB. Euler proved that even perfect numbers are of form (2^(p-1)) x (2^p - 1)
NB. if both p and 2^p - 1 are prime, ie if 2^p - 1 is a
"Mersenne Prime"
N =. y
l =. '' NB. seed for list of perfect numbers
p =. 1
NB. loop over all primes until N < 2^(p-1) x (2^p - 1)
while. N >: cand =. <.@*/ n =. 0 _1 + 2x^ _1 0 + p =. 4 p: p do.
if. 1 p: {: n do. NB. ok if 2^p - 1 is prime
l =. l, cand
end.
end.
l
)
NB.This is fast already, but I did go on and solve the quadratic in 2^(p-1) for maximum candidate
NB. prime p, so the following should do very slightly less work than fast_perfect, but seems to take more space!
fp =: {{
N =. y
maxp =. 1 + 2 ^. 4 %~1 + %:1 + 8x*N NB. Max required p -
not necessarily integer or prime yet
p =. (p: @: i. @(p:inv + 1&p:)) maxp NB. generate all
candidate primes
perfectnos =. (* -:@>:) mersennes =. (#~1&p:) <:
twotop =. 2x ^ p NB. Needs extended for large perfect nos.
}}
(fast_perfect -: fp)1e128
1
ts'fp 1e128'
0.0057702 19272
Mike
To unsubscribe from this group and stop receiving emails from it, send an email to forum+un...@jsoftware.com.
