Prev: what are some of the uses for ada?
Next: Library_Symbol_File in Stand-alone Library Projects
From: Nasser M. Abbasi on 9 Jun 2010 06:49
I never used complex variables before in Ada, so for the last 3 hrs I
was learning how to do it. I wrote a simple program, to find the DFT of
an array of 3 elements {1,2,3} (DFT=discrete Fourier transform).

The definition of DFT is one equation with summation, here it is, first
equation you'll see:

http://en.wikipedia.org/wiki/Discrete_Fourier_transform

Since I have not programmed in Ada nor in FORTRAN for a looong time, I
am very rusty, I program in Mathematica and sometimes in matlab these
days, but I wanted to try Ada on complex numbers.

I also wrote a FORTRAN equivalent of the small Ada function. Here is
below the Ada code, and the FORTRAN code. Again, do not scream too much
if it is not good code, I just learned this now, I am sure this can be
improved a lot.

And for comparing, I show the Matlab and the Mathematica output just to
make sure.


====== START ADA ============
--
-- dtf.adb, compiled with GNAT 4.3.4 20090804 (release) 1
-- under CYGWIN 1.7.5
-- gnatmake dft.adb
--
-- ./dft.exe
-- ( 6.00000E+00, 0.00000E+00)
-- (-1.50000E+00, 8.66026E-01)
-- (-1.50000E+00,-8.66025E-01)
-- $


with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;

with Ada.Numerics; use Ada.Numerics;

with Ada.Numerics.Complex_Elementary_Functions;
use Ada.Numerics.Complex_Elementary_Functions;

with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;

procedure dft is
N : positive := 3;
J : constant complex :=(0.0,1.0); -- SQRT(-1)
X : array(0 .. N-1) of Complex := (others=>(0.0,0.0));
data : array(0 .. N-1) of float :=(1.0,2.0,3.0);

begin
FOR k in X'range LOOP
FOR m in data'range LOOP
X(k) := X(k) + data(m) * exp(- J*(2.0*Pi)/float(N) *
float(m*k) );
END LOOP;
put(X(k)); new_line;
END LOOP;

end dft;
================== END ADA ==============

======= FORTRAN code ===========
! dtf.f90, compiled with GCC 4.3.4
! under CYGWIN 1.7.5
! gfortran -Wall dft.f90
! ./a.exe
! ( 6.0000000 , 0.0000000 )
! ( -1.4999999 , 0.86602557 )
! ( -1.5000005 ,-0.86602497 )
!

PROGRAM dft

IMPLICIT NONE

INTEGER, PARAMETER :: N = 3
COMPLEX, parameter :: J =(0,1)

REAL, parameter :: Pi = ACOS(-1.0)
INTEGER :: k,m
COMPLEX, dimension(N) :: X
REAL, dimension(N) :: data=(/1.0,2.0,3.0/)

DO k=1,N
X(k)=(0,0)
DO m=1,N
X(k) = X(k) + data(m) * EXP(-1.0*J*2.0*Pi/N *(m-1)*(k-1) )
END DO
print *,X(k)

END DO

END PROGRAM dft
==================================

==== Matlab code ====
EDU>> fft([1,2,3])'

ans =

6.0000
-1.5000 - 0.8660i
-1.5000 + 0.8660i
===============================

=== Mathematica ====
In[5]:= Chop[Fourier[{1, 2, 3}, FourierParameters -> {1, -1}]]

Out[5]= {6., -1.5 + 0.8660254037844386*I, -1.5 - 0.8660254037844386*I}
=========================

Conclusion:
I actually liked the Ada implementation more than FORTRAN because:

1. In Ada, I did not have to change the index of m and k in the
summation to reflect the 1-off per the definition of DFT.
DFT starts from 0 to N-1. In Ada, using 'Range and defining the arrays
to go from 0 .. N-1 solved the problem.

2. In Ada, the compiler complained more times more about types being
mixed up. I placed float() around the places it complained about.

3. It actually took me less time to do the Ada function than the FORTRAN
one, even though I am equally not familiar with both at this time :)

ok, this was a fun learning exercise


--Nasser
From: Ludovic Brenta on 9 Jun 2010 07:26
On Jun 9, 12:49 pm, "Nasser M. Abbasi" <n...(a)12000.org> wrote:
> I never used complex variables before in Ada, so for the last 3 hrs I
> was learning how to do it. I wrote a simple program, to find the DFT of
> an array of 3 elements {1,2,3} (DFT=discrete Fourier transform).
>
> The definition of DFT is one equation with summation, here it is, first
> equation you'll see:
>
> http://en.wikipedia.org/wiki/Discrete_Fourier_transform
>
> Since I have not programmed in Ada nor in FORTRAN for a looong time, I
> am very rusty, I program in Mathematica and sometimes in matlab these
> days, but I wanted to try Ada on complex numbers.
>
> I also wrote a FORTRAN equivalent of the small Ada function.  Here is
> below the Ada code, and the FORTRAN code.  Again, do not scream too much
> if it is not good code, I just learned this now, I am sure this can be
> improved a lot.
>
> And for comparing, I show the Matlab and the Mathematica output just to
> make sure.
>
> ====== START ADA ============
> --
> -- dtf.adb, compiled with GNAT 4.3.4 20090804 (release) 1
> -- under CYGWIN 1.7.5
> -- gnatmake dft.adb
> --
> -- ./dft.exe
> -- ( 6.00000E+00, 0.00000E+00)
> -- (-1.50000E+00, 8.66026E-01)
> -- (-1.50000E+00,-8.66025E-01)
> -- $
>
> with Ada.Text_IO; use Ada.Text_IO;
> with Ada.Numerics.Complex_Types; use  Ada.Numerics.Complex_Types;
>
> with Ada.Numerics; use  Ada.Numerics;
>
> with Ada.Numerics.Complex_Elementary_Functions;
> use  Ada.Numerics.Complex_Elementary_Functions;
>
> with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;
>
> procedure dft is
>      N : positive := 3;
>      J : constant complex :=(0.0,1.0);  -- SQRT(-1)
>      X : array(0 .. N-1) of Complex  := (others=>(0.0,0.0));
>      data : array(0 .. N-1) of float :=(1.0,2.0,3.0);
>
> begin
>      FOR k in X'range LOOP
>          FOR m in data'range LOOP
>              X(k) := X(k) + data(m) * exp(- J*(2.0*Pi)/float(N) *
> float(m*k) );
>          END LOOP;
>          put(X(k)); new_line;
>      END LOOP;
>
> end dft;
> ================== END ADA ==============
>
> ======= FORTRAN code ===========
> ! dtf.f90, compiled with GCC 4.3.4
> ! under CYGWIN 1.7.5
> ! gfortran -Wall dft.f90
> ! ./a.exe
> ! (  6.0000000    ,  0.0000000    )
> ! ( -1.4999999    , 0.86602557    )
> ! ( -1.5000005    ,-0.86602497    )
> !
>
> PROGRAM dft
>
>    IMPLICIT NONE
>
>    INTEGER, PARAMETER :: N = 3
>    COMPLEX, parameter :: J =(0,1)
>
>    REAL, parameter :: Pi = ACOS(-1.0)
>    INTEGER :: k,m
>    COMPLEX, dimension(N) :: X
>    REAL, dimension(N) :: data=(/1.0,2.0,3.0/)
>
> DO k=1,N
>     X(k)=(0,0)
>     DO m=1,N
>        X(k) = X(k) + data(m) * EXP(-1.0*J*2.0*Pi/N *(m-1)*(k-1) )
>     END DO
>     print *,X(k)
>
> END DO
>
> END PROGRAM dft
> ==================================
>
> ==== Matlab code ====
> EDU>> fft([1,2,3])'
>
> ans =
>
>     6.0000
>    -1.5000 - 0.8660i
>    -1.5000 + 0.8660i
> ===============================
>
> === Mathematica ====
> In[5]:= Chop[Fourier[{1, 2, 3}, FourierParameters -> {1, -1}]]
>
> Out[5]= {6., -1.5 + 0.8660254037844386*I,  -1.5 - 0.8660254037844386*I}
> =========================
>
> Conclusion:
> I actually liked the Ada implementation more than FORTRAN because:
>
> 1. In Ada, I did not have to change the index of m and k in the
> summation to reflect the 1-off per the definition of DFT.
> DFT starts from 0 to N-1. In Ada, using 'Range and defining the arrays
> to go from 0 .. N-1 solved the problem.
>
> 2. In Ada, the compiler complained more times more about types being
> mixed up. I placed float() around the places it complained about.
>
> 3. It actually took me less time to do the Ada function than the FORTRAN
> one, even though I am equally not familiar with both at this time :)
>
> ok, this was a fun learning exercise

You should use constants and named numbers instead of variables
wherever possible; this simplifies your Ada program. Also, i and j are
predefined so you do not need to declare a "new" value for J:

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;

with Ada.Numerics; use Ada.Numerics;

with Ada.Numerics.Complex_Elementary_Functions;
use Ada.Numerics.Complex_Elementary_Functions;

with Ada.Complex_Text_IO; use Ada.Complex_Text_IO;

procedure dft is
N : constant := 3; -- named number, no conversion to Float needed
X : array(0 .. N-1) of Complex := (others=>(0.0,0.0));
data : constant array(0 .. N-1) of float :=(1.0,2.0,3.0);
Two_Pi_Over_N : constant := 2 * Pi / N;
-- named number, outside the loop, like in ARM 3.3.2(9)
begin
FOR k in X'range LOOP
FOR m in data'range LOOP
X(k) := X(k) + data(m) * exp(- J*Two_Pi_Over_N *
float(m*k) );
END LOOP;
put(X(k)); new_line;
END LOOP;
end dft;

--
Ludovic Brenta.
From: Niklas Holsti on 9 Jun 2010 08:43
Nasser M. Abbasi wrote [much snipped]:
> I wrote a simple program, to find the DFT of
> an array of 3 elements {1,2,3} (DFT=discrete Fourier transform).
> I also wrote a FORTRAN equivalent of the small Ada function. Here is
> below the Ada code, and the FORTRAN code.
>
> ====== START ADA ============
> X(k) := X(k) + data(m) * exp(- J*(2.0*Pi)/float(N) *
> float(m*k) );

> ======= FORTRAN code ===========
> X(k) = X(k) + data(m) * EXP(-1.0*J*2.0*Pi/N *(m-1)*(k-1) )

There is an obscure difference in these computations, which can cause
trouble when the data vector is much longer than the N = 3 in the
example. The Ada version may then fail because the product m*k may
overflow the Integer range. It would be safer to write float(m)*float(k).

In the Fortran version, as far as I understand the Fortran rules, the
compiler may choose to multiply the integer subexpressions (m-1) and
(k-1) as integers, and then convert the product to REAL (as is done in
Nasser's Ada version), or the compiler may choose to convert separately
(m-1) to REAL and (k-1) to REAL and use REAL multiplication, as I
suggest above for the Ada code. IIUC it would be safer in Fortran, too,
to prevent integer overflow by explicit conversions of (m-1) and (k-1)
to REAL before multiplication.

--
Niklas Holsti
Tidorum Ltd
niklas holsti tidorum fi
. @ .
From: Jerry on 9 Jun 2010 19:50
On Jun 9, 4:26 am, Ludovic Brenta <ludo...(a)ludovic-brenta.org> wrote:
snip
>      Two_Pi_Over_N : constant := 2 * Pi / N;
>       -- named number, outside the loop, like in ARM 3.3.2(9)
snip
> Ludovic Brenta.

So 2 and N are universal_integer and Pi is probably a universal_real.
I didn't know that one could ever do mixed-mode arithmetic in Ada, but
this is obviously allowed here. By the way, I don't see where in ARM
3.3.2 this mixed operation is described. Line 9 on that ARM page uses
2.0*Ada.Numerics.Pi which is not mixed mode. Is this described
somewhere else?

Jerry
From: Jeffrey R. Carter on 9 Jun 2010 21:03
Jerry wrote:
>
> So 2 and N are universal_integer and Pi is probably a universal_real.
> I didn't know that one could ever do mixed-mode arithmetic in Ada, but
> this is obviously allowed here. By the way, I don't see where in ARM
> 3.3.2 this mixed operation is described. Line 9 on that ARM page uses
> 2.0*Ada.Numerics.Pi which is not mixed mode. Is this described
> somewhere else?

In ARM 4.5.5 you'll find

function "*"(Left : root_real; Right : root_integer) return root_real
function "*"(Left : root_integer; Right : root_real) return root_real
function "/"(Left : root_real; Right : root_integer) return root_real

--
Jeff Carter
"I've got to stay here, but there's no reason
why you folks shouldn't go out into the lobby
until this thing blows over."
Horse Feathers
50
 |  Next  |  Last
Pages: 1 2 3 4 5 6 7
Prev: what are some of the uses for ada?
Next: Library_Symbol_File in Stand-alone Library Projects