Need help with code for number theory problem

I'm completely new to Mathematica (used previously only for very simple cases). I need to write a quite complex function. The function must do the following:

1. Input consists of two numbers: a and b.

2. Factorization of expression a * (27*b^2 + a^6) must be performed, i.e. FactorInteger

3. For each prime factor powered by its exponent from the factorization of prevoius step the following procedure must be performed:

Reduce[{(x + a)^3 == -b, x^3 == -b}, x, Modulus -> f],


where f is that powered prime factor.

Basically, it solves the system of two cubic congruences by modulo prime power.

4. That solutions for prime factors must be placed in ChineseRemainder and output of it should be returned.

For example, for a = 17 and b = 697 function should return 5675. I can do it step by step but can't combine it in function.

EDIT. Pseudocode (Python-like + Mathematica):

def MyFunc(a, b):
T = []
F = FactorInteger[a * (27*b^2 + a^6)]
# assuming F is of form [(p1,e1), (p2,e2),...]
for f in F:
p = f^f
# assuming Reduce returns 'x -> <num_value>'
Reduce[{(x + a)^3 == -b, x^3 == -b}, x, Modulus -> p]
T.append(tuple(<numvalue>,p))
return ChineseRemainder[T]

• Please provide all steps, since you have already written them. Give us a script that works, and we can help you turn it into a function. – C. E. May 14 '15 at 14:09
• Your question reads as if you are looking for someone to do your homework for you. We don't do that. Especially as you don't provide complete specification. We are willing to help with specific coding issues expressed in Mathematica code – m_goldberg May 14 '15 at 14:11
• It is not a homework. OK, if so, I will try it on my own. I just need two constructs - 'for' loop for primes in factorization and a way of getting numerical value from Reduce. – user144765 May 14 '15 at 14:15
• I cannot get your example to work out. For a=17 and b=697, a*(27*b^2 + a^6) is 633325004. This has 4 as one factor. Solve[{(x + a)^3 == -b, x^3 == -b} /. {a -> 17, b -> 697}, x, Modulus -> 4] claims there is no solution. – Daniel Lichtblau May 14 '15 at 16:28
• Oh man! I wasted so much time trying to figure out why my function wasn't getting the result in your example, and then you had it wrong! )-: – MarcoB May 14 '15 at 16:42

Walk through this a step at a time with each line of input followed by the line of (* output *)

a = 17; b = 697; F = FactorInteger[a*(27*b^2 + a^6)]
(* {{2, 2}, {13, 1}, {17, 3}, {37, 1}, {67, 1}} *)

v = Map[First[#]^Last[#] &, F]
(* {4, 13, 4913, 37, 67} *)

(* When it returns {ReplaceAll[x],p} there was no solution *)
{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x, Modulus -> 4]], 4}
(* {ReplaceAll[x], 4} *)

(* When it returns {x,p} there was a solution *)
{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x, Modulus -> 13]], 13}
(* {7, 13} *)

{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x, Modulus -> 4913]], 4913}
(* {ReplaceAll[x], 4913} *)

{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x, Modulus -> 37]], 37}
(* {14, 37} *)

{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x, Modulus -> 67]], 67}
(* {47, 67} *)

(* Now do all those steps at once *)
Map[{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x, Modulus->#]],#}&, v]
(* {{ReplaceAll[x],4}, {7,13}, {ReplaceAll[x],4913}, {14,37}, {47,67}} *)

(* Eliminate the cases where there was no solution *)
r = DeleteCases[Map[{x /. ToRules[Reduce[{(x + a)^3 == -b, x^3 == -b}, x,
Modulus -> #]], #} &, v], {ReplaceAll[x], _}]
(* {{7, 13}, {14, 37}, {47, 67}} *)

(* Regroup to create ChineseRemainder input *)
Transpose[r]
(* {{7, 14, 47}, {13, 37, 67}} *)

ChineseRemainder[{7, 14, 47}, {13, 37, 67}]
(* 5675 BUT That doesn't match what you say it should be *)
(* If you can point out the error then perhaps it can be fixed *)

(* Now wrap it all up as a single function *)
MyFunc[a_, b_] := Module[{F, v, r},
F = FactorInteger[a*(27*b^2 + a^6)];
v = Map[First[#]^Last[#] &, F];
r = DeleteCases[Map[{x/.ToRules[Reduce[{(x+a)^3==-b, x^3==-b}, x,
Modulus -> #]], #} &, v], {ReplaceAll[x], _}];
ChineseRemainder @@ Transpose[r]
];
MyFunc[17, 697]
(* 5675 *)

• Beat me to it. But I got the same error in the output. – N.J.Evans May 14 '15 at 16:03
• Thanks! This is what I wanted! Your code is fine, the error was due my miscalculations. – user144765 May 14 '15 at 16:30
• Yep, was working on the same thing but couldn't figure out why I didn't get the same value that the OP posted, and it turns out that it was a mistake on his part... Oh well. I think I had come up with something along the same lines, but my solution below uses Solve rather than Reduce, which may allow one to skip a few steps. – MarcoB May 14 '15 at 16:46

This seems to accomplish what you wanted, but it's not returning your result. I have commented the code a bit more heavily than usual to show how I put it together and what it accomplishes: hopefully this may be helpful to you in constructing similar functions in the future.

function[a_, b_] :=
Module[
{factorlist, listofsolutions, toCR},

(* Generate the list of prime factors *)
factorlist = FactorInteger[a (27 b^2 + a^6)];
Print["List of prime factors: ", factorlist]; (* Debug: remove if not needed *)

(* Solve the modulo equations, and discard entries with no solution *)
(* Also keep track of the prime factor corresponding to each solution *)
listofsolutions =
Select[
{x /. Flatten@ Solve[{(x + a)^3 == -b, x^3 == -b}, x, Modulus -> (#[]^#[]) ],
#[]^#[] }& /@ factorlist,
NumberQ[ #[] ]& (* Select only those entries for which Solve had returned a number*)
];
Print["List of solutions: ", listofsolutions]; (* Debug: remove if not needed *)

(* Transpose the list and feed each sublist to ChineseRemainder *)
(* Since this is the last expression in the Module,             *)
(* its values will be returned by the function                  *)
toCR = Transpose@listofsolutions;
ChineseRemainder[ #[], #[] ]& @toCR
]


I have left in the Print statements that I used for debugging as I was writing the code, to see what was going on internally; if you don't need that kind of output, you can just delete those lines from the code.