Static Array Initialization

The Fortran DATA Statement

In C, if we initialize a static array, the initialization occurs before the program begins running.

static double   list[] = { 2, 4, 10 };
            

Without the static modifier, the array would be initialized at run time, when the block containing the variable definition is entered. With the static modifier, the initialization takes place at compile time, so the numbers are already in the array as soon as the process begins executing.

Another way to look at it is that variables are always initialized when they are instantiated, i.e. when their memory is allocated. Static variables are instantiated at compiled time and auto variables at run time, when they are needed.

Another way to initialize a Fortran array is using a data statement.

integer :: list(3)

data list / 2, 4, 10 /
            

This differs from assigning the array an array constant in that the data statement is processed at compile-time, whereas the assignment is done at run-time. When the program begins execution, the values from a data statement will already be in the array. Since the data statement is not performed at run-time, it actually makes no difference where it is placed within the subprogram. It is processed and removed by the compiler, not translated to machine code.

You may need to assign values to an array at run-time, or even assign the array several different sets of values at different times. In this case, an assignment statement is necessary.

However, if the contents of an array will not change during program execution, then using a data statement will save execution time, since the data was loaded into the array at the same time the program itself was loaded into memory.

Look-up Tables

A look-up table is a list of precomputed values. If a value is expensive to compute (i.e. requires a loop) and frequently needed, it may be quicker to store precomputed values in an array and look them up instead. This only involves indexing the array (computing the address of an element) which is trivial (a few machine instructions).

Ideal candidates for look-up tables are mathematical functions with a small domain that require a loop to compute. Since we are storing precomputed results for all cases, a look-up table will generally require more memory than the code to compute results on-the-fly.

The factorial is one such function. Only a small number of factorials are within the range of typical integer or floating point types, and they must be computed iteratively.

To achieve a performance gain using a look-up table, we must use a Fortran data statement or a C static array. Assigning an array constant at run-time or initializing an automatic array is costly, since it is itself an implied loop.

uint64_t   fastfact(unsigned int n)

{
    static uint64_t table[] = { 1ul, 1ul, 2ul, 6ul, 24ul, 120ul, 720ul,
        5040ul, 40320ul, 362880ul, 3628800ul, 39916800ul, 479001600ul,
        6227020800ul, 87178291200ul, 1307674368000ul, 20922789888000ul,
        355687428096000ul, 6402373705728000ul, 121645100408832000ul };
    
    return table[n];
}
function fastfact(n)

    ! Disable implicit declarations (i-n rule)
    implicit none
    
    ! Function type
    integer(8) :: fastfact
    
    ! Argument variables
    integer, intent(in) :: n
    
    ! Local variables
    integer(8) :: table(0:19)
    
    fastfact = table(n)
    
    data table / 1_8, 1_8, 2_8, 6_8, 24_8, 120_8, 720_8, 5040_8, 40320_8, &
        362880_8, 3628800_8, 39916800_8, 479001600_8, 6227020800_8, &
        87178291200_8, 1307674368000_8, 20922789888000_8, &
        355687428096000_8, 6402373705728000_8, 121645100408832000_8 /
end function

To generate the factorial data for the lookup table requires a program. Integers are limited to about 19!. Double precision floating point can represent up to 170!, but there's another problem: Even 64-bit floating point systems are limited to about 16 significant figures. As soon as the factorial value reaches 17 significant figures, it gets rounded. All later factorials are computed from this rounded value, and so the round-off error grows rapidly from this point on. 170! has 305 significant digits, but at most 16 of them will be accurate if the value is computed using double precision. Output from a C program using double precision follows. Note that 171! is reported as infinity, since 64-bit IEEE format cannot represent it.

    169! =
    4.26906800900470559560e+304
    
    170! =
    7.25741561530800402460e+306
    
    171! =
    inf
    

What we need to compute large factorials precisely is an arbitrary-precision integer package. Such a package works by concatenating integer values end-to-end, and processing these large integers in chunks (e.g. 64 bits at a time on a 64 bit computer). There are libraries available for C, C++, Java, etc. However, the Unix bc is quite convenient for this purpose, and includes a C-like scripting language. Below is a bc script that prints factorials. You can run it by saving it to a file (e.g. fact.bc) and running bc fact.bc at the Unix prompt.

f = 1
for (c = 1; c <= 170; ++c) {
    f
    f = f * c
}
quit

    169! =
    42690680090047052749392518888995665380688186360567361038491634111179\
    77554942180092854323970142716152653818230139905012221568248567907501\
    77960523574894559464847084134121076211998036035274015123788150487897\
    50405684196703601544535852628274771797464002689372589486243840000000\
    000000000000000000000000000000000
    
    170! =
    72574156153079989673967282111292631147169916812964513765435777989005\
    61843401706157852350749242617459511490991237838520776666022565442753\
    02532890077320751090240043028005829560396661259965825710439855829425\
    75689663134396122625710949468067112055688804571933402126614528000000\
    00000000000000000000000000000000000
    

A modified version of the bc script can generate output that can be copied and pasted directly into a program. Note that this version requires GNU bc, which is not standard on all systems.

fact = 1
for (c = 1; c <= 170; ++c) {
    print fact, ",\n";
    fact = fact * c;
}
quit