Discussion between Mihir Singhal and user51819

Stack Exchange

15:41

A: $\binom{n}{k}$ modulo prime power for large $n$ and small $k$

Try using $\binom{n}{k} = \frac{n-k+1}{k} \binom{n}{k-1}$. To avoid issues with dividing by $k$, in case $k$ has an exponent of $p$, just keep track of the exponent of $p$ separately. For example have a cumulative tracker called pow. Now, if $p \mid n$, then replace $n$ with $n/p$ in the numerato...

user51819

As mentioned in the original post, iterative identities like these are prone to coprimality issues when the denominator $k$ is not invertible mod $p^a$. I am not sure what you mean by tracking exponents separately.

Mihir Singhal

That's why I said to keep track of the exponent of $p$ separately. I've elaborated on my answer.

user51819

Although I thought it was $\binom{n}{k} = \binom{n}{k-1}\frac{n-k+1}{k}$ ?

Mihir Singhal

That's also a true identity, which could also yield a valid solution.

user51819

I am sorry but I still do not understand the answer. Can you provide an example? (just to be clear, this is for prime powers, not just primes. I can use Lucas Theorem if the modulus is a single prime $p$)

I think this still might be too slow since I am computing the coefficients over a range of $k$. Performing this function for each coefficient does not do better than the Granville method.

Mihir Singhal

15:41

I have added pseudocode. Note that you can always divide by a number that is not a multiple of p (mod p^a). Hopefully you understand now.

This solution is O(k log(n + p^a)). With your bounds on n, k, p^a mentioned (10^8, 10^5, 10^8), this should easily run in under a second.

user51819

I have to compute the coefficients over many $k$ (so for example for $k=1$ to $10^5$, for some unchanging $n$), and so if max $k$ is $10^5$ then that's $O(k^2)$ runtime for $10^5$. I also tested it for speed.

Mihir Singhal

Oh, I didn't understand what you were asking for. I've fixed the code to use your identity instead; it should still work in approximately O(k) time.

user51819

I tested it and it does not seem to work (although the error may be mine). For the line "currentans * numer / denom (mod p)" I tried "currentans * numer * denom^(p-2) mod p"

Mihir Singhal

Oops. It should have said mod $p^a$. To find 1/denom efficiently, try using the Extended Euclidean Algorithm.

user51819

So with the edit to p^a I now use inverse mod (extended gcd), answers still do not seem to match the correct results (I am comparing against a slow but correct n choose k mod p^e implementation)

Mihir Singhal

15:41

I moved this discussion to chat

but

i made another mistake

should say n-i+1 instead of n-i-1

Try it now.

user51819

Same error

Mihir Singhal

what do you mean?

is there an error, is it giving the wrong answer, is it taking too long, or what?

user51819

(I am testing 41 choose k mod 3^3 over the range 0<=k<=41)

results should be 1,14, 10, 22, 20, 13, 24, 12, 24, 16, 8, 25, 22, 20, 13, 18, 9, 18, 15, 21, 15, 15, 21, 15, 18, 9, 18, 13, 20, 22, 25, 8, 16, 24, 12, 24, 13, 20, 22, 10, 14, 1

Mihir Singhal

what did you get

user51819

but i get 1, 14, 3, 19, 11, 24, 16, 20, 0, 12, 6, 9, 12, 24, 18, 3, 24, 0, 24, 12, 18, 24, 21, 9, 6, 21, 0, 22, 11, 12, 13, 26, 15, 1, 8, 0, 15, 21, 18, 15, 3, 9

Mihir Singhal

15:44

wait, are you still doing denom^(p - 2)

user51819

Mihir Singhal

because that won't work any more

what are you doing

user51819

currentAns = currentAns * numer * inverse(denom, p**a)

inverse is egcd algorithm

@ErickWong Hello, what brings you here

Mihir Singhal

can you send your code or something

Erick Wong

Just helping to debug :)

Mihir Singhal

15:45

because doing my algorithm by hand, it should work for i=2 at least

put it on pastebin

user51819

pastebin.com/raw.php?i=LEh9WvgL

Erick Wong

Why is there a k-=1, n-=1 at the end?

Mihir Singhal

First of all, in the fifth to last line, you should have currentAns = (currentAns * numer * inverse(denom,p**a)) % p ** a

that's the problem

yes, there should not be an n -= 1

user51819

output matches without n-=1

wow!

Mihir Singhal

but for efficiency, also add the % p ** a at the end of the fifth to last line

otherwise you'll be dealing with very large numbers

Erick Wong

15:49

the k-=1 doesn't make sense either (but is harmless because k isn't used anywhere inside the loop)

Mihir Singhal

yeah, those are remnants from my previous answer that didn't work

erick is right

user51819

i wish i understood how this worked! this is a nice approach

accepted/upvoted answer

may be a small problem when n=1

because it will get stuck on the numer mod p = 0 line

Mihir Singhal

When n=1, then k is at most 1 as well (and thus i is at most 1), so the lowest numer will be is 1, so that should not be a problem. Try it for yourself.

user51819

16:04

nevermind it seems the real problem is when k>n and numer goes negative i think

i can just add early break

Mihir Singhal

16:15

Yes, I assumed that k ≤ n in my solution. It's easy to fix nonetheless.

last day (15 days later) »

all rooms