on Project Euler & BigInts

About a year ago, I discovered Project Euler and started blowing through problems. The first two dozen or so are actually pretty simple and shouldn't take much more than 50 lines of code to do (most require less). I even wrote a MODULE that contained a fair number of repeated functions (like a prime finder and a digit sum) to make life more simple.

But some of the posted solutions for Problem 16 bother me. The question asks,

2¹⁵ = 32768 and the sum of its digits is 3 + 2 + 7 + 6 + 8 = 26. What is the sum of the digits of the number 2¹⁰⁰⁰?

The digit sum is a fairly simple thing to do:

digit = 0
do while(n /= 0)
    digit = digit + mod(n,10)
    n = n/10
enddo

But the way I solved this, this part isn't even necessary (though it's still useful to know).
The problem you'll run into is that 2¹⁰⁰⁰ contains 302 digits and even using quadruple precision (via the ISO_FORTRAN_ENV module and assuming you have the library for it) the maximum value is a mere 19 digit number: 9223372036854775807. Obviously way too small.
So some people, namely C/C++ and Java programmers, utilized the BigInt package that extends integers to arbitrary positions.

I personally think this is cheating. I think this problem is to be solved without resorting to external libraries that extend precision to arbitrary limits. And this problem can be done in any language without doing so too.

All you have to remember is this: arrays. If you think of any n-digit number as an n-celled array, then a 302-digit number is a 302-celled array.

2 → (2)
22 → (2,2)
222 → (2,2,2)


and so on. If you add 22+222, you get 244:

222+22 → (2,2,2) + (2,2) = (2,4,4)


Adding 66, you have to carry your 10 over:

244+66 → (2,4,4) + (6,6) = (2,10,10) = (3,0,10) = (3,1,0)


With multiplication, you need to multiply each cell by the value & carry the tens multiple times:

222*22 → (2,2,2) * 22 = (2*22,2*22,2*22) = (44,44,44) = (4,8,8,4)


If you understand this, then 2¹⁰⁰⁰ is painlessly easy.

I clocked this method with the Fortran intrinsic subroutine cpu_time, which is accurate to the microsecond, and kept getting zero's across the board (i.e., 0.000000E+00).

On Shaping Arrays

The best thing about Fortran is its handling of arrays, in particular the multi-dimensional ones. As someone who models astrophysical fluid flows for a living, I tend use 4-dimensional arrays (3 for space and 1 for variables) in most of my relevant codes.

In defining arrays, there are three methods by which we can use to define them:
   (1) nested do loops
   (2) forall assignment
   (3) nested implied-do loops

I presume every Fortran user knows of the first of these methods very well, as it is the most used and most versatile method of defining an array. I think most people have heard of the forall assignment, which was once part of the now-defunct High Performance Fortran but added to Fortran 95 update. It is used as follows,

forall(i=1:iMax, j=1:jMax, k=1:kMax) A(i,j,k)=i*k+j/jMax

I can pass the i,j,k indices to a function that returns a value, so it can work much like the nested do loop. The original point of the forall assigment was for optimization in defining arrays, but since most compilers are pretty good at optimizing loops for memory purposes, the use of this command is declining. For an array of size (750,750,750), it takes about half a second to define the array with the above definition using either of the two methods.

The last method, the nested implied do loop, is probably about as useful as the forall assignment, but it introduces the command reshape, which has some useful applications for clusters and file I/O. Suppose we have a 2D array defined in the nested do loop,

do j=1,jMax
   do i=1,iMax
     B(i,j) = i*j
   enddo
enddo

we can collapse the inner loop as,

do j=1,jMax
   B(1:iMax,j) = (/ (i*j, i=1,iMax) /)
enddo

But in collapsing this other do loop,

B =(/ ((/ (i*j, i=1,iMax) /), j=1,jMax) /)

we would get an "incompatible ranks" error. This is because B is defined as a 2-dimensional array, B(iMax,jMax), and the above doubly-collapsed loop is a 1-dimensional array with dimensions of (iMax*jMax). To get into the 2D array from the 1D array, we can use the reshape command. For ease, I have written the (iMax*jMax) array as a separate line:

tempArray=(/ ((/ (i*j, i=1,iMax) /), j=1,jMax) /)
B =reshape( tempArray, (/iMax,jMax/) )

Extending this to a 3D array should be straight-forward.

Nested do loops should be the primary method for defining your arrays as it is the most clear method (you just have to remember that Fortran is column-major). The forall assignment, while sometimes convenient, just is not as useful as it was thought to be, probably should be dropped from codes that currently employ it, a should be marked as obsolescent in future updates to Fortran. The nested implied do loops probably should not be used for defining arrays, but it is useful for understanding how reshape works for use in other areas of codes.

Comments are always welcome!