Integrate the product of four infinite series functions

Question

I wrote the following code to integrate the product of four functions.

Two of the functions are hypergeometric functions and the other two are incomplete gamma functions. I have written each function as an infinite series, one per nested loop. The first loop (variable k) is for time.

y1 = 0;
for k = 1:L-1
    for i = 0:N
        total1 = i*log(x) + 2*(n-m-i+k)*log(sig);
        for j = 0:N
            total2 = j*log(y) + 2*(n-m-j+k)*log(sig);
            for p = 0:100
                total3 = 0;
                for q = 0:100
                    temp = (n-m+i+k+p)*log(a);
                    temp = temp - 2*(2*n - 2*m + i + j + 2*k + p + q)*log(sig);
                    temp = temp - func1(k, j+1, q+1) + func2(i+1, j+1, p+1, q+1, k);
                    temp = temp + (i + j + 2*k + p + q)*log(sig^2/(1+a));
                    temp = temp + log(gammainc(b*(1 + x)/sig^2,i + j + p + q));
                    total3 = total3 + exp(temp);
                end
                total4 = func3(k, i+1) + (n - m + i + k + p)*log(a/sig^2);
                total4 = total4 + func5(k, i+1, p+1);
                total4 = total4 + (n - m + i + k + p + 2)*log(sig^2/a);
                total4 = total4 + log(gammainc(a*b/sig^2, i + k + p)));
                total4 = exp(total4) - total3;
                y1 = y1 + exp(total1 + total2 - func5(k, i+1, p+1) + log(total4));
            end
        end
    end
end

Thereafter, I removed the innermost loop the following way to improve speed:

y1 = 0;
for k = 1:L-1
    for i = 0:N
        total1 = i*log(x) + 2*(n-m-i+k)*log(sig);
        for j = 0:N
            total2 = j*log(y) + 2*(n-m-j+k)*log(sig);
            for p = 0:100
                q = 0:100;
                temp = (n-m+i+k+p)*log(a);
                temp = temp - 2*(2*n - 2*m + i + j + 2*k + p + q)*log(sig);
                temp = temp - reshape(func1(k, j+1, q+1),1,101)
                temp = temp + reshape(func2(i+1, j+1, p+1, q+1, k), 1, 101) ;
                temp = temp + (i + j + 2*k + p + q)*log(sig^2/(1+a));
                temp = temp + log(gammainc(b*(1 + x)/sig^2,i + j + p + q));
                total3 = total3 + exp(temp);
                total4 = func3(k, i+1) + (n - m + i + k + p)*log(a/sig^2);
                total4 = total4 + reshape(func5(k, i+1, p+1), 1, 101);
                total4 = total4 + (n - m + i + k + p + 2)*log(sig^2/a);
                total4 = total4 + log(gammainc(a*b/sig^2,i + k + p)));
                total4 = exp(total4) - total3;
                y1 = y1 + exp(total1 + total2 - reshape(func5(k, i+1, p+1),1,101) + log(total4));
            end
        end
    end
end

Can you suggest further changes to improve speed?

Basically this is the function I am trying to integrate:

\begin{eqnarray} F &=& \sum\limits_{k=1}^{L-1}\sum\limits_{i=1}^{\infty}\sum\limits_{j=1}^{\infty}\sum\limits_{p=1}^{\infty}\left(\frac{x^i\sigma^{2\left(n-m-i+k\right)}}{\left(n-m+1\right)_ii!}\frac{y^j\sigma^{2\left(n-m-j+k\right)}}{\left(n-m+1\right)_jj!}\frac{1}{s\left(s+1\right)\cdots\left(s+p\right)}\right)\nonumber\\ &&\left[\left(\Gamma\left(n-m+j+k\right)\left(\frac{a}{\sigma^2}\right)^{n-m+i+k+p}\int_0^{b*}e^{-\frac{ab}{\sigma^2}}b^{n-m+i+k+p+1}db\right)\right]\nonumber\\ &-&\left(\sum\limits_{q=0}^{\infty}\frac{a^{n-m+i+k+p}}{\sigma^{2\left(2n-2m+i+j+2k+p+q\right)}}\frac{1}{s(s+1)\cdots(s+q)}\int_0^{b*}e^{-\frac{ab+b}{\sigma^2}b^{2n-2m+i+j+p+q+2k+2}}db\right) \end{eqnarray}

Answer 1

Caveat: I am not a Matlab expert, nor am I competent at an advanced undergraduate level in mathematics.

Efficiency

Removing data structure impedance and redundant or unnecessary calculations is likely to improve the running time of the code under review.

1. func5 [whatever the func that is] is called twice with exactly the same arguments.

2. Each time it is wrapped in reshape which adds to the cost.

3. Indeed func1, func2, and func5 are defined each returning a data structure that requires reshaping. This can be avoided by either working with their data structure directly in the code under review or redefining the functions to return the data structure used by the code.

4. (n - m + i + k + p)*log(a/sig^2) is calculated twice. Sure in one location it needs to be doubled.

5. Each time through the inner loop (n - m + i + k + p) increases by one. There's no need to recalculate it in full. It can be passed as a parameter from the enclosing function.

6. a [defined elsewhere?] is a constant within the code under review. Therefore log(a/sig^2) is also a constant. Yet is being calculated twice in the inner loop.

7. (log sig) [sig defined elsewhere?] is calculated for each iteration at the each level of the loop. It's a constant.

8. There are probably other places where redundant computations can be avoided.

Code Clarity

The code is hard to read.

The more that the names of items in the code reflect the way we would label the pieces of the problem domain [mathematics here], the easier it becomes to reason about the way the code implements the problem.

sigma_k better expresses the mathematics than total1. And though I know temp = temp + ... does something, but it's not obvious what. function1_result = start_value + ... is better. partialgamma1_result = start_value + ... is better.

Final thoughts

Giving things meaningful names helps a lot. logsig = log(sig) seems like more work because it is more typing [by 8 more keystrokes]. But it's a false economy because debugging and refactoring code for speed is orders of magnitude more difficult than banging on the keyboard.