1
$\begingroup$

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[
   {i, answer},
   answer = {};
   For[i = 1, i <= n, i++,
    If[Apply[Times, Divisors[i]] == i^3, Append[answer, i]]];
   Return[answer];
   ];

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?

screenshot

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – LLlAMnYP Feb 22 '17 at 22:28
  • $\begingroup$ Can someone come up with a better title for this post? I couldn't, but the current title seems misleading. $\endgroup$ – Chris K Feb 23 '17 at 14:37
  • $\begingroup$ @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? $\endgroup$ – LLlAMnYP Feb 24 '17 at 10:32
3
$\begingroup$

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.

Please try:

ClearAll[findDivisors]

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

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}, # &]

enter image description here

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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) $\endgroup$ – Hasek Feb 22 '17 at 12:13
  • $\begingroup$ Thank you very much! Second version of code ("more idiomatic way") works as desired. $\endgroup$ – Hasek Feb 22 '17 at 12:17
  • 1
    $\begingroup$ @Hasek Did you restart Mathematica completely? If not please try that now. $\endgroup$ – Mr.Wizard Feb 22 '17 at 12:18
  • $\begingroup$ You're right, after I restart it both versions run without errors. $\endgroup$ – Hasek Feb 22 '17 at 12:22
  • 1
    $\begingroup$ @Mr.W I'm really surprised you used conditional patterns, an If statement, but not Select[Times @@ Divisors[#] == #^3 &][Range[n]] $\endgroup$ – LLlAMnYP Feb 23 '17 at 11:46
1
$\begingroup$

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}}]

enter image description here

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.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.