Assignment 2.1 - PPDI 水仙花数

功能

求出所有k位水仙花数,按照从小到大的顺序进行排序。

支持范围:$1\leq k \leq 38$

输入 输出

输入:$k$(水仙花数位数)

输出:所有$k$位的水仙花数,以空格分割,从小到大输出。

实现

通过dfs进行递归搜索,递归所有“数码从左到右依次从小到大排列”的数,

计算这些数的水仙花和(即$F_{10}(n)$),判断水仙花和是否为水仙花数

若是水仙花数,则存储至数列中,最终经过排序后输出

优化

  • uint128类型计算,最大支持到38位水仙花数计算,在3.0GHz下38位水仙花数需要计算45分钟。(39位可能需要高精度)

  • 组合数学优化,对于由同样数字组成的数,只遍历一次(即其数字从小到大排列的情况)。

  • 快速幂优化,使用快速幂计算数字的高次幂,并将计算结果存储在数列中,后次计算时直接读取

  • 归并排序,将最终求得的水仙花数数列进行归并排序后输出

生成

文件结构如下:

$./ppdi
├── CMakeLists.txt
├── README.md
├── src
│   ├── main_big.c      # 支持至k=38
└── └── main_small.c    # 支持至k=19

可以通过如下方式生成并运行:

gcc ./src/main_big.c -o ppdi
./ppdi

亦可通过CMake进行生成:

mkdir build
cd build
cmake ..
make
./ppdi

代码

#include <stdio.h>

// PluPerfect Digital Invariant (PPDI)

__uint128_t powlist[10];

__uint128_t ppdi_list[1000] = {0};
__uint128_t ppdi_sorttemp[1000] = {0};
__uint128_t ppdi_counter = 0;
const __uint128_t max64 = 18446744073709551615U;
const __uint128_t splitter = 10000000000000000000U;

// 快速幂 - 求x的k次方

__uint128_t qpow(int x, int k)
{
    // 0次方
    if (k == 0)
        return 1;
    // 奇次方
    if (k & 1)
        return (x * qpow(x, k - 1));
    // 偶次方
    __uint128_t t = qpow(x, k / 2);
    return (t * t);
}

// 求一位数的x次方列表 存储在array中

void gen_powlist(int k)
{
    for (int i = 0; i < 10; i++)
        powlist[i] = qpow(i, k);
}

// 判断n(在10进制下有k位)是否为水仙花数

char isPPDI(__uint128_t n, int k)
{
    if (n < qpow(10, k - 1) || n >= qpow(10, k))
        return 0;
    __uint128_t bak = n; // 预存n值
    __uint128_t F = 0;   // PPDI函数Sum值
    while (n != 0)
    {
        F += powlist[n % 10];
        n /= 10;
    }
    return F == bak;
}

// 求n(k位数)的F值(水仙花和)

__uint128_t sumPPDI(__uint128_t n, int k)
{
    __uint128_t bak = n; // 预存n值
    __uint128_t F = 0;   // PPDI函数Sum值
    while (n != 0)
    {
        F += powlist[n % 10];
        n /= 10;
    }
    return F;
}

// 获得长为length的num的第digit位(从左起为第1位)
// "12345" - len=5
// get_digit(12345, 5, 4) -> 4
int get_digit(__uint128_t num, int length, int digit)
{
    if (digit < 1)
        return 0;
    return (num / qpow(10, length - digit)) % 10;
}

// DFS 根据数字组合搜索水仙花数 并存入ppdi_list

int searchPPDI(__uint128_t num, int length, int index) //index:12345 6(search)
{
    if (index != length + 1)
        for (int i = get_digit(num, length, index - 1); i < 10; i++)
        {
            searchPPDI(num + i * qpow(10, length - index), length, index + 1);
        }
    else
    {
        if (isPPDI(sumPPDI(num, length), length))
        {
            ppdi_list[ppdi_counter] = sumPPDI(num, length);
            ppdi_counter++;
        }
    }
    return 0;
}

// 归并 将array[l, r-1]进行排序

void mergesort(__uint128_t *array, __uint128_t *temp, int l, int r)
{
    if (r - l <= 1)
        return; // 单元素排序
    int mid = (l + r) >> 1;

    mergesort(array, temp, l, mid);
    mergesort(array, temp, mid, r);

    // 此时[l, mid)与[mid, r)都已经排序完毕,接下来需要将两个子数列进行合并
    // 双指针扫一遍 每次取较小元素放入temp 随后将temp怼入array

    int i = l, j = mid, head = l; // i: 左半部分指针,j: 右半部分指针
    while (i < mid && j < r)      // 双指针扫完半边之前 取较小的元素放入临时数组中
        temp[head++] = array[i] < array[j] ? array[i++] : array[j++];
    while (i < mid) // 将左半剩余元素怼到temp里
        temp[head++] = array[i++];
    while (j < r) // 将右半剩余元素怼到temp里
        temp[head++] = array[j++];

    for (i = l; i < r; i++) // 将temp整个怼到array里
        array[i] = temp[i];

    return;
}

// 输入:PPDI数的位数 - 0 < k <= 19
// 输出:k位的所有PPDI,按从小到大排序

int main()
{
    int k;

    printf("Input k:\n> ");
    scanf("%d", &k);

    // 输入数据处理 不支持的直接扔掉结束
    if (k <= 0 || k >= 20) // 用了uint128_t 如果再要更高的位数的话需要高精度 不过再高的话太慢了也没什么意义
    {
        printf("Using 128bit calculation, maybe slow\n");
        //return -1;
    }

    // 0是水仙花数(特例)
    if (k == 1)
        printf("0 ");

    // 存储1~9的k次方 优化速度
    gen_powlist(k);

    // 搜索所有的水仙花数 直接输出
    searchPPDI(0, k, 1);

    // 对搜到的PPDI进行排序
    mergesort(ppdi_list, ppdi_sorttemp, 0, ppdi_counter);

    // 输出
    if (!ppdi_counter)
        printf("Not Found");
    else
        for (int i = 0; i < ppdi_counter; i++)
        {
            if (ppdi_list[i] <= max64)
                printf("%llu ", (unsigned long long)ppdi_list[i]);
            else
            {
                __uint128_t ppdi1 = (ppdi_list[i] / splitter);
                __uint128_t ppdi2 = (ppdi_list[i] - (__uint128_t)(ppdi_list[i] / splitter) * splitter);
                printf("%llu%llu ",
                       (unsigned long long)ppdi1,
                       (unsigned long long)ppdi2);
            }
        }
    printf("\n");

    return 0;
}
最后修改:2021 年 10 月 22 日 07 : 10 PM
如果觉得这篇文章对你有用,请随意赞赏~