Problem:
3D Printing We are designing an installation art piece consisting of 
a number of cubes with 3D printing technology for submitting one to
Installation art Contest with Printed Cubes (ICPC).
At this time, we are trying to model a piece consisting of exactly k cubes 
of the same size facing the same direction. First, using a CAD system, 
we prepare n ( n ≥ k ) positions as candidates
in the 3D space where cubes can be placed.

When cubes would be placed at all the candidate positions, 
the following three conditions are satisfied.  Each cube may overlap zero, one or two other cubes, but not three or more. When a cube overlap two other cubes, those two cubes do not overlap. Two non-overlapping cubes do not touch at their surfaces, edges or corners.  Second, choosing appropriate k different positions from n candidates
and placing cubes there, 
we obtain a connected polyhedron as a union of the k cubes. 

When we use a 3D printer,
we usually print only the thin surface of a 3D object.

In order to save the amount of filament material for the 3D printer, 
we want to find the polyhedron with the minimal surface area. Your job is to find the polyhedron with the minimal surface area
consisting of k connected
cubes placed at k selected positions among n given ones.   Figure E1. A polyhedron formed with connected identical cubes.  Input The input consists of multiple datasets. 
The number of datasets is at most 100.
Each dataset is in the following format.  n  k  s  x 1  y 1  z 1 ...  x n  y n  z n  In the first line of a dataset, n is the number of the candidate positions, k is the number of the cubes to form the connected polyhedron,
and s is the edge length of cubes. n, k and s are integers separated by a space. 
The following n lines specify the n candidate positions.
In the i -th line, there are three integers x i , y i and z i that specify the coordinates of a position,  
where the corner of the cube with the smallest coordinate values
may be placed.
Edges of the cubes are to be aligned with either of three axes.
All the values of coordinates are integers separated by a space.

The three conditions on the candidate positions mentioned above are satisfied. The parameters satisfy the following conditions:
1 ≤ k ≤ n ≤ 2000, 3 ≤ s ≤ 100,  and
-4×10 7 ≤ x i , y i , z i ≤ 4×10 7 . The end of the input is indicated by a line containing three zeros
separated by a space. Output For each dataset, output a single line containing one integer
indicating the surface area of the connected polyhedron with 
the minimal surface area.

When no k cubes form a connected polyhedron, output -1. Sample Input 1 1 100
100 100 100
6 4 10
100 100 100
106 102 102
112 110 104
104 116 102
100 114 104
92 107 100
10 4 10
-100 101 100
-108 102 120
-116 103 100
-124 100 100
-132 99 100
-92 98 100
-84 100 140
-76 103 100
-68 102 100
-60 101 100
10 4 10
100 100 100
108 101 100
116 102 100
124 100 100
132 102 100
200 100 103
192 100 102
184 100 101
176 100 100
168 100 103
4 4 10
100 100 100
108 94 100
116 100 100
108 106 100
23 6 10
100 100 100
96 109 100
100 118 100
109 126 100
118 126 100
127 118 98
127 109 104
127 100 97
118 91 102
109 91 100
111 102 100
111 102 109
111 102 118
111 102 91
111 102 82
111 114 96
111 114 105
102 114 114
93 114 114
84 114 105
84 114 96
93 114 87
102 114 87
10 3 10
100 100 100
116 116 102
132 132 104
148 148 106
164 164 108
108 108 108
124 124 106
140 140 104
156 156 102
172 172 100
0 0 0 Output for the Sample Input 60000
1856
-1
1632
1856
2796
1640

Solution (Python):
class Cube:
    def __init__(self, x, y, z, s):
        self.x, self.y, self.z = x, y, z
        self.s = s

    def is_in_cube(self, x, y, z):
        return self.x <= x <= self.x + self.s and self.y <= y <= self.y + self.s and self.z <= z <= self.z + self.s

    def intersect(self, C):
        dxyz = [(0, 0, 0),
                (C.s, 0, 0), (0, C.s, 0), (0, 0, C.s),
                (C.s, C.s, 0), (C.s, 0, C.s), (0, C.s, C.s),
                (C.s, C.s, C.s)]
        for dx1, dy1, dz1 in dxyz:
            nx1, ny1, nz1 = C.x + dx1, C.y + dy1, C.z + dz1
            if self.is_in_cube(nx1, ny1, nz1):
                for dx2, dy2, dz2 in dxyz:
                    nx2, ny2, nz2 = self.x + dx2, self.y + dy2, self.z + dz2
                    if C.is_in_cube(nx2, ny2, nz2):
                        a, b, c = abs(nx1 - nx2), abs(ny1 - ny2), abs(nz1 - nz2)
                        if a * b * c == 0:
                            continue
                        
                        return 2 * (a * b + b * c + c * a)
        return 0


edges = list()
inters = dict()


def calc_overlap(vs):
    ret = sum(inters.get((vs[i], vs[i + 1]), 0) for i in range(len(vs) - 1))
    if len(vs) > 2:
        ret += inters.get((vs[-1], vs[0]), 0)
    return ret


def dfs(v, par, vs, res):
    if res == 0:
        return calc_overlap(vs)

    ret = -1

    for e in edges[v]:
        if e != par:
            vs.append(e)
            ret = max(ret, dfs(e, v, vs, res - 1))
            vs.pop()
    return ret


INF = 10 ** 9
while True:
    N, K, S = map(int, input().split())
    
    if not (N | K | S):
        break
    cubes = []
    for _ in range(N):
        x, y, z = map(int, input().split())
        cubes.append(Cube(x, y, z, S))
    
    edges = [[] for _ in range(N)]
    inters = dict()
    for i in range(N):
        for j in range(i + 1, N):
            sur = cubes[i].intersect(cubes[j])
            if sur > 0:
                
                inters[i, j] = inters[j, i] = sur
                edges[i].append(j)
                edges[j].append(i)
    

    ans = -1
    for i in range(N):
        ans = max(ans, dfs(i, -1, [i], K - 1))

    print(-1 if ans == -1 else S * S * 6 * K - ans)

Solution (Python):
class Cube:
    def __init__(self, x, y, z, s):
        self.x, self.y, self.z = x, y, z
        self.s = s

    def is_in_cube(self, x, y, z):
        return self.x <= x <= self.x + self.s and self.y <= y <= self.y + self.s and self.z <= z <= self.z + self.s

    def intersect(self, C):
        dxyz = [(0, 0, 0),
                (C.s, 0, 0), (0, C.s, 0), (0, 0, C.s),
                (C.s, C.s, 0), (C.s, 0, C.s), (0, C.s, C.s),
                (C.s, C.s, C.s)]
        for dx1, dy1, dz1 in dxyz:
            nx1, ny1, nz1 = C.x + dx1, C.y + dy1, C.z + dz1
            if self.is_in_cube(nx1, ny1, nz1):
                for dx2, dy2, dz2 in dxyz:
                    nx2, ny2, nz2 = self.x + dx2, self.y + dy2, self.z + dz2
                    if C.is_in_cube(nx2, ny2, nz2):
                        a, b, c = abs(nx1 - nx2), abs(ny1 - ny2), abs(nz1 - nz2)
                        if a * b * c == 0:
                            continue
                        
                        return 2 * (a * b + b * c + c * a)
        return 0


edges = list()
inters = dict()


def calc_overlap(vs):
    ret = sum(inters.get((vs[i], vs[i + 1]), 0) for i in range(len(vs) - 1))
    if len(vs) > 2:
        ret += inters.get((vs[-1], vs[0]), 0)
    return ret


def dfs(v, par, vs, res):
    if res == 0:
        return calc_overlap(vs)

    ret = -1

    for e in edges[v]:
        if e != par:
            vs.append(e)
            ret = max(ret, dfs(e, v, vs, res - 1))
            vs.pop()
    return ret


while True:
    N, K, S = map(int, input().split())
    
    if not (N | K | S):
        break
    cubes = []
    for _ in range(N):
        x, y, z = map(int, input().split())
        cubes.append(Cube(x, y, z, S))
    
    edges = [[] for _ in range(N)]
    inters = dict()
    for i in range(N):
        for j in range(i + 1, N):
            sur = cubes[i].intersect(cubes[j])
            if sur > 0:
                
                inters[i, j] = inters[j, i] = sur
                edges[i].append(j)
                edges[j].append(i)
    

    ans = -1
    for i in range(N):
        ans = max(ans, dfs(i, -1, [i], K - 1))

    print(-1 if ans == -1 else S * S * 6 * K - ans)

Solution (Python):
from collections import deque
while 1:
    n, k, s = map(int, raw_input().split())
    if n==0:break
    ps = [map(int, raw_input().split()) for i in xrange(n)]

    G = [[] for i in xrange(n)]
    for i in xrange(n):
        xi, yi, zi = ps[i]
        for j in xrange(i+1, n):
            if i==j: continue
            xj, yj, zj = ps[j]

            dx = abs(xi-xj)
            dy = abs(yi-yj)
            dz = abs(zi-zj)

            if dx < s and dy < s and dz < s:
                cost = 2*((s-dx)*(s-dy)+(s-dy)*(s-dz)+(s-dz)*(s-dx))
                G[i].append((j, cost))
                G[j].append((i, cost))
    if k == 1:
        print 6*s*s
        continue

    ans = -1
    leaf = set()
    used = [0]*n
    for i in xrange(n):
        if len(G[i]) == 0:
            used[i] = 1
        elif len(G[i]) == 1:
            leaf.add(i)
    for v in leaf:
        if used[v]:
            continue
        used[v] = 1
        prev = t = None
        deq = deque()
        su = 0
        while 1:
            if prev is not None and len(G[v]) == 1:
                break
            for t, cost in G[v]:
                if t == prev:
                    continue
                used[t] = 1
                if len(deq) < k-1:
                    deq.append(cost)
                    su += cost
                else:
                    su -= deq.popleft()
                    deq.append(cost)
                    su += cost
                break
            if len(deq) == k-1:
                ans = max(ans, su)
            v, prev = t, v
    for v in xrange(n):
        if used[v]:
            continue
        prev = t = None
        used[v] = 1
        u = set([v])
        while v is not None and used[v] < 3:
            for t, cost in G[v]:
                if t == prev:
                    continue
                used[t] += 1
                u.add(t)
                break
            v, prev = t, v
        cont = k if len(u) == k else k-1
        prev = t = None
        deq = deque(); su = 0
        for i in xrange(2*len(u)):
            for t, cost in G[v]:
                if t == prev:
                    continue
                if len(deq) < cont:
                    deq.append(cost)
                    su += cost
                else:
                    su -= deq.popleft()
                    su += cost
                    deq.append(cost)
                break
            if len(deq) == cont:
                ans = max(ans, su)
            v, prev = t, v

    if ans == -1:
        print -1
    else:
        print 6*k*s*s - ans

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools
input = sys.stdin.readline
INF = 10**20
EPS = 1.0 / 10**10
MOD = 10**9 + 7
def LI(): return [int(x) for x in input().split()]
def LF(): return [float(x) for x in input().split()]
def LS(): return input().split()
def I(): return int(input())
def F(): return float(input())
def S(): return input()

class Graph:
  def __init__(self, _n):
    self.n = _n
    self.m = 0
    self.g = [[] for i in range(_n)]
  def add_edge(self, s, t, c = 1):
    self.g[s].append((t,c))
    self.g[t].append((s,c))
  def __getitem__(self, v):
    return self.g[v]

def edge_cost(i,j):
  L = [s - abs(ps[i][k]-ps[j][k]) for k in range(3)]
  if len([x for x in L if x <= 0]) > 0:
    return -1
  return 2*(L[0]*L[1]+L[1]*L[2]+L[2]*L[0])

def dfs(v, pv, k):
  if k == 0:
    for e in G[v]:
      if e[0] != pv and used[e[0]]:
        return e[1]
    return 0
  if used[v]:
    return INF
  used[v] = True
  res = INF
  for e in G[v]:
    if e[0] == pv:
      continue
    dd = dfs(e[0],v,k-1)
    if dd < INF:
      res = min(res,dd+e[1])
  used[v] = False
  return res


if __name__ == '__main__':
  while True:
    n,k,s = LI()
    if n == 0:
      break
    ps = [LI() for i in range(n)]

    G = Graph(n)
    for i in range(n-1):
      for j in range(i+1,n):
        c = edge_cost(i,j)
        if c > 0:
          G.add_edge(i,j,-c)

    used = [False] * n
    sub = INF
    for i in range(n):
      sub = min(sub,dfs(i,-1,k-1))
    if sub == INF:
      print(-1)
    else:
      print(6*k*s*s+sub)

Solution (Python):
import math,string,itertools,fractions,heapq,collections,re,array,bisect,sys,random,time,copy,functools
input = sys.stdin.readline
INF = 10**20
EPS = 1.0 / 10**10
MOD = 10**9 + 7
def LI(): return [int(x) for x in input().split()]
def LF(): return [float(x) for x in input().split()]
def LS(): return input().split()
def I(): return int(input())
def F(): return float(input())
def S(): return input()

class Graph:
  def __init__(self, _n):
    self.n = _n
    self.m = 0
    self.g = [[] for i in range(_n)]
  def add_edge(self, s, t, c = 1):
    self.g[s].append((t,c))
    self.g[t].append((s,c))
  def __getitem__(self, v):
    return self.g[v]

def edge_cost(i,j):
  L = [s - abs(ps[i][k]-ps[j][k]) for k in range(3)]
  if len([x for x in L if x <= 0]) > 0:
    return -1
  return 2*(L[0]*L[1]+L[1]*L[2]+L[2]*L[0])

def dfs(v, pv, k):
  if k == 0:
    for e in G[v]:
      if e[0] != pv and used[e[0]]:
        return e[1]
    return 0
  if used[v]:
    return INF
  used[v] = True
  res = INF
  for e in G[v]:
    if e[0] == pv:
      continue
    dd = dfs(e[0],v,k-1)
    if dd < INF:
      res = min(res,dd+e[1])
  used[v] = False
  return res


if __name__ == '__main__':
  while True:
    n,k,s = LI()
    if n == 0:
      break
    ps = [LI() for i in range(n)]

    G = Graph(n)
    for i in range(n-1):
      for j in range(i+1,n):
        c = edge_cost(i,j)
        if c > 0:
          G.add_edge(i,j,-c)

    used = [False] * n
    sub = INF
    for i in range(n):
      sub = min(sub,dfs(i,-1,k-1))
    if sub == INF:
      print(-1)
    else:
      print(6*k*s*s+sub)