基于OpenMP实现利用泰勒级数计算正弦函数sin(x)的并行程序

LT_real

已于 2024-03-23 20:30:59 修改

阅读量336

点赞数 6

文章标签：算法 c语言

于 2024-03-23 20:27:55 首次发布

本文链接：https://blog.csdn.net/qq_47814951/article/details/136974283

版权

一、正弦函数的泰勒级数

$sin(x) = x - \frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\cdots=\sum_{i=0}^{\infty}(-1)^i\frac{x^{2i+1}}{(2i+1)!}$

基于此泰勒级数，完成正弦函数sin(x)计算程序，并在多核系统中，比较不同线程个数与串行程序的运行时间。

二、程序设计

1.串行程序

并行程序的基础是串行程序，我们先实现一个简单的利用泰勒级数计算正弦函数的程序：

#include <stdio.h>
#include <math.h>
#define pi 3.1415926

int main(){

    int i=1,count=0;
    double x=pi/3;
    double sum,term;
    
    // printf("Input x: ");
    // scanf("%f",&x);
    sum=x;   //赋初始值
    term=x;

    for(i;i<10;i+=2){
        term =-term*x*x/((i+1)*(i+2));//计算相应项
        // printf("%f\n",term);
        sum+=term;
        count++;
    }
    printf("sin(%f)=%f,count=%d\n",x,sum,count);
    return 0;
}

程序结果为：

基于此程序，我们才可以开展并行程序设计。

2.并行程序

观察泰勒级数，发现每一次计算第i个数时，我们都是对第i-1个数直接进行操作以简化阶乘、次方的运算，即后续计算结果强依赖于前面的计算结果，这与阶乘的并行化有异曲同工之妙。

那么现在我们只需固定计算总次数（for循环次数）为MAXC，按线程数量NUM_THREADS对MAXC进行切分，让每一个线程都对不同的i从1开始计算，最终我们会得到线程数量个数的计算偏移量（offset1、offset2、...），将这些偏移量按顺序相乘，便可以得到该线程对应第i块的级数：

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <time.h>
#include <math.h>
#define pi 3.1415926
#define MAXC 1000000000
#define NUM_THREADS 16

int main(){
    double x;
    x=pi/3;
    // printf("Input x: ");
    // scanf("%lf",&x);

    while(x>pi){
        x-=pi;
    }
    while(x<-pi){
        x+=pi;
    }
    int i=1;
    double sum=0,term;
    double mid[NUM_THREADS+1];//offset for every threads
    mid[0]=x;
    
    // Start measuring time
    clock_t start = clock();
    #pragma omp parallel num_threads(NUM_THREADS) shared(mid,sum) private(term)
    {
        int itd = omp_get_thread_num();
        int ntd = omp_get_num_threads();
        int jud=1;
        double midsum=0;
        double offset=1.0;
        term=1;

        #pragma omp for schedule(static) private(i)
        for(i=1;i<MAXC;i++){
            jud=(i%2==0)?1:-1;
            term=(term<0)?-term:term;
            term =(double)jud*term*x*x/(double)((2*i)*(2*i+1));
            midsum+=term;
        }
        mid[itd+1]=term;

        #pragma omp barrier
        for(i=0;i<(itd+1);i++){
            offset*=mid[i];
        }
        midsum*=offset;

        #pragma omp critial
        sum+=midsum;
    }
    // Stop measuring time and calculate the elapsed time
    clock_t end = clock();
    double elapsed = (double)(end - start)/CLOCKS_PER_SEC;
    printf("sin(%f)=%f\n",x,sum+x);
    printf("用时: %.3f 秒\n", elapsed);

    return 0;
}