desdif.cpp

// g++ -O3 desdif.cpp -fPIC -shared -o desdif.so -Wall -Wextra -g3

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <limits>
#include <map>
#include <queue>
#include <set>
#include <unordered_map>
#include <unordered_set>
#include <vector>

#include "prev_table.h"
#include <chrono>
#include <memory>
#include <random>

//std::mt19937_64 rnd(std::chrono::system_clock::now().time_since_epoch().count());
std::mt19937_64 rnd(42);

#define eprintf(format, ...)                             \
    fprintf(stderr, "\x1b[3%dm[%s:%d(%s)]\x1b[m" format, \
        __LINE__ % 6 + 1, __FILE__, __LINE__,            \
        __PRETTY_FUNCTION__, ##__VA_ARGS__)

#define ewatchi(EXPR) eprintf(#EXPR " = %lld\n", (long long)(EXPR));
#define ewatchx(EXPR) eprintf(#EXPR " = %llx\n", (long long)(EXPR));
#define ewatchf(EXPR) eprintf(#EXPR " = %Lf\n", (long double)(EXPR));
#define ewatcha_gen(EXPR, NNN, TYPE, TYPEFIX)                                 \
    {                                                                         \
        eprintf(#EXPR " = [%lld] {\n", (long long)(NNN));                     \
        for (int __i__ = 0; __i__ < NNN; __i__ += 8) {                        \
            for (int __j__ = 0; __j__ < 8 && __i__ + __j__ < NNN; __j__++)    \
                fprintf(stderr, " %" TYPEFIX, (TYPE)((EXPR)[__i__ + __j__])); \
            fprintf(stderr, "\n");                                            \
        }                                                                     \
        fprintf(stderr, "}\n");                                               \
    }
#define ewatcha(EXPR, NNN) ewatcha_gen(EXPR, NNN, long long, "5lld")
#define ewatchax(EXPR, NNN) ewatcha_gen(EXPR, NNN, long long, "17llx")
#define ewatchaf(EXPR, NNN) ewatcha_gen(EXPR, NNN, long double, "3.6Lf")
#define ewatchs(EXPR) eprintf(#EXPR " = %s\n", (const char*)(EXPR));

using std::array;
using std::map;
using std::pair;
using std::unordered_map;
using std::vector;

static constexpr uint64_t INVALID_KEY = 1;

static constexpr int ROUNDS = 8;
static constexpr uint8_t IP[64] = { 62, 54, 46, 38, 30, 22, 14, 6, 60, 52, 44, 36, 28, 20, 12, 4, 58, 50, 42, 34, 26, 18, 10, 2, 56, 48, 40, 32, 24, 16, 8, 0, 63, 55, 47, 39, 31, 23, 15, 7, 61, 53, 45, 37, 29, 21, 13, 5, 59, 51, 43, 35, 27, 19, 11, 3, 57, 49, 41, 33, 25, 17, 9, 1 };
static constexpr uint8_t IP_I[64] = { 31, 63, 23, 55, 15, 47, 7, 39, 30, 62, 22, 54, 14, 46, 6, 38, 29, 61, 21, 53, 13, 45, 5, 37, 28, 60, 20, 52, 12, 44, 4, 36, 27, 59, 19, 51, 11, 43, 3, 35, 26, 58, 18, 50, 10, 42, 2, 34, 25, 57, 17, 49, 9, 41, 1, 33, 24, 56, 16, 48, 8, 40, 0, 32 };
static constexpr uint8_t SBOX[8][4][16] = {
    {
        { 7, 12, 4, 10, 11, 6, 13, 0, 2, 5, 15, 9, 8, 3, 1, 14 },
        { 0, 5, 7, 9, 14, 3, 11, 12, 15, 6, 4, 10, 2, 13, 8, 1 },
        { 2, 15, 11, 12, 7, 9, 4, 10, 8, 3, 6, 5, 1, 14, 13, 0 },
        { 15, 10, 2, 5, 1, 12, 8, 6, 3, 13, 9, 0, 4, 7, 14, 11 },
    },
    {
        { 15, 9, 6, 3, 1, 4, 12, 10, 8, 14, 13, 0, 7, 11, 2, 5 },
        { 12, 3, 15, 6, 2, 8, 1, 13, 11, 0, 4, 9, 14, 5, 7, 10 },
        { 0, 10, 5, 9, 14, 3, 11, 4, 7, 1, 2, 12, 13, 6, 8, 15 },
        { 11, 13, 12, 0, 5, 14, 2, 7, 1, 6, 15, 10, 8, 3, 4, 9 },
    },
    {
        { 5, 8, 6, 13, 9, 3, 15, 4, 0, 11, 12, 2, 7, 14, 10, 1 },
        { 11, 4, 12, 3, 0, 10, 6, 15, 14, 1, 2, 13, 9, 7, 5, 8 },
        { 11, 13, 1, 10, 2, 4, 12, 7, 6, 8, 15, 5, 9, 3, 0, 14 },
        { 8, 2, 6, 13, 11, 7, 1, 4, 5, 15, 9, 10, 0, 12, 14, 3 },
    },
    {
        { 14, 8, 0, 13, 7, 1, 9, 2, 11, 4, 6, 3, 12, 10, 5, 15 },
        { 11, 2, 6, 8, 13, 4, 0, 7, 1, 14, 15, 5, 10, 3, 12, 9 },
        { 5, 15, 3, 10, 9, 12, 14, 1, 6, 8, 13, 4, 0, 7, 11, 2 },
        { 12, 9, 5, 3, 0, 10, 11, 4, 15, 2, 8, 14, 6, 13, 1, 7 },
    },
    {
        { 4, 1, 14, 11, 2, 12, 13, 7, 3, 10, 5, 0, 8, 15, 6, 9 },
        { 7, 10, 2, 12, 4, 15, 11, 1, 13, 0, 14, 9, 3, 5, 8, 6 },
        { 2, 15, 5, 6, 8, 3, 14, 0, 4, 9, 11, 12, 13, 10, 1, 7 },
        { 13, 6, 8, 5, 3, 0, 4, 10, 1, 15, 7, 2, 14, 9, 11, 12 },
    },
    {
        { 3, 0, 9, 7, 5, 12, 6, 10, 8, 11, 4, 14, 15, 2, 1, 13 },
        { 5, 6, 14, 0, 2, 11, 9, 12, 15, 8, 3, 13, 4, 7, 10, 1 },
        { 9, 14, 4, 8, 15, 2, 3, 13, 7, 0, 1, 11, 10, 5, 12, 6 },
        { 2, 13, 9, 6, 4, 8, 15, 1, 12, 7, 10, 0, 3, 14, 5, 11 },
    },
    {
        { 2, 12, 15, 10, 4, 9, 1, 6, 13, 3, 0, 5, 7, 14, 11, 8 },
        { 11, 7, 2, 4, 13, 10, 8, 1, 0, 12, 9, 15, 14, 3, 5, 6 },
        { 8, 5, 3, 0, 13, 6, 14, 9, 2, 15, 12, 10, 11, 1, 7, 4 },
        { 6, 9, 8, 7, 11, 0, 5, 12, 13, 10, 2, 4, 1, 15, 14, 3 },
    },
    {
        { 11, 5, 6, 10, 1, 12, 13, 3, 4, 9, 15, 0, 2, 7, 8, 14 },
        { 8, 3, 5, 0, 11, 6, 14, 9, 15, 10, 12, 7, 1, 13, 2, 4 },
        { 14, 0, 9, 15, 2, 5, 7, 10, 13, 6, 3, 12, 8, 11, 4, 1 },
        { 4, 15, 2, 12, 7, 9, 1, 6, 8, 3, 5, 10, 14, 0, 11, 13 },
    },
};
static constexpr uint8_t P[32] = { 15, 6, 19, 20, 28, 11, 27, 16, 0, 14, 22, 25, 4, 17, 30, 9, 1, 7, 23, 13, 31, 26, 2, 8, 18, 12, 29, 5, 21, 10, 3, 24 };
static constexpr uint8_t P_I[32] = { 8, 16, 22, 30, 12, 27, 1, 17, 23, 15, 29, 5, 25, 19, 9, 0, 7, 13, 24, 2, 3, 28, 10, 18, 31, 11, 21, 6, 4, 26, 14, 20 };
static constexpr uint8_t KEY_PERM[8][48] = {
    { 14, 53, 38, 60, 55, 23, 39, 63, 6, 15, 21, 46, 5, 37, 30, 31, 44, 62, 61, 7, 36, 29, 22, 47, 18, 28, 33, 50, 35, 4, 41, 26, 3, 51, 17, 27, 59, 19, 34, 57, 9, 20, 43, 10, 11, 58, 49, 25 },
    { 6, 45, 30, 52, 47, 15, 31, 55, 61, 7, 13, 38, 60, 29, 22, 23, 36, 54, 53, 62, 63, 21, 14, 39, 10, 20, 25, 42, 27, 57, 33, 18, 28, 43, 9, 19, 51, 11, 26, 49, 1, 12, 35, 2, 3, 50, 41, 17 },
    { 53, 29, 14, 36, 31, 62, 15, 39, 45, 54, 60, 22, 44, 13, 6, 7, 55, 38, 37, 46, 47, 5, 61, 23, 59, 4, 9, 26, 11, 41, 17, 2, 12, 27, 58, 3, 35, 28, 10, 33, 50, 57, 19, 51, 20, 34, 25, 1 },
    { 37, 13, 61, 55, 15, 46, 62, 23, 29, 38, 44, 6, 63, 60, 53, 54, 39, 22, 21, 30, 31, 52, 45, 7, 43, 49, 58, 10, 28, 25, 1, 51, 57, 11, 42, 20, 19, 12, 59, 17, 34, 41, 3, 35, 4, 18, 9, 50 },
    { 21, 60, 45, 39, 62, 30, 46, 7, 13, 22, 63, 53, 47, 44, 37, 38, 23, 6, 5, 14, 15, 36, 29, 54, 27, 33, 42, 59, 12, 9, 50, 35, 41, 28, 26, 4, 3, 57, 43, 1, 18, 25, 20, 19, 49, 2, 58, 34 },
    { 5, 44, 29, 23, 46, 14, 30, 54, 60, 6, 47, 37, 31, 63, 21, 22, 7, 53, 52, 61, 62, 55, 13, 38, 11, 17, 26, 43, 57, 58, 34, 19, 25, 12, 10, 49, 20, 41, 27, 50, 2, 9, 4, 3, 33, 51, 42, 18 },
    { 52, 63, 13, 7, 30, 61, 14, 38, 44, 53, 31, 21, 15, 47, 5, 6, 54, 37, 36, 45, 46, 39, 60, 22, 28, 1, 10, 27, 41, 42, 18, 3, 9, 57, 59, 33, 4, 25, 11, 34, 51, 58, 49, 20, 17, 35, 26, 2 },
    { 36, 47, 60, 54, 14, 45, 61, 22, 63, 37, 15, 5, 62, 31, 52, 53, 38, 21, 55, 29, 30, 23, 44, 6, 12, 50, 59, 11, 25, 26, 2, 20, 58, 41, 43, 17, 49, 9, 28, 18, 35, 42, 33, 4, 1, 19, 10, 51 },
};

static constexpr uint16_t CROSS_AV[8][64] = {
    { 1, 65128, 60136, 61438, 65256, 38909, 61438, 64510, 65256, 40687, 27646, 55574, 48872, 65531, 23454, 27134, 65256, 55295, 30718, 32767, 38910, 65531, 53246, 32703, 65406, 32508, 65490, 31743, 40822, 65019, 44990, 65495, 65256, 38910, 32119, 64991, 34815, 49017, 32731, 32255, 57341, 57343, 65494, 11191, 65397, 57215, 27071, 65430, 24575, 49018, 63999, 65023, 49135, 63485, 65535, 61435, 65531, 45023, 28607, 65527, 65279, 65535, 57303, 64511 },
    { 1, 64712, 65248, 63334, 63144, 40575, 57246, 46590, 61096, 59307, 64510, 63998, 65192, 63475, 46070, 16382, 64744, 48951, 65502, 20479, 65534, 65279, 65534, 57087, 30590, 65023, 63486, 60413, 49018, 49135, 48990, 30695, 65256, 65406, 31615, 52734, 57311, 63452, 65527, 59886, 64767, 46837, 61375, 48470, 65469, 64443, 56310, 32214, 65407, 65519, 44765, 31679, 61407, 65531, 65530, 65527, 65407, 63215, 65471, 65535, 59390, 32725, 63359, 65257 },
    { 1, 61160, 65128, 29694, 63208, 61375, 61950, 61302, 65256, 30719, 63358, 57214, 65224, 65527, 5630, 64958, 65192, 65405, 32190, 65151, 34814, 65519, 64478, 65279, 23390, 39791, 49150, 65523, 65534, 63227, 64958, 49151, 65224, 26613, 46079, 53247, 65485, 65466, 65407, 60927, 55295, 61119, 49151, 63951, 45051, 64991, 64878, 32767, 57239, 49131, 61439, 63295, 8191, 40959, 61439, 61307, 65263, 65524, 28670, 47103, 65527, 61423, 59262, 61307 },
    { 1, 30856, 61152, 32190, 32488, 28663, 38206, 61430, 64200, 32190, 33622, 4382, 32488, 31711, 38206, 64478, 32488, 28663, 32488, 31711, 28662, 32766, 31710, 28662, 38206, 61430, 38206, 64478, 31710, 33825, 28662, 32766, 1632, 28150, 64959, 31646, 63967, 32702, 60919, 32702, 64959, 31646, 60855, 28150, 64415, 64959, 61367, 64959, 60919, 64959, 61367, 32254, 14602, 50421, 32190, 59799, 63967, 64959, 64415, 32254, 32190, 33345, 20636, 47883 },
    { 1, 65256, 61160, 32702, 65256, 40959, 61246, 56214, 65128, 65535, 31742, 32694, 65256, 47071, 32734, 64190, 65256, 63459, 32662, 15855, 46906, 65005, 64918, 42934, 55294, 56937, 64470, 32767, 63486, 32379, 64446, 28639, 65256, 63206, 31711, 40862, 65387, 62973, 32766, 61398, 64867, 65263, 64495, 61334, 49017, 65483, 60927, 57278, 63487, 65133, 64190, 31143, 65535, 64377, 59347, 32766, 63231, 32747, 65438, 65471, 65451, 65517, 65430, 28671 },
    { 1, 61160, 65192, 49050, 65256, 65531, 65534, 53182, 47848, 21831, 57278, 65342, 32488, 61435, 63230, 28638, 65256, 49143, 49146, 32765, 65530, 65231, 59390, 55215, 65022, 24559, 56702, 49142, 59390, 65535, 65470, 64511, 65256, 32727, 30719, 65015, 65019, 28639, 65519, 61375, 65367, 48634, 63487, 57343, 56959, 65535, 57135, 65534, 32746, 65518, 61335, 48119, 61311, 28411, 65534, 32732, 59391, 65517, 61339, 62975, 65535, 65273, 57342, 40958 },
    { 1, 65128, 65256, 47094, 64104, 47103, 37886, 65438, 65256, 40943, 65374, 65398, 57064, 65469, 65398, 61150, 61128, 46447, 30454, 61423, 46974, 61417, 65430, 58357, 59390, 65515, 32758, 32511, 63474, 61309, 48894, 31735, 65256, 49151, 64959, 60894, 16247, 65515, 65006, 32755, 56045, 65143, 64383, 60878, 32631, 55295, 61435, 65343, 46975, 65404, 32671, 31734, 65022, 32767, 63486, 65259, 46969, 44799, 65463, 32127, 42735, 65403, 32607, 51919 },
    { 1, 46760, 63720, 57310, 65256, 38903, 47998, 65534, 32488, 57079, 31230, 64958, 32488, 40959, 32254, 31678, 65256, 52989, 32750, 65469, 46718, 16235, 56606, 63949, 65278, 65007, 48062, 32759, 65214, 64895, 63486, 61231, 40168, 40700, 65526, 28607, 65387, 57071, 59901, 65467, 44923, 49150, 60382, 59807, 64255, 48831, 27567, 60373, 46847, 63339, 65535, 64287, 65423, 53231, 65527, 16383, 63319, 61423, 60926, 32735, 65271, 65407, 49151, 28660 },
};

static uint32_t debug_pin = 0;

static constexpr uint8_t SUBKEY_I67[48] = { 14, 8, 48, 48, 20, 6,
    4, 16, 22, 15, 13, 17,
    10, 1, 11, 23, 3, 9,
    0, 5, 48, 48, 2, 7,
    38, 44, 46, 48, 33, 41,
    39, 48, 37, 48, 26, 42,
    43, 28, 27, 48, 47, 32,
    36, 31, 35, 40, 29, 30 };

// the characteristic used in the article

static constexpr uint32_t target_character[8] = {
    /*
    0x405c0000,
    0x04000000,
    0x00540000,
    0,
    0x00540000,
    0x04000000,
    0x405c0000,
    */
    // bswap
    0x00003a02,
    0x00000020,
    0x00002a00,
    0,
    0x00002a00,
    0x00000020,
    0x00003a02,
};

static void pretty_diff64(uint64_t x, uint64_t y)
{
    fprintf(stderr, "%016lx:%016lx:%016lx\n", x, y, x ^ y);
}

static void pretty_diff48(uint64_t x, uint64_t y)
{
    fprintf(stderr, "%012lx:%012lx:%012lx\n", x, y, x ^ y);
}

static void pretty_diff32(uint32_t x, uint32_t y)
{
    fprintf(stderr, "%08x:%08x:%08x\n", x, y, x ^ y);
}
#define ewatchdiff(EXPR, DIGIT)            \
    {                                      \
        eprintf(#EXPR " => ");             \
        pretty_diff##DIGIT(EXPR, EXPR##_); \
    }

static constexpr uint64_t bvx(uint64_t i, uint64_t l, uint64_t w, uint64_t p)
{
    return ((i >> l) & ((1ULL << w) - 1)) << p;
}

static constexpr uint64_t apply_ip(uint64_t input)
{
    uint64_t v1 = 0;
    for (int i = 0; i < 64; i++) {
        v1 |= bvx(input, IP[i], 1, i);
    }
    return v1;
}

static constexpr uint64_t apply_ipi(uint64_t input)
{
    uint64_t v1 = 0;
    for (int i = 0; i < 64; i++) {
        v1 |= bvx(input, i, 1, IP[i]);
    }
    return v1;
}

static constexpr uint64_t input_diff(uint32_t ldiff, uint32_t rdiff)
{
    uint64_t input = bvx(ldiff, 0, 32, 0) | bvx(rdiff, 0, 32, 32);
    return apply_ipi(input);
}

static constexpr bool lead_to(uint8_t sbox_i, uint8_t input_xor, uint8_t output_xor)
{
    return (CROSS_AV[sbox_i][input_xor] >> output_xor) & 1;
}

static constexpr uint64_t input_diff1 = input_diff(target_character[0], target_character[1]);

static constexpr size_t DEFAULT_NPAIRS = 1 << 19; //1 << 17; //32768;
static constexpr size_t NSUGGEST_THRESHOLD = 2048;
static constexpr uint64_t RAND_STRIDE = input_diff1 | 0x1000000000000001;

static constexpr uint64_t plaintext_i(uint64_t i)
{
    // cannot return i directly - if plaintext is small, then they may be all right pairs or wrong pairs, without variation
    return i * RAND_STRIDE;
}

static void gen_pairs(uint64_t* blocks, size_t npairs)
{
    for (size_t i = 0; i < npairs; i++) {
        uint64_t blk = plaintext_i(i);
        blocks[i * 2] = blk;
        blocks[i * 2 + 1] = blk ^ input_diff1;
    }
}

static void gen_subkey(uint64_t key, uint64_t subkeys[ROUNDS])
{
    // TODO: many places bit endians!
    for (int i = 0; i < ROUNDS; i++) {
        uint64_t sk = 0;
        for (int j = 0; j < 48; j++) {
            sk |= bvx(key, KEY_PERM[i][j], 1, j);
        }
        subkeys[i] = sk;
    }
}

static constexpr uint64_t des_f_e(uint32_t input)
{
    uint64_t input_r = bvx(input, 0, 32, 0) | bvx(input, 0, 32, 32);
    uint64_t e = 0;
    for (int i = 0; i < 8; i++) {
        e |= bvx(input_r, i == 0 ? 32 - 1 : i * 4 - 1, 6, i * 6);
    }
    return e;
}

static constexpr uint32_t des_f_p(uint32_t s)
{
    uint64_t p = 0;
    for (int i = 0; i < 32; i++) {
        p |= bvx(s, P[i], 1, i);
    }
    return p;
}

static constexpr uint32_t des_f_pi(uint32_t s)
{
    uint64_t p = 0;
    for (int i = 0; i < 32; i++) {
        p |= bvx(s, i, 1, P[i]);
    }
    return p;
}

static constexpr uint8_t des_sbox(int sbox_i, uint64_t e)
{
    int row = bvx(e, sbox_i * 6, 1, 1) | bvx(e, sbox_i * 6 + 5, 1, 0);
    int col = bvx(e, sbox_i * 6 + 1, 4, 0);
    return SBOX[sbox_i][row][col];
}

static constexpr uint32_t des_f(uint32_t input, uint64_t subkey)
{
    uint64_t e = des_f_e(input);
    e ^= subkey;
    uint64_t s = 0;
    for (int i = 0; i < 8; i++) {
        s |= bvx(des_sbox(i, e), 0, 4, i * 4);
    }
    return des_f_p(s);
}

static uint64_t des_block(uint64_t input, uint64_t subkeys[ROUNDS], int dec = 0)
{
    uint64_t v1 = apply_ip(input);
    for (int h = 0; h < ROUNDS; h++) {
        uint32_t left = bvx(v1, 0, 32, 0), right = bvx(v1, 32, 32, 0);
        left ^= des_f(right, subkeys[dec ? ROUNDS - h - 1 : h]);
        v1 = bvx(left, 0, 32, 32) | bvx(right, 0, 32, 0);
    }
    v1 = bvx(v1, 0, 32, 32) | bvx(v1, 32, 32, 0);
    uint64_t v2 = apply_ipi(v1);
    return v2;
}

static bool right_pair(uint64_t base, uint64_t subkeys[ROUNDS])
{
    uint64_t base_ = input_diff1 ^ base;
    uint64_t v1 = apply_ip(base), v1_ = apply_ip(base_);
    uint32_t left = bvx(v1, 0, 32, 0), left_ = bvx(v1_, 0, 32, 0);
    uint32_t right = bvx(v1, 32, 32, 0), right_ = bvx(v1_, 32, 32, 0);
    assert((left ^ left_) == target_character[0]);
    assert((right ^ right_) == target_character[1]);
    for (int i = 2; i < 7; i++) {
        uint32_t tmp = left ^ des_f(right, subkeys[i - 2]);
        left = right;
        right = tmp;
        uint32_t tmp_ = left_ ^ des_f(right_, subkeys[i - 2]);
        left_ = right_;
        right_ = tmp_;
        if ((right ^ right_) != target_character[i]) {
            return false;
        }
    }
    return true;
}

static bool check_lead_to_packed(uint64_t inp_d, uint32_t f_d, uint8_t sbox_i)
{
    return lead_to(sbox_i, bvx(inp_d, sbox_i * 6, 6, 0), bvx(f_d, sbox_i * 4, 4, 0));
}

static uint64_t find_inv(uint64_t inp, uint64_t inp_d, uint32_t f_d, uint8_t sbox_i)
{
    uint8_t inp_p = bvx(inp, sbox_i * 6, 6, 0);
    uint8_t inp_pd = bvx(inp_d, sbox_i * 6, 6, 0);
    uint8_t f_pd = bvx(f_d, sbox_i * 4, 4, 0);
    uint64_t orig_mask = CROSS_AV_DET[sbox_i][inp_pd][f_pd];
    uint64_t new_mask = 0;
    for (int i = 0; i < 64; i++) {
        if (bvx(orig_mask, i, 1, 0)) {
            new_mask |= 1ull << (i ^ inp_p);
        }
    }
    return new_mask;
}

size_t count_mask5_n(uint8_t mask)
{
    return __builtin_popcount(mask);
}

static void apply_sum(uint32_t* counter, array<uint64_t, 6> term, uint8_t sum_mask)
{
    // find first bit
    assert(sum_mask);
    int first_bit = -1, nbit = 0;
    for (int i = 4; i >= 0; i--) {
        if (bvx(sum_mask, i, 1, 0)) {
            if (first_bit == -1)
                first_bit = i;
            nbit++;
        }
    }
    uint8_t next_mask = sum_mask & ~(1ull << first_bit);
    for (int i = 0; i < 64; i++) {
        if (bvx(term[first_bit], i, 1, 0)) {
            if (nbit == 1) {
                counter[i]++;
            } else {
                apply_sum(counter + (i << (6 * (nbit - 1))), term, next_mask);
            }
        }
    }
}

static vector<uint32_t> count_crit(
    const vector<array<uint64_t, 6>>& terms,
    uint32_t cond,
    uint8_t cond_mask,
    uint8_t sum_mask)
{
    size_t sum_n = count_mask5_n(sum_mask);
    vector<uint32_t> ret(1ull << (sum_n * 6));
    for (size_t i = 0, s = terms.size(); i < s; i++) {
        bool skip = false;
        for (int j = 0; j < 5; j++) {
            if (bvx(cond_mask, j, 1, 0)) {
                if (!bvx(terms[i][j], bvx(cond, j * 6, 6, 0), 1, 0)) {
                    skip = true;
                    break;
                }
            }
        }
        if (skip)
            continue;
        apply_sum(&ret[0], terms[i], sum_mask);
    }
    return ret;
}

static int term_suggest_strength(uint64_t mask)
{
    int ret = __builtin_popcountl(mask);
    //if (ret > 16)
    //    ret = 16;
    //if (ret < 4)
    //    ret = 4;
    return ret;
}

static uint64_t crack_post_process(vector<array<uint64_t, 4>> right_pairs, uint32_t part_key)
{
    // the most stable method in the case where part_key is correct should be to bruteforce
    // first determine wave8 sbox1+sbox4 key
    // G'=f'^h'=f'^Tr', useful G' should be f_pi(G') to invert perm
    // g=H^Tl=f(h,cand_key)^Tl=f(Tr,cand_key)^Tl
    // key7 should be found via certain bits in key8
    // find out whether corresponding key bit can lead to sbox'(g^k) = G'
    // sboxi=2, 3rd box
    // structure:
    // <IPed>Tl,<IPed>Tr,<IPed>Tl*,<IPed>Tr*,f_pi(G')
    bool debug = debug_pin && part_key == debug_pin;
    vector<array<uint64_t, 3>> preproc;
    for (size_t i = 0, s = right_pairs.size(); i < s; i++) {
        const auto& rp = right_pairs[i];
        uint64_t T = apply_ip(rp[2]), T_ = apply_ip(rp[3]);
        uint64_t G_p = target_character[6] ^ bvx(T, 32, 32, 0) ^ bvx(T_, 32, 32, 0);
        uint64_t G_pp = des_f_pi(G_p);
        preproc.push_back({ T, T_, G_pp });
    }
    uint64_t subkey_h_bg = bvx(part_key, 0, 6, 6) | bvx(part_key, 6, 24, 24);
    if (debug) {
        ewatchx(subkey_h_bg);
    }
    //ewatchi(right_pairs.size());
    //ewatchx(subkey_h_bg);
    {
        // sbox 0+3 are found with sbox 2
        uint32_t best_key = 0;
        size_t best_count = 0;
        for (uint32_t i = 0; i < (1 << 12); i++) {
            uint64_t subkey_h = subkey_h_bg | bvx(i, 0, 6, 0) | bvx(i, 6, 6, 18);
            uint64_t subkey_g = 0;
            for (int j = 0; j < 48; j++) {
                subkey_g |= bvx(subkey_h, SUBKEY_I67[j], 1, j);
            }
            size_t count = 0;
            for (size_t j = 0, s = right_pairs.size(); j < s; j++) {
                if (count + (s - j) <= best_count)
                    break;
                uint32_t Tl = bvx(preproc[j][0], 0, 32, 0), Tr = bvx(preproc[j][0], 32, 32, 0);
                uint32_t Tl_ = bvx(preproc[j][1], 0, 32, 0), Tr_ = bvx(preproc[j][1], 32, 32, 0);
                uint32_t g = Tl ^ des_f(Tr, subkey_h);
                uint32_t g_ = Tl_ ^ des_f(Tr_, subkey_h);
                uint32_t G_p = preproc[j][2];
                //if()
                uint32_t so_g = des_sbox(2, des_f_e(g) ^ subkey_g);
                uint32_t so_g_ = des_sbox(2, des_f_e(g_) ^ subkey_g);
                if (bvx(G_p, 2 * 4, 4, 0) == (so_g ^ so_g_)) {
                    count += 1;
                }
            }
            if (count > best_count) {
                best_count = count;
                best_key = i;
            }
        }
        subkey_h_bg |= bvx(best_key, 0, 6, 0) | bvx(best_key, 6, 6, 18);

        if (debug) {
            ewatchx(best_key);
            ewatchi(best_count);
        }
    }
    if (debug) {
        ewatchx(subkey_h_bg);
    }
    {
        // sbox 2 is found with sbox 1
        uint32_t best_key = 0;
        size_t best_count = 0;
        for (uint32_t i = 0; i < (1 << 6); i++) {
            uint64_t subkey_h = subkey_h_bg | bvx(i, 0, 6, 12);
            uint64_t subkey_g = 0;
            for (int j = 0; j < 48; j++) {
                subkey_g |= bvx(subkey_h, SUBKEY_I67[j], 1, j);
            }
            size_t count = 0;
            for (size_t j = 0, s = right_pairs.size(); j < s; j++) {
                if (count + (s - j) <= best_count)
                    break;
                uint32_t Tl = bvx(preproc[j][0], 0, 32, 0), Tr = bvx(preproc[j][0], 32, 32, 0);
                uint32_t Tl_ = bvx(preproc[j][1], 0, 32, 0), Tr_ = bvx(preproc[j][1], 32, 32, 0);
                uint32_t g = Tl ^ des_f(Tr, subkey_h);
                uint32_t g_ = Tl_ ^ des_f(Tr_, subkey_h);
                uint32_t G_p = preproc[j][2];
                //if()
                uint32_t so_g = des_sbox(1, des_f_e(g) ^ subkey_g);
                uint32_t so_g_ = des_sbox(1, des_f_e(g_) ^ subkey_g);
                if (bvx(G_p, 1 * 4, 4, 0) == (so_g ^ so_g_)) {
                    count += 1;
                }
            }
            if (count > best_count) {
                best_count = count;
                best_key = i;
            }
        }
        subkey_h_bg |= bvx(best_key, 0, 6, 12);
        if (debug) {
            ewatchx(best_key);
            ewatchi(best_count);
        }
    }
    if (debug) {
        ewatchx(subkey_h_bg);
    }
    uint64_t des_key = 0;
    for (int i = 0; i < 48; i++) {
        des_key |= bvx(subkey_h_bg, i, 1, KEY_PERM[7][i]);
    }
    uint64_t subkey_g = 0;
    for (int i = 0; i < 48; i++) {
        subkey_g |= bvx(subkey_h_bg, SUBKEY_I67[i], 1, i);
    }
    {
        // remaining bit that the last wave does not cover
        for (int h = 0; h < 8; h++) {
            uint8_t perm[6] = {}, bit_n = 0;
            for (int j = 0; j < 6; j++) {
                if (SUBKEY_I67[h * 6 + j] == 48) {
                    perm[bit_n++] = j;
                }
            }
            if (bit_n == 0)
                continue;
            uint32_t best_key = 0;
            size_t best_count = 0;
            for (uint32_t i = 0; i < (1u << bit_n); i++) {
                uint64_t subkey_g_test = subkey_g;
                for (int j = 0; j < bit_n; j++) {
                    subkey_g_test |= bvx(i, j, 1, 6 * h + perm[j]);
                }
                size_t count = 0;
                for (size_t j = 0, s = right_pairs.size(); j < s; j++) {
                    if (count + (s - j) <= best_count)
                        break;
                    uint32_t Tl = bvx(preproc[j][0], 0, 32, 0), Tr = bvx(preproc[j][0], 32, 32, 0);
                    uint32_t Tl_ = bvx(preproc[j][1], 0, 32, 0), Tr_ = bvx(preproc[j][1], 32, 32, 0);
                    uint32_t g = Tl ^ des_f(Tr, subkey_h_bg);
                    uint32_t g_ = Tl_ ^ des_f(Tr_, subkey_h_bg);
                    uint32_t G_p = preproc[j][2];
                    //if()
                    uint32_t so_g = des_sbox(h, des_f_e(g) ^ subkey_g_test);
                    uint32_t so_g_ = des_sbox(h, des_f_e(g_) ^ subkey_g_test);
                    if (bvx(G_p, h * 4, 4, 0) == (so_g ^ so_g_)) {
                        count += 1;
                    }
                }
                if (count > best_count) {
                    best_count = count;
                    best_key = i;
                }
            }
            if (debug) {
                ewatchx(h);
                ewatchx(best_key);
                ewatchi(best_count);
            }
            for (int i = 0; i < bit_n; i++) {
                des_key |= bvx(best_key, i, 1, KEY_PERM[6][6 * h + perm[i]]);
            }
        }
    }
    return des_key;
}

static uint64_t try_crack_with_subkey(const vector<array<uint64_t, 6>>& cand_support, const uint64_t* encblocks, uint32_t cand_subkey)
{
    vector<array<uint64_t, 4>> right_pairs;

    for (size_t i = 0, s = cand_support.size(); i < s; i++) {
        size_t nsuggest = 1;
        for (int j = 0; j < 5; j++) {
            nsuggest *= term_suggest_strength(cand_support[i][j]);
        }
        bool good = true;
        for (int j = 0; j < 5; j++) {
            if (!bvx(cand_support[i][j], bvx(cand_subkey, j * 6, 6, 0), 1, 0)) {
                good = false;
                break;
            }
        }
        if (good) {
            //is_right_pair[i] = 1;
            int idx = cand_support[i][5];
            right_pairs.push_back({
                plaintext_i(idx),
                plaintext_i(idx) ^ input_diff1,
                encblocks[idx * 2],
                encblocks[idx * 2 + 1],
            });
        }
    }

    if (debug_pin && cand_subkey == debug_pin) {
        ewatchi(right_pairs.size());
    }

    if (right_pairs.size() < 1) {
        return INVALID_KEY;
    }

    uint64_t des_key = crack_post_process(right_pairs, cand_subkey);

    if (debug_pin && cand_subkey == debug_pin) {
        ewatchx(des_key);
    }
    {
        uint64_t subkeys[ROUNDS];
        gen_subkey(des_key, subkeys);
        for (size_t i = 0, s = cand_support.size(); i < s && i < 16; i++) {
            int idx = cand_support[i][5];
            if (des_block(plaintext_i(idx), subkeys) != encblocks[idx * 2] || des_block(plaintext_i(idx) ^ input_diff1, subkeys) != encblocks[idx * 2 + 1]) {
                return INVALID_KEY;
            }
        }
    }
    return des_key;
}

static uint64_t max_path(const vector<array<uint64_t, 6>>& terms, const uint64_t* encblocks)
{
    // filter
    vector<array<uint64_t, 6>> filtered_terms;
    {
        vector<size_t> count;
        for (size_t i = 0, s = terms.size(); i < s; i++) {
            size_t nsuggest = 1;
            for (int j = 0; j < 5; j++) {
                nsuggest *= term_suggest_strength(terms[i][j]);
            }
            count.push_back(nsuggest);
        }
        vector<size_t> count_sorted = count;
        std::sort(count_sorted.begin(), count_sorted.end());
        size_t threshold = count_sorted[(count_sorted.size() - 1) * 3 / 4]; //* 15 / 16];
        ewatchi(threshold);
        for (size_t i = 0, s = terms.size(); i < s; i++) {
            if (count[i] >= threshold) {
                filtered_terms.push_back(terms[i]);
            }
        }
        ewatchi(terms.size());
        ewatchi(filtered_terms.size());
    }
    vector<uint32_t> root_counter = count_crit(filtered_terms, 0, 0, 0b11100);
    vector<pair<uint32_t, uint32_t>> root_index(root_counter.size());
    for (size_t i = 0, s = root_index.size(); i < s; i++) {
        root_index[i] = { root_counter[i], i };
    }
    std::sort(root_index.begin(), root_index.end());
    std::reverse(root_index.begin(), root_index.end());
    size_t max_occ = 0;
    uint32_t max_cand = 0;
    for (size_t i = 0, s = root_index.size() // very unstable heuristic that seems to work with high prob - d&c
         ;
         i < s; i++) {
        if (i % 4096 == 0 || i < 256) {
            eprintf("i = %ld / %ld, max_occ = %ld, cur_occ = %u, max_cand = %08x\n", i, s, max_occ, root_index[i].first, max_cand);
        }
        // unstable heuristic

        //if (root_index[i].first <= 2 * max_occ) {
        //    break;
        //}

        // this is stable
        if (root_index[i].first <= max_occ) {
            break;
        }
        uint32_t known = root_index[i].second << 12;
        vector<uint32_t> branch_counter = count_crit(filtered_terms, known, 0b11100, 0b00011);
        for (size_t j = 0, t = branch_counter.size(); j < t; j++) {
            if (branch_counter[j] > max_occ) {
                max_occ = branch_counter[j];
                max_cand = j | known;
            }
        }
        for (size_t j = 0, t = branch_counter.size(); j < t; j++) {
            if (branch_counter[j] >= max_occ || branch_counter[j] > 5) { //(branch_counter[j] * 2 > max_occ) {
                // try out
                uint64_t cand_key = try_crack_with_subkey(terms, encblocks, j | known);
                if (cand_key != INVALID_KEY) {
                    ewatchi(i);
                    ewatchx((j | known));
                    return cand_key;
                }
            }
        }
    }
    return INVALID_KEY;
}

static uint64_t do_crack(const uint64_t* encblocks, size_t npairs)
{
    // plaintext is assumed to be generated by gen_pairs
    // 2678 are empty
    //static uint8_t counter[1 << 24];
    //memset(counter, 0, sizeof(counter));
    //unordered_map<uint32_t, size_t> counter;
    vector<array<uint64_t, 6>> cand_support;
    cand_support.reserve(npairs);
    size_t right_pair_predict = 0;
    for (size_t i = 0; i < npairs; i++) {
        uint64_t enc1 = apply_ip(encblocks[i * 2]), enc2 = apply_ip(encblocks[i * 2 + 1]);
        uint64_t enc_d = enc1 ^ enc2;
        uint64_t inp_h_d = des_f_e(bvx(enc_d, 32, 32, 0));
        uint64_t inp_h = des_f_e(bvx(enc1, 32, 32, 0));
        uint32_t fh_d = des_f_pi(bvx(enc_d, 0, 32, 0) ^ target_character[5]);
        // 1,4,5,6,7
        if (
            !check_lead_to_packed(inp_h_d, fh_d, 1) || !check_lead_to_packed(inp_h_d, fh_d, 4) || !check_lead_to_packed(inp_h_d, fh_d, 5) || !check_lead_to_packed(inp_h_d, fh_d, 6) || !check_lead_to_packed(inp_h_d, fh_d, 7)) {
            //wrong pair
            continue;
        }
        array<uint64_t, 6> cand = {
            find_inv(inp_h, inp_h_d, fh_d, 1),
            find_inv(inp_h, inp_h_d, fh_d, 4),
            find_inv(inp_h, inp_h_d, fh_d, 5),
            find_inv(inp_h, inp_h_d, fh_d, 6),
            find_inv(inp_h, inp_h_d, fh_d, 7),
            i
        };
        if (debug_pin) {
            bool good = true;
            for (int j = 0; j < 5; j++) {
                if (!bvx(cand[j], bvx(debug_pin, j * 6, 6, 0), 1, 0)) {
                    good = false;
                    break;
                }
            }
            if (good)
                right_pair_predict++;
        }
        cand_support.push_back(cand);
        for (int j = 0; j < 5; j++) {
            assert(cand[j]);
        }
    }
    if (debug_pin) {
        ewatchi(right_pair_predict);
    }
    //uint32_t cand_subkey = max_path(cand_support);
    return max_path(cand_support, encblocks);
}

extern "C" {
void des_encrypt(const uint64_t* input, uint64_t* output, const uint64_t* key, size_t nblock);
void des_decrypt(const uint64_t* input, uint64_t* output, const uint64_t* key, size_t nblock);
void crack(const uint64_t* encblocks, size_t npairs, uint64_t* key);
void pair_stream(uint64_t* blocks, size_t npairs);
}

void des_encrypt(const uint64_t* input, uint64_t* output, const uint64_t* key, size_t nblock)
{
    uint64_t subkeys[ROUNDS];
    gen_subkey(*key, subkeys);
    for (size_t i = 0; i != nblock; i++) {
        output[i] = des_block(input[i], subkeys, 0);
    }
}
void des_decrypt(const uint64_t* input, uint64_t* output, const uint64_t* key, size_t nblock)
{
    uint64_t subkeys[ROUNDS];
    gen_subkey(*key, subkeys);
    for (size_t i = 0; i != nblock; i++) {
        output[i] = des_block(input[i], subkeys, 1);
    }
}

void crack(const uint64_t* encblocks, size_t npairs, uint64_t* key)
{
    *key = do_crack(encblocks, npairs);
}

void pair_stream(uint64_t* blocks, size_t npairs)
{
    return gen_pairs(blocks, npairs);
}

static void pretty_des_f(uint32_t input, uint32_t input_, uint64_t subkey, uint32_t& fv, uint32_t& fv_)
{
    uint64_t e = des_f_e(input);
    uint64_t e_ = des_f_e(input_);
    eprintf("e return: ");
    pretty_diff48(e, e_);
    e ^= subkey;
    e_ ^= subkey;
    eprintf("xord key: ");
    pretty_diff48(e, e_);
    uint64_t s = 0;
    uint64_t s_ = 0;
    for (int i = 0; i < 8; i++) {
        s |= bvx(des_sbox(i, e), 0, 4, i * 4);
        s_ |= bvx(des_sbox(i, e_), 0, 4, i * 4);
    }
    eprintf("sbox    : ");
    pretty_diff32(s, s_);
    fv = des_f_p(s);
    fv_ = des_f_p(s_);
}

static uint64_t pretty_des_block(uint64_t input, uint64_t input_, uint64_t subkeys[ROUNDS])
{
    eprintf("Input: ");
    pretty_diff64(input, input_);

    uint64_t v1 = apply_ip(input);
    uint64_t v1_ = apply_ip(input_);
    eprintf("   ip: ");
    pretty_diff64(v1, v1_);
    for (int h = 0; h < ROUNDS; h++) {
        eprintf("Wave %d - key %012lx\n", h, subkeys[h]);
        uint32_t left = bvx(v1, 0, 32, 0), right = bvx(v1, 32, 32, 0);
        uint32_t left_ = bvx(v1_, 0, 32, 0), right_ = bvx(v1_, 32, 32, 0);
        eprintf("left    : ");
        pretty_diff32(left, left_);
        eprintf("right   : ");
        pretty_diff32(right, right_);
        uint32_t fv, fv_;
        pretty_des_f(right, right_, subkeys[h], fv, fv_);
        eprintf("f return: ");
        pretty_diff32(fv, fv_);
        left ^= fv;
        left_ ^= fv_;
        eprintf("new left: ");
        pretty_diff32(left, left_);
        v1 = bvx(left, 0, 32, 32) | bvx(right, 0, 32, 0);
        v1_ = bvx(left_, 0, 32, 32) | bvx(right_, 0, 32, 0);
        eprintf("   block: ");
        pretty_diff64(v1, v1_);
    }
    v1 = bvx(v1, 0, 32, 32) | bvx(v1, 32, 32, 0);
    v1_ = bvx(v1_, 0, 32, 32) | bvx(v1_, 32, 32, 0);
    eprintf("end - shift back: ");
    pretty_diff64(v1, v1_);
    uint64_t v2 = apply_ipi(v1);
    uint64_t v2_ = apply_ipi(v1_);
    eprintf("ipi             : ");
    pretty_diff64(v2, v2_);
    return v2 ^ v2_;
}

int main()
{
    uint64_t key = 0x123456789abcdef0;
    uint64_t subkeys[ROUNDS];
    uint64_t* pairs = new uint64_t[DEFAULT_NPAIRS * 2];
    uint64_t* pairs_res = new uint64_t[DEFAULT_NPAIRS * 2];
    for (int t = 0; t < 1024; t++) {
        ewatchi(t);
        key = rnd() & 0xfefefefefefefefeUL;
        ewatchx(key);
        gen_subkey(key, subkeys);
        uint32_t cand = bvx(subkeys[7], 6, 6, 0) | bvx(subkeys[7], 24, 24, 6);
        debug_pin = cand;
        ewatchx(cand);
        ewatchx(subkeys[7]);
        gen_pairs(pairs, DEFAULT_NPAIRS);
        eprintf("encrypt\n");
        des_encrypt(pairs, pairs_res, &key, DEFAULT_NPAIRS * 2);
        eprintf("crack\n");
        uint64_t cracked = do_crack(pairs_res, DEFAULT_NPAIRS);
        if (cracked != key) {
            for (size_t i = 0; i < DEFAULT_NPAIRS; i++) {
                //ewatchi(i)
                uint64_t cand = plaintext_i(i);
                if (right_pair(cand, subkeys)) {
                    printf("Right pair: %ld %lx\n", i, cand);
                    pretty_des_block(cand, cand ^ input_diff1, subkeys);
                    break;
                }
            }
            eprintf("WRONG cracked: key = %016lx, cracked key = %016lx : diff %016lx %016lx\n", key, cracked, key ^ cracked, cracked & ~key);
            char buf[64];
            fgets(buf, 64, stdin);
        }
    }
    delete[] pairs;
    delete[] pairs_res;
    return 0;
}