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

Reduce[A*x^2+B*x+C==M*p^k,{x,M},Integers]

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 \mod p^k$?

Thanks!

-
Could try to rewrite it in a way that uses PowerMod[...,1/2,p^k]. –  Daniel Lichtblau Dec 20 '13 at 16:53
So there's nothing that implements Hensels Lemma a.k.a lifting? I guess I'll have to write my own... –  Aeryk Dec 20 '13 at 19:00
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. –  Daniel Lichtblau Dec 20 '13 at 19:54

EDIT

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

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

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


yields:

x == 103 || x == 144


and

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


yields:

{{x -> 103}, {x -> 144}}

-
For large numbers you get Reduce::munsp: Reduce was unable to solve the system. Increasing the value of "MaxModularPoints" in SystemOptions["ReduceOptions"] may help. >> –  belisarius Dec 20 '13 at 4:25
But Reduce[Mod[a*x^2 + b*x + c, p^k] == 0, x, Integers] works (although not very fast) –  belisarius Dec 20 '13 at 4:27
@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... –  ubpdqn Dec 20 '13 at 4:36
But I haven't tried to increase "MaxModularPoints". Perhaps you should, at least enough to cope with the figures mentioned by the OP –  belisarius Dec 20 '13 at 4:39

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.

-