Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
    
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
add comment

2 Answers

EDIT

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]

yields:

x == 103 || x == 144

and

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

yields:

{{x -> 103}, {x -> 144}}
share|improve this answer
    
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
1  
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
1  
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
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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