Pull to refresh

Comments 28

А можно ли сгенерировать один квадрат (для нечетной стороны это почти тривиально, известным алгоритмом), а потом тасовать строки и столбцы всеми возможными способами?

Мне кажется — так мы получим не все возможные квадраты.
Спасибо за идею, как уже сказали выше, сумма диагоналей не совпадет.
Спасибо за статью, хорошую задачу подобрали.
1) ИМХО Вы очень зря используете 16 процессоров на 16 потоков. Во-первых, у процессора 32 потока, так что половина из них просто стоят, во-вторых, если у Вас десктопная видеокарта, то процессоров там явно больше тысячи. Странно сравнивать CPU и GPU, задействуя при этом все потоки процессора, и игнорируя большую часть ресурсов видеокарты.
2) Думаю, если Вы откажетесь в CPU реализации от использования std::set, она взлетит как гордый орел.
3) Еще на счет CPU, не думали задействовать векторные инструкции? Все таки без них потенциал не полностью раскрыт.
И укажите, пожалуйста, железо, на котором проводились тесты.
Спасибо за советы. Железо нормальное: i7 + GeForce 1060, разумеется, выжать можно больше, осталось придумать как :)

Разбиение на 16 потоков и процессов я взял потому, что это легко «ложится» на алгоритм — индекс потока/процесса равен ячейке матрицы. Так-то да, можно больше задействовать.

Без std:set попробую, спасибо, еще кстати некий прирост скорости получился когда заменил int на char, возможно были лишние преобразования типов, или просто кеш эффективнее используется.
И еще раз спасибо :) Заменил вычисление через 16 блоков и 16 процессов на 256 блоков по 1 процессу, время стало 1.2 вместо 2.6, еще в 2 раза меньше.

Добавил такой вариант в текст.

Я кстати не нашел ответа о максимально допустимом числе cuda-блоков на видеокарте — оно же «железом» должно определяться?
Вы можете узнать эти ограничения через cudaGetDeviceProperties. Но они совершенно неадекватные, просто для галочки. На 940MX они составляют 2147483647 x 65536 x 65536, видимо, только для того, что бы какие-то счетчики не переполнились. Так что можно запускать на сколько угодно блоков, просто одновременно будут исполняться только 1280 (на Geforce 1060).
Кстати, о самом то важном и не сказал, измерять время в CUDA-реализации с помощью clock() — не очень правильно, в CUDA для этого есть cudaEvent_t.
Если написать универсальную функцию для любого N с обращением к матрице N х N, замедлит ли это программу (м.б. наоборот ускорит)?
Имхо да, «разворачивание» цикла до простых линейных операций делает код быстрее, так что универсальный вариант будет работать медленнее.
Мне кажется, что с учетом симметрий, число уникальных квадратов должно быть 220 (7040 ÷ 32).
А почему 32? Квадратов 3х3 всего 8, но все они отражения или повороты одного, так что возможно надо на 8 делить.
Нет, 3 на 3 здесь не причем. 32 = 8×4. 8 — за счет поворотов и отражений, 4 — за счет двух внутренних симметрий.
Да, внутренние симметрии позволяют из квадрата с номерами ячеек

1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16

создать квадрат

6 5 8 7
2 1 4 3
14 13 16 15
10 9 12 11

и квадрат

1 3 2 4
9 11 10 12
5 7 6 8
13 15 14 16

Ну и комбинация двух этих преобразований дает еще один вариант.

Я ничего не понимаю в С++, но немного переработал вашу программу, получив ускорение в 300 раз. В одном потоке, на CPU.


Ключевая идея заключается в том, что в магическом квадрате 4х4 помимо вертикалей, горизонталей и диагоналей магическую сумму дает также сумма угловых ячеек и сумма центральных ячеек. Это позволяет намного раньше останавливать перебор вариантов.


Скрытый текст
#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;
}

Мне кажется, что если вместо Set использовать просто битовую маску, то перебор ускорится еще больше. А потом уже можно и многопоточность, и GPU прикрутить.

Спасибо, интересно. А как доказывается это свойство, и распространяется ли оно на квадраты большей размерности?

Идея отказаться от Set была, да. Суть использования Set в том, чтобы заведомо не перебирать и не проверять уже повторяющиеся числа, например 1 с 1.

Во-первых, x22 + x23 + x32 + x33 = 2S — (x21 + x24 + x31 + x34) = 2S — (2S — (x11 + x14 + x41 + x44)) = x11 + x14 + x41 + x44. Но (x22 + x23 + x32 + x33) + (x11 + x14 + x41 + x44) = 2S как сумма двух диагоналей. Отсюда x22 + x23 + x32 + x33 = x11 + x14 + x41 + x44 = S. Насколько я помню, непосредственно это свойство на квадраты высших порядков не обобщается, но я не специалист.

Спасибо, добавил ваше описание в текст. Если аналогичный подход сработает для 5х5 и выше, это заметно ускорит поиск.

Замена Set на битовую маску дала ускорение еще в 20 раз, программа отрабатывает за пару сотых секунды.


Скрытый текст
#include <set>
#include <stdio.h>
#include <ctime>
// #include "stdafx.h"

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 mask = 0;
  for (int x11 = 1; x11 <= MAX; x11++) {
    mask += 1 << x11;
    for (int x14 = 1; x14 <= MAX; x14++) {
      if (mask & (1 << x14)) continue;
      mask += 1 << x14;
      for (int x41 = 1; x41 <= MAX; x41++) {
        if (mask & (1 << x41)) continue;
        int x44 = S - x11 - x14 - x41;
        if (x44 < 1 || x44 > MAX) continue;
        if (mask & (1 << x44) || x44 == x41) continue;
        mask += 1 << x41;
        mask += 1 << x44;
        for (int x22 = 1; x22 <= MAX; x22++){
          if (mask & (1 << x22)) continue;
          int x33 = S - x11 - x22 - x44;
          if (x33 < 1 || x33 > MAX) continue;
          if (mask & (1 << x33) || x33 == x22) continue;
          mask += 1 << x22;
          mask += 1 << x33;
          for (int x23 = 1; x23 <= MAX; x23++){
            if (mask & (1 << x23)) continue;
            int x32 = S - x14 - x23 - x41;
            if (x32 < 1 || x32 > MAX) continue;
            if (mask & (1 << x32) || x32 == x23) continue;
            mask += 1 << x23;
            mask += 1 << x32;
            for (int x12 = 1; x12 <= MAX; x12++){
              if (mask & (1 << x12)) continue;
              int x13 = S - x11 - x12 - x14;
              if (x13 < 1 || x13 > MAX) continue;
              if (mask & (1 << x13) || x12 == x13) continue;
              int x42 = S - x12 - x22 - x32;
              if (x42 < 1 || x42 > MAX) continue;
              if (mask & (1 << x42) || x42 == x12 || x42 == x13) continue;
              int x43 = S - x13 - x23 - x33;
              if (x43 < 1 || x43 > MAX) continue;
              if (mask & (1 << x43) || x43 == x12 || x43 == x13 || x43 == x42) continue;
              mask += 1 << x12;
              mask += 1 << x13;
              mask += 1 << x42;
              mask += 1 << x43;
              for (int x21 = 1; x21 <= MAX; x21++){
                if (mask & (1 << x21)) continue;
                int x31 = S - x11 - x21 - x41;
                if (x31 < 1 || x31 > MAX) continue;
                if (mask & (1 << x31) || x31 == x21) continue;
                int x24 = S - x21 - x22 - x23;
                if (x24 < 1 || x24 > MAX) continue;
                if (mask & (1 << x24) || x24 == x21 || x24 == x31) continue;
                int x34 = S - x31 - x32 - x33;
                if (x34 < 1 || x34 > MAX) continue;
                if (mask & (1 << x34) || x34 == x21 || x34 == x24 || x34 == x31) 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++;
                }
              mask -= 1 << x12;
              mask -= 1 << x13;
              mask -= 1 << x42;
              mask -= 1 << x43;
              }
            mask -= 1 << x23;
            mask -= 1 << x32;
            }
          mask -= 1 << x22;
          mask -= 1 << x33;
          }
        mask -= 1 << x41;
        mask -= 1 << x44;
        }
      mask -= 1 << x14;
      }
    mask -= 1 << x11;
    }

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

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

  return 0;
}
А говорите, не знаете C ;)

Респект, отличный результат.

Хе-хех… Есть чрезвычайно простой алгоритм для создания "магических квадратов". Никакого перебора в нём вообще нет. Помню, в детстве у меня была детская энциклопедия по математике (от Аванта+), и в ней был описан алгоритм построения ("Меланхолия и магия в квадрате"). Заболев на 2-м курсе воспалением лёгких, лежал в больничке и, помимо решения японских кроссвордов, в тетрадке писал код на паскале, строящий "магические квадраты" заданной размерности. Неплохо развлёкся :)
Впрочем, переборный способ, наверное, неплох для демонстрации многопоточных техник и вычислений на GPU.

Да, методы создания квадратов разных размерностей известны еще со средних веков, но здесь-то вопрос был в нахождении всех вариантов.
Кстати, подсчитать общее количество вариантов для квадрата 4x4 можно вообще без компьютера — с помощью головы, бумаги и карандаша где-то за половину часа.

Спасибо, за задачу я попробую сам сделать алгоритм для магического квадрата, а затем прочитаю вашу статью))
Взял старенький ноут с Core 2 Duo 1.66GHz десятилетней давности, скомпилировал самое первое решение автора под ubuntu 18.04 32-bit gcc 5.4, выкинув одну строчку с #include «stdafx.h», ибо мешалась. Получил T = 309, поставил ключик -O2 для gcc получил уже T~ 85-90. Так что не очень понятно почему i7 так плохо выступил.

Лет десять не программировал, но написал своё решение. Может где ошибся конечно, но без разворачивания циклов ручками и прочей ручной оптимизации, получил время T между 18-19 на всё том же стареньком ноуте, а с ключиком -O2 время было T ~ 4-5. Что я не так понял для однопоточного случая?
своё решение:
#include <stdio.h>
#include #define N 4
#define MAX (N*N)
#define S ((N)*((N)*(N)+1)/2)

int main(int argc,char *argv[]) {
const clock_t begin_time = clock();
int squares = 0;
struct Ts{
int pole[MAX]; // поле с цифирьками
int numb[MAX]; // где стоит каждое число
int* ocher[MAX*2]; // какие изменения делались
int o; // сколько чисел поменяли
int num; // текущее число
int last[MAX]; // список предположений
int l; // указатель последнего предположения
Ts():o(0),num(0),l(0){
for (int i=0;i<MAX;i++) { //начальные значения всё в ноль
pole[i]=0;
numb[i]=0;
ocher[i]=NULL;
ocher[i+MAX]=NULL;
last[i]=0;
}
}
inline bool set(int m, int n) {
ocher[o++] = &numb[n-1]; numb[n-1]=m+1; // запоминаем что модифицировали, между корректными точками
ocher[o++] = &pole[m]; pole[m]=n;
return check(m); // возможна рекурсия...
}
inline bool add(int m, int n) {
last[l++]=n-1; // записываем куда вернёмся в случае roolback
return set(m,n);
}
inline bool test(int m,int s){
if (m==-1) return s==S; // если пустых ячеек нет, то сумма всех должна равняться заранее известному числу
s = S - s;
if (s>MAX || s<=0 || numb[s-1]) return false;
return set(m,s); // если единственная оставшаяся ячейка проходит проверки записываем в неё правильное значение вне очереди
}
bool diag1() { // проверяем главную диагональ
int m=-1,sum=0;
for (int i=0;i<N;i++) {
int n=i+i*N;
if (pole[n]) sum+=pole[n];
else {
if(m==-1) m=n;
else return true; // если две пустых ячейки ничего дальше считать не удем, рано.
}
}
return test(m,sum);
}
bool diag2() { // проверяем другую диагональ
int m=-1,sum=0;
for (int i=0;i<N;i++) {
int n=N-1-i+i*N;
if(! pole[n]) {
if(m==-1) m=n;
else return true;
} else sum+=pole[n];
}
return test(m,sum);
}
bool checkS(int s){ // проверяем строку
int m=-1,sum=0;
for (int i=0;i<N;i++) {
int n=i+s*N;
if(! pole[n]) {
if(m==-1) m=n;
else return true;
} else sum+=pole[n];
}
return test(m,sum);
}
bool checkH(int h){ // проверяем столец
int m=-1,sum=0;
for (int i=0;i<N;i++) {
int n=h+i*N;
if(! pole[n]) {
if(m==-1) m=n;
else return true;
} else sum+=pole[n];
}
return test(m,sum);
}
inline bool check(int m) { // решаем что удем проверять после очередного предположения
int s = m/N;
int h = m%N;
return checkS(s) && checkH(h) && ((h==s)?diag1():((h==N-1-s)?diag2():true));
}
int roolback() { // откат предыдущего преположения
int lnum=last[--l];
if (num<MAX) numb[num]=0; // текущее значение уже не интересно
for (o--;ocher[o]!=&numb[lnum];o--)
*ocher[o]=0;
//*ocher[o]=0;
numb[lnum]--;
return lnum;
}
inline void next() { // следующая позиция цифры
numb[num]++;
}
inline void print() { // печать поля
printf("[");
for (int i=0;i<MAX-1;i++)
printf((i%N==N-1)?"%d ":"%d ",pole[i]);
printf("%d]\n",pole[MAX-1]);
};
} p;
for (;;p.next()) { // не очень красиво, вечный цикл
int m = p.numb[p.num];
if (m>=MAX)
if (p.num==0) break; // всё конец, все варианты перерали
else p.num=p.roolback(); // не смогли для текущего числа найти позицию
else if (p.pole[m]==0) {
if (p.add(m,++p.num)) { // ставим число в новую позицию
while (p.num<MAX && p.numb[p.num]>0) p.num++; // а вдруг мы уже всё поставили?
if (p.num==MAX) {
p.print();squares++; // новый квадрат
p.num=p.roolback(); // ищем следующий
} else p.numb[p.num]=-1; // не очень красиво, коменсируем следующий next
}
else p.num=p.roolback(); // не получилось
}
}

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

float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
printf("T = %.2fs\n", diff_t);
return 0;
}
Спасибо. Я тоже пробовал, заметно быстрее без std::set, все же медленный он для таких целей.

С первым вариантом, я нашел у себя ошибку, где-то в цикле было char а где-то int, часть времени «съелось» на преобразовании типов. После исправления стало 60с.
Кстати, перепроверил сейчас еще раз — все еще более запутанно :)

Если в качестве базовой переменной использовать int, время выполнения 100с и 19с в одно и многопоточном варианте. Если заменить int на char, то однопоточный вариант дает 60с, а многопоточный наоборот ухудшается до 34с.
(Компилятор: Visual Studio 2015)

Думаю, в однопоточном варианте программа с char лучше помещается в кеш, и работает быстрее.
Можно очень просто ускорить вычисление, просто выводя все данные в .txt файл:

#include <set>
#include <stdio.h>
#include <ctime>
#include <fstream>
#include "stdafx.h"

...

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

        char *fileName = "magicSquares.txt"
        std::ofstream magicSquares(fileName);
        magicSquares.open(fileName);
        freopen(fileName, "w", stdout);

	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++)
       {
	       ...
       }
	
	printf("CNT: %d\n", squares);

	float diff_t = float(clock() - begin_time)/CLOCKS_PER_SEC;
	printf("T = %.2fs\n", diff_t);
        
        magicSquares.close();
       
	return 0;
}
Only those users with full accounts are able to leave comments. Log in, please.