What an annoying task. There are multiple valid solutions
The serial key is taken as `argv[1]` is of the format `dddd_dddd_dddd_dddd`
Here's my z3 script

```
from z3 import *
from hashlib import md5

"""
`get_models` returns atmost M solutions
https://stackoverflow.com/questions/11867611/z3py-checking-all-solutions-for-equation
"""
def get_models(F, M):
s = Solver()
count = 0
s.add(F)
while count < M and s.check() == sat:
m = s.model()
count += 1
yield m
# Create a new constraint the blocks the current model
block = []
for d in m:
# d is a declaration
if d.arity() > 0:
raise Z3Exception("uninterpreted functions are not supported")
# create a constant from declaration
c = d()
if is_array(c) or c.sort().kind() == Z3_UNINTERPRETED_SORT:
raise Z3Exception("arrays and uninterpreted sorts are not supported")
block.append(c != m[d])
s.add(Or(block))

d = [BitVec('i%d' % i, 16) for i in xrange(19)]
#s = Solver()
s = []

for i in xrange(4):
s.append(And(d[i] > 0, d[i] <= 9))
s.append(And(d[i+5] > 1, d[i+5] < 9))
s.append(And(d[i+10] >= 0, d[i+10] <= 9))
s.append(And(d[i+15] > 2, d[i+15] < 8))

s.append(d[4+5*0] == ord('_')-0x30)
s.append(d[4+5*1] == ord('_')-0x30)
s.append(d[4+5*2] == ord('_')-0x30)

r11 = sum(d[0:4])
r13 = sum(d[5:9])
esi = sum(d[10:14])
r8 = sum(d[15:])

s.append((r11+r13+esi)/3 == r8)
s.append(esi/3 == r11)
s.append(r11 != r13)

l = []
for i in xrange(4):
l.append(d[i] != d[i+10])

s.append(And(l))

l = []
for i in xrange(4):
l.append(d[i+5] != d[i+15])

s.append(And(l))

for A in get_models(s, 10):
ans = []
for i in d:
ans.append(chr(A[i].as_long()+0x30))

x = "".join(ans)
print x
```

Here's the output

```
1612_8456_8769_6537
2262_8738_9999_6477
2262_8738_9999_6576
2262_2728_9999_4576
2262_2728_9999_4477
2262_2728_9999_6637
2262_2728_9999_6655
2262_5728_9999_6647
2262_3388_9999_6647
2222_7375_5388_5643
```

All of them are valid !

There was a hint that -
`The correct md5 starts with one of the two options: 0c79 or 98fc`
and that took 20 points !! How annoying

So, now the modified script is

```
...
for A in get_models(s, 10000000):
ans = []
for i in d:
ans.append(chr(A[i].as_long()+0x30))
x = "".join(ans)
m = md5(x).hexdigest()
if m.startswith('0c79') or m.startswith('98fc'):
print "[*] Serial: %s" % x
print "[*] MD5: %s" % m
```

The output
```
[*] Serial: 8111_2562_9997_6653
[*] MD5: 98fcc746a6ab0b2d776c331875f245b9
[*] Serial: 1118_3538_9897_5673
[*] MD5: 0c79badeb7544ce59bbb958a812e11e8
[*] Serial: 2144_2362_9789_3457
[*] MD5: 0c79842c0d7ecd5cab431da828c19134
[*] Serial: 1241_6888_9790_5646
[*] MD5: 0c79a44c4104531c569399f4811ac4a5
...
```