如何并行斐波那契数列直到 10^5 项

Question

我需要从这个顺序代码中并行（使用 openmp）斐波那契数列来计算该数列的第 10⁵ 项，但我已经被卡住了 3 周没有任何好主意，有人有任何想法或好方法吗？

C中的顺序代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAX 100010
#define LEN 25001

char seq[MAX][LEN];

void add(int a, int b) {
    int i, aux, s;

    for (i = 0, aux = 0; seq[a][i] != '[=10=]' && seq[b][i] != '[=10=]'; i++) {
        s = seq[a][i] + seq[b][i] + aux - '0' - '0';
        aux = s / 10;
        seq[a + 1][i] = s % 10 + '0';
    }

    while (seq[a][i] != '[=10=]') {
        s = seq[a][i] + aux - '0';
        aux = s / 10;
        seq[a + 1][i] = s % 10 + '0';
        i++;
    }

    while (seq[b][i] != '[=10=]') {
        s = seq[b][i] + aux - '0';
        aux = s / 10;
        seq[a + 1][i] = s % 10 + '0';
        i++;
    }

    if (aux != 0)
        seq[a + 1][i++] = aux + '0';

    seq[a + 1][i] = '[=10=]';
}

int main() {
    int n, i, len;

    seq[0][0] = '0';
    seq[0][1] = '[=10=]';
    seq[1][0] = '1';
    seq[1][1] = '[=10=]';

    for (i = 2; i < MAX; i++)
        add(i - 1, i - 2);

    scanf("%d", &n);

    len = strlen(seq[n]);
    for (i = 0; i <= len - 1; i++)
        printf("%c", seq[n][len - 1 - i]);
    printf("\n");
    fflush(stdout);

    return 0;
}

Answer 1

如您所知，F(100000) 是一个天文数字。要计算该值，您必须将另外两个非常大的数字 F(99999) 和 F(99998) 相加。

这是我的提示：

你有两个数字，有很多（千）位长和 N 个处理器。您可以将加法拆分到多个线程中。例如：

F(167) == 35600075545958458963222876581316753
F(168) == 57602132235424755886206198685365216

要计算 F(169)，您必须将上面的两个数字相加。但让我们将其视为 4 个单独的加法，每个加法为 9 位数字。

          (A)         (B)         (C)        (D)
F(167) == 035600075   545958458   963222876  581316753
F(168) == 057602132   235424755   886206198  685365216
======================================================= + 
           93202207   781383213  1849429074 1266681969
                                 ^          ^

所以现在我们得到了 4 个总和。其中两个 C 和 D 具有进位操作。所以我们只需要将每个左边的结果调整+1

           93202207   781383213   849429074  266681969
                             +1          +1
======================================================= + 
F(169) ==  93202207   781383214   849429075 266681969

因此，您计算 Fib(100000) 的算法类似于以下内容。其中 BigNumber 是一个表示您的数字的结构。您使用的是字符数组，这也是可以接受的。

void fib(int stop)
{
    BigNumber  f0 = 0;
    BigNumber  f1 = 1;
    BigNumber  f2 = 1;

    BigNumber* pF0 = &f0;
    BigNumber* pF1 = &f1;
    BigNumber* pF2 = &f2;

    for (int i = 2; i <= stop; i++)
    {
        ParallelAdd(pF0, pF1, pF2); // *pF2 = *pF1 + *pF0

        // shift via pointers
        pF0 = pF1;
        pF1 = pF2;
        pF2 = pF0;
    }
    Print(f2);
}

main()
{
    fib(100000);
}

您的 ParallelAdd 会将通过指针 f1 和 f0 传递的数字分成 N 组，每组 K 位，其中 N 是您可用或想要使用的处理器数量。然后每个处理器计算每个的加法 - 使用您已有的代码。在这 N 个操作完成后，扫描结果集以查看哪些添加包含 K+1 位长的结果，然后按上述 +1 逻辑进行调整。合并成一个字符串并推回 pF2.

引用的指针地址

Answer 2

与其尝试并行化 bignum 加法，这很棘手，您可以尝试并行计算多个项：

F(n+1) = F(n) + F(n-1)
F(n+2) = F(n+1) + F(n) = 2*F(n) + F(n-1)
F(n+3) = F(n+2) + F(n+1) = 2*F(n+1) + F(n) = 2*F(n) + 2*F(n-1) + F(n) = 3*F(n) + 2*F(n-1)
...

另请注意，您应该一次计算数字块：可以使用 32 位数组元素计算 8 或 9 个以 10 为基数的数字。

这是一个有多项改进的修改版本：

它一次计算 8 位数字的块
它可以接受命令行参数
它使用更少的内存
它可以处理更大的值
效率更高 (20x)

您应该能够轻松地将其并行化。

/* Parallelisable bignum Fibonacci computation by chqrlie */
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#include <stdlib.h>
#include <string.h>

#if 0

/* 2500ms for fib(100000) */
#define DIGIT  10
#define NDIGIT 1
#define FMT "d"
typedef unsigned char digit_t;

#elif 1

/* 279ms for fib(100000) */
#define DIGIT  100000000
#define NDIGIT 8
#define FMT PRIu32
typedef uint32_t digit_t;

#else

/* 720ms for fib(100000) */
#define DIGIT  1000000000000000000
#define NDIGIT 18
#define FMT PRIu64
typedef uint64_t digit_t;

#endif

int add1(digit_t *c, const digit_t *a, int alen, const digit_t *b, int blen) {
    digit_t aux, s;
    int i;

    /* assuming alen >= blen */
    for (i = 0, aux = 0; i < alen; i++) {
        s = a[i] + b[i] + aux;
        aux = s / DIGIT;
        c[i] = s % DIGIT;
    }
    if (aux != 0) {
        c[i++] = (digit_t)aux;
    }
    return i;
}

int add2(digit_t *c, const digit_t *a, int alen, const digit_t *b, int blen) {
    digit_t aux, s;
    int i;

    /* assuming alen >= blen */
    for (i = 0, aux = 0; i < alen; i++) {
        s = 2 * a[i] + b[i] + aux;
        aux = s / DIGIT;
        c[i] = s % DIGIT;
    }
    if (aux != 0) {
        c[i++] = (digit_t)aux;
    }
    return i;
}

int add3(digit_t *c, const digit_t *a, int alen, const digit_t *b, int blen) {
    digit_t aux, s;
    int i;

    /* assuming alen >= blen */
    for (i = 0, aux = 0; i < alen; i++) {
        s = 3 * a[i] + 2 * b[i] + aux;
        aux = s / DIGIT;
        c[i] = s % DIGIT;
    }
    if (aux != 0) {
        c[i++] = (digit_t)aux;
    }
    return i;
}

int add4(digit_t *c, const digit_t *a, int alen, const digit_t *b, int blen) {
    digit_t aux, s;
    int i;

    /* assuming alen >= blen */
    for (i = 0, aux = 0; i < alen; i++) {
        s = 5 * a[i] + 3 * b[i] + aux;
        aux = s / DIGIT;
        c[i] = s % DIGIT;
    }
    if (aux != 0) {
        c[i++] = (digit_t)aux;
    }
    return i;
}

void printseq(const digit_t *s, int len) {
    printf("%"FMT, s[len - 1]);
    for (int i = 1; i < len; i++)
        printf("%.*"FMT, NDIGIT, s[len - 1 - i]);
    printf("\n");
}

int main(int argc, char *argv[]) {
    int MIN, i, LEN, MAX;

    if (argc > 1) {
        MAX = MIN = strtol(argv[1], NULL, 0);
        if (argc > 2)
            MAX = strtol(argv[2], NULL, 0);
    } else {
        scanf("%d", &MIN);
        MAX = MIN;
    }

    /* length if fib(n) is less than n*log10(phi)+2 */
    LEN = (MAX * 20910ULL) / 100000 / NDIGIT + 2;
    /* allocate 8 bignums */
    int *slen = calloc(sizeof(*slen), 8);
    digit_t (*seq)[LEN] = calloc(sizeof(*seq), 8);

    if (slen == NULL || seq == NULL) {
        fprintf(stderr, "memory allocation error\n");
        return 1;
    }

    seq[0][0] = 0;
    slen[0] = 1;
    if (0 >= MIN) printseq(seq[0], slen[0]);
    seq[1][0] = 1;
    slen[1] = 1;
    if (1 >= MIN) printseq(seq[1], slen[1]);

    for (i = 2; i <= MAX && (MAX + 1 - i) % 4 != 0; i++) {
        slen[i] = add1(seq[i], seq[i - 1], slen[i - 1], seq[i - 2], slen[i - 2]);
        if (i >= MIN) printseq(seq[i], slen[i]);
    }
    for (; i <= MAX; i += 4) {
        int im2 = (i - 2) & 7;
        int im1 = (i - 1) & 7;
        int i0 = (i + 0) & 7;
        int i1 = (i + 1) & 7;
        int i2 = (i + 2) & 7;
        int i3 = (i + 3) & 7;
        /* the next 4 calls can be parallelised */
        slen[i0] = add1(seq[i0], seq[im1], slen[im1], seq[im2], slen[im2]);
        slen[i1] = add2(seq[i1], seq[im1], slen[im1], seq[im2], slen[im2]);
        slen[i2] = add3(seq[i2], seq[im1], slen[im1], seq[im2], slen[im2]);
        slen[i3] = add4(seq[i3], seq[im1], slen[im1], seq[im2], slen[im2]);
        /* the print calls must be called sequentially */
        if (i + 0 >= MIN) printseq(seq[i0], slen[i0]);
        if (i + 1 >= MIN) printseq(seq[i1], slen[i1]);
        if (i + 2 >= MIN) printseq(seq[i2], slen[i2]);
        if (i + 3 >= MIN) printseq(seq[i3], slen[i3]);
    }
    free(slen);
    free(seq);
    return 0;
}

如何并行斐波那契数列直到 10^5 项

How to parallel Fibonacci sequence until 10^5 term

c

parallel-processing

multithreading

openmp

fibonacci