c - 查找具有一定步数/计数的最小起始值的数字

Closed. This question needs to be more focused. It is not currently accepting answers. Learn more。
想改进这个问题吗？更新问题，使其只关注一个问题editing this post。
去年关门了。
我正在研究一个Collatz猜想程序。我有一个遵循简单规则的序列：
如果当前学期是偶数：下一学期是当前学期的一半。
如果当前项是奇数：下一个项是当前项的三倍，加1。
我们继续这个，直到我们达到一个和计数每个使用的术语。
例如：对于数字17，我们有：
17 52 26 13 40 20 10 5 16 8 4 1
因此，计数是13。
我特别需要帮助的是，我需要找到计数超过1234的最小起始值（数字）。
即：
1有计数1
2有计数2
3有计数8
4有计数3
5有6
六有九
......
18有21
......
97有119
等
这是我认为有效的代码，但需要非常长的时间来处理。我怎样才能优化这个以便它能更快地找到计数？
有人建议我进行二进制搜索，但计数不是线性的，所以不起作用。。

#include <stdio.h>

int main (void) {


    unsigned long long int num1 = 1;
    while (num1 <= 99999999999999999) {

    unsigned long long int num = num1;
    int count = 1;
    while (num != 1) {

        if (num % 2 == 0) {
            num = (num/2);
            count++;
        }

        else {
            num = ((num*3) + 1);
            count++;
        }

    }

    if (count > 1234) {
        printf("%llu\n", num1);
        break;
    }

    num1++;
    }

    return 0;
}

最佳答案

产生大于1234步的Collatz序列的最小根是133561134663。它达到最大值319497287463520，这是一个49位的值。因此，许多C库可以使用uint64_t中提供的<stdint.h>来评估这个序列（通常通过<inttypes.h>）。
不幸的是，有许多从低位开始的序列，需要多达71位整数才能正确解析，比如从110243094271开始的573步序列。
最烦人的是，如果您使用64位无符号整数实现Collatz序列，而不检查溢出，那么您将为这些序列计算错误的长度；但是，由于它们恰好具有无趣的长度，因此该错误被屏蔽，您仍然可以找到上述解决方案！
本质上，如果这是一个C练习，那就是错误：即使是一个可怕的错误和有限的实现也能找到正确的答案。
我碰巧在64位Linux上使用GCC-5.4.0（在笔记本电脑上的i5-7200U核心CPU上运行），所以我可以使用128位无符号整数类型unsigned __int128。我还有16gib的RAM，其中我使用了12gib来存储序列长度高达130亿（13×109）的奇数索引。
当你发现序列的长度从某个奇数i开始时，从2i开始的序列的长度就变长了一步（但在其他方面是一样的）。事实上，有一个从2ki开始的序列，正好比k步长。如果只为某个最小长度的序列寻找最小的起始值，则可以忽略所有偶数个起始点，直到k足够小（基本上是将需要查找的起始值范围括起来）。
考虑到即使是起始值也会有显著的加速（20%或更多），但我在这里画了一条“线”，而是检查每个起始值。
为了加快计算速度，我使用了GCC内置的（64位的__builtin_ctzll()）来查找连续的最低有效零位的数目。这样我可以一次性处理所有连续的“偶数”步骤，但是正确地计算步骤数。
我还将序列步进器分成两个并行的步进器，一个处理状态适合64位变量的情况，另一个处理需要完全128位的情况。
在64位部分，如果当前状态的序列长度是可缓存的，我会查找它。如果它不是零，我把它加到我们已经完成的步骤数上，就得到了结果。否则，我会将该缓存条目索引和根目录中的步骤数添加到一个小的更新缓存中。这样，当我得到结果时，我不需要重复序列，只需在一个紧密的循环中应用更新。在乘三加一之前，我会检查这个值是否溢出（6148914691236517205溢出，到264）；如果溢出，我会切换到128位部分，然后在那里进行乘法相加。
在128位部分，我假设我们没有那么多内存（在exabyte范围内），所以我根本不需要担心缓存部分。在我乘以3再加1（我使用unsigned long long）之前，我会检查状态是否溢出，只是确定一下。
在配备英特尔酷睿i5-7200U的笔记本电脑上，使用单核，计算从1到133561134663的所有Collatz序列的长度需要大约4253秒的CPU时间（大约一小时十分钟）。
下一个更长的序列从158294678119开始，是1243步长；为了达到这个目的，我的笔记本电脑消耗了一个半小时的CPU时间。
然而，为了达到这一点，代码已经变得不像OP，而且相当可怕。特别是，为了混合64位/128位循环，我不得不求助于curr128 += curr128 + curr128 + 1；整个实现中到处都是GCC内置代码。我认为它几乎只写代码；也就是说，如果我怀疑它的任何部分，我会从头重写它。当优化优先于可维护性和健壮性时，这是典型的最终结果。
无论如何，为了让其他在x86-64上的Linux上使用GCC的用户验证上述内容，下面是可怕的代码collatz.c：

#define  _POSIX_C_SOURCE  200809L
#define  _GNU_SOURCE
#include <stdlib.h>
#include <unistd.h>
#include <inttypes.h>
#include <string.h>
#include <stdarg.h>
#include <stdio.h>
#include <errno.h>
#include <time.h>

#define  MAX_SEQ_LEN  2000

#define  OVER64       UINT64_C(6148914691236517205)

typedef  unsigned __int128   u128;
typedef  uint64_t            u64;

static inline u64  u128_lo(const u128  u) { return (u64)u;         }
static inline u64  u128_hi(const u128  u) { return (u64)(u >> 64); }

static inline u64 highest_bit_set_64(const u64  u)
{
    if (sizeof u == sizeof (unsigned int))
        return 63 - __builtin_clz(u);
    else
    if (sizeof u == sizeof (unsigned long))
        return 63 - __builtin_clzl(u);
    else
    if (sizeof u == sizeof (unsigned long long))
        return 63 - __builtin_clzll(u);
    else
        exit(EXIT_FAILURE);
}

static unsigned int highest_bit_set_128(u128 u)
{
    u64  hi = u128_hi(u);
    if (hi)
        return 64 + highest_bit_set_64(hi);
    else
        return highest_bit_set_64(u128_lo(u));
}

static inline unsigned int  ctz64(const u64  u)
{
    if (sizeof u <= sizeof (unsigned int))
        return __builtin_ctz(u);
    else
    if (sizeof u <= sizeof (unsigned long))
        return __builtin_ctzl(u);
    else
    if (sizeof u <= sizeof (unsigned long long))
        return __builtin_ctzll(u);
    else
        exit(EXIT_FAILURE);
}

static inline unsigned int  ctz128(const u128  u)
{
    if (sizeof u == sizeof (unsigned long long))
        return __builtin_ctzll(u);
    else {
        const u64  lo = u128_lo(u);
        if (lo)
            return ctz64(u);
        else
            return 64 + ctz64(u128_hi(u));
    }
}

static const char *s128(u128 u)
{
    static char  buffer[40];
    char        *p = buffer + sizeof buffer;

    *(--p) = '\0';
    do {
        *(--p) = '0' + (u % 10);
        u /= 10;
    } while (u);

    return p;
}

static struct timespec  wall_started, wall_now;
static inline void      wall_start(void)
{
    clock_gettime(CLOCK_MONOTONIC, &wall_started);
}
static inline double    wall_seconds(void)
{
    clock_gettime(CLOCK_MONOTONIC, &wall_now);
    return (double)(wall_now.tv_sec - wall_started.tv_sec)
         + (double)(wall_now.tv_nsec - wall_started.tv_nsec) / 1000000000.0;
}

static struct timespec  cpu_elapsed;
static inline double    cpu_seconds(void)
{
    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &cpu_elapsed);
    return (double)cpu_elapsed.tv_sec + (double)cpu_elapsed.tv_nsec / 1000000000.0;
}

static int out_and_err = 0;
static void print(const char *fmt, ...)
{
    va_list  args;

    va_start(args, fmt);
    vdprintf(STDOUT_FILENO, fmt, args);
    va_end(args);

    if (out_and_err) {
        va_start(args, fmt);
        vdprintf(STDERR_FILENO, fmt, args);
        va_end(args);
    }
}

static size_t        cache_index[MAX_SEQ_LEN];
static unsigned int  cache_depth[MAX_SEQ_LEN];

static u64           seq_max = 0;  /* 2*seq_num */
static size_t        seq_num = 0;
static uint16_t     *seq_len = NULL;

static unsigned int  tip_bits = 0;
static u128          tip_max  = 0;

static inline unsigned int  collatz_length(const u64 root, u128 *maxto)
{
    u128          curr128, max128;
    u64           curr64, max64, lo;
    size_t        cached = 0, i;
    unsigned int  steps = 1, n;

    curr128 = max128 = root;
    curr64 = max64 = root;

any64bit:

    if (!(curr64 & 1)) {
        n = ctz64(curr64);
        curr64 >>= n;
        steps += n;
    }

    if (curr64 >= OVER64) {
        curr128 = curr64;
        goto odd128bit;
    }

odd64bit:

    if (curr64 <= 1)
        goto done;

    if (curr64 < seq_max) {
        i = curr64 >> 1;
        if (seq_len[i]) {
            steps += seq_len[i];
            goto done;
        }

        cache_index[cached] = i;
        cache_depth[cached] = steps;
        cached++;
    }

    curr64 += curr64 + curr64 + 1;
    steps++;

    if (max64 < curr64)
        max64 = curr64;

    goto any64bit;


any128bit:

    if (!(curr128 & 1)) {
        n = ctz128(curr128);
        curr128 >>= n;
        steps += n;

        if (!u128_hi(curr128)) {
            lo = u128_lo(curr128);
            if (lo <= 1)
                goto done;
            if (lo < OVER64) {
                curr64 = lo;
                goto odd64bit;
            }
        }
    }

odd128bit:

    if (u128_hi(curr128) >= OVER64) {
        print("Overflow at root %" PRIu64 ".\n", root);
        exit(EXIT_FAILURE);
    }

    curr128 += curr128 + curr128 + 1;
    steps++;

    if (max128 < curr128)
        max128 = curr128;

    goto any128bit;

done:

    if (cached >= MAX_SEQ_LEN) {
        print("Update cache overrun.\n");
        exit(EXIT_FAILURE);
    }

    while (cached-->0)
        seq_len[ cache_index[cached] ] = steps - cache_depth[cached];

    if (max128 < (u128)max64)
        max128 = max64;

    if (maxto)
        *maxto = max128;

    if (tip_max <= max128) {
        const unsigned int  maxbits = highest_bit_set_128(max128) + 1;
        tip_max = max128;
        if (tip_bits <= maxbits) {
            tip_bits = maxbits;
            print("%" PRIu64 " length %u (reaches %s - %u bits).\n", root, steps, s128(max128), maxbits);
        }
    }

    return steps;
}

int main(void)
{
    unsigned int  n, nmax = 0;
    u128          m;
    uint64_t      i = 1;

    wall_start();

    /* If standard output is redirected to a file, print everything to standard error also. */
    out_and_err = (isatty(STDERR_FILENO) && !isatty(STDOUT_FILENO));

    /* Try allocating up to 16 GiB of cache. */
    seq_num = (size_t)1024 * 1024 * 1024 * 16 / sizeof seq_len[0];
    while (1) {
        seq_len = malloc(seq_num * sizeof seq_len[0]);
        if (seq_len)
            break;
        seq_num = ( seq_num * 7 ) / 8;
    }
    seq_max = 2 * (uint64_t)seq_num;
    memset(seq_len,~0, seq_num * sizeof seq_len[0]);
    memset(seq_len, 0, seq_num * sizeof seq_len[0]);

    print("Allocated %zu entries (%.3f GiB)\n", seq_num, (double)(seq_num * sizeof seq_len[0]) / 1073741824.0);

    do {
        n = collatz_length(i, &m);
        if (n >= nmax) {
            const double cs = cpu_seconds();
            const double ws = wall_seconds();
            const char  *s = s128(m);
            nmax = n;
            print("%" PRIu64 " length %u (reaches %s) [%.3f seconds elapsed, %.3f seconds CPU time]\n", i, n, s, ws, cs);
        }

        i++;
    } while (nmax < MAX_SEQ_LEN);

    return EXIT_SUCCESS;
}

我使用goto编译并运行它。（例如，优化可以确保在编译时解析gcc -Wall -O2 collatz.c -lrt -o collatz && ./collatz > results.txtif子句，并使用条件移动来生成最少数量的慢速条件跳转。因此，上面的程序设计为至少使用sizeof进行编译。）它包含额外的代码，例如显示所用的实时和CPU时间。因为它们与手头的问题无关，所以它们确实可以帮助助教轻松检测是否有人试图将此代码作为自己的作业提交。
当它工作时，我可以在另一个终端中使用-O2来查看到目前为止获得的结果，按增加的序列长度排序；或者使用grep -gk 3 results.txt来查看根据序列中的峰值（序列中的最大值）排序的结果。