c - 如何用C语言编写求和函数的代码？

Question

编辑:我添加了 main、factorial 和 trapGamma 函数来提供完整的图片，但我特别讨论的是 I 函数中 iSum 的 for 循环。

基本上，我已经没有想法了，并且在我所知道的所有地方都竭尽全力寻找这个问题的答案。我需要编写一个程序来计算代表 M/M/1 队列的复杂函数。

该函数包含计算 Gamma 函数积分、计算阶乘等子函数。我已经编写了所有计算代码，但我的总和给了我巨大的数字，而我预计不会高于 0.35

#include <math.h>
#include <stdio.h>

double I(int k, double t);
double trapGamma(double z);
unsigned long long int factorial(unsigned int n);

int main()
{
    int k;
    int i = 0;
    double dt = 0.1;

    printf("Ikx = [ \n");

        for (t = 14.0 ; t <= 15.0; t += dt)
        {
            printf("%f " , t);

            for (k = 1 ; k <= 10 ; k++)
            {
                I(k, t);
                printf("%f " , I(k, t));
            }
            printf("\n");

        }

        printf(" ];\n");

    return (0);
}


double I(int k, double t)
{
    unsigned long long int x;
    unsigned int n = 20;
    double numerator, y, pow1, c;
    double iSum;
    double Ix;
    int i = 0;

    iSum = 0.0;
    Ix = 0.0;

    a = .25 * pow(t , 2);
    b = pow(a, i);
    x = factorial(n);
    y = trapGamma(k + i + 1);

    iSum = (b / (x * y));

    //This is the sum loop that I'm having trouble with, I've broke the iSum equation down for my own readability while coding right above this comment

    for (i = 0; i <= 100 ; i++)
        {
        iSum += i;
        }

    Ix = (pow((.5 * t), k) ) * iSum;
    return Ix;
}


/*
I've checked both the factorial and trapGamma functions and they are giving me the expected results. 
*/

unsigned long long int factorial(unsigned int n)
{
    if(n <= 1)
        return 1;
    else
        return (n * factorial(n - 1));
}


double trapGamma (double z)
{

    int i , N = 100;
    double gamma;
    double a = 0.0;
    double b = 15.0;
    double x1, x2, y1, y2;
    double areai;
    double w = (b - a) / N; 
    gamma = 0.0;

    for (i = 1; i < N; i++)
        {
            x1 = a + ((i - 1) * w); //the left bound point
            x2 = a + (i*w); //the right bound point
            y1 = pow(x1,z - 1)*exp(-x1); //the height of our left bound
            y2 = pow(x2, z - 1)*exp(-x2); //the height of our right bound
            areai = ((y1 + y2) / 2.0) * (x2 - x1);
            gamma += areai;
        }
    return gamma;
}

这是建立在另一个项目的基础上的，我在该项目中使用贝塞尔函数在 60 秒的时间内创建了 M/M/1 队列，这样我就可以看到这个队列应该是什么。我已经检查了 trapGamma 和阶乘函数的结果，它们都按预期工作。

求和应该如何编码？

Answer 1

如果发布代码的目的是计算 modified Bessel function I ，有一些陷阱和有用的简化需要注意。鉴于

尝试分别计算阶乘、Gamma 函数的值、它们的乘积以及总和的每一项的幂迟早会导致整数溢出。

最好更新总和的每个加数的值。

此外，鉴于 k 是一个整体，我们有 Γ(n) = (n - 1)!

由于 double 类型的精度有限，加数越来越小，经过一些迭代后，加数太小而无法添加到总和中。

// Evaluates x^k / k! trying not to overflow
double power_over_factorial(double x, int k)
{
    double result = 1.0;
    for ( int i = 1; i <= k; ++i )
    {
        result *= x / i;
    }
    return result;
}

#define MAX_ITERS 20

double modified_Bessel_I(int k, double x)
{
    x /= 2;
    const double xx = x * x;
    double partial = power_over_factorial(x, k);
    double old_sum, sum = partial;
    int m = 1;
    do
    {
        old_sum = sum;
        partial *= xx / ((m + k) * m);
        sum += partial;   
    }
    while ( old_sum != sum  &&  ++m < MAX_ITERS );

    return sum;
}

可测试here .