//------------------------------------------------------------------------------------------------ // // ファイル名 // plarail.c // // 作った人 // Makoto Kamada // // 更新した日 // 2019-06-15 // // ウェブサイト // https://stdkmd.net/plarail/ // //------------------------------------------------------------------------------------------------ // // タイトル // プラレールの閉路の数え上げ // // 問題 // プラレールの直線レールと曲線レールを合わせて N 本使ってできる閉路は何通りあるか。 // // 条件 // 曲線レールは直線レールの長さを半径とする中心角 45 度の円弧とする。 // ある閉路とそれを回転または裏返してできる閉路は同じと見なす。 // // 参考 // 三谷先生のツイート // https://twitter.com/jmitani/status/880070117635235840 // https://twitter.com/jmitani/status/1136228098314125313 // //------------------------------------------------------------------------------------------------ #include #include #include #include #include //レールの座標 // 右が+x、下が+yの平面座標を用いる // ┌→+x // ↓ // +y // 直線レールはS、曲線レールは進行方向に対して右に曲がるものをR、左に曲がるものをLと書く // Sは長さが2の線分、RとLは半径が2で中心角が45度の円弧とする // 右を向いているとき // Sは右に2進む // Rは右にsqrt(2)、下に2-sqrt(2)進み、右に45度方向転換する // Lは右にsqrt(2)、上に2-sqrt(2)進み、左に45度方向転換する // 　　　　　　　　　　　　　　　　　　　　　　　　○ // 　　　　　　　　　　　　　　　　　　　　　　　／│ // 　　　　　　　　　　　　　　　　　　　　　　／　│ // 　　　　　　　　　　　　　　　　　　　　　／　　│ // ●ＳＳＳＳＳＳＳＳＳ●ＲＲ────┐　　／　　　│ // 　　　　　　　　　　│　　ＲＲ　　│　／　　　　│ // 　　　　　　　　　　│　　　　ＲＲ│／　　　　　│ // 　　　　　　　　　　│　　　　　　●　　　　　　│ // 　　　　　　　　　　│　　　　　／　ＬＬ　　　　│ // 　　　　　　　　　　│　　　　／　　　＼ＬＬ　　│ // 　　　　　　　　　　│　　　／　　　　　＼　ＬＬ● // 　　　　　　　　　　│　　／　　　　　　　＼　／　 // 　　　　　　　　　　│　／　　　　　　　　　　　　 // 　　　　　　　　　　│／　　　　　　　　　　　　　 // 　　　　　　　　　　○　　　　　　　　　　　　　　 // +-----+---------------------------+-------------------------------+-------------------------------+ // | | S=0 | R=1 | L=2 | // | h +----------+----------+-----+------------+------------+-----+------------+------------+-----+ // | | dx | dy | dh | dx | dy | dh | dx | dy | dh | // +-----+----------+----------+-----+------------+------------+-----+------------+------------+-----+ // | 0 | 2 | 0 | 0 | sqrt(2) | 2-sqrt(2) | 45 | sqrt(2) | -2+sqrt(2) | -45 | // | 45 | sqrt(2) | sqrt(2) | 0 | 2-sqrt(2) | sqrt(2) | 45 | sqrt(2) | 2-sqrt(2) | -45 | // | 90 | 0 | 2 | 0 | -2+sqrt(2) | sqrt(2) | 45 | 2-sqrt(2) | sqrt(2) | -45 | // | 135 | -sqrt(2) | sqrt(2) | 0 | -sqrt(2) | 2-sqrt(2) | 45 | -2+sqrt(2) | sqrt(2) | -45 | // | 180 | -2 | 0 | 0 | -sqrt(2) | -2+sqrt(2) | 45 | -sqrt(2) | 2-sqrt(2) | -45 | // | 225 | -sqrt(2) | -sqrt(2) | 0 | -2+sqrt(2) | -sqrt(2) | 45 | -sqrt(2) | -2+sqrt(2) | -45 | // | 270 | 0 | -2 | 0 | 2-sqrt(2) | -sqrt(2) | 45 | -2+sqrt(2) | -sqrt(2) | -45 | // | 315 | sqrt(2) | -sqrt(2) | 0 | sqrt(2) | -2+sqrt(2) | 45 | 2-sqrt(2) | -sqrt(2) | -45 | // +-----+----------+----------+-----+------------+------------+-----+------------+------------+-----+ // 座標を2*a+sqrt(2)*bの形で表現する static const int coeff[3][8][5] = { { //S { 1, 0, 0, 0, 0 }, // 0 { 0, 1, 0, 1, 0 }, // 45 { 0, 0, 1, 0, 0 }, // 90 { 0, -1, 0, 1, 0 }, //135 { -1, 0, 0, 0, 0 }, //180 { 0, -1, 0, -1, 0 }, //225 { 0, 0, -1, 0, 0 }, //270 { 0, 1, 0, -1, 0 }, //315 }, { //R { 0, 1, 1, -1, 1 }, // 0 { 1, -1, 0, 1, 1 }, // 45 { -1, 1, 0, 1, 1 }, // 90 { 0, -1, 1, -1, 1 }, //135 { 0, -1, -1, 1, 1 }, //180 { -1, 1, 0, -1, 1 }, //225 { 1, -1, 0, -1, 1 }, //270 { 0, 1, -1, 1, 1 }, //315 }, { //L { 0, 1, -1, 1, -1 }, // 0 { 0, 1, 1, -1, -1 }, // 45 { 1, -1, 0, 1, -1 }, // 90 { -1, 1, 0, 1, -1 }, //135 { 0, -1, 1, -1, -1 }, //180 { 0, -1, -1, 1, -1 }, //225 { -1, 1, 0, -1, -1 }, //270 { 1, -1, 0, -1, -1 }, //315 }, }; //文字 // S=0,R=1,L=2 static const char chr[4] = { 'S', 'R', 'L', 'X' }; // メイン int main (int argc, char *argv[]) { //パラメータ int numofRails = 8; //レールの数 char outputFileName[1024] = ""; //出力ファイル名 switch (argc) { case 3: if (sscanf (argv[2], "%1023[-().0-9A-Z_a-z~]", outputFileName) != 1) { goto help; } case 2: if (sscanf (argv[1], "%d", &numofRails) != 1 || numofRails < 2 || 32 < numofRails) { //レールの数が2～32でない goto help; } break; default: help: printf ("usage: plarail \n"); return 1; } //計測開始 struct timespec startTime; clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &startTime); //最上位の3進数のリストを作る int (*ternaryList1)[8192] = calloc (8192 * (8 + 1), sizeof (ternaryList1[0][0])); //最上位の3進数のリスト ternaryList1[0][0] = 0; //0桁 ternaryList1[0][1] = -1; for (int d = 1; d <= 8; d++) { //最上位の桁数 int w = d << 1; //ビット数 int m = 0xffff >> (16 - w); //下位wビットが1 int *p = ternaryList1[d]; //d桁の3進数のリストの先頭 for (int n = 0; n <= m; n++) { //nはd桁の4進数 int x = (n >> 1) & n & 0x5555; //3を抽出する if (x) { //3があるとき n |= (x & -x) - 1; //一番右にある3の右側をすべて3にする continue; } //nはd桁の3進数 for (int i = w - 2; 0 < i; i -= 2) { if ((n & (m >> i)) < (n >> i)) { //下位w-iビットが上位w-iビットより若い goto nextN; } } x = ((n >> 1) & 0x5555) | ((n & 0x5555) << 1); //RとLを入れ替える for (int i = w - 2; 0 <= i; i -= 2) { if ((x & (m >> i)) < (n >> i)) { //RとLを入れ替えた下位w-iビットが元の上位w-iビットより若い goto nextN; } } //逆順にする x = ((n & 0x3333) << 2) | ((n >> 2) & 0x3333); x = ((x & 0x0f0f) << 4) | ((x >> 4) & 0x0f0f); x = ((x << 8) | (x >> 8)) & 0xffff; x >>= 16 - w; for (int i = w - 2; 0 <= i; i -= 2) { if ((x & (m >> i)) < (n >> i)) { //逆順にした下位w-iビットが元の上位w-iビットより若い goto nextN; } } x = ((x >> 1) & 0x5555) | ((x & 0x5555) << 1); //RとLを入れ替える for (int i = w - 2; 0 <= i; i -= 2) { if ((x & (m >> i)) < (n >> i)) { //逆順にしてRとLを入れ替えた下位w-iビットが元の上位w-iビットより若い goto nextN; } } *p++ = n; nextN: ; } *p = -1; //番兵 } //最上位以外の3進数のリストを作る int *ternaryList2 = calloc (8192, sizeof (ternaryList2[0])); //最上位以外の3進数のリスト { int *p = ternaryList2; //8桁の3進数のリストの先頭 for (int n = 0x0000; n <= 0xffff; n++) { //nは8桁の4進数 int x = (n >> 1) & n & 0x5555; //3を抽出する if (x) { //3があるとき n |= (x & -x) - 1; //一番右にある3の右側をすべて3にする continue; } //nは8桁の3進数 *p++ = n; } *p = -1; //番兵 } //8桁の3進数の閉路条件のハッシュテーブルを作る // 下位8桁の3進数の方向条件と閉路条件の符号を反転する // 下位8桁が0の3進数の閉路条件(xa,xb,ya,yb)のハッシュ関数を(((xa*7+xb)*7+ya)*7+yb)&(512-1)とする // 下位8桁が0の3進数の方向条件(h0)毎に8個のハッシュテーブルを作る // 下位8桁が0の3進数の方向条件と閉路条件に合う下位8桁を連結することで方向条件と閉路条件を満たす3進数を構成する #define CIRCUIT_HASH_SIZE 512 #define CIRCUIT_LIST_SIZE 256 // circuitHashTable[h0][g][5*i]={n,-xa,-xb,-ya,-yb} // [CIRCUIT_LIST_SIZE-1]はカウンタ。ハッシュテーブルを作るときだけ使うので邪魔にならないところに置く int (*circuitHashTable)[CIRCUIT_HASH_SIZE][CIRCUIT_LIST_SIZE] = calloc (CIRCUIT_LIST_SIZE * CIRCUIT_HASH_SIZE * 8, sizeof (circuitHashTable[0][0][0])); for (int h0 = 0; h0 < 8; h0++) { for (int g = 0; g < CIRCUIT_HASH_SIZE; g++) { circuitHashTable[h0][g][CIRCUIT_LIST_SIZE - 1] = 0; //カウンタ } } for (int n = 0x0000; n <= 0xffff; n++) { //nは8桁の4進数 int x = (n >> 1) & n & 0x5555; //3を抽出する if (x) { //3があるとき n |= (x & -x) - 1; //一番右にある3の右側をすべて3にする continue; } //nは8桁の3進数 for (int h0 = 0; h0 < 8; h0++) { //方向条件と閉路条件を計算する int xa = 0; int xb = 0; int ya = 0; int yb = 0; int h = h0; for (int i = 2 * 8 - 2; 0 <= i; i -= 2) { const int *c = coeff[(int) (n >> i) & 3][h]; xa += c[0]; xb += c[1]; ya += c[2]; yb += c[3]; h = (h + c[4]) & 7; } if (h == 0) { //方向条件 xa = -xa; xb = -xb; ya = -ya; yb = -yb; //ハッシュテーブルを更新する int g = (((xa * 7 + xb) * 7 + ya) * 7 + yb) & (CIRCUIT_HASH_SIZE - 1); //ハッシュコード int *p = circuitHashTable[h0][g]; int i = p[CIRCUIT_LIST_SIZE - 1]; p[i ] = n; p[i + 1] = xa; p[i + 2] = xb; p[i + 3] = ya; p[i + 4] = yb; p[CIRCUIT_LIST_SIZE - 1] = i + 5; } //if 方向条件 } //for h0 } //for n for (int h0 = 0; h0 < 8; h0++) { for (int g = 0; g < CIRCUIT_HASH_SIZE; g++) { circuitHashTable[h0][g][circuitHashTable[h0][g][CIRCUIT_LIST_SIZE - 1]] = -1; //番兵 } } //閉路の配列を確保する // 8 10 12 14 16 18 20 22 24 26 28 30 32 レールの数 n // 1 1 4 9 42 161 847 4739 29983 198683 1375928 9786630 35658701? 閉路の数 // 1 2 6 18 61 238 1058 5380 31012 201439 1465610 11879054 106724526 確保する数 floor(pow(0.03*n+0.8,n+0.7)) unsigned long long *uniqueCircuitArray = calloc ((size_t) floor (pow (0.03 * (double) numofRails + 0.8, (double) numofRails + 0.7)), sizeof (unsigned long long)); //閉路の配列 //閉路を作る int circuitWidth = numofRails << 1; //閉路のビット数 unsigned long long circuitMask = 0xffffffffffffffffULL >> (64 - circuitWidth); //閉路のマスク。(1ULL<> 1) | x) & 0x5555555555555555ULL) != 0) { //00がある continue; } } //先頭がSSでないときSSを含むものを除く if (circuitSRSS <= n00) { unsigned long long x = n00 | invertedCircuitMask2; if ((~((x >> 3) | (x >> 2) | (x >> 1) | x) & 0x5555555555555555ULL) != 0) { //0000がある continue; } } //先頭がSSSでないときSSSを含むものを除く if (circuitSSRSS <= n00) { unsigned long long x = n00 | invertedCircuitMask2; if ((~((x >> 5) | (x >> 4) | (x >> 3) | (x >> 2) | (x >> 1) | x) & 0x5555555555555555ULL) != 0) { //000000がある continue; } } while (0 <= *tp1) { unsigned long long n0 = n00 | ((unsigned long long) *tp1++ << 16); //ビット31～16を連結する //n0はnumofRails桁の先頭が2ではない下位8桁が0の3進数 //先頭がSでないときSを含むものを除く if (circuitRSS <= n0) { unsigned long long x = n0 | invertedCircuitMask1; if ((~((x >> 1) | x) & 0x5555555555555555ULL) != 0) { //00がある continue; } } //先頭がSSでないときSSを含むものを除く if (circuitSRSS <= n0) { unsigned long long x = n0 | invertedCircuitMask1; if ((~((x >> 3) | (x >> 2) | (x >> 1) | x) & 0x5555555555555555ULL) != 0) { //0000がある continue; } } //先頭がSSSでないときSSSを含むものを除く if (circuitSSRSS <= n0) { unsigned long long x = n0 | invertedCircuitMask1; if ((~((x >> 5) | (x >> 4) | (x >> 3) | (x >> 2) | (x >> 1) | x) & 0x5555555555555555ULL) != 0) { //000000がある continue; } } //方向条件と閉路条件を計算する int xa = 0; int xb = 0; int ya = 0; int yb = 0; int h = 0; for (int i = circuitWidth - 2; 2 * 8 <= i; i -= 2) { //下位8桁はここにはない const int *c = coeff[(int) (n0 >> i) & 3][h]; xa += c[0]; xb += c[1]; ya += c[2]; yb += c[3]; h = (h + c[4]) & 7; } if (xa < -8 || 8 < xa || xb < -8 || 8 < xb || ya < -8 || 8 < ya || yb < -8 || 8 < yb) { //閉路条件 continue; } for (int *p = circuitHashTable[h][(((xa * 7 + xb) * 7 + ya) * 7 + yb) & (CIRCUIT_HASH_SIZE - 1)]; 0 <= *p; p += 5) { if (p[1] == xa && p[2] == xb && p[3] == ya && p[4] == yb) { //閉路条件 unsigned long long circuit = n0 | (unsigned long long) *p; //ビット15～0を連結する //circuitはnumofRails桁の3進数で方向条件と閉路条件を満たしている //先頭がSでないときSを含むものを除く if (circuitRSS <= circuit) { unsigned long long x = circuit | invertedCircuitMask; x = ~((x >> 1) | x) & 0x5555555555555555ULL; //00を抽出する if (x != 0) { //00がある circuit |= (x & -x) - 1ULL; //一番右にある00の右側をすべて1にする goto nextCircuit; } } //先頭がSSでないときSSを含むものを除く if (circuitSRSS <= circuit) { unsigned long long x = circuit | invertedCircuitMask; x = ~((x >> 3) | (x >> 2) | (x >> 1) | x) & 0x5555555555555555ULL; //0000を抽出する if (x != 0) { //0000がある circuit |= (x & -x) - 1ULL; //一番右にある0000の右側をすべて1にする goto nextCircuit; } } //先頭がSSSでないときSSSを含むものを除く if (circuitSSRSS <= circuit) { unsigned long long x = circuit | invertedCircuitMask; x = ~((x >> 5) | (x >> 4) | (x >> 3) | (x >> 2) | (x >> 1) | x) & 0x5555555555555555ULL; //000000を抽出する if (x != 0) { //000000がある circuit |= (x & -x) - 1ULL; //一番右にある000000の右側をすべて1にする goto nextCircuit; } } //ローテートすると若くなるものを除く { unsigned long long x = circuit; for (int i = circuitWidth - 2; 0 < i; i -= 2) { if ((unsigned long long) (((x << (circuitWidth - i)) | (x >> i)) & circuitMask) < circuit) { goto nextCircuit; } } } //RとLを入れ替えてローテートすると若くなるものを除く { unsigned long long x = circuit; x = ((x & 0x5555555555555555ULL) << 1) | ((x >> 1) & 0x5555555555555555ULL); //RとLを入れ替える if (x < circuit) { goto nextCircuit; } for (int i = circuitWidth - 2; 0 < i; i -= 2) { if ((unsigned long long) (((x << (circuitWidth - i)) | (x >> i)) & circuitMask) < circuit) { goto nextCircuit; } } } //逆順に繋いだ閉路を作る unsigned long long reversedCircuit; { //16ビットずつテーブルを使って反転するより速い unsigned long long x = circuit; //x = ((x & 0x5555555555555555ULL) << 1) | ((x >> 1) & 0x5555555555555555ULL); x = ((x & 0x3333333333333333ULL) << 2) | ((x >> 2) & 0x3333333333333333ULL); x = ((x & 0x0f0f0f0f0f0f0f0fULL) << 4) | ((x >> 4) & 0x0f0f0f0f0f0f0f0fULL); x = ((x & 0x00ff00ff00ff00ffULL) << 8) | ((x >> 8) & 0x00ff00ff00ff00ffULL); x = ((x & 0x0000ffff0000ffffULL) << 16) | ((x >> 16) & 0x0000ffff0000ffffULL); x = (x << 32) | (x >> 32); reversedCircuit = x >> (64 - circuitWidth); } //逆順に繋いでローテートすると若くなるものを除く { unsigned long long x = reversedCircuit; if (x < circuit) { goto nextCircuit; } for (int i = circuitWidth - 2; 0 < i; i -= 2) { if ((unsigned long long) (((x << (circuitWidth - i)) | (x >> i)) & circuitMask) < circuit) { goto nextCircuit; } } } //逆順に繋いでRとLを入れ替えてローテートすると若くなるものを除く { unsigned long long x = reversedCircuit; x = ((x & 0x5555555555555555ULL) << 1) | ((x >> 1) & 0x5555555555555555ULL); //RとLを入れ替える if (x < circuit) { goto nextCircuit; } for (int i = circuitWidth - 2; 0 < i; i -= 2) { if ((unsigned long long) (((x << (circuitWidth - i)) | (x >> i)) & circuitMask) < circuit) { goto nextCircuit; } } } //重複を除いた閉路を記録する uniqueCircuitArray[numofUniqueCircuits++] = circuit; nextCircuit: ; } //if 閉路条件 } //for p } //while tp1 tp1 = ternaryList2; } //while tp2 tp2 = ternaryList2; } //while tp3 //計測終了 struct timespec endTime; clock_gettime (CLOCK_PROCESS_CPUTIME_ID, &endTime); time_t sec = endTime.tv_sec - startTime.tv_sec; long nsec = endTime.tv_nsec - startTime.tv_nsec; if (nsec < 0) { sec--; nsec += 1000000000L; } double elapsedTime = (double) sec + (double) nsec / 1000000000.0; //閉路を出力する // 閉路は昇順に見つかるのでソートする必要はない if (outputFileName[0] != '\0') { FILE *fp; if (strcmp (outputFileName, "-") == 0) { fp = stdout; } else { fp = fopen (outputFileName, "wb"); if (fp == NULL) { printf ("%s cannot be opened for writing.\n", outputFileName); } } if (fp != NULL) { for (size_t index = 0; index < numofUniqueCircuits; index++) { unsigned long long circuit = uniqueCircuitArray[index]; char text[256]; char *p = text; for (int i = circuitWidth - 2; 0 <= i; i -= 2) { *p++ = chr[(int) (circuit >> i) & 3]; } *p = '\0'; fprintf (fp, "%s\n", text); //比較するとき邪魔なので個々の閉路に番号を振らない } if (fp != stdout) { fclose (fp); printf ("%s has been updated.\n", outputFileName); } } } //閉路の数と所要時間を出力する printf ("It took %.3lf seconds to find %zd circuits made by connecting %d rails.\n", elapsedTime, numofUniqueCircuits, numofRails); return 0; }