Central Delannoy numbers

I am solving this problem on SPOJ Upper Right King (Hard):

There is a king in the lower left corner of the $n \times n$ chess board. The king can move one step right, one step up or one step up-right. How many ways are there for him to reach the upper right corner of the board? Give the answer modulo 1000003.

The numbers to be output are the central Delannoy numbers.

I've written a solution in Python with a time complexity of $O(n(\log m)^2)$. But it's showing "time limit exceeded". Can anyone help me out with this?

import math
MAX = 10000002
MOD = 1000003

# Extended euclidean algorithm for modular inverse

def egcd(a, b):
if a == 0:
return (b, 0, 1)
else:
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

def modinv(a, m):
g, x, y = egcd(a, m)
if g != 1:
return 0
else:
return x % m

t = raw_input()
t = long(t)

fact = []

for i in range(0,MAX):
fact.append(0)
fact=1
fact=1
for i in range(2,MAX):
fact[i] = (i%MOD * fact[i-1]%MOD)%MOD

while t:
n = raw_input()
n = long(n)
n = n-1
Sum = 1
prod = ((n+1)%MOD * (n)%MOD)%MOD
increment = n+1
decrement = n
Sum = (Sum%MOD + prod%MOD)%MOD

for k in range(0,n-1):
temp_prod = (prod%MOD * (increment+1)%MOD * (decrement-1)%MOD)%MOD
prod = temp_prod
fct = fact[k+2]
fat2 = (fct%MOD * fct%MOD)%MOD
#   print fat2,modinv(fat2,MOD)
result = (temp_prod%MOD * modinv(fat2,MOD))%MOD
Sum = (Sum%MOD + result%MOD)%MOD
increment = increment + 1
decrement = decrement - 1

print "%ld" %(Sum)
t = t-1

How can I reduce the complexity of this problem in order to have it accepted?

• See this code by Nick Hobson. – Gareth Rees Oct 31 '13 at 14:00
• The above code won't work for large constraints mentioned in problem statement. any ways, i got it accepted after so many optimizations.thanks for your time. – rohspeed Nov 1 '13 at 6:26

Flatten your method body

def egcd(a, b):
if a == 0:
return (b, 0, 1)
else:
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

can be re-written as:

def egcd(a, b):
if a == 0:
return (b, 0, 1)
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

because returning halts the function.

Use in-line 'ternary' to reduce verbosity

def modinv(a, m):
g, x, y = egcd(a, m)
if g != 1:
return 0
else:
return x % m

should become

def modinv(a, m):
g, x, y = egcd(a, m)
return 0 if g != 1 else x % m

Say something meaningful when asking for input

t = raw_input()

The user sees nothing, only an hanging interpreter.

Use a main function

1. Python code runs faster in functions
2. You can use an if name ... block to avoid execution when the module is imported.

Convert directly (minor)

Do not waste lines:

n = raw_input()
n = long(n)

should become

n = long(raw_input()) # <- also add a meaningful prompt

Follow naming conventions (minor)

You follow the naming conventions pretty well, just one imprecision:

Sum = 1

should be:

sum = 1

Use proper spacing (minor)

For example this line:

prod = ((n+1)%MOD * (n)%MOD)%MOD

is hard to read because all the symbols are squashed together.

Modularize more (important)

You defined two functions, that is very good, but you could do more to reduce the 35 lines of top level imperative code that is too much.

Usually each function (including main) should be no more than 10-15 lines of code.

Use doctests to for automatic testing

In Math code, automatically testing your functions each time you run the programme can be invaluable, just use something like:

import doctest

def egcd(a, b):
"""
Meaningful description here...

>>> egcd(4132, 434)
(2, -48, 457)
"""

if a == 0:
return (b, 0, 1)

g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

doctest.testmod()

and the code will be tested automatically.

Tell why

I can see many calculations in your code, but why? Following which algorithm, please write this inside the docstrings of your functions.

• A couple of points: 1. Adding a prompt is wrong here: SPOJ makes strict requirements on the program's output, so adding the prompt will cause the program to fail. 2. "Python code runs faster in functions" needs an explanation, and anyway it's only true of CPython, and for this problem we'll probably use PyPy instead. – Gareth Rees Apr 14 '15 at 17:09

1. Review

I won't do a detailed review of your code, because as you'll see this answer is quite long enough already, but I'll just note that I found it very hard to figure out what it was doing. I think that the formula you are using is $$D(n) = \sum_{0 ≤ k ≤ n} { n \choose k } { n + k \choose k} = \sum_{0 ≤ k ≤ n} { (n + k)(n + k - 1)\dotsm(n - k + 1) \over (k!)^2 }.$$ But this is not remotely obvious from the code. A few comments would have made a big difference to the understandability.

2. A better algorithm

The Wikipedia article on the central Delannoy numbers gives a three-term recurrence relationship: $$nD(n) = 3(2n - 1)D(n - 1) - (n - 1)D(n - 2).$$ The means that it's easy to write a function that efficiently generates the numbers:

from itertools import count

def delannoy():
"""Generate the central Delannoy numbers."""
x, y = 1, 3
for n in count(2):
yield x
x, y = y, (((6 * n - 3) * y - (n - 1) * x) // n)

However, for the SPOJ problem we need the numbers modulo $m = 1000003$, and to do this we need to be able to do division modulo $m$, which requires finding the multiplicative inverse of $n$ modulo $m$. The fastest way to do this is via the extended Euclidean algorithm, and I see that you used the implementation from the "Algorithm Implementation" wikibook. That leads to this solution:

from itertools import count, islice

def egcd(a, b):
"""Return g, x, y, where g is is greatest common divisor of a and b,
and x and y are integers such that ax + by = 1.

"""
if a == 0:
return (b, 0, 1)
else:
g, y, x = egcd(b % a, a)
return (g, x - (b // a) * y, y)

def modinv(a, m):
"""Return the multiplicative inverse of a modulo m, if it exists."""
gcd, x, y = egcd(a, m)
if gcd != 1:
raise ValueError("no inverse of {} modulo {}".format(a, m))
else:
return x % m

def delannoy(m):
"""Generate the central Delannoy numbers modulo m."""
x, y = 1, 3
for n in count(2):
yield x
x, y = y, (((6 * n - 3) * y - (n - 1) * x) * modinv(n, m)) % m

def main():
d = list(islice(delannoy(1000003), 1000000))
T = int(input())
for _ in range(T):
n = int(input())
print(d[n - 1])

if __name__ == '__main__':
main()

Unfortunately this is not fast enough: with the biggest possible test case ($T = 10000$), this takes about 1.9 seconds on my machine under PyPy, and the time limit at SPOJ is 1.183 seconds.

The only place where we can clearly save time is in the computation of the modular inverses. There are several observations that will help here:

1. We need inverses modulo $m$ for all $1 ≤ i ≤ n$, so we might as well compute them all at once.

2. Function calls are quite expensive in Python, so we should switch to the iterative version of the extended Euclidean algorithm, and inline it.

3. It happens to be the case that $m = 1000003$ is prime, so we know that modinv(a, m) will always succeed for $a < m$. So the computation of the GCD and the associated test can be omitted.

4. We can avoid computation of three-quarters of the modular inverses. That's because if we have computed $i^{-1} ≡ j \mod m$ then we immediately know \eqalign{ j^{-1} &≡ i &\pmod m \\ (-i)^{-1} &≡ -j &\pmod m \\ (-j)^{-1} &≡ -i &\pmod m } too.

Applying all of these improvements leads to this implementation:

from itertools import count, islice

def invmod(p):
"""Return a list of inverses modulo the prime number p."""
# Cache function to avoid global lookup
_divmod = divmod
inv =  * p
for i in range(1, p // 2 + 1):
if inv[i]:
continue

# Specialization of Extended Euclidean algorithm for the case
# where we know that p is prime and i < p, so the GCD is 1.
j, k, a, b = 0, 1, p, i
while b:
a, (q, b) = b, _divmod(a, b)
j, k = k, j - q * k

# Having found that j is the inverse of i, we immediately know
# the inverses of j, -i and -j.
inv[i] = j
inv[j] = i
inv[-i] = -j
inv[-j] = -i
return inv

def delannoy(n, m):
"""Return a list of the first n central Delannoy numbers, modulo m."""
inv = invmod(m)
def d():
x, y = 1, 3
for i in count(2):
yield x
x, y = y, (((6 * i - 3) * y - (i - 1) * x) * inv[i]) % m
return list(islice(d(), n))

def main():
d = delannoy(1000000, 1000003)
T = int(input())
for _ in range(T):
n = int(input())
print(d[n - 1])

if __name__ == '__main__':
main()

And this runs in just over 1 second under PyPy, so it's good to go.