Fortran Function explanation

Question

I have this function in Fortran and i'm trying to recode it in C#

C **************************************************************** 
C    FUNCTION   POLY 
C***************************************************************** 
      FUNCTION POLY(N,A,X) 
      DIMENSION A(N) 
C
      POLY    = 0. 
      L       = N 
      DO 1 K  = 1,N 
      POLY    = POLY*X + A(L) 
1     L       = L-1 
      RETURN 
      END 
C*****************************************************************

I found out that DIMENSION A(N) creates a vector of N values, but we already have a variable A from the function parameters, does this means that the array values all equal to A? If so then what is the use of A(N). By the way can anyone just explain what does this function do so i can re-implement it in C#

score 5 · Accepted Answer · edited Apr 16 '12 at 14:15

5

      FUNCTION POLY(N,A,X)      ! implicitly real (float) function poly(int n,int a,real x)
      DIMENSION A(N)            ! shape A as 1d array of n values in this scope
C                               ! say nothing (blank comment)
      POLY    = 0.              ! initialise return variable to float value 0
      L       = N               ! set L (implicitly integer) to n
      DO 1 K  = 1,N             ! for(int k=1; k<=n; ++k)
      POLY    = POLY*X + A(L)   !    update return variable
1     L       = L-1             !    decrement L
      RETURN                    ! return current value for poly
      END

so in c-like syntax:

float poly(int n, int a, float x) {
    // redim a(n)
    float result = 0;
    int l = n;
    for(int k=1; k <= n; ++k) {
        result = result*x + a(l);
        --l;
    }
    return result;
}

The bit that doesn't translate is redimensioning A as an array. In C you would pass a pointer and use it as an array, and in C++/C# you'd probably pass a vector-like structure with its own length property.

In C#, using a list:

float poly(List<float> coeffs, float x) {
    float result = 0;
    for(int i=coeffs.Count-1; i >= 0; --i) {
        result = result*x + coeff[i];
    }
    return result;
}

edited Apr 16 '12 at 14:15

David Heffernan

601,492
42
1,072
1,490

answered Apr 16 '12 at 12:54

Phil H

19,928
7
68
105

1

The commenting of the FORTRAN code is nice, but the C# and Python code you present are incorrect. – David Heffernan Apr 16 '12 at 13:46
@DavidHeffernan: ...in what way? – Phil H Apr 16 '12 at 13:55
In the way that they give the wrong answers. In the C# code you mean `=` rather than `+=`. And you iterate over the coefficients in the wrong order. You want to `coeff[n-1]`, `coeff[n-2]` etc. And the Python code bears no relation to polynomial evalutation. You can fix it by replacing the `+` with a `*` but then it's not Horner's method and it inefficient. It's pretty enough, it just gives the wrong answers. The only code in your answer that is correct is the C pseudo code. I'm sure you can correct it easily enough, but as it stands this answer is down-vote bait. – David Heffernan Apr 16 '12 at 13:57
@DavidHeffernan: Thanks, I spotted the order mistake after your first comment, but not the +=. Removed the python code because the point was more a sensible way to achieve the result in python, rather than replicate the algorithm, since python has its overheads. – Phil H Apr 16 '12 at 14:06
OK, it was getting hard to do this in comments. I corrected your code properly now. – David Heffernan Apr 16 '12 at 14:16
shouldn't that be `int *a` and `a[l]`? Or perhaps `float *a` (as it would be implicitly that in the fortran code)? – laxxy Apr 16 '12 at 22:15
Thank you guys, you really helped me a lot, I'm now moving toward [A C# Class for Complex Polynomials - CodeProject®](http://www.codeproject.com/Articles/19446/A-C-Class-for-Complex-Polynomials) because I found out that I'll be using polynomials in other parts of the code so for the sake of consistency I'll re-implement the entire data processing section – workoverflow Apr 17 '12 at 02:49

David Heffernan · Answer 2 · 2012-04-16T12:37:41.013

It evaluates a polynomial in x of the form:

a[1] + a[2]x + a[3]x^2 + ... a[N]x^(N-1)

Remember that Fortran uses 1-based array indices and I have followed that convention in this equation.

You could write it in C# like this:

double EvaluatePolynomial(double[] a, double x)
{
    double result = 0.0;
    int i = a.Length;
    while (i>0)
    {
        i--;
        result = result*x + a[i];
    }
    return result;
}

Here we are using 0-based array indexing appropriate to C#. And so this function evalutates a polynomial in x of the form:

a[0] + a[1]x + a[2]x^2 + ... a[N-1]x^(N-1)

score 2 · Answer 3 · answered Apr 16 '12 at 12:38

guessing slightly I think this is specifying that the parameter A itself is an array of N elements.

So for the C# equivalent you wouldn't need a separate N parameter; you would only need to pass A as a double[], as in .NET arrays can tell you their .Length.

The function evaluates polynomials using Horner's method.

score 2 · Answer 4 · answered Apr 16 '12 at 12:51

2

The line DIMENSION A(N) just declares details of the A dummy argument (PARAMETERs are something very different in Fortran), i.e. it says it is an array from 1 to N. The other ones are not declared this way, because the function uses implicit typing.

answered Apr 16 '12 at 12:51

Vladimir F Героям слава

57,977
4
76
119

1

Yes. It is not creating the variable A, it is declaring that dummy argument to describe the actual argument that is being passed to this function. If the calling routine hasn't created the storage there is a problem. – M. S. B. Apr 16 '12 at 15:44
It's pretty much equivalent to K&R prototypes in C – David Heffernan Apr 16 '12 at 15:55

Fortran Function explanation

4 Answers4

Linked