c - 如何并行斐波那契数列直到 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] != '\0' && seq[b][i] != '\0'; 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] != '\0') {
        s = seq[a][i] + aux - '0';
        aux = s / 10;
        seq[a + 1][i] = s % 10 + '0';
        i++;
    }

    while (seq[b][i] != '\0') {
        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] = '\0';
}

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

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

    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

与其尝试并行化 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)
...

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

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

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

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

/* 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;
}