24 lines
514 B
Mathematica
24 lines
514 B
Mathematica
N = 2e4
|
|
X = [(1:N)' zeros(N,1)]
|
|
|
|
X(find(members(X(:,1), primes(N)) == 1), :) = []
|
|
|
|
for i = 1:size(X,1)
|
|
x = X(i,1)
|
|
p = 1:sqrt(x)
|
|
q = x./p
|
|
p = unique([p(find(q == round(q)))' q(find(q == round(q)))'])
|
|
X(i,2) = sum(p(p<x))
|
|
end
|
|
|
|
X(find(sum(X,2) == 2*X(:,1)), :) = [] // delete perfect numbers
|
|
|
|
for i = 1:size(X,1)
|
|
for j = find(X(:,2) == X(i,1))
|
|
if X(i,2) == X(j,1) then
|
|
amicables($+1,:) = X(i,:)
|
|
end
|
|
end
|
|
end
|
|
|
|
amicables = unique(gsort(amicables, 'c', 'i'), 'r')
|