DouKaku?

import copy
import time
def mahobunkatsu(N):
    m = set()
    for i in xrange(1,N**2+1):
        m.add(i)
    ave = N*(N**2+1) / 2
    
    def createPair(restNumber,numberSum, pair, count, pairList):
        for i in restNumber:
            if count and (count+1)%N == 0:
                if numberSum + i == ave:
                    if count+1 == N**2:
                        pair.append(i)
                        pairList.append(copy.copy(pair))
                        pair.pop()
                        continue
                else:
                    continue
            if numberSum + i > ave:
                continue
            restNumber.remove(i)
            pair.append(i)
            createPair(restNumber, (numberSum + i) % ave, pair, count + 1, pairList)
            restNumber.add(i)
            pair.pop()
        return

    pairList = []
    createPair(m, 0, [], 0, pairList)

    def listToSortedPairs(pairList):
        temp = []
        for l in pairList:
            temp.append([])
            for i in range(N):
                temp[-1].append(l[i*N:(i+1)*N])
                temp[-1][-1].sort()
                temp[-1].sort()
        return temp

    sortedList = listToSortedPairs(pairList)

    uniqueList = set()
    for x in sortedList:
        uniqueList.add(str(x))

    print uniqueList
    return uniqueList


t = time.time()
for x in range(2,6):
    n = len(mahobunkatsu(x))
    print "Size=",x,"Mahozin_num=",n,"time=",time.time()-t
    t = time.time()

f[n_] := f[Range[n^2], n, Total@Range[n^2]/n]

f[in_, n_, s_] :=
  If[Length@in == n, 1,
    Total[f[Complement[Rest@in, #], n, s] & /@ (
    Select[Subsets[Rest@in, {n - 1}], First@in + Total@# == s &])]]

module NS = Set.Make(
struct
   type t = int
   let compare = Pervasives.compare
end)

module SS = Set.Make(
struct
   type t = NS.t
   let compare = NS.compare
end)

module LS = Set.Make(
struct
   type t = SS.t
   let compare = SS.compare
end)

let make_set_with_sum n =
   let rec loop set sum = function
      | 0  -> (set, sum)
      | n' -> loop (NS.add n' set) (sum + n') (n' - 1)
   in
   loop NS.empty 0 n

let ns_map f ns = NS.fold (fun s acc -> NS.add (f s) acc) ns NS.empty
let ss_map f ss = SS.fold (fun s acc -> SS.add (f s) acc) ss SS.empty
let ls_map f ls = LS.fold (fun s acc -> LS.add (f s) acc) ls LS.empty

let rec take_subsets set num limit =
   if NS.is_empty set || NS.cardinal set < num || NS.choose set > limit
   then SS.empty
   else begin
      match num with
      | 1 ->
           if NS.mem limit set
           then SS.singleton (NS.singleton limit)
           else SS.empty
      | n when n > 1 ->
           let result = ref SS.empty in
           NS.iter begin fun i ->
              let set' = NS.remove i set in
              let num' = num - 1 in
              let limit' = limit - i in
              begin match take_subsets set' num' limit' with
              | s when SS.is_empty s -> ()
              | ss -> result := SS.union !result (ss_map (NS.add i) ss)
              end
           end set;
           !result
      | _ -> invalid_arg "num is required positive number."
   end

let make_sets n =
   if n < 2 then invalid_arg "required greater than or equal to 2.";
   let full, max = make_set_with_sum (n * n) in
   let heads, _  = make_set_with_sum n in
   let diff  = NS.diff full heads in
   let limit = max / n in
   let subsets =
      NS.fold begin fun i acc ->
         LS.add
            (ss_map (NS.add i) (take_subsets diff (n - 1) (limit - i)))
            acc
      end heads LS.empty
   in
   let inters =
      let hd = LS.choose subsets in
      let tl = LS.remove hd subsets in
      LS.fold begin fun ss acc ->
         let result = ref LS.empty in
         LS.iter begin fun ss' ->
            SS.iter begin fun ns ->
               if SS.for_all (fun s -> NS.is_empty (NS.inter s ns)) ss'
               then result := LS.add (SS.add ns ss') !result
            end ss
         end acc;
         !result
      end tl (SS.fold (fun e acc -> LS.add (SS.singleton e) acc) hd LS.empty)
   in
   inters

let ns_print ns =
   print_string "{ ";
   NS.iter (Printf.printf "%d, ") ns;
   print_string "}"
let ss_print ss =
   SS.iter begin fun s ->
      ns_print s;
      print_string ", "
   end ss;
   print_newline ()
let ls_print ls = LS.iter ss_print ls

let exam n = LS.cardinal (make_sets n)

let main () =
   let num =
      match Sys.argv with
      | [|_; n |] -> int_of_string n
      | _ -> 3
   in
   Printf.printf "n = %d => %d patterns\n" num (exam num)
let () = if not !Sys.interactive then main ()

def gen_comb(list,n)
  if n==0
    yield([])
  else
    (list.size-n+1).times do |i|
      gen_comb(list[i+1..-1],n-1) do|ls|
        yield([list[i]]+ls)
      end
    end
  end
end

def comb_num(list,n,sum)
  ret=[]
  gen_comb(list,n) do |a|
    ret<<=a if (a.inject(0){|r,e| r+=e})==sum
  end
  ret
end

def comb_array(list,n)
  ret=[]
  gen_comb(list,n) do |a|
    ret<<=a if a.inject{|r,e| break unless (r&e).empty?;r+=e}
  end
  ret
end

def maho(n)
  m=n**2
  sum=m*(m+1)/2/n
  comb_array(comb_num((1..m).to_a,n,sum),n)
end

start_time = Time.now
puts maho(4).size
puts Time.now-start_time

(* list solutions of
     x1+...+xn=m
     1<=x1<...<xn<=l *)
let list_solutions m n k =
  let rec loop m n lbd tail =
    if m < n || m < lbd then []
    else if n = 1 then
      if lbd <= m && m <= k then [m::tail] else []
    else
      let acc = ref [] in
        for i = lbd to (min (m-n+1) k) do
          acc :=
            List.map (fun sols -> i::sols)
              (loop (m-i) (n-1) (i+1) tail) :: !acc
        done;
        List.concat !acc
  in loop m n 1

let rec disjoint xs ys =
  match xs, ys with
    | [], _ | _, [] -> true
    | x::xs', y::ys' ->
        if x > y then disjoint xs ys'
        else if x < y then disjoint xs' ys
        else false

let (@<) = List.merge compare

let rec count_choices n lists ex =
  if n = 1 then
    List.fold_left (fun i x -> if disjoint x ex then i+1 else i) 0 lists
  else
    match lists with
      | [] -> 0
      | x::rest ->
          let b = count_choices n rest ex in
            if disjoint x ex then
              b + count_choices (n-1) rest (x@<ex)
            else
              b

let count_partitions n =
  count_choices n (list_solutions (n * (n*n + 1) / 2) n (n*n) []) []

let _ = print_int (count_partitions 5); exit 0

(use srfi-1)
(use util.combinations)
(define (maho n)
  (define (maho-in m l)
    (if (= m 0) 1
      (letrec
       ((cn (/ (* n (+ (* n n) 1)) 2))
        (mylist
         (map (lambda (li) (cons (car l) li))
              (filter
               (lambda (li) (= (apply + li) (- cn (car l))))
               (combinations (cdr l) (- n 1))))))
       (apply + (map (lambda (a) (maho-in (- m n) (lset-difference equal? l a))) mylist)))))
  (maho-in (* n n) (iota (* n n) 1)))

(define (maho-by-enm n)
  (define (center n) (/ (* n (+ (* n n) 1)) 2))
  (define (flatten2 l c)
    (define (flatten1 l c)
      (if (null? l) c
        (cons (car l) (flatten1 (cdr l) c))))
    (if (null? l) c
      (flatten2 (cdr l) (flatten1 (car l) c))))
  (define (my-equal? l1 l2)
    (null? (lset-xor eq? (flatten2 l1 '()) l2)))
  (define (enm-n n)
    (filter (lambda (l) (= (apply + l) (center n)))
            (combinations (iota (* n n) 1) n)))
  (filter (lambda (l) (my-equal? l (iota (* n n) 1)))
          (combinations (enm-n n) n)))

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class Answer108 {
    private final int size_;
    private final int maxNumber_;
    private final int average_;

    private int count_ = 0;

    public Answer108(int size) {
        size_ = size;
        maxNumber_ = size * size;
        average_ = size * (maxNumber_ + 1) / 2;

        countPair();
    }

    private void countPair() {
        if (size_ <= 1) return;
        List<List<Integer>> array = new ArrayList<List<Integer>>();
        for (int index = 0; index < size_; index++) {
            array.add(new ArrayList<Integer>());
        }
        array.get(0).add(1);
        countPair(array, 2);
    }
    private void countPair(List<List<Integer>> array, int nextNumber) {
        if (nextNumber <= maxNumber_) {
            for (int index = 0; index < size_; index++) {
                List<Integer> list = array.get(index);
                int size = list.size();
                if (size == size_) continue;
                if (size == size_ - 1) {
                    if (sum(list) + nextNumber != average_) continue;
                } else {
                    int rest = 0;
                    for (int lastIndex = 0; lastIndex < size_ - size - 1; lastIndex++) {
                        rest += maxNumber_ - lastIndex;
                    }
                    if (sum(list) + nextNumber + rest < average_) continue;
                }

                if (index >= nextNumber) continue;
                list.add(nextNumber);
                countPair(array, nextNumber + 1);
                list.remove(Integer.valueOf(nextNumber));
                if (list.size() == 0) break;
            }
        } else {
            //System.out.println(toString(array));
            count_++;
        }
    }
    private int sum(List<Integer> array) {
        int sum = 0;
        for (int num: array) {
            sum += num;
        }
        return sum;
    }

    public int getCount() {
        return count_;
    }


    public static String toString(List<List<Integer>> array) {
        String[] strs = new String[array.size()];
        for (int index = 0; index < strs.length; index++) {
            strs[index] = array.get(index).toString();
        }
        return Arrays.toString(strs);
    }

    public static void main(String[] args) {
        long start = System.currentTimeMillis();
        Answer108 ans = new Answer108(5);
        System.out.println(ans.getCount());
        long end = System.currentTimeMillis();
        System.out.println("elapse: " + (end - start) + "(ms)");
    }
}

#include <iostream>
#include <vector>
typedef unsigned int bit_t;

std::vector<bit_t> bits;
int n, n2, cnt = 0;

void comb(int a, int k, int s, bit_t b) {
  if (k < n-1)
    for (int i = a; i < n2; ++i)
      comb(i+1, k+1, s-i, b | 1 << (i-1));
  else
    if (a <= s && s <= n2)
      bits.push_back(b | 1 << (s-1));
}

void calc(int s, int k, bit_t b) {
  if (k == n) { ++cnt; return; }

  for (int i = s; i < (int)bits.size(); ++i)
    if (!(b & bits[i]))
      calc(i+1, k+1, b | bits[i]);
}

int main() {
  n = 5;
  n2 = n * n;
  int m = n * (n2+1) / 2;

  comb(1, 0, m, 0);
  calc(0, 0, 0);

  std::cout << cnt << std::endl;
}