# Tag List is Protected Error in simple script

I want to have a simple script that gives me all numbers $k \in \mathbb{N}$ in some range such that product of all their divisors is equal to $k^3$. Here is the code that I wrote:

  findDivisors[n_] := Module[
For[i = 1, i <= n, i++,
If[Apply[Times, Divisors[i]] == i^3, Append[answer, i]]];
];


So it must check all integers from $1$ to $n$ and gives me the list "answer" which elements satisfy the desired property. But when I run it the output is

SetDelayed::write: Tag List in {}[n_] is Protected.


So where is my mistake?

• For another approach I can hint: you're looking for numbers with exactly 6 divisors. These are, for instance, the fifth power of any prime. Also the square of one prime times another one. Since the prime factorization of integers is unique, I strongly suspect that this exhausts the possibilities. – LLlAMnYP Feb 22 '17 at 22:28
• Can someone come up with a better title for this post? I couldn't, but the current title seems misleading. – Chris K Feb 23 '17 at 14:37
• @Chris you see, there's hardly anything that can be done here. OP ran into an error very specific to his problem, that was fixable by a reset (the real answer not posted is "restart your kernel") although his real problem was the difference between Append and AppendTo. And Mr.W and I went off on a tangent to solve his actual task. What do you think the question is? – LLlAMnYP Feb 24 '17 at 10:32

I get no such error from your code but I also get no useful output because you used Append where you needed AppendTo. See Assuming commands will have side effects when they don't.

ClearAll[findDivisors]

findDivisors[n_] :=
For[i = 1, i <= n, i++,
If[Apply[Times, Divisors[i]] == i^3, AppendTo[answer, i]]];
]

findDivisors[100]

{1, 12, 18, 20, 28, 32, 44, 45, 50, 52, 63, 68, 75, 76, 92, 98, 99}


A more idiomatic way to write that is:

findDivisors2[n_] :=
Table[
If[Times @@ Divisors[i] == i^3, i, ## &[]],
{i, n}
]


Or taking additional liberties with the translation:

findDivisors3[n_] := Cases[Range @ n, i_ /; Times @@ Divisors[i] == i^3]


Reference:

## Performance

Using LLlAMnYP's observation that each output number besides one has exactly six divisors here is a function targeting performance making use of DivisorSigma:

findDivisors4[n_] :=
Pick[#, DivisorSigma[0, #], 6] & @ Range @ n // Prepend[1]

Needs["GeneralUtilities"]

BenchmarkPlot[{findDivisors, findDivisors2, findDivisors3, findDivisors4}, # &]


Though this does not improve on the computational order of the other methods is it constantly three times faster.

• I even create new notebook document to test it but unfortunately it gives me the same strange error (link to screenshot: savepic.ru/13002986.png) – Hasek Feb 22 '17 at 12:13
• Thank you very much! Second version of code ("more idiomatic way") works as desired. – Hasek Feb 22 '17 at 12:17
• @Hasek Did you restart Mathematica completely? If not please try that now. – Mr.Wizard Feb 22 '17 at 12:18
• You're right, after I restart it both versions run without errors. – Hasek Feb 22 '17 at 12:22
• @Mr.W I'm really surprised you used conditional patterns, an If statement, but not Select[Times @@ Divisors[#] == #^3 &][Range[n]] – LLlAMnYP Feb 23 '17 at 11:46

Like I commented, every single number required by the OP has exactly six divisors (including one and itself, with the trivial exception of 1). Consequentially, there are only two possible options: either the number is a fifth power of a prime: $p^5$ and has divisors $p^{0..5}$ whose product is, of course, $p^{15} = {(p^5)}^3$; otherwise it is the product $p_i^2 p_j$ and has divisors $1, p_i, p_j, p_i p_j, p_i^2, p_i^2 p_j$. Their product is of course $p_i^6 p_j^3 = {(p_i^2 p_j)}^3$.

I now want to generate all possible tuples of the form {i, i, j} for i<=n and j<=n (and i != j). These kind of tasks are best met by built-in functions like Subsets, but since the first number is repeated, a customized compiled function to suit my needs should also work very well.

Code dump:

getTuples =
Compile[{{n, _Integer}},
Block[{out = Table[0, {n * (n - 1)}, {3}], i = 1, j = 1, k = 1},
For[i = 1, i <= n, i++,
For[j = 1, j < i, j++; k++,
out[[k, {1, 2}]] = Table[i, {2}]; out[[k, 3]] = j];
For[j = i + 1, j <= n, j++; k++,
out[[k, {1, 2}]] = Table[i, {2}]; out[[k, 3]] = j]
]; out]]


to actually get the required numbers we run the following:

getTuples[10] // Prime // Times @@@ # &


We lack here the fifth powers of primes though. To remedy, I define a little helper function:

getNumbers[n_] :=
Sort[{1}~Join~((Prime[Range[n]])^5)~
Join~(getTuples[n] // Prime // Times @@@ # &)]


It gets me 1000001 of such numbers in 3 seconds, however there is no guarantee that they are the first 1000001; since $p_i^2 p_j$ for $j = 1, i = 1000$ is present, but the likely much smaller case of, say, $j = 1001, i = 2$ is not.

Obligatory sample output:

getNumbers[10]
(* {1, 12, 18, 20, 28, 32, 44, 45, 50, 52, 63, 68, 75, 76, 92, 98, 99, \
116, 117, 147, 153, 171, 175, 207, 242, 243, 245, 261, 275, 325, 338, \
363, 425, 475, 507, 539, 575, 578, 605, 637, 722, 725, 833, 845, 847, \
867, 931, 1058, 1083, 1127, 1183, 1421, 1445, 1573, 1587, 1682, 1805, \
1859, 2023, 2057, 2299, 2523, 2527, 2645, 2783, 2873, 3125, 3179, \
3211, 3509, 3703, 3757, 3887, 3971, 4205, 4693, 4901, 5491, 5819, \
5887, 6137, 6647, 6877, 8303, 8381, 8993, 9251, 10051, 10469, 10933, \
14297, 15341, 15979, 16807, 19343, 161051, 371293, 1419857, 2476099, \
6436343, 20511149} *)


Observe how it fails though:

ListPlot[{getNumbers[10], getNumbers[100][[;; 101]]}, PlotRange -> {Full, {0, 1000}}]


The output of the first 19 primes matches for both lists, but then getNumbers[10]` shoots up, missing several suitable values. There may be ways to improve this.