Вычисляем «магические квадраты» с помощью GPU

#include <set>
#include <stdio.h>
#include <ctime>

typedef std::set<int> Set;
typedef Set::iterator SetIterator;

#define N 4
#define MAX (N*N)
#define S 34

int main(int argc, char *argv[])
{
  // x11 x12 x13 x14
  // x21 x22 x23 x24
  // x31 x32 x33 x34
  // x41 x42 x43 x44

  const clock_t begin_time = clock();

  int squares = 0;
  int digits[] = { 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 };
  Set mset(digits, digits + N*N);
  for (int x11 = 1; x11 <= MAX; x11++) {
    Set set14(mset); set14.erase(x11);
    for (SetIterator it14 = set14.begin(); it14 != set14.end(); it14++) {
      int x14 = *it14;
      Set set41(set14); set41.erase(x14);
      for (SetIterator it41 = set41.begin(); it41 != set41.end(); it41++) {
        int x41 = *it41;
        int x44 = S - x11 - x14 - x41;
        if (x44 < 1 || x44 > MAX) continue;
        if (x44 == x11 || x44 == x14 || x44 == x41) continue;
        Set set22(set41); set22.erase(x41); set22.erase(x44);
        for (SetIterator it22 = set22.begin(); it22 != set22.end(); it22++){
          int x22 = *it22;
          int x33 = S - x11 - x22 - x44;
          if (x33 < 1 || x33 > MAX) continue;
          if (x33 == x11 || x33 == x14 || x33 == x41 || x33 == x44 || x33 == x22) continue;
          Set set23(set22); set23.erase(x22); set23.erase(x33);
          for (SetIterator it23 = set23.begin(); it23 != set23.end(); it23++){
            int x23 = *it23;
            int x32 = S - x14 - x23 - x41;
            if (x32 < 1 || x32 > MAX) continue;
            if (x32 == x11 || x32 == x14 || x32 == x41 || x32 == x44 || x32 == x22 || x32 == x33 || x32 == x23) continue;
            Set set12(set23); set12.erase(x23); set12.erase(x32);
            for (SetIterator it12 = set12.begin(); it12 != set12.end(); it12++){
              int x12 = *it12;
              int x13 = S - x11 - x12 - x14;
              if (x13 < 1 || x13 > MAX) continue;
              if (x13 == x11 || x13 == x14 || x13 == x41 || x13 == x44 || x13 == x22 || x13 == x33 || x13 == x23 || x13 == x32 || x13 == x12) continue;
              int x42 = S - x12 - x22 - x32;
              if (x42 < 1 || x42 > MAX) continue;
              if (x42 == x11 || x42 == x14 || x42 == x41 || x42 == x44 || x42 == x22 || x42 == x33 || x42 == x23 || x42 == x32 || x42 == x12 || x42 == x13) continue;
              int x43 = S - x13 - x23 - x33;
              if (x43 < 1 || x43 > MAX) continue;
              if (x43 == x11 || x43 == x14 || x43 == x41 || x43 == x44 || x43 == x22 || x43 == x33 || x43 == x23 || x43 == x32 || x43 == x12 || x43 == x13 || x43 == x42) continue;
              Set set21(set12); set21.erase(x12); set21.erase(x13); set21.erase(x42); set21.erase(x43);
              for (SetIterator it21 = set21.begin(); it21 != set21.end(); it21++){
                int x21 = *it21;
                int x31 = S - x11 - x21 - x41;
                if (x31 < 1 || x31 > MAX) continue;
                if (x31 == x11 || x31 == x14 || x31 == x41 || x31 == x44 || x31 == x22 || x31 == x33 || x31 == x23 || x31 == x32 || x31 == x12 || x31 == x13 || x31 == x42 || x31 == x43 || x31 == x21) continue;
                int x24 = S - x21 - x22 - x23;
                if (x24 < 1 || x24 > MAX) continue;
                if (x24 == x11 || x24 == x14 || x24 == x41 || x24 == x44 || x24 == x22 || x24 == x33 || x24 == x23 || x24 == x32 || x24 == x12 || x24 == x13 || x24 == x42 || x24 == x43 || x24 == x21 || x24 == x31) continue;
                int x34 = S - x31 - x32 - x33;
                if (x34 < 1 || x34 > MAX) continue;
                if (x34 == x11 || x34 == x14 || x34 == x41 || x34 == x44 || x34 == x22 || x34 == x33 || x34 == x23 || x34 == x32 || x34 == x12 || x34 == x13 || x34 == x42 || x34 == x43 || x34 == x21 || x34 == x31 || x34 == x24) continue;

                printf("%d %d %d %d  %d %d %d %d  %d %d %d %d  %d %d %d %d\n", x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44);
                squares++;
                }
              }
            }
          }
        }
      }
    }

  printf("CNT: %d\n", squares);

  float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
  printf("T = %.2fs\n", diff_t);

  return 0;
}