Find all integers $$k\le100$$, so that there exists an integer $$n$$, satisfying $k\mid3n^6+26n^4+33n^2+1.$

By number theory knowledge it suffices to check $$n\in[1,k]$$, but we'll do $$[1,100]$$ for convenience.

My code, using mma, is

Select[Range[100],
Or @@ (Divisible[3 #^6 + 26 #^4 + 33 #^2 + 1, #] & /@ Range[100]) &]


But Mathematica output all numbers in Range[100]. What is wrong here?

• Note. This was an awful problem from a Beijing high school contest a few days ago. Commented May 22, 2023 at 13:31
• 1. "By number theory knowledge it suffices to check n∈[1,k]" Are you sure? With[{n = 2, k = 1}, Divisible[3 n^6 + 26 n^4 + 33 n^2 + 1, k]] is True. 2. Why do you think Select can be used in this manner? Commented May 22, 2023 at 14:36
• @xzczd when k=1, n=1 also works. I mean all n works Commented May 22, 2023 at 15:36
• @OkkesDulgerci OK, I see. Numeric experiment for $k \le 100$ also shows the statement seems to be correct. Commented May 23, 2023 at 2:41

Re: What is wrong here?

Puzzling outcome you are seeing is caused by using a nested pure function (function with anonymous arguments) as the selector function in Select.

yourSelectorFunction =
Or @@ (Divisible[3 #^6 + 26 #^4 + 33 #^2 + 1, #] & /@ Range[100]) &

 Or @@ (Divisible[3 #1^6 + 26 #1^4 + 33 #1^2 + 1, #1] &) /@ Range[100] &


This function returns True for any input passed to it:

yourSelectorFunction /@ {12, blahblah, ☺, xyz}

 {True, True, True, True}


The # in the second argument of the inner function Divisible[...] is filled with a value before the outer function is called to work. (Slots (#) are consumed by the nearest &.)

Simple rule of thumb to avoid such problems is to never use a nested pure function (an expression that contains more than a single &) in the second argument of Select. Instead, use functions with named arguments (Function[x, body] or x |-> body) instead.

Some examples:

Select[Range[100],
Function[u, Or @@ (Function[b, Divisible[3 b^6 + 26 b^4 + 33 b^2 + 1, u]] /@
Range[100])]]

{1, 3, 7, 9, 13, 19, 21, 27, 29, 31, 39, 49, 57, 59, 63, 71,
81, 83, 87, 91, 93}

Select[Range[100],
Or @@ (Function[b, Divisible[3 b^6 + 26 b^4 + 33 b^2 + 1, #]] /@
Range[100]) &]

{1, 3, 7, 9, 13, 19, 21, 27, 29, 31, 39, 49, 57, 59, 63, 71,
81, 83, 87, 91, 93}

Select[Range[100],
Function[a,
Or @@ (Divisible[3 #^6 + 26 #^4 + 33 #^2 + 1, a] & /@ Range[100])]]

{1, 3, 7, 9, 13, 19, 21, 27, 29, 31, 39, 49, 57, 59, 63, 71,
81, 83, 87, 91, 93}


An aside: Since Divisible is Listable, we can pass Range[100] into the first argument of Divisible and leave # in the second argument for the outer function to consume:

Select[Range[100],
Or @@ (Divisible[3 #^6 + 26 #^4 + 33 #^2 + 1 & @ Range[100], #]) &]

{1, 3, 7, 9, 13, 19, 21, 27, 29, 31, 39, 49, 57, 59, 63, 71,
81, 83, 87, 91, 93}

• Thank you very much. By the way, are these methods faster compared to Do[] loops? Commented May 22, 2023 at 15:58
• @youthdoo you can always write the code yourself using a Do loop and using Timing functions to test it. But as a rule of thumb: Built in loops (For, While, Do) in mathematica are much much slower than designated functions (that in practice also loop behind the scenes). Commented May 23, 2023 at 6:28

Perhaps in two steps:

knTF = Table[{k, n, Divisible[3 n^6 + 26 n^4 + 33 n^2 + 1, k]}, {k, 1,100}, {n, 1, k}] // Flatten[#, 1] &
Select[knTF, #[[-1]] &]

(*{{1, 1, True}, {3, 1, True}, {3, 2, True}, {7, 1, True}, {7, 6, True},...} *)


This can be done in 13.2.1 on Windows 10 directly.

Resolve[Exists[{n, s},  s*k == 3 n^6 + 26 n^4 + 33 n^2 + 1 && k <= 100], PositiveIntegers]


k==1||k==3||k==7||k==9||k==13||k==19||k==21||k==27||k==29||k==31||k==39||k==49||k==57||k==59||k==63||k==71||k==81||k==83||k==87||k==91||k==93

I post this answer as a hybrid between "pen and paper" and using Mathematica to solve. It uses various Mathematica functions.

Firstly, some observations:

• 1 is trivially a solution

• from inspection no even numbers can be solutions

• the solutions must be derived from primes (primes, powers of primes, or products of the preceding)

• to further constrain, note 3 x 5 x 7 =105>1, so apart from powers we need only look at products of pairs of relatively prime numbers

Set up

f[n_] := 3 n^6 + 26 n^4 + 33 n^2 + 1;
p = Prime[Range[PrimePi[100]]];
pm = PolynomialMod[f[n], #] & /@ p;


where p is prime numbers < 100 and pm reduces of polynomial.

Finding prime number solutions

atom = Cases[
MapThread[{#2, n /. Solve[#1 == 0, n, Modulus -> #2]} &, {pm, p}],
Except[{_, n}]];
basicatom = atom[[All, 1]]


The prime solutions are: {3, 7, 13, 19, 29, 31, 59, 71, 83}

Prime power solutions

Hensels lemma allows lifting of solutions $$\mod p^k$$ to $$\mod p^{k+m}$$ where $$1\leq m\leq k$$.

    hensel[{{x_, y_, z_}, c_}] := Module[{d = f'[n], bp = x^(1/c)},
{{bp^(c + 1), n /. Solve[f[n] == 0, n, Modulus -> bp^(c + 1)],
Mod[d /. n -> # & /@ y, bp]}, c + 1}
]
prep[x_] := {{x[[1]], x[[2]],
Mod[f'[n] /. n -> # & /@ x[[2]], x[[1]]]}, 1}
nwf[x_] := Most@NestWhileList[hensel, prep[x], #[[1, 1]] < 100 &]


hensel input number, solutions mod number and value of derivative of solution mod base prime and outputs lifted solution (prime power, new solution, value derivative of previous solution mod p). This is more illustrative. All that was needed was to find prime power solutions. It was really just an aside.

prep just reformats the atom entries to allow processing. 

nwf is just iterating hensel till prime power exceeds 100.

As 3 and 7 are the only primes who powers>1 remain less than 100 they are the only ones to consider. So finding, prime powers of relevance:

powers = Rest /@ (nwf /@ atom[[{1, 2}]]);


.

Combining all generators of solutions

The Chinese remainder theorem means that distinct generators where (p,q)=1 imply p q is a solution

fullatom =
Sort@Join[basicatom, powers[[1, All, 1, 1]], powers[[2, All, 1, 1]]];


Looking at prducts of pairs (as discussed in preamble):

comb = Times @@@
Select[GCD @@ # == 1 && Times @@ # < 100 &][Subsets[fullatom, {2}]]


and all solutions:

ans = Sort@Join[{1},fullatom, comb]


yields

{1,3, 7, 9, 13, 19, 21, 27, 29, 31, 39, 49, 57, 59, 63, 71, 81, 83, 87,
91, 93}

This is consistent with the solution from the other answers (which I have voted for).

I note reality checking the lack of need to check products of 3 or more generators(given constraint) as argued above:

test = Times @@@
Select[GCD @@ # == 1 && Times @@ # < 100 &][Subsets[fullatom, {3}]]


yields {}.

Why bother with this tripe? It's a mess.

I understand this and I upfront apologize for errors.

Some reasons:

• by their nature Olympiad and similar puzzles and problems are more than "get the answer" but about how to get the answer. So, in this spirit I post. I was trying to illustrate constraining search space. Of course, Mathematica can do this in almost one line.
• it was a build up approach on a smaller search space and allowed illustration of various Mathemaitca capabilities, e.g. PolynomialMod, Solve[...,Modulus->...]` etc. Some already illustrated.
• it was a nice example to combine some theory and Mathematica

If this answer, offends, is conceptually or otherwise wrong. I would delete after understanding issue.