text:
Paint by Numbers
You're given a true colour image. Your task is to generate a version of this image, which looks like it was painted using paint-by-numbers (the children's activity, not nonograms). Along with the image, you're given two parameters: P, the maximum size of the colour palette (i.e. the maximum number of distinct colours to use), and N, the maximum number of cells to use. Your algorithm does not have to use all P colours and N cells, but it must not use more than that. The output image should have the same dimensions as the input.
A cell is defined as a contiguous area of pixels which all have the same colour. Pixels touching only at a corner are not considered contiguous. Cells may have holes.
In short, you are to approximate the input image with only N flat-shaded/solid-colour areas and P different colours.
Just to visualise the parameters, here is a very simple example (for no particular input image; showing off my mad Paint skills). The following image has P = 6 and N = 11:

Here are a few images to test your algorithm on (mostly our usual suspects). Click the pictures for larger versions.

Please include a handful of results for different parameters. If you want to show a large number of results, you can create a gallery over on imgur.com, to keep the size of the answers reasonable. Alternatively, put thumbnails in your post and make them links to larger images, like I did above. Also, feel free to use other test images, if you find something nice.
I assume that parameters around N ≥ 500, P ~ 30 would be similar to real paint-by-number templates.
This is a popularity contest, so the answer with the most net votes wins. Voters are encouraged to judge answers by

how well the original images are approximated.
how well the algorithm works on different kinds of images (paintings are probably generally easier than photographs).
how well the algorithm works with very restrictive parameters.
how organic/smooth the cell shapes look.

I will use the following Mathematica script, to validate results:
image = <pastedimagehere> // ImageData;
palette = Union[Join @@ image];
Print["P = ", Length@palette];
grid = GridGraph[Reverse@Most@Dimensions@image];
image = Flatten[image /. Thread[palette -> Range@Length@palette]];
Print["N = ",
 Length@ConnectedComponents[
   Graph[Cases[EdgeList[grid],
     m_ <-> n_ /; image[[m]] == image[[n]]]]]]

Sp3000 was kind enough to write a verifier in Python 2 using PIL, which you find at this pastebin.

Not the most efficient thing, but here's a Python 2 PIL verifier.

What a lovely question but I was hoping we would also see the proper "paint by numbers " version. That is with numbers in place so I could use the answers :)

@Lembik I originally wanted to include that, but I felt that it distracted from the interesting part of the question. It shouldn't be too hard to take the output of one of the submissions and convert into a template, though.

This is a fascinating post. Has anyone gone the extra step of adding the color numbers like an actual Paint by Number?

Python 2 with PIL (Gallery)
from __future__ import division
from PIL import Image
import random, math, time
from collections import Counter, defaultdict, namedtuple

"""
Configure settings here
"""

INFILE = "spheres.png"
OUTFILE_STEM = "out"
P = 30
N = 300
OUTPUT_ALL = True # Whether to output the image at each step

FLOOD_FILL_TOLERANCE = 10
CLOSE_CELL_TOLERANCE = 5
SMALL_CELL_THRESHOLD = 10
FIRST_PASS_N_RATIO = 1.5
K_MEANS_TRIALS = 30
BLUR_RADIUS = 2
BLUR_RUNS = 3

"""
Color conversion functions
"""

X = xrange

# http://www.easyrgb.com/?X=MATH
def rgb2xyz(rgb):
 r,g,b=rgb;r/=255;g/=255;b/=255;r=((r+0.055)/1.055)**2.4 if r>0.04045 else r/12.92
 g=((g+0.055)/1.055)**2.4 if g>0.04045 else g/12.92;b=((b+0.055)/1.055)**2.4 if b>0.04045 else b/12.92
 r*=100;g*=100;b*=100;x=r*0.4124+g*0.3576+b*0.1805;y=r*0.2126+g*0.7152+b*0.0722
 z=r*0.0193+g*0.1192+b*0.9505;return(x,y,z)
def xyz2lab(xyz):
 x,y,z=xyz;x/=95.047;y/=100;z/=108.883;x=x**(1/3)if x>0.008856 else 7.787*x+16/116
 y=y**(1/3)if y>0.008856 else 7.787*y+16/116;z=z**(1/3)if z>0.008856 else 7.787*z + 16/116
 L=116*y-16;a=500*(x-y);b=200*(y-z);return(L,a,b)
def rgb2lab(rgb):return xyz2lab(rgb2xyz(rgb))
def lab2xyz(lab):
 L,a,b=lab;y=(L+16)/116;x=a/500+y;z=y-b/200;y=y**3 if y**3>0.008856 else(y-16/116)/7.787
 x=x**3 if x**3>0.008856 else (x-16/116)/7.787;z=z**3 if z**3>0.008856 else(z-16/116)/7.787
 x*=95.047;y*=100;z*=108.883;return(x,y,z)
def xyz2rgb(xyz):
 x,y,z=xyz;x/=100;y/=100;z/=100;r=x*3.2406+y*-1.5372+z*-0.4986
 g=x*-0.9689+y*1.8758+z*0.0415;b=x*0.0557+y*-0.2040+z*1.0570
 r=1.055*(r**(1/2.4))-0.055 if r>0.0031308 else 12.92*r;g=1.055*(g**(1/2.4))-0.055 if g>0.0031308 else 12.92*g
 b=1.055*(b**(1/2.4))-0.055 if b>0.0031308 else 12.92*b;r*=255;g*=255;b*=255;return(r,g,b)
def lab2rgb(lab):rgb=xyz2rgb(lab2xyz(lab));return tuple([int(round(x))for x in rgb])

"""
Stage 1: Read in image and convert to CIELAB
"""

total_time = time.time()

im = Image.open(INFILE)
width, height = im.size

if OUTPUT_ALL:
  im.save(OUTFILE_STEM + "0.png")
  print "Saved image %s0.png" % OUTFILE_STEM

def make_pixlab_map(im):
  width, height = im.size
  pixlab_map = {}

  for i in X(width):
    for j in X(height):
      pixlab_map[(i, j)] = rgb2lab(im.getpixel((i, j)))

  return pixlab_map

pixlab_map = make_pixlab_map(im)

print "Stage 1: CIELAB conversion complete"

"""
Stage 2: Partitioning the image into like-colored cells using flood fill
"""

def d(color1, color2):
  return (abs(color1[0]-color2[0])**2 + abs(color1[1]-color2[1])**2 + abs(color1[2]-color2[2])**2)**.5

def neighbours(pixel):
  results = []

  for neighbour in [(pixel[0]+1, pixel[1]), (pixel[0]-1, pixel[1]),
            (pixel[0], pixel[1]+1), (pixel[0], pixel[1]-1)]:

    if 0 <= neighbour[0] < width and 0 <= neighbour[1] < height:
      results.append(neighbour)

  return results

def flood_fill(start_pixel):
  to_search = {start_pixel}
  cell = set()
  searched = set()
  start_color = pixlab_map[start_pixel]

  while to_search:
    pixel = to_search.pop()

    if d(start_color, pixlab_map[pixel]) < FLOOD_FILL_TOLERANCE:
      cell.add(pixel)
      unplaced_pixels.remove(pixel)

      for n in neighbours(pixel):
        if n in unplaced_pixels and n not in cell and n not in searched:
          to_search.add(n)

    else:
      searched.add(pixel)

  return cell

# These two maps are inverses, pixel/s <-> number of cell containing pixel
cell_sets = {}
pixcell_map = {}
unplaced_pixels = {(i, j) for i in X(width) for j in X(height)}

while unplaced_pixels:
  start_pixel = unplaced_pixels.pop()
  unplaced_pixels.add(start_pixel)
  cell = flood_fill(start_pixel)

  cellnum = len(cell_sets)
  cell_sets[cellnum] = cell

  for pixel in cell:
    pixcell_map[pixel] = cellnum

print "Stage 2: Flood fill partitioning complete, %d cells" % len(cell_sets)

"""
Stage 3: Merge cells with less than a specified threshold amount of pixels to reduce the number of cells
     Also good for getting rid of some noise
"""

def mean_color(cell, color_map):
  L_sum = 0
  a_sum = 0
  b_sum = 0

  for pixel in cell:
    L, a, b = color_map[pixel]
    L_sum += L
    a_sum += a
    b_sum += b

  return L_sum/len(cell), a_sum/len(cell), b_sum/len(cell)

def remove_small(cell_size):
  if len(cell_sets) <= N:
    return

  small_cells = []

  for cellnum in cell_sets:
    if len(cell_sets[cellnum]) <= cell_size:
      small_cells.append(cellnum)

  for cellnum in small_cells:
    neighbour_cells = []

    for cell in cell_sets[cellnum]:
      for n in neighbours(cell):
        neighbour_reg = pixcell_map[n]

        if neighbour_reg != cellnum:
          neighbour_cells.append(neighbour_reg)

    closest_cell = max(neighbour_cells, key=neighbour_cells.count)

    for cell in cell_sets[cellnum]:
      pixcell_map[cell] = closest_cell

    if len(cell_sets[closest_cell]) <= cell_size:
      small_cells.remove(closest_cell)

    cell_sets[closest_cell] |= cell_sets[cellnum]
    del cell_sets[cellnum]

    if len(cell_sets) <= N:
      return

for cell_size in X(1, SMALL_CELL_THRESHOLD):
  remove_small(cell_size)

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for cellnum in cell_sets:
    cell_color = mean_color(cell_sets[cellnum], pixlab_map)

    for pixel in cell_sets[cellnum]:
      frame_im.putpixel(pixel, lab2rgb(cell_color))

  frame_im.save(OUTFILE_STEM + "1.png")
  print "Saved image %s1.png" % OUTFILE_STEM

print "Stage 3: Small cell merging complete, %d cells" % len(cell_sets)

"""
Stage 4: Close color merging
"""

cell_means = {}

for cellnum in cell_sets:
  cell_means[cellnum] = mean_color(cell_sets[cellnum], pixlab_map)

n_graph = defaultdict(set)

for i in X(width):
  for j in X(height):
    pixel = (i, j)
    cell = pixcell_map[pixel]

    for n in neighbours(pixel):
      neighbour_cell = pixcell_map[n]

      if neighbour_cell != cell:
        n_graph[cell].add(neighbour_cell)
        n_graph[neighbour_cell].add(cell)

def merge_cells(merge_from, merge_to):
  merge_from_cell = cell_sets[merge_from]

  for pixel in merge_from_cell:
    pixcell_map[pixel] = merge_to

  del cell_sets[merge_from]
  del cell_means[merge_from]

  n_graph[merge_to] |= n_graph[merge_from]
  n_graph[merge_to].remove(merge_to)

  for n in n_graph[merge_from]:
    n_graph[n].remove(merge_from)

    if n != merge_to:
      n_graph[n].add(merge_to)

  del n_graph[merge_from]

  cell_sets[merge_to] |= merge_from_cell
  cell_means[merge_to] = mean_color(cell_sets[merge_to], pixlab_map)

# Go through the cells from largest to smallest. Keep replenishing the list while we can still merge.
last_time = time.time()
to_search = sorted(cell_sets.keys(), key=lambda x:len(cell_sets[x]), reverse=True)
full_list = True

while len(cell_sets) > N and to_search:
  if time.time() - last_time > 15:
    last_time = time.time()
    print "Close color merging... (%d cells remaining)" % len(cell_sets)

  while to_search:
    cellnum = to_search.pop()
    close_cells = []

    for neighbour_cellnum in n_graph[cellnum]:
      if d(cell_means[cellnum], cell_means[neighbour_cellnum]) < CLOSE_CELL_TOLERANCE:
        close_cells.append(neighbour_cellnum)

    if close_cells:
      for neighbour_cellnum in close_cells:
        merge_cells(neighbour_cellnum, cellnum)

        if neighbour_cellnum in to_search:
          to_search.remove(neighbour_cellnum)

      break

  if full_list == True:
    if to_search:
      full_list = False

  else:
    if not to_search:
      to_search = sorted(cell_sets.keys(), key=lambda x:len(cell_sets[x]), reverse=True)
      full_list = True

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for cellnum in cell_sets:
    cell_color = cell_means[cellnum]

    for pixel in cell_sets[cellnum]:
      frame_im.putpixel(pixel, lab2rgb(cell_color))

  frame_im.save(OUTFILE_STEM + "2.png")
  print "Saved image %s2.png" % OUTFILE_STEM

print "Stage 4: Close color merging complete, %d cells" % len(cell_sets)

"""
Stage 5: N-merging - merge until <= N cells
     Want to merge either 1) small cells or 2) cells close in color
"""

# Weight score between neighbouring cells by 1) size of cell and 2) color difference
def score(cell1, cell2):
  return d(cell_means[cell1], cell_means[cell2]) * len(cell_sets[cell1])**.5

n_scores = {}

for cellnum in cell_sets:
  for n in n_graph[cellnum]:
    n_scores[(n, cellnum)] = score(n, cellnum)

last_time = time.time()

while len(cell_sets) > N * FIRST_PASS_N_RATIO:
  if time.time() - last_time > 15:
    last_time = time.time()
    print "N-merging... (%d cells remaining)" % len(cell_sets)

  merge_from, merge_to = min(n_scores, key=lambda x: n_scores[x])

  for n in n_graph[merge_from]:
    del n_scores[(merge_from, n)]
    del n_scores[(n, merge_from)]

  merge_cells(merge_from, merge_to)

  for n in n_graph[merge_to]:
    n_scores[(n, merge_to)] = score(n, merge_to)
    n_scores[(merge_to, n)] = score(merge_to, n)

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for cellnum in cell_sets:
    cell_color = cell_means[cellnum]

    for pixel in cell_sets[cellnum]:
      frame_im.putpixel(pixel, lab2rgb(cell_color))

  frame_im.save(OUTFILE_STEM + "3.png")
  print "Saved image %s3.png" % OUTFILE_STEM

del n_graph, n_scores

print "Stage 5: N-merging complete, %d cells" % len(cell_sets)

"""
Stage 6: P merging - use k-means
"""

def form_clusters(centroids):
  clusters = defaultdict(set)

  for cellnum in cell_sets:
    # Add cell to closest centroid.
    scores = []

    for centroid in centroids:
      scores.append((d(centroid, cell_means[cellnum]), centroid))

    scores.sort()
    clusters[scores[0][1]].add(cellnum)

  return clusters

def calculate_centroid(cluster):
  L_sum = 0
  a_sum = 0
  b_sum = 0

  weighting = 0

  for cellnum in cluster:
    # Weight based on cell size
    color = cell_means[cellnum]
    cell_weight = len(cell_sets[cellnum])**.5

    L_sum += color[0]*cell_weight
    a_sum += color[1]*cell_weight
    b_sum += color[2]*cell_weight

    weighting += cell_weight

  return (L_sum/weighting, a_sum/weighting, b_sum/weighting)

def db_index(clusters):
  # Davies-Bouldin index
  scatter = {}

  for centroid, cluster in clusters.items():
    scatter_score = 0

    for cellnum in cluster:
      scatter_score += d(cell_means[cellnum], centroid) * len(cell_sets[cellnum])**.5

    scatter_score /= len(cluster)
    scatter[centroid] = scatter_score**2 # Mean squared distance

  index = 0

  for ci, cluster in clusters.items():
    dist_scores = []

    for cj in clusters:
      if ci != cj:
        dist_scores.append((scatter[ci] + scatter[cj])/d(ci, cj))

    index += max(dist_scores)

  return index

best_clusters = None
best_index = None

for i in X(K_MEANS_TRIALS):
  centroids = {cell_means[cellnum] for cellnum in random.sample(cell_sets, P)}
  converged = False

  while not converged:
    clusters = form_clusters(centroids)
    new_centroids = {calculate_centroid(cluster) for cluster in clusters.values()}

    if centroids == new_centroids:
      converged = True

    centroids = new_centroids

  index = db_index(clusters)

  if best_index is None or index < best_index:
    best_index = index
    best_clusters = clusters

del cell_means
newpix_map = {}

for centroid, cluster in best_clusters.items():
  for cellnum in cluster:
    for pixel in cell_sets[cellnum]:
      newpix_map[pixel] = centroid

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for pixel in newpix_map:
    frame_im.putpixel(pixel, lab2rgb(newpix_map[pixel]))

  frame_im.save(OUTFILE_STEM + "4.png")
  print "Saved image %s4.png" % OUTFILE_STEM

print "Stage 6: P-merging complete"

"""
Stage 7: Approximate Gaussian smoothing
     See http://blog.ivank.net/fastest-gaussian-blur.html
"""

# Hindsight tells me I should have used a class. I hate hindsight.
def vec_sum(vectors):
  assert(vectors and all(len(v) == len(vectors[0]) for v in vectors))
  return tuple(sum(x[i] for x in vectors) for i in X(len(vectors[0])))

def linear_blur(color_list):
  # Can be made faster with an accumulator
  output = []

  for i in X(len(color_list)):
    relevant_pixels = color_list[max(i-BLUR_RADIUS+1, 0):i+BLUR_RADIUS]
    pixsum = vec_sum(relevant_pixels)
    output.append(tuple(pixsum[i]/len(relevant_pixels) for i in X(3)))

  return output

def horizontal_blur():
  for row in X(height):
    colors = [blurpix_map[(i, row)] for i in X(width)]
    colors = linear_blur(colors)

    for i in X(width):
      blurpix_map[(i, row)] = colors[i]

def vertical_blur():
  for column in X(width):
    colors = [blurpix_map[(column, j)] for j in X(height)]
    colors = linear_blur(colors)

    for j in X(height):
      blurpix_map[(column, j)] = colors[j]

blurpix_map = {}

for i in X(width):
  for j in X(height):
    blurpix_map[(i, j)] = newpix_map[(i, j)]

for i in X(BLUR_RUNS):
  vertical_blur()
  horizontal_blur()

# Pixel : color of smoothed image
smoothpix_map = {}

for i in X(width):
  for j in X(height):
    pixel = (i, j)
    blur_color = blurpix_map[pixel]
    nearby_colors = {newpix_map[pixel]}

    for n in neighbours(pixel):
      nearby_colors.add(newpix_map[n])

    smoothpix_map[pixel] = min(nearby_colors, key=lambda x: d(x, blur_color))

del newpix_map, blurpix_map

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for pixel in smoothpix_map:
    frame_im.putpixel(pixel, lab2rgb(smoothpix_map[pixel]))

  frame_im.save(OUTFILE_STEM + "5.png")
  print "Saved image %s5.png" % OUTFILE_STEM

print "Stage 7: Smoothing complete"

"""
Stage 8: Flood fill pass 2
     Code copy-and-paste because I'm lazy
"""

def flood_fill(start_pixel):
  to_search = {start_pixel}
  cell = set()
  searched = set()
  start_color = smoothpix_map[start_pixel]

  while to_search:
    pixel = to_search.pop()

    if start_color == smoothpix_map[pixel]:
      cell.add(pixel)
      unplaced_pixels.remove(pixel)

      for n in neighbours(pixel):
        if n in unplaced_pixels and n not in cell and n not in searched:
          to_search.add(n)

    else:
      searched.add(pixel)

  return cell

cell_sets = {}
pixcell_map = {}
unplaced_pixels = {(i, j) for i in X(width) for j in X(height)}

while unplaced_pixels:
  start_pixel = unplaced_pixels.pop()
  unplaced_pixels.add(start_pixel)
  cell = flood_fill(start_pixel)

  cellnum = len(cell_sets)
  cell_sets[cellnum] = cell

  for pixel in cell:
    pixcell_map[pixel] = cellnum

cell_colors = {}

for cellnum in cell_sets:
  cell_colors[cellnum] = smoothpix_map[next(iter(cell_sets[cellnum]))]

print "Stage 8: Flood fill pass 2 complete, %d cells" % len(cell_sets)

"""
Stage 9: Small cell removal pass 2
"""

def score(cell1, cell2):
  return d(cell_colors[cell1], cell_colors[cell2]) * len(cell_sets[cell1])**.5

def remove_small(cell_size):
  small_cells = []

  for cellnum in cell_sets:
    if len(cell_sets[cellnum]) <= cell_size:
      small_cells.append(cellnum)

  for cellnum in small_cells:
    neighbour_cells = []

    for cell in cell_sets[cellnum]:
      for n in neighbours(cell):
        neighbour_reg = pixcell_map[n]

        if neighbour_reg != cellnum:
          neighbour_cells.append(neighbour_reg)

    closest_cell = max(neighbour_cells, key=neighbour_cells.count)

    for cell in cell_sets[cellnum]:
      pixcell_map[cell] = closest_cell

    if len(cell_sets[closest_cell]) <= cell_size:
      small_cells.remove(closest_cell)

    cell_color = cell_colors[closest_cell]

    for pixel in cell_sets[cellnum]:
      smoothpix_map[pixel] = cell_color

    cell_sets[closest_cell] |= cell_sets[cellnum]
    del cell_sets[cellnum]
    del cell_colors[cellnum]

for cell_size in X(1, SMALL_CELL_THRESHOLD):
  remove_small(cell_size)

if OUTPUT_ALL:
  frame_im = Image.new("RGB", im.size)

  for pixel in smoothpix_map:
    frame_im.putpixel(pixel, lab2rgb(smoothpix_map[pixel]))

  frame_im.save(OUTFILE_STEM + "6.png")
  print "Saved image %s6.png" % OUTFILE_STEM

print "Stage 9: Small cell removal pass 2 complete, %d cells" % len(cell_sets)

"""
Stage 10: N-merging pass 2
     Necessary as stage 7 might generate *more* cells
"""

def merge_cells(merge_from, merge_to):
  merge_from_cell = cell_sets[merge_from]

  for pixel in merge_from_cell:
    pixcell_map[pixel] = merge_to

  del cell_sets[merge_from]
  del cell_colors[merge_from]

  n_graph[merge_to] |= n_graph[merge_from]
  n_graph[merge_to].remove(merge_to)

  for n in n_graph[merge_from]:
    n_graph[n].remove(merge_from)

    if n != merge_to:
      n_graph[n].add(merge_to)

  del n_graph[merge_from]

  cell_color = cell_colors[merge_to]

  for pixel in merge_from_cell:
    smoothpix_map[pixel] = cell_color

  cell_sets[merge_to] |= merge_from_cell

n_graph = defaultdict(set)

for i in X(width):
  for j in X(height):
    pixel = (i, j)
    cell = pixcell_map[pixel]

    for n in neighbours(pixel):
      neighbour_cell = pixcell_map[n]

      if neighbour_cell != cell:
        n_graph[cell].add(neighbour_cell)
        n_graph[neighbour_cell].add(cell)

n_scores = {}

for cellnum in cell_sets:
  for n in n_graph[cellnum]:
    n_scores[(n, cellnum)] = score(n, cellnum)

last_time = time.time()

while len(cell_sets) > N:
  if time.time() - last_time > 15:
    last_time = time.time()
    print "N-merging (pass 2)... (%d cells remaining)" % len(cell_sets)

  merge_from, merge_to = min(n_scores, key=lambda x: n_scores[x])

  for n in n_graph[merge_from]:
    del n_scores[(merge_from, n)]
    del n_scores[(n, merge_from)]

  merge_cells(merge_from, merge_to)

  for n in n_graph[merge_to]:
    n_scores[(n, merge_to)] = score(n, merge_to)
    n_scores[(merge_to, n)] = score(merge_to, n)

print "Stage 10: N-merging pass 2 complete, %d cells" % len(cell_sets)

"""
Stage last: Output the image!
"""

test_im = Image.new("RGB", im.size)

for i in X(width):
  for j in X(height):
    test_im.putpixel((i, j), lab2rgb(smoothpix_map[(i, j)]))

if OUTPUT_ALL:
  test_im.save(OUTFILE_STEM + "7.png")
else:
  test_im.save(OUTFILE_STEM + ".png")

print "Done! (Time taken: {})".format(time.time() - total_time)

Update time! This update features a simple smoothing algorithm to make images look less fuzzy. If I update again I'll have to revamp a fair chunk of my code though, because it's getting messy & I hd 2 glf a fw thngs 2 mke t char lim.
I've also made k-means weight colours based on cell sizes, which loses some details for more restrictive parameters (e.g. the centre of the nebula and American Gothic's pitchfork) but makes the overall colour choice sharper and nicer. Interestingly, it loses the whole background for raytraced spheres for P = 5.
Algorithm summary:

Convert the pixels to the CIELAB colour space: CIELAB approximates human vision better than RGB. Originally I used HSL (hue, saturation, lightness) but this had two problems — hue of white/grey/black is undefined, and hue is measured in degrees which wrap around, making k-means difficult to use.
Divide the image into like-coloured cells using flood fill: Pick a pixel not in a cell and do a flood fill using a specified tolerance. To measure the distance between two colours I use the standard Euclidean norm. More complicated formulae are available on this wiki article.
Merge together small cells with their neighbours: The flood fill generates a lot of 1 or 2 pixel cells — merge cells less than a specified size with the neighbouring cell with the most adjacent pixels. This considerably reduces the number of cells, improving running time for later steps.
Merge together similarly-coloured regions: Go through the cells in order of decreasing size. If any neighbouring cell has mean colour less than a certain distance away, merge the cells. Keep going through the cells until no more can be merged.
Merge until we have less than 1.5N cells (N-merging): Merge cells together, using a scoring based on cell size and colour difference, until we have at most 1.5N cells. We allow a bit of leeway as we'll merge again later.
Merge until we have less than P colours, using k-means (P-merging): Use the k-means clustering algorithm some specified number of times to generate clusterings of cell colours, weighting based on cell size. Score each clustering based on a variation of the Davies-Bouldin index and pick the best clustering to use.
Approximate Gaussian smoothing: Use several linear blurs to approximate Gaussian blurring (details here). Then for each pixel, pick the colour out of itself and its neighbours in the pre-blurred image closest to its colour in the blurred image. This part can be optimised more time-wise if necessary as I have yet to implement the optimal algorithm.
Do another flood fill pass to work out the new regions: This is necessary as the previous step may actually generate more cells.
Do another small cell merging pass
Do another N-merging pass: This time we go down to N cells rather than 1.5N.

The processing time for each image highly depends on its size and complexity, with times ranging from 20 seconds to 7 minutes for the test images.
Because the algorithm uses randomisation (e.g. merging, k-means), you can get different results on different runs. Here's a comparison of two runs for the bear image, with N = 50 and P = 10:

Note: All images below are links. Most of these images are straight from the first run, but if I didn't like the output I allowed myself up to three attempts to be fair.
N = 50, P = 10

N = 500, P = 30

But I'm pretty lazy when it comes to paint by colours, so just for fun...
N = 20, P = 5

Additionally, it's amusing to see what happens when you try to squeeze 1 million colors into N = 500, P = 30:

Here's a step-by-step walkthrough of the algorithm for the underwater image with N = 500 and P = 30, in animated GIF form:

I've also made a gallery for the previous versions of the algorithm here. Here's some of my favourites from the last version (from when the nebula had more stars and the bear looked furrier):

If anyone's getting an Exception when the program's trying to unpack pixels, it looks like im = im.convert("RGB") is needed for some pics. I'll put that in after I restructure the code a bit.

Python 2 with PIL
Also a Python solution and probably very much a work in progress:
from PIL import Image, ImageFilter
import random

def draw(file_name, P, N, M=3):
    img = Image.open(file_name, 'r')
    pixels = img.load()
    size_x, size_y = img.size

    def dist(c1, c2):
        return (c1[0]-c2[0])**2+(c1[1]-c2[1])**2+(c1[2]-c2[2])**2

    def mean(colours):
        n = len(colours)
        r = sum(c[0] for c in colours)//n
        g = sum(c[1] for c in colours)//n
        b = sum(c[2] for c in colours)//n
        return (r,g,b)

    def colourize(colour, palette):
        return min(palette, key=lambda c: dist(c, colour))

    def cluster(colours, k, max_n=10000, max_i=10):
        colours = random.sample(colours, max_n)
        centroids = random.sample(colours, k)
        i = 0
        old_centroids = None
        while not(i>max_i or centroids==old_centroids):
            old_centroids = centroids
            i += 1
            labels = [colourize(c, centroids) for c in colours]
            centroids = [mean([c for c,l in zip(colours, labels)
                               if l is cen]) for cen in centroids]
        return centroids

    all_coords = [(x,y) for x in xrange(size_x) for y in xrange(size_y)]
    all_colours = [pixels[x,y] for x,y in all_coords]
    palette = cluster(all_colours, P)
    print 'clustered'

    for x,y in all_coords:
        pixels[x,y] = colourize(pixels[x,y], palette)
    print 'colourized'

    median_filter = ImageFilter.MedianFilter(size=M)
    img = img.filter(median_filter)
    pixels = img.load()
    for x,y in all_coords:
        pixels[x,y] = colourize(pixels[x,y], palette)
    print 'median filtered'

    def neighbours(edge, outer, colour=None):
        return set((x+a,y+b) for x,y in edge
                   for a,b in ((1,0), (-1,0), (0,1), (0,-1))
                   if (x+a,y+b) in outer
                   and (colour==None or pixels[(x+a,y+b)]==colour))

    def cell(centre, rest):
        colour = pixels[centre]
        edge = set([centre])
        region = set()
        while edge:
            region |= edge
            rest = rest-edge
            edge = set(n for n in neighbours(edge, rest, colour))
        return region, rest

    print 'start segmentation:'
    rest = set(all_coords)
    cells = []
    while rest:
        centre = random.sample(rest, 1)[0]
        region, rest = cell(centre, rest-set(centre))
        cells += [region]
        print '%d pixels remaining'%len(rest)
    cells = sorted(cells, key=len, reverse=True)
    print 'segmented (%d segments)'%len(cells)

    print 'start merging:'
    while len(cells)>N:
        small_cell = cells.pop()
        n = neighbours(small_cell, set(all_coords)-small_cell)
        for big_cell in cells:
            if big_cell & n:
                big_cell |= small_cell
                break
        print '%d segments remaining'%len(cells)
    print 'merged'

    for cell in cells:
        colour = colourize(mean([pixels[x,y] for x,y in cell]), palette)
        for x,y in cell:
            pixels[x,y] = colour
    print 'colorized again'

    img.save('P%d N%d '%(P,N)+file_name)
    print 'saved'

draw('a.png', 11, 500, 1)

The algorithm follows a different approach than SP3000's, starting with colours first:

Find a colour palette of P colours by k-means clustering and paint the image in this reduced palette.
Apply a slight median filter to get rid of some noise.
Make a list of all monochromatic cells and sort it by size.
Merge the smallest cells with their respective largest neighbour until there are only N cells left.

There is quite some room for improvement, both in terms of speed and quality of the results. Especially the cell merging step can take up to many minutes and gives far from optimal results.

P = 5, N = 45

P = 10, N = 50

P = 4, N = 250

P = 11, N = 500

I first tried to use about the same approach (tried to do it in Javascript on a canvs) but eventaully gave up because it was taking way too long, so it's really nice to see what it could have looked like, nice work!

Very nice work. I loved the bear with 20 cells.

Mathematica
At the moment, this takes the number of colors and the Gaussian radius to be used in the Gaussian filter.  The larger the radius, the greater the blurring and merging of colors.
Because it does not allow for input of the number of cells, it doesn't meet one of the basic requirements of the challenge.
Output includes the number of cells for each color and also the total number of cells.
quantImg[img_,nColours_,gaussR_]:=ColorQuantize[GaussianFilter[img,gaussR],nColours,
Dithering-> False]

colours[qImg_]:=Union[Flatten[ImageData[qImg],1]]

showColors[image_,nColors_,gaussR_]:=
   Module[{qImg,colors,ca,nCells},
   qImg=quantImg[image,nColors,gaussR];
   colors=colours[qImg];
   ca=ConstantArray[0,Reverse@ImageDimensions[image]];
   nCells[qImgg_,color_]:=
   Module[{r},
   r=ReplacePart[ca,Position[ImageData@qImg,color]/.{a_,b_}:> ({a,b}->1)];
   (*ArrayPlot[r,ColorRules->{1\[Rule]RGBColor[color],0\[Rule]White}];*)
   m=MorphologicalComponents[r];
   {RGBColor@color,Max[Union@Flatten[m,1]]}];
   s=nCells[qImg,#]&/@colors;
   Grid[{
    {Row[{s}]}, {Row[{"cells:\t\t",Tr[s[[All,2]]]}]},{Row[{"colors:\t\t",nColors}]},
    {Row[{"Gauss. Radius: ", gaussR}]}},Alignment->Left]]

Update
quantImage2 allows to specify the desired number of cells as input.  It determines the a best Gaussian Radius by looping through scenarios with greater radii until a close match is found.
quantImage2 outputs (picture, requested cells, used cells, error, gaussian Radius used).
It is, however, very slow.  To save time, you may start with an initial radius, the default value of which is 0.
gaussianRadius[img_,nCol_,nCells_,initialRadius_:0]:=
Module[{radius=initialRadius,nc=10^6,results={},r},
While[nc>nCells,(nc=numberOfCells[ape,nColors,radius]);
results=AppendTo[results,{nColors,radius,nc}];radius++];
r=results[[{-2,-1}]];
Nearest[r[[All,3]],200][[1]];
Cases[r,{_,_,Nearest[r[[All,3]],nCells][[1]]}][[1,2]]
]

quantImg2[img_,nColours_,nCells1_,initialRadius_:0]:={ColorQuantize[GaussianFilter[img,
g=gaussianRadius[img,nColours,nCells1,initialRadius]],nColours,Dithering->False],
nCells1,nn=numberOfCells[img,nColours,g],N[(nn-nCells1)/nCells1],g}

Example for which we specify the number of cells desired in the output.
Example requesting 90 cells with 25 colors. Solution returns 88 cells, 2% error.  The function chose the Gaussian radius of 55. (Lots of distortion).

Examples for which the input includes the Gaussian radius, but not the number of cells.
25 Colors, Gaussian radius of 5 pixels
nColors = 25;
gR = 5;
quantImg[balls, nColors, gR]

Three Colors, radius of 17 pixels
nColors=3;gaussianRadius=17;
showColors[wave,nColors,gaussianRadius]
quantImg[wave,nColors,gaussianRadius]

Twenty colors, radius of 17 pixels
We increased the number of colors but not the focus.  Note the increase in the number of cells.

Six Colors, radius of 4 pixels
nColors=6;gaussianRadius=4;
showColors[wave,nColors,gaussianRadius]
quantImg[wave,nColors,gaussianRadius]

nColors = 6; gaussianRadius = 17;
showColors[ape, nColors, gaussianRadius]
quantImg[ape, nColors, gaussianRadius]

nColors = 6; gaussianRadius = 3;
showColors[ape, nColors, gaussianRadius]
quantImg[ape, nColors, gaussianRadius]

Starry Night
With only 6 colors and 60 cells.
There is a color mismatch in the colors the showColors claims it uses.
(Yellow doesn't appear among the 5 colors but it is used in the drawing.)
I'll see if I can figure this out.

This is absolutely gorgeous, and works really well for restrictive parameters. Any chance of turning the number of cells into a parameter? (I suppose you could always find some estimate for the radius from the number of cells, apply that, and then merge small cells until you're below the limit.)

It is possible to make a Table of showColors, looping through a range of numbers of colors and radii and choosing the combination that comes closest to the desired number of cells.  Not sure if I have the gas to do that at the moment. Perhaps later.

Sure, let me know if you do. (I'd also love to see some more results for the other images. :))

@MartinBüttner,  I implemented a brute force method to find the best Gaussian radius for n Colors.  It's really slow, even if you suggest a starting radius value.

That's fine. Thanks for playing by the rules. ;)

I'm liking the spheres! They're nice and round

I love this answer, enough to give it a bounty, but I'm reluctant because it doesn't fit the criteria...

No problem. It was fun to make.  The first example appears to satisfy the criteria (desired number of cells is  input) but it takes way too long to run unless I start with an reasonably large Gaussian radius.  There are quicker ways, e.g. approaching through a binary tree but I'm reluctant to invest the time needed.

Martin, Thanks for the award.  I didn't know there were going to be two winners.  It was a lovely challenge.

Python 2 with PIL
This is still somewhat a work in progress. Also, the code below is a horrible mess of spaghetti, and should not be used as an inspiration. :)
from PIL import Image, ImageFilter
from math import sqrt
from copy import copy
from random import shuffle, choice, seed

IN_FILE = "input.png"
OUT_FILE = "output.png"

LOGGING = True
GRAPHICAL_LOGGING = False
LOG_FILE_PREFIX = "out"
LOG_FILE_SUFFIX = ".png"
LOG_ROUND_INTERVAL = 150
LOG_FLIP_INTERVAL = 40000

N = 500
P = 30
BLUR_RADIUS = 3
FILAMENT_ROUND_INTERVAL = 5
seed(0) # Random seed

print("Opening input file...")

image = Image.open(IN_FILE).filter(ImageFilter.GaussianBlur(BLUR_RADIUS))
pixels = {}
width, height = image.size

for i in range(width):
    for j in range(height):
        pixels[(i, j)] = image.getpixel((i, j))

def dist_rgb((a,b,c), (d,e,f)):
    return (a-d)**2 + (b-e)**2 + (c-f)**2

def nbors((x,y)):
    if 0 < x:
        if 0 < y:
            yield (x-1,y-1)
        if y < height-1:
            yield (x-1,y+1)
    if x < width - 1:
        if 0 < y:
            yield (x+1,y-1)
        if y < height-1:
            yield (x+1,y+1)

def full_circ((x,y)):
    return ((x+1,y), (x+1,y+1), (x,y+1), (x-1,y+1), (x-1,y), (x-1,y-1), (x,y-1), (x+1,y-1))

class Region:

    def __init__(self):
        self.points = set()
        self.size = 0
        self.sum = (0,0,0)

    def flip_point(self, point):
        sum_r, sum_g, sum_b = self.sum
        r, g, b = pixels[point]
        if point in self.points:
            self.sum = (sum_r - r, sum_g - g, sum_b - b)
            self.size -= 1
            self.points.remove(point)
        else:
            self.sum = (sum_r + r, sum_g + g, sum_b + b)
            self.size += 1
            self.points.add(point)

    def mean_with(self, color):
        if color is None:
            s = float(self.size)
            r, g, b = self.sum
        else:
            s = float(self.size + 1)
            r, g, b = map(lambda a,b: a+b, self.sum, color)
        return (r/s, g/s, b/s)

print("Initializing regions...")

aspect_ratio = width / float(height)
a = int(sqrt(N)*aspect_ratio)
b = int(sqrt(N)/aspect_ratio)

num_components = a*b
owners = {}
regions = [Region() for i in range(P)]
borders = set()

nodes = [(i,j) for i in range(a) for j in range(b)]
shuffle(nodes)
node_values = {(i,j):None for i in range(a) for j in range(b)}

for i in range(P):
    node_values[nodes[i]] = regions[i]

for (i,j) in nodes[P:]:
    forbiddens = set()
    for node in (i,j-1), (i,j+1), (i-1,j), (i+1,j):
        if node in node_values and node_values[node] is not None:
            forbiddens.add(node_values[node])
    node_values[(i,j)] = choice(list(set(regions) - forbiddens))

for (i,j) in nodes:
    for x in range((width*i)/a, (width*(i+1))/a):
        for y in range((height*j)/b, (height*(j+1))/b):
            owner = node_values[(i,j)]
            owner.flip_point((x,y))
            owners[(x,y)] = owner

def recalc_borders(point = None):
    global borders
    if point is None:
        borders = set()
        for i in range(width):
            for j in range(height):
                if (i,j) not in borders:
                    owner = owner_of((i,j))
                    for pt in nbors((i,j)):
                        if owner_of(pt) != owner:
                            borders.add((i,j))
                            borders.add(pt)
                            break
    else:
        for pt in nbors(point):
            owner = owner_of(pt)
            for pt2 in nbors(pt):
                if owner_of(pt2) != owner:
                    borders.add(pt)
                    break
            else:
                borders.discard(pt)

def owner_of(point):
    if 0 <= point[0] < width and 0 <= point[1] < height:
        return owners[point]
    else:
        return None

# Status codes for analysis
SINGLETON = 0
FILAMENT = 1
SWAPPABLE = 2
NOT_SWAPPABLE = 3

def analyze_nbors(point):
    owner = owner_of(point)
    circ = a,b,c,d,e,f,g,h = full_circ(point)
    oa,ob,oc,od,oe,of,og,oh = map(owner_of, circ)
    nbor_owners = set([oa,oc,oe,og])
    if owner not in nbor_owners:
        return SINGLETON, owner, nbor_owners - set([None])
    if oc != oe == owner == oa != og != oc:
        return FILAMENT, owner, set([og, oc]) - set([None])
    if oe != oc == owner == og != oa != oe:
        return FILAMENT, owner, set([oe, oa]) - set([None])
    last_owner = oa
    flips = {last_owner:0}
    for (corner, side, corner_owner, side_owner) in (b,c,ob,oc), (d,e,od,oe), (f,g,of,og), (h,a,oh,oa):
        if side_owner not in flips:
            flips[side_owner] = 0
        if side_owner != corner_owner or side_owner != last_owner:
            flips[side_owner] += 1
            flips[last_owner] += 1
        last_owner = side_owner
    candidates = set(own for own in flips if flips[own] == 2 and own is not None)
    if owner in candidates:
        return SWAPPABLE, owner, candidates - set([owner])
    return NOT_SWAPPABLE, None, None

print("Calculating borders...")

recalc_borders()

print("Deforming regions...")

def assign_colors():
    used_colors = {}
    for region in regions:
        r, g, b = region.mean_with(None)
        r, g, b = int(round(r)), int(round(g)), int(round(b))
        if (r,g,b) in used_colors:
            for color in sorted([(r2, g2, b2) for r2 in range(256) for g2 in range(256) for b2 in range(256)], key=lambda color: dist_rgb(color, (r,g,b))):
                if color not in used_colors:
                    used_colors[color] = region.points
                    break
        else:
            used_colors[(r,g,b)] = region.points
    return used_colors

def make_image(colors):
    img = Image.new("RGB", image.size)
    for color in colors:
        for point in colors[color]:
            img.putpixel(point, color)
    return img

# Round status labels
FULL_ROUND = 0
NEIGHBOR_ROUND = 1
FILAMENT_ROUND = 2

max_filament = None
next_search = set()
rounds = 0
points_flipped = 0
singletons = 0
filaments = 0
flip_milestone = 0
logs = 0

while True:
    if LOGGING and (rounds % LOG_ROUND_INTERVAL == 0 or points_flipped >= flip_milestone):
        print("Round %d of deformation:\n %d edit(s) so far, of which %d singleton removal(s) and %d filament cut(s)."%(rounds, points_flipped, singletons, filaments))
        while points_flipped >= flip_milestone: flip_milestone += LOG_FLIP_INTERVAL
        if GRAPHICAL_LOGGING:
            make_image(assign_colors()).save(LOG_FILE_PREFIX + str(logs) + LOG_FILE_SUFFIX)
            logs += 1
    if max_filament is None or (round_status == NEIGHBOR_ROUND and rounds%FILAMENT_ROUND_INTERVAL != 0):
        search_space, round_status = (next_search & borders, NEIGHBOR_ROUND) if next_search else (copy(borders), FULL_ROUND)
        next_search = set()
        max_filament = None
    else:
        round_status = FILAMENT_ROUND
        search_space = set([max_filament[0]]) & borders
    search_space = list(search_space)
    shuffle(search_space)
    for point in search_space:
        status, owner, takers = analyze_nbors(point)
        if (status == FILAMENT and num_components < N) or status in (SINGLETON, SWAPPABLE):
            color = pixels[point]
            takers_list = list(takers)
            shuffle(takers_list)
            for taker in takers_list:
                dist = dist_rgb(color, owner.mean_with(None)) - dist_rgb(color, taker.mean_with(color))
                if dist > 0:
                    if status != FILAMENT or round_status == FILAMENT_ROUND:
                        found = True
                        owner.flip_point(point)
                        taker.flip_point(point)
                        owners[point] = taker
                        recalc_borders(point)
                        next_search.add(point)
                        for nbor in full_circ(point):
                            next_search.add(nbor)
                        points_flipped += 1
                    if status == FILAMENT:
                        if round_status == FILAMENT_ROUND:
                            num_components += 1
                            filaments += 1
                        elif max_filament is None or max_filament[1] < dist:
                            max_filament = (point, dist)
                    if status == SINGLETON:
                        num_components -= 1
                        singletons += 1
                    break
    rounds += 1
    if round_status == FILAMENT_ROUND:
        max_filament = None
    if round_status == FULL_ROUND and max_filament is None and not next_search:
        break

print("Deformation completed after %d rounds:\n %d edit(s), of which %d singleton removal(s) and %d filament cut(s)."%(rounds, points_flipped, singletons, filaments))

print("Assigning colors...")

used_colors = assign_colors()

print("Producing output...")

make_image(used_colors).save(OUT_FILE)

print("Done!")

How it works
The program divides the canvas into P regions, each of which consists of some number of cells without holes. Initially, the canvas is just divided into approximate squares, which are randomly assigned to the regions. Then, these regions are "deformed" in an iterative process, where a given pixel can change its region if

the change would decrease the pixel's RGB distance from the average color of the region that contains it, and
it does not break or merge cells or introduce holes in them.

The latter condition can be enforced locally, so the process is a bit like a cellular automaton. This way, we don't have to do any pathfinding or such, which speeds the process up greatly. However, since the cells can't be broken up, some of them end up as long "filaments" that border other cells and inhibit their growth. To fix this, there is a process called "filament cut", which occasionally breaks a filament-shaped cell in two, if there are less than N cells at that time. Cells can also disappear if their size is 1, and this makes room for the filaments cuts.
The process ends when no pixel has the incentive to switch regions, and after that, each region is simply colored by its average color. Usually there will be some filaments remaining in the output, as can be seen in the examples below, especially in the nebula.
P = 30, N = 500

More pictures later.
Some interesting properties of my program are that it is probabilistic, so that the results may vary between different runs, unless you use the same pseudorandom seed of course. The randomness is not essential, though, I just wanted to avoid any accidental artifacts that may result from the particular way Python traverses a set of coordinates or something similar. The program tends to use all P colors and almost all N cells, and the cells never contain holes by design. Also, the deformation process is quite slow. The colored balls took almost 15 minutes to produce on my machine. On the upside, it you turn on the GRAPHICAL_LOGGING option, you'll get a cool series of pictures of the deformation process. I made the Mona Lisa ones into a GIF animation (shrunk 50 % to reduce the file size). If you look closely at her face and hair, you can see the filament cutting process in action.

Wow, these results look really beautiful (although not quite like it's painted by numbers, but still very nice :) ).

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Implement t9-like functionality
Your challenge today is to implement a t9-like functionality.
You will implement a function that will only have 2 parameters.
You will receive 1 phone number in a string and the content of a text file with a list of words (don't assume a specific newline style).
You can use the link https://raw.githubusercontent.com/eneko/data-repository/master/data/words.txt to test the functionality, or use /usr/share/dict/words (check A text file with a list of words [closed] for more information).
You can assume that you will always receive at least 2 numbers.
Given the number, you will read from a list of words and returns the words starting with the letters mapping to those words. This means that the input should be only numbers from 2 to 9.
You can do whatever you want if you receive invalid input.
If no match is found, you can return an empty list, null/nil or 0.
Remember that the cellphone keys are mapped to their equivalent chars:

0 and 1 are invalid
2 matches [abc]
3 matched [def]
4 matches [ghi]
5 matches [jkl]
6 matches [mno]
7 matches [pqrs]
8 matches [tuv]
and 9 matches [wxyz]

Examples:
f('52726')
//returns ["Japan","japan","Japanee","Japanese","Japanesque"...,"larbowlines"]

f('552')
//returns ["Kjeldahl","kjeldahlization","kjeldahlize"...,"Lleu","Llew"]

f('1234')
//makes demons fly out your nose or divide by 0

f('9999')
//returns ["Zyzzogeton"]

f('999999')
//returns [] or null/nil or 0

After you run your function, you can print it in any way you wish.
Rules:

Standard loopholes are INVALID
You must return something, even if it is null/nil
Javascript will return undefined if you don't return something, therefore this rule.
You cannot use or re-implement other's answers or copy my implementation.
You can assume, for Javascript, that the browser will be already opened and that the innerText/textContent of the automatic element will be passed as the 2nd parameter
For compiled languages, you cannot pass special arguments to the compiler
You can receive the file name over compiler arguments
Variables, macros, global variables, constants, non-standard classes and all the sort passing other values inside the function will be considered invalid.
In Javascript, variables without the keyword var render your code invalid
Your function will be named f
You can only and only have 2 arguments on your function
Try to keep your code under 500 seconds to run.
You don't have to worry about whitespace
You must use only ASCII printable characters.
Exceptions are languages that only use non-printable characters (APL and whitespace are 2 examples).

Scoring:

Lowest number of bytes win
Having invalid ASCII printable characters in your answer, will count as the answer being encoded in UTF-32
The exception to the encoding will make your answer be count by characters.
Only the function body counts, don't count anything else you do outside it
Bonus of -30% if you make a prediction system based on the neighborhood or most common words
Bonus of -20% in size if you only return the first 5 matches for each letter corresponding to the first number (e.g.: 245 would returns 5 words starting with 'a', 5 starting with 'b' and 5 starting with 'c').

Here is an example of an implementation, using Javascript:
function f(phone, words)
{
    var keypad=['','','abc','def','ghi','jkl','mno','pqrs','tuv','wxyz'];
    var regex='';

    for(var i=0,l=phone.length;i<l;i++)
    {
        regex+='['+keypad[phone[i]]+']';
    }

    var regexp=new RegExp('\\s('+regex+'[a-z]*)\\s','gi');

    return words.match(regexp);
}

To run it, open the list link and run, for example:
f('9999',document.getElementsByTagName('pre')[0].innerText);
//returns [" Zyzzogeton "]

This example was tested and works under Opera 12.17 64bits on Windows 7 Home Edition 64bits.

Is the second argument to the program a file name containing the words or the list of words itself ?

@MartinBüttner UTF-8 isn't unfair (it still counts ASCII chars as being 1 byte), but I changed the rule.

@Optimizer The 2nd argument is a list of the words. You can pass the filename over a compiler argument and read the file, if you want to. But the only thing that counts is the function body.

@MartinBüttner By counting as ASCII, it is being counted as bytes. You want me to say that APL code will have 1 byte having the size of 8 bits?

@MartinBüttner I think it is more fair now. The required encoding is specified (printable ASCII characters) but APL answers will be count by characters.

@MartinBüttner Now it is more fair: All languages that only use non-printable ascii chars are counted by chars. All the others are counted by bytes.

Can you explain your first bonus a bit more please ?

I'm pretty sure space, tab and newline (which are the only characters Whitespace uses) are printable. Also, correct me if I'm wrong, but by "You don't have to worry about whitespace" I assume you mean whitespace doesn't count. In that case, a Whitespace answer would instantly win with a score of 0.

@Optimizer The t9 function tries to predict the most used words. You can use this to, for example, list the most commonly used words in english. An example would be to only return ["Japan", "japan", "Japanese"] with the example number "52726".

@nyuszika7h The whitespace being returned don't count. you can return ["a"] or [" a "]. And you made a nice point about whitespace. In this case, what I am referring to is characters that show when you print. A space won't be printed, unless you put a background from a different color than the paper (some printers optimize its printing based on the paper color), but you don't have to go on this extend.

-1 for inappropriate restrictions

@AJMansfield "-1 for inappropriate restrictions" --> please, develop more about this. This type of comment is as useful as going to a concert with an air guitar.

@IsmaelMiguel Specifically, rules 2, 3, 5, 6, 7, 8, and 9 are all inappropriate. Additionally, rules 4, 12, S2, and S3 are also poor.

@AJMansfield rule 2: if it doesn't return something, it isn't a function. 3: you want to copy people's answers? 5: the goal is to create a function, not a full-blown program. 6: that is to enforce the scoring, to avoid cheaters. 7: refer to rule 6. 8: this is to avoid multiple function names. 9: check the implementations: you only need 2 arguments. 'poor rules': 4: indeed the wording is poor, but check my javascript implementation for the reason why. 12: indeed it is badly worded. s2: refer to rule 13. s3: the goal is to make a function, this is FAIR

@IsmaelMiguel Which leads to the heart of the issue. Limiting allowed programs to just functions is unreasonable, as that makes a large class of possible solutions 'invalid', such as shell scripts or command-line applications. Those rules get in the way of the fun of making super-optimized programs.

So 2 is invalid because the goal shouldn't be to make a function. 3 is invalid because it means that I can't post a super-tightly optimized version of another solution. 5 is simply wrong. 6 is already covered by standard loopholes, and your version of it is bad. Same for 7, plus the 'function' thing. 8 is overly restrictive. 9 but that may preclude certain types of optimization, such as multi-argument recursion. 4 is just an excuse to let you shorten your code. You should just use one of the regular scoring rules; S2 is worse than any of them. S3 is bad because of the function thing.

@AJMansfield 2 is perfectly fine. Look up at the tags. 3 is perfectly fine too. The goal is exactly that. If you want to answer, make your own answer. 5 is fine: no need for compiler arguments when all you have to do is a function. Without the rule 8, I would have to count the whole code and making a function wouldn't make any sense. 9: if you want recursion, use eval or lambda functions. 4 is not an excuse but a facility. If you post an answer in javascript, you can count that the page will already be opened on the file. s3 is perfect. you can make shorter functions with this.

CJam, 28 bytes
q~{el{'h-_9/-D+3/}%s1$#!},p;

Takes input in the form of "<number>" [<list of words>]
Example:
"52726" ["Japan" "japan" "Japanee" "Japanese" "Japanesque" "larbowlines" "ablution" "ablutionary" "abluvion" "ably" "abmho" "Abnaki" "abnegate"]

Output:
["Japan" "japan" "Japanee" "Japanese" "Japanesque" "larbowlines"]

Not going for any bonus for now.
Try the code online here but for actual time measurements, run it on the Java compiler
Note that CJam represents empty lists like ""
To convert the raw wordlist into CJam list use the following code with the wordlist as input:
qN/p

"You will receive 1 phone number in a string and the content of a text file with a list of words" -> can you implement, on a different block, the needed code to read the file into an usable list?

@IsmaelMiguel You mean not the part of this code, but just a helper code to convert the list to correct format ?

Exactly. Which your code isn't enough to prove it can use a list of words as the provided examples. But I upvoted anyway, I just wanted that helper code.

Can you add it to the answer? As an edit, in a different part

Exactly. That's what I'm talking about! Lets see if you can optimize it even further

@IsmaelMiguel note that the online compiler is not the official one and is at least 10 times slower than the official command line Java compiler.

Don't worry, timing isn't critical. As stated on the rules: "Try to keep your code under 500 seconds to run". Nobody wants to wait 300 years for a word list while writing an important SMS for your boss.

If you want to go for the 20% bonus, check my answer. You have an example of how to implement it.

Java: 395
This forms a regex pattern based on the letters that are allowed for each number, and then tacks on a .* at the end to match any following characters.
Here is the golfed version:
static ArrayList<String> f(String n,ArrayList<String> d){String[] k={"","","([A-Ca-c])","([D-Fd-f])","([G-Ig-i])","([J-Lj-l])","([M-Om-o])","([P-Sp-s])","([T-Vt-v])","([W-Zw-z])"};String r="";for(int i=0;i<n.length();++i)r+=k[n.charAt(i)-'0'];r += ".*";Pattern p=Pattern.compile(r);ArrayList<String> a=new ArrayList<String>();for(String w:dictionary)if(p.matcher(w).matches())a.add(w);return a;}

And here is the ungolfed version for read-ability
public static ArrayList<String> f(String phoneNumber, ArrayList<String> dictionary) {

    String[] KEY_VALUES = {"", "", "([A-Ca-c])", "([D-Fd-f])", "([G-Ig-i])",
                                            "([J-Lj-l])", "([M-Om-o])", "([P-Sp-s])",
                                            "([T-Vt-v])", "([W-Zw-z])"};

    String regex = "";
    for (int i = 0; i < phoneNumber.length(); ++i) {
        regex += KEY_VALUES[phoneNumber.charAt(i) - '0'];
    }
    regex += ".*";
    Pattern p = Pattern.compile(regex);
    ArrayList<String> answers = new ArrayList<String>();
    for (String word : dictionary) {
        if (p.matcher(word).matches()) {
            answers.add(word);
        }
    }
    return answers;
}

Your code goes against rule number 7: "Variables, macros, global variables, constants, non-standard classes and all the sort passing other values inside the function will be considered invalid." and it kinda goes against rule number 3: "You cannot use or re-implement other's answers or copy my implementation.", but on your code it is kinda debatable. And it also goes against the rule 9: "Your function will be named f".

@IsmaelMiguel Oops.  Rule 7 can be easily fixed by moving the constant inside the function.  I was just pulling it outside the function for better programming style.  Rule 9 is also an easy fix.  I confess I didn't read your answer, so I didn't intentionally try to copy it.  I can remove my answer if you think it is too close for the contest.

Your answer is alright. You have a bug on your code. On the last constant (([W-Zw-z)]) it should be ([W-Zw-z]). And on Code-golf you don't have to worry about programming styles and good practices: your code must simply do it's thing withing the required parameters. If you check my answer, you will see this line: $s=[2=>abc,def,ghi,jkl,mno,pqrs,tuv,wxyz];. This is an awful 'crime' in PHP. Basically I am forcing PHP to convert non-existing constants into strings. This is perfectly acceptable. You will also see that I'm not even setting the variable $t to an array before using it as such

@IsmaelMiguel Good catch on the regex error.  Thanks for pointing it out.  I'll try to actually golf it tomorrow; maybe find some Java examples on this site.

I'm not a java programmer, but I tell you some things. You can check http://codegolf.stackexchange.com/questions/6671/tips-for-golfing-in-java to get some tips. General tips include removing useless whitespace (newlines, spaces, tabs), one-letter-long variable names and do all that you can to reduce the code size as much as possible.

Nice!!! I don't know if it is possible, but maybe you can store (inside the function) a reference to ArrayList and Pattern? (I don't know if this is possible: I'm not a Java programmer) But really nice looking code! You really golfed the heck out of it!

@IsmaelMiguel when you say store a reference to, do you mean as a way to shorten the names?

Exactly, as long as it is inside the body. You have a few repeated data types (I think I can call them that) that you can try to shorten.

 C# .NET 4.5 235
This should work:
IEnumerable<string>F(string n,string d){IEnumerable<string>w=d.Split(null).ToList();string[]a={"","","abc","def","ghi", "jkl","mno","pqrs","tuv","wxyz"};foreach(var i in n){w=w.Where(x=>x.IndexOfAny(a[i-'0'].ToArray())>0);}return w;}

Welcome to PPCG. Your code will work, but you still need to reduce it a lot more. By removing all useless whitespace (spaces, tabs, newlines) I have managed to reduce your code to 167 bytes. This code can be reduced a lot more, I'm sure of it. I recommend reading http://codegolf.stackexchange.com/questions/173/tips-for-code-golfing-in-c to shorten even more your code. To help you a little, the word list is a string separated by newlines, and you seem to assume it is already possible to use a foreach in it. If you expect it to already be IEnumerable, include the code used outside

@IsmaelMiguel TY I will look at it. List are a IEnumerable there is no code outside what I posted.

If you look at the specification of the function, you will see the 2nd parameter is also a string. (Quoting: "You will receive 1 phone number in a string and the content of a text file with a list of words (don't assume a specific newline style).") And you have 1 useless whitespace on your a var.

I've noticed the improvements on your question, and I gave you an upvote. But you still can save a byte on your a var. But I really see noticeable improvements! Keep up the good work.

Python 2 (155 bytes)
Should also work in Python 3 with the appropriate replacements (string->bytes, b prefix on strings, etc.).
I wasn't sure if having the maketrans call outside the function is considered "fair"; if not, the function is 134 bytes with it moved inside.
EDIT: Dropped one byte from a stupid oversight.
With prepared maketrans, 67 bytes:
from string import maketrans
t=maketrans('abcdefghijklmnopqrstuvwxyz','22233344455566677778889999')

def f(n,w):
    return[x for x in w.split()if x.lower().translate(t).startswith(n)]

With maketrans in body, 134 bytes:
from string import maketrans

def f(n,w):
    return[x for x in w.split()if x.lower().translate(maketrans('abcdefghijklmnopqrstuvwxyz','22233344455566677778889999')).startswith(n)]

With import and maketrans in body, 155 bytes:
def f(n,w):
    return[x for x in w.split()if x.lower().translate(__import__('string').maketrans('abcdefghijklmnopqrstuvwxyz','22233344455566677778889999')).startswith(n)]

Test call:
print f('9999',open('words.txt','rt').read())

The maketrans is part of the function body. You should move it. I don't know if it is even possible, but you can try to directly use the import. I think I saw it somewhere... But your code is really nice!

Do you mean to move the import and call into the body? Yes, I think that can be done, too.

I was thinking about t=(from stirng import maketrans)([...]). I have no idea if it is even possible. But maybe you can use from string import as x  t=x([...]) which I'm not sure if it is possible too :/

The correct version is the last one. But the answer as-is is acceptable in my opinion. +1 for __import__('string').maketran.

OK, thanks. I've removed the invalid answers.

Those were alright, you could keep them anyway

PHP 5.4+ (171 186-20% = 148.8 bytes):
Well, this is quite a huge answer, but well.
I hope this brings more people to answer.
This function expects the raw content being read.
Here is the code:
function f($_,$a){$s=[2=>abc,def,ghi,jkl,mno,pqrs,tuv,wxyz];$a=preg_split('@\r\n|\r|\n@',$a);for($i=0;$c=$_[$i];++$i)foreach($a as$k=>$v)if(!strpos(1..$s[$c],$v[$i])||$t[$v[0]]++>4)unset($a[$k]);return$a;}

This works by verifying that the letter is in the list of allowed letters.
Example: the input 36 would make to check if 1abc has the first letter of the word and that 1def has the 2nd letter.
I append 1 so it doesn't check if the letter is in the 1st position (which would return 0 and that would evaluate to false). if(!strpos(1..$s[$c],$v[$i])) or if(!strpos($c.$s[$c],$v[$i])) will have the same effect, but the 1st confuses more and I like it.
Failing to do so, will remove the word.
With no words left, it returns an empty array.
To test this online, go to http://writecodeonline.com/php/ and create a simple variable with a word for line.
A testable example:
function f($_,$a)
{
    $s=array(2=>abc,def,ghi,jkl,mno,pqrs,tuv,wxyz);
    $a=preg_split('@\r\n|\r|\n@',$a);

    for($i=0;$c=$_[$i];++$i)
        foreach($a as$k=>$v)
            if(!strpos(1..$s[$c],$v[$i]) || $t[$v[0]]++>4)
                unset($a[$k]);
    return$a;
}

$lines=<<<WORDS
one
two
three
four
five
six
seven
eight
nine
ten
WORDS;

var_dump(f('36',$lines));

This should output:
array(1) {
    [3]=>
      string(4) "four"
}

To work on older php versions, replace $s=[2=>abc,def,ghi,jkl,mno,pqrs,tuv,wxyz]; by $s=array(2=>abc,def,ghi,jkl,mno,pqrs,tuv,wxyz);

For the 20% bonus:
To reduce the code I simply added ||$t[$v[0]]++>4, which checks how many times the first letter was used.
In php, $t doesn't need to be defined, helping to reduce a big chunk of 37.2 bytes.
To see this effect, use the following variable as the 2nd argument:
$lines=<<<WORDS
one
two
three
four
five
six
seven
eight
nine
ten
twelve
time
tutor
test
truth
WORDS;

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Coding a recursive function for highest possible input
Challenge
You are given the following function:-  which is the same as:-

with the base cases q(r, b, L) = 1 whenever r ≤ L, q(r, 0, L) = 0, if r > L and  q(r, 0, L) = 1, if r ≤ L.
Your task is to code a program that takes r as input, and outputs the value of q(r, r - L, L) for all L taking integer values from 1 to (r-1), where r is any nonnegative integer.
Example 1
Input

Enter the value of r: 2

Output

q(2,1,1) = 0.3333333333

Example 2
Input

Enter the value of r: 3

Output

q(3,2,1) = 0.1
q(3,1,2) = 0.5

Winning criterion
The code that can correctly output q(r, r-L, L) for all L taking integer values from 1 to (r-1), for the highest value of r, in less than 300 seconds. In case of a tie, the code with lesser runtime will be considered. As this is a runtime-based challenge, I shall test all submissions on my machine.

Your sum/product notations both use the letter k as their variable. I think this is technically valid, but it might be easier to read if you used two different letters.

@stokastic: Edited, thanks.

Which of the two base cases takes priority when r < L and b = 0? And what value of b will be used in the winning criterion?

@PeterTaylor: Thanks for pointing this out. The first case (r ≤ L) takes priority. Edited accordingly. Also, the highest input value of b for which the code produces a correct output in less than 300 seconds, will be used.

@PeterTaylor let me guess: are you working on a closed-form?

That doesn't really answer my second question. The current statement of the winning criterion says that it's the highest value of r with tie-breaking on the highest value of b, but since the program always takes both values it doesn't make any sense to talk about the highest value of r unless you define a function from the input value of r in the non-tiebreak case to the value of b. In other words, you need to specify a sequence of (r, b) pairs for the winning criterion to make any sense.

I'm sorry, but I don't understand what you're trying to say. Why doesn't it make sense to talk about the highest value of r? To compare two submissions, first the highest values of r that can be evaluated by both the codes are looked at, regardless of the value of b. In the case of a tie-breaker, the highest value of b that can be evaluated by both codes, for that tied value of r, is looked at. Could you provide an example in terms of two submissions to show how this current winning criterion breaks down?

One way of looking at it is that the recursion is in b as well as r, so for every r there will be a different b for which a submission can meet the time constraint. If that function from r to b has different asymptotics for two submissions, they're incomparable. Another way is that if we take b = 0, the limit on how far a submission can go is I/O-bound; if we take b = 1000, it's CPU-bound.

Oh, I think I somewhat understand the problem now. I can't think of a way to get around it though. Do you have suggestions for an alternate winning criterion?

I presume you came across the function in some context which gives you more intuition for how it should behave than I have, but the easiest fix would be to say the highest value of x such that input r=x, b=x completes within the time constraint.

Yes, I do know the context behind the function, but sadly, it doesn't give me much insight into how the function behaves, other than the fact that the number of recursions increases quite rapidly with increasing values of r and b. It's one of the reasons I started this challenge; to see if others could come up with better, possibly less recursive approaches to understand this function.

That fix sounds good enough. Let's go with it.

@JobinIdiculla One problem with using r = b is that the values tend towards 1 very quickly: http://codepad.org/RT33DjpB

@primo: Thank you, working on a fix.

@primo: Could you run your code for r = 50, b = 22 please?

@JobinIdiculla slightly more interesting values: http://codepad.org/2IWYwqSj Although I think something like b = r-L might be better: http://codepad.org/5ZbrN1ai

@primo: Hmm, I see. Could you try the code for r = 50 and b = 25 as well, please? It's just that my code isn't efficient enough, for values beyond r = 24.

@JobinIdiculla Also including r = 100; b = 50: http://codepad.org/QWb9v2uz

@primo: Thank you! I think your idea of b = r - L works much better. I think it'll be better to go with that. :)

Though this would mean the challenge statement would change (r becomes the only input). Is it alright if I change the challenge statement, given that other coders might already be working on the original statement?

@JobinIdiculla I don't think it would be a problem, because it only affects the input, rather than the function definition. It does affect the winning criterion, but they're aren't any submissions yet anyway.

@primo: I see. It's done.

q(r, r-1, 1) == 2(r!)^2/(2r)! (http://oeis.org/A001700) and r(r, 1, r-1) == (r-1)/(r+1) (trivial). Haven't checked the others yet.

In the interests of clarity, q(r,r-1,1) is the reciprocal of A001700.

Java
I've applied some algebraic transformation and dynamic programming.
public class CodeGolf42234 {
    public static void main(String[] args) {
        int r = Integer.parseInt(args[0]);
        for (int L = 1; L < r; L++) System.out.println(qDouble5(r, r-L, L));
    }

    private static double qDouble5(int _r, int _b, int L) {
        double[][] q = new double[_r+1][_b+1];
        for (int r = 0; r <= L; r++) {
            for (int b = 0; b < q[r].length; b++) q[r][b] = 1;
        }
        for (int r = L + 1; r < q.length; r++) {
            for (int b = 1 + (r - L - 1)/L; b < q[r].length; b++) {
                double sum = 0, m = 1;
                for (int k = 0; k <= L; k++) {
                    sum += m * q[r-k][b-1];
                    m = m * (r - k) / (r + b - 1 - k);
                }
                q[r][b] = sum * b / (r + b);
            }
        }
        return q[_r][_b];
    }
}

Thank you for the submission. On my system (Windows, Netbeans IDE), the highest value of r that can be evaluated using this code, in under 300 seconds (297 seconds) is 731. Good job, Peter! :)

Java
Here's a recursive solution. Just like you wanted. You can always uncomment the output and display it if you want.
import java.util.ArrayList;
import java.util.Scanner;

public class HeaderFinder
{
static ArrayList<Double> answers = new ArrayList<Double>();
public static void main(String[]args)
{
    Scanner scanner = new Scanner(System.in);
    int r = scanner.nextInt();
    for(int b=r-1; b!=0; b--)
        for(int L=r-1; L!=0; L--)
            answers.add(formulate(r,b,L));

//  int i=-1;
//      for(int b=r-1; b!=0; b--)
//          for(int L=r-1; L!=0; L--)
//              System.out.println("q("+r+", "+b+", "+L+") = "+answers.get(++i));

}

public static double formulate(double r, double b, double L)
{
    if(b==0)
        if(r>L)
            return 0.0;
        else
            return 1;
    if(r<=L)
        return 1;
    double partone = (b/(r+b));
    double result = 0.0;
    for(int k=1; k<=L; k++)
    {
        result+=summation(r, b, k) * (b/(r+b-k)) * formulate(r-k, b-1, L);
    }
    return (partone * formulate(r, b-1, L)) + result;
}

public static double summation(double r, double b, int k)
{
    double rb = r+b;
    double mul=r/rb;
    for(int i=1; i<k; i++)
    {
        mul*=(r-i)/(rb-i);
    }
    return mul;
}
}