As part of a larger program, I need to solve $A x^2+B x+C \equiv 0 \pmod {p^k}$ for prime $p$. Right now I'm doing this by calling


and then parsing this to find what I need (i.e. usually just assigning C[1]->0 and ignoring M results).

But when $A$, $B$, $C$, and especially $p^k$ start to get large (e.g. $A,B,C \approx 10^{30}$ and $p^k \approx 10^{70}$), there is quite a lag.

Is there a faster way to find solutions for $x \bmod p^k$?


  • $\begingroup$ Could try to rewrite it in a way that uses PowerMod[...,1/2,p^k]. $\endgroup$ – Daniel Lichtblau Dec 20 '13 at 16:53
  • $\begingroup$ So there's nothing that implements Hensels Lemma a.k.a lifting? I guess I'll have to write my own... $\endgroup$ – Aeryk Dec 20 '13 at 19:00
  • 1
    $\begingroup$ If I'm not mistaken PowerMod for prime power moduli will do lifting. Also Roots (and thus Solve and Reduce) will use Hensel lifting. The reason to prefer PowerMod is that if I recall correctly it has special case code for square roots. The general case uses degree 1 of the distinct degree factorization algorithm and I'm guessing that's the expensive step. $\endgroup$ – Daniel Lichtblau Dec 20 '13 at 19:54

When a is invertible and there are two solutions (the "nice" case), the code below will find one root. It seems to be around 3x faster than `Reduce[...,x, Modulus->p^k]. It does not check for various error conditions though.

solveQuadratic[a_, b_, c_, p_, k_] := Module[
  {x, ainv, cp, corr, pk = p^k, rhs},
  ainv = PowerMod[a, -1, pk];
  corr = Mod[ainv*b*PowerMod[2, -1, pk], pk];
  cp = Mod[-c + a*PowerMod[corr, 2, pk], pk];
  rhs = Mod[ainv*cp, pk];
  Mod[PowerMod[rhs, 1/2, pk] - corr, pk]

One can negate this mod p^k to get the other root.

When there are more than 2 roots this still seems to find one.

  • 2
    $\begingroup$ It appears that PowerMod[x,1/2,p] uses the Tonelli-Shanks modular square root algorithm (or a variant of this) for p prime. $\endgroup$ – John McGee Feb 27 '16 at 15:15


Please see comments by belisarius below in relation to the large $p^k$.

My original answer (below):

You can add Modulus option to Solve or Reduce. Just a simple example:

Reduce[x^2 - 4 x + 9 == 0, x, Modulus -> 3^5]


x == 103 || x == 144


Solve[x^2 - 4 x + 9 == 0, x, Modulus -> 3^5]


{{x -> 103}, {x -> 144}}
  • $\begingroup$ For large numbers you get Reduce::munsp: Reduce was unable to solve the system. Increasing the value of "MaxModularPoints" in SystemOptions["ReduceOptions"] may help. >> $\endgroup$ – Dr. belisarius Dec 20 '13 at 4:25
  • 1
    $\begingroup$ But Reduce[Mod[a*x^2 + b*x + c, p^k] == 0, x, Integers] works (although not very fast) $\endgroup$ – Dr. belisarius Dec 20 '13 at 4:27
  • $\begingroup$ @belisarius thank you for the feedback. I just wished to express note this option. I completely missed the point of behaviour for large prime multiples. So, have learned a lot from your comments... $\endgroup$ – ubpdqn Dec 20 '13 at 4:36
  • 1
    $\begingroup$ But I haven't tried to increase "MaxModularPoints". Perhaps you should, at least enough to cope with the figures mentioned by the OP $\endgroup$ – Dr. belisarius Dec 20 '13 at 4:39

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.