73 lines
1.3 KiB
Text
73 lines
1.3 KiB
Text
function factor(_n) {
|
|
n = _n
|
|
a = J(14, 2, .)
|
|
i = 0
|
|
if (mod(n, 2)==0) {
|
|
j = 0
|
|
while (mod(n, 2)==0) {
|
|
j++
|
|
n = n/2
|
|
}
|
|
i++
|
|
a[i,1] = 2
|
|
a[i,2] = j
|
|
}
|
|
for (k=3; k*k<=n; k=k+2) {
|
|
if (mod(n, k)==0) {
|
|
j = 0
|
|
while (mod(n, k)==0) {
|
|
j++
|
|
n = n/k
|
|
}
|
|
i++
|
|
a[i,1] = k
|
|
a[i,2] = j
|
|
}
|
|
}
|
|
if (n>1) {
|
|
i++
|
|
a[i,1] = n
|
|
a[i,2] = 1
|
|
}
|
|
return(a[1::i,.])
|
|
}
|
|
|
|
function sumdigits(_n) {
|
|
n = _n
|
|
for (s=0; n>0; n=floor(n/10)) s = s+mod(n,10)
|
|
return(s)
|
|
}
|
|
|
|
function smith(n) {
|
|
a = J(n, 1, .)
|
|
i = 0
|
|
for (j=2; j<=n; j++) {
|
|
f = factor(j)
|
|
m = rows(f)
|
|
if (m>1 | f[1,2]>1) {
|
|
s = 0
|
|
for (k=1; k<=m; k++) s = s+sumdigits(f[k,1])*f[k,2]
|
|
if (s==sumdigits(j)) a[++i] = j
|
|
}
|
|
}
|
|
return(a[1::i])
|
|
}
|
|
|
|
a = smith(10000)
|
|
n = rows(a)
|
|
n
|
|
376
|
|
|
|
a[1::10]'
|
|
|
|
1 2 3 4 5 6 7 8 9 10
|
|
+-------------------------------------------------------------+
|
|
1 | 4 22 27 58 85 94 121 166 202 265 |
|
|
+-------------------------------------------------------------+
|
|
|
|
a[n-9::n]'
|
|
|
|
1 2 3 4 5 6 7 8 9 10
|
|
+-----------------------------------------------------------------------+
|
|
1 | 9843 9849 9861 9880 9895 9924 9942 9968 9975 9985 |
|
|
+-----------------------------------------------------------------------+
|