def run(n_target):
 global save,solsave, MAX,UPPERBOUND, opt,sol_opt, count
## INITIALIZE PARAMETERS:
 MAX = n_target**3 # == "infinity" = upper bound on cost
 UPPERBOUND = int(n_target**2.5/3-n_target**2/13) ## empirical!
 count = 0 # number of calls to "store"
## INITIALIZE TABLES:
 save = {} # initialize data
 solsave = {} # remember for backtracking
 for h in range(0,n_target+2): # h = height
   save[h]={}     # table for storing optimal values
   solsave[h]={}  # table for storing optimal solutions
   for w in range(n_target+1): # w = width
     for n in range(n_target+1): # n = number of cells (redundant?)
       save[h][w,n]={}
       solsave[h][w,n]={}
 opt = {}
 sol_opt = {}

## RUN THE RECURSION
# import math
# store (h=0,w=2*int((math.sqrt(n_target))), ...
 store (h=0, w=n_target,w_prev=0, n=0, ind=(0,0,0,0,0,0),ind_prev=0,cost=0)
 for h in range(0,n_target+1):
   for w in range(n_target,0,-1):
     for n in range(n_target+1):
       process(h,w,n,n_target)
   save[h]={} # save storage: discard previous entries
   print "h =",h,"...", # output while user is waiting
   print ("+"+str(count-count_prev)+" =" if h else ""), count, "steps"
   count_prev = count

## PRINT THE RESULTS
 for n in range(1,n_target+1):
   if opt[n]>UPPERBOUND:
     error("UPPERBOUND was set too small.",UPPERBOUND)
     abort
   else:
     print "n="+str(n)+': cost='+str(opt[n])+',',
     printsolutions(n)

def printsolutions(n):
  sols = []
  for (h, w, ind) in sol_opt[n]:
     retrieve_all_solutions([w],h-1,w,ind,sols)
  nsols = 0
  for w in sols: # count solutions
    if len(w)<w[-1] or (len(w)==w[-1] and w >= dualpartition(w)):
       # symmetric solutions are discarded
       nsols +=1
  print nsols, "solution"+("s." if nsols!=1 else ".")
  for w in sols:
    if len(w)<w[-1] or (len(w)==w[-1] and w <= dualpartition(w)):
      print ' ',str(n)+' = ' + '+'.join(map(str,w))
      displaysolution(w,'    ')

def retrieve_all_solutions(partial,h,w0,ind0,sols):
  if h==1:
     sols.append(partial)
  else:
     (s_left_down, s_left_up, s_right_down, s_right_up,
          n_left, n_right) = ind0
     n = n_left+n_right + h*w0
     for (w,ind) in solsave[h][w0, n][ind0]:
        retrieve_all_solutions(partial+[w],h-1,w,ind,sols)

def displaysolution(w, prefix='', nix = '.', foll= 'X'):
  for i in range(0,len(w),2)+range(len(w)/2*2-1,0,-2):
    le = w[i]
    print prefix,
    print nix*(w[-1]/2-le/2)+ foll*le +nix*((w[-1]+1)/2-(le+1)/2)

def dualpartition(w):
  "spap horizontal and vertical direction"
  wd = [0]*w[-1]
  for le in w:
    for i in range(w[-1]/2-le/2, w[-1]/2-le/2+le):   wd[i] += 1
  wd.sort()
  return wd
  
def process(h,w0,n,n_target):
  "process one index triple h,w0,n of the table"
  global UPPERBOUND
  for (ind_prev, cost) in save[h][w0,n].items():
    (s_left_up, s_left_down, s_right_up, s_right_down,
          n_left, n_right) = ind_prev
    s_left_up  += n_left  # add 1 (vertical) to all up-distances
    s_right_up += n_right

    w = w0
    while w>=0:
        n_total = n_left+n_right + (h+1)*w # h == old value of h
        if n_total <= n_target:
           new_cost = cost + ( (s_left_up + s_right_up) * w +
                               (n_left + n_right) * w*(w-1)/2 +
                               (w+1)*w*(w-1)/6 * (2*h+1) +
                                w*w * h*(h+1)/2 )
           if new_cost <= UPPERBOUND:
	     # empirical cut-off value for pruning unnecessary branches
               ## potential additional speed-up by taking into account
               ## s_left, n_left,h,w etc.
             store (h+1, w, w0, n_total,
              (s_left_down, s_left_up, s_right_down, s_right_up,
                                                n_left, n_right), 
              ind_prev, new_cost) # exchange up and down when storing

        if (w%2)==1: # remove element from the left of the top row
            n_left += h
            s_left_up   += n_left + h*(h+1)/2
            s_left_down += n_left + h*(h-1)/2
        else: # remove one element from the right of the top row
            n_right += h
            s_right_up   += n_right + h*(h+1)/2
            s_right_down += n_right + h*(h-1)/2
        w -= 1

def store(h, w, w_prev, n, ind, ind_prev, cost):
  global save,solsave, MAX ,opt,sol_opt, count
  count +=1
  if w==0: # width is reduced to 0: solution is complete.
    if cost < opt.get(n, MAX):
        opt[n] = cost
        sol_opt[n] = [(h, w_prev, ind_prev)]
    elif cost == opt[n]: # store alternative optimal solutions
        sol_opt[n].append((h, w_prev, ind_prev))
  # otherwise store it in the array if it is better than previous:
  elif save[h][w, n].get(ind, MAX) > cost:
    save   [h][w, n][ind] = cost
    solsave[h][w, n][ind] = [(w_prev,ind_prev)]
  elif save[h][w, n].get(ind, MAX) == cost:
    solsave[h][w, n][ind].append((w_prev,ind_prev))

run(80)
