diff: extras/complex_number changes

Fergus Henderson fjh at cs.mu.oz.au
Wed Jul 16 17:52:51 AEST 1997
Previous message: diff: changes to CVSROOT admin files
Next message: diff: type_name for higher order, round 2
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
extras/complex_numbers/*:
	Add an Mmakefile to build the complex number library, using the
	new support for user-defined libraries.
	Move the test cases into a new subdirectory `tests'.
	Add a new subdirectory `samples', containing an example
	module `fft.m' that uses the library.
	Add a README file.

extras/complex_numbers/complex.m:
	Add new functions `complex/1' (converts float to complex),
	`cis' (cos + i * sin), and `conj'.
	Rename `norm' as `abs2' (square of absolute value).

cvs diff: Diffing .
Index: README
===================================================================
RCS file: README
diff -N README
--- /dev/null	Wed Jul 16 17:48:20 1997
+++ README	Tue Jul 15 20:12:38 1997
@@ -0,0 +1,6 @@
+This directory contains a library for handling complex numbers.
+
+The subdirectory `tests' contains a test case.
+
+The subdirectory `samples' contains an example module using this library.
+
Index: complex.m
===================================================================
RCS file: /home/staff/zs/imp/mercury/extras/complex_numbers/complex.m,v
retrieving revision 1.1
diff -u -r1.1 complex.m
--- complex.m	1997/07/10 09:08:44	1.1
+++ complex.m	1997/07/16 05:42:54
@@ -31,12 +31,18 @@
 % both constructing and pattern matching; with intermodule optimization
 % enabled, the compiler should generate equally good code for it.
 
+	% convert float to complex
+:- func complex(float) = complex.
+
 	% extract real part
 :- func real(complex) = float.
 
 	% extract imaginary part
 :- func imag(complex) = float.
 
+	% square of absolute value
+:- func abs2(complex) = float.
+
 	% absolute value (a.k.a. modulus)
 :- func abs(complex) = float.
 
@@ -45,9 +51,8 @@
 	% for all Z, -pi < arg(Z) and arg(Z) =< pi.
 :- func arg(complex) = float.
 
-	% norm (square of absolute value)
-	% XXX is `norm' the right terminology?
-:- func norm(complex) = float.
+	% complex conjugate
+:- func conj(complex) = complex.
 
 	% addition
 :- func complex + complex = complex.
@@ -85,6 +90,9 @@
 :- func sqrt(complex) = complex.
 :- mode sqrt(in) = out is det.
 
+	% cis(Theta) = cos(Theta) + i * sin(Theta)
+:- func cis(float) = complex.
+
 	% polar_to_complex(R, Theta):
 	% conversion from polar coordinates
 :- func polar_to_complex(float, float) = complex.
@@ -100,6 +108,8 @@
 :- implementation.
 :- import_module float, math.
 
+complex(Real) = cmplx(Real, 0.0).
+
 real(cmplx(Real, _Imag)) = Real.
 imag(cmplx(_Real, Imag)) = Imag.
 
@@ -123,12 +133,14 @@
 + cmplx(R, I) = cmplx(+ R, + I).
 - cmplx(R, I) = cmplx(- R, - I).
 
-norm(cmplx(R, I)) = R*R + I*I.
+abs2(cmplx(R, I)) = R*R + I*I.
 
-abs(Z) = sqrt(norm(Z)).
+abs(Z) = sqrt(abs2(Z)).
 
 arg(cmplx(R, I)) = atan2(I, R).
 
+conj(cmplx(R, I)) = cmplx(R, -I).
+
 sqr(cmplx(Re0, Im0)) = cmplx(Re, Im) :-
 	Re = Re0 * Re0 - Im0 * Im0,
 	Im = 2.0 * Re0 * Im0.
@@ -144,5 +156,7 @@
 polar_to_complex(Magnitude, Theta) = cmplx(Real, Imag) :-
 	Real = Magnitude * cos(Theta),
 	Imag = Magnitude * sin(Theta).
+
+cis(Theta) = cmplx(cos(Theta), sin(Theta)).
 
 %------------------------------------------------------------------------------%
Index: complex_lib.m
===================================================================
RCS file: complex_lib.m
diff -N complex_lib.m
--- /dev/null	Wed Jul 16 17:48:20 1997
+++ complex_lib.m	Mon Jul 14 20:55:29 1997
@@ -0,0 +1,5 @@
+:- module complex_lib.
+:- import_module complex, float, imag.
+:- import_module complex_imag, imag_complex.
+:- import_module float_imag, imag_float.
cvs diff: complex_test.exp was removed, no comparison available
cvs diff: complex_test.m was removed, no comparison available
cvs diff: Diffing tests
Index: tests/Mmakefile
===================================================================
RCS file: Mmakefile
diff -N Mmakefile
--- /dev/null	Wed Jul 16 17:51:42 1997
+++ Mmakefile	Mon Jul 14 20:58:35 1997
@@ -0,0 +1,21 @@
+# COMPLEX_NUMBER_DIR specifies the location of the complex number library,
+# e.g. `/foo/bar/mercury-<version>/extra/complex_numbers'.
+COMPLEX_NUMBER_DIR = ..
+
+# The following stuff tells Mmake to use the complex number library
+VPATH = $(COMPLEX_NUMBER_DIR):$(MMAKE_VPATH)
+MCFLAGS = -I$(COMPLEX_NUMBER_DIR) $(EXTRA_MCFLAGS)
+MLFLAGS = -R$(COMPLEX_NUMBER_DIR) -L$(COMPLEX_NUMBER_DIR) $(EXTRA_MLFLAGS)
+MLLIBS = -lcomplex_lib $(EXTRA_MLLIBS)
+C2INITFLAGS = $(COMPLEX_NUMBER_DIR)/complex_lib.init
+
+MAIN_TARGET = complex_test.res
+
+depend: complex_test.depend
+
+complex_test.out: complex_test
+	./complex_test > complex_test.out
+
+complex_test.res: complex_test.out complex_test.exp
+	diff -c complex_test.out complex_test.exp > complex_test.res
+
cvs diff: tests/complex_test.exp is a new entry, no comparison available
cvs diff: tests/complex_test.m is a new entry, no comparison available
cvs diff: Diffing samples
Index: samples/Mmakefile
===================================================================
RCS file: Mmakefile
diff -N Mmakefile
--- /dev/null	Wed Jul 16 17:51:42 1997
+++ Mmakefile	Mon Jul 14 18:52:17 1997
@@ -0,0 +1,12 @@
+# COMPLEX_NUMBER_DIR specifies the location of the complex number library,
+# e.g. `/foo/bar/mercury-<version>/extra/complex_numbers'.
+COMPLEX_NUMBER_DIR = ..
+
+# The following stuff tells Mmake to use the complex number library
+VPATH = $(COMPLEX_NUMBER_DIR):$(MMAKE_VPATH)
+MCFLAGS = -I$(COMPLEX_NUMBER_DIR)
+MLFLAGS = -R$(COMPLEX_NUMBER_DIR) -L$(COMPLEX_NUMBER_DIR)
+MLLIBS = -lmercomplex
+
+MAIN_TARGET = fft
+depend : fft.depend
Index: samples/fft.m
===================================================================
RCS file: fft.m
diff -N fft.m
--- /dev/null	Wed Jul 16 17:51:42 1997
+++ fft.m	Mon Jul 14 16:16:23 1997
@@ -0,0 +1,127 @@
+%------------------------------------------------------------------------------%
+%
+% file: fft.m
+% main author: conway.
+% data: August 1996
+%
+% This module provides a predicate for performing the Fast Fourier Transfrom
+% (fft) on a list of complex numbers. The list must be a power of 2 in length
+% (ie 4, 8, 16, ... elements).
+%
+% The algorithm used here is derived from a few sources:
+%	- "Numerical Recipes in C": Press, et al
+%	- "Elements of Computer Music": Moore
+%	- Lecture notes from 433-325 "Mathemematical Software B": Dr Rex Harris
+%
+% It actually took the combination of these sources for me to uncover the
+% algorithm for combining the component transforms. This code is not maximally
+% efficient, but rather is intended as a clear presentation of the FFT
+% algorithm.
+%
+%------------------------------------------------------------------------------%
+:- module fft.
+%------------------------------------------------------------------------------%
+:- interface.
+
+:- import_module list, complex.
+
+:- pred fft(list(complex), list(complex)).
+:- mode fft(in, out) is det.
+
+%------------------------------------------------------------------------------%
+:- implementation.
+
+:- import_module float, int, math, require.
+
+fft(Ins, Outs) :-
+		% First put the list into bit-reversed order.
+	bit_rev(Ins, Shuffle),
+	list__length(Shuffle, NInt),
+	int__log2(NInt, R),
+	int__to_float(NInt, N),
+		% Now recombine the component transforms
+	combine(N, 1.0, R, Shuffle, Outs).
+
+:- pred bit_rev(list(T), list(T)).
+:- mode bit_rev(in, out) is det.
+
+bit_rev([], []).
+bit_rev([X], [X]).
+bit_rev([X1,X2|Xs], List) :-
+	split([X1,X2|Xs], List1, List2),
+	bit_rev(List1, List3),
+	bit_rev(List2, List4),
+	list__append(List3, List4, List).
+
+:- pred split(list(T), list(T), list(T)).
+:- mode split(in, out, out) is det.
+
+split([], [], []).
+split([_], _, _) :-
+	error("Input sequence not a power of two in length").
+split([X1, X2|Xs], [X1|Ys], [X2|Zs]) :-
+	split(Xs, Ys, Zs).
+
+:- pred combine(float, float, int, list(complex), list(complex)).
+:- mode combine(in, in, in, in, out) is det.
+
+combine(N, K, R, Ins, Outs) :-
+	(
+		R = 0
+	->
+		Outs = Ins
+	;
+		R1 = R - 1,
+		L = 1 << R1,
+			 % Split the list in half
+		divide(L, Ins, Fs0, Ss0),
+		K2 = K * 2.0,
+			% Now combine the transforms in the respective halves
+		combine(N, K2, R1, Fs0, Fs),
+		combine(N, K2, R1, Ss0, Ss),
+			% Now perform the 'butterfly'
+		xform(Fs, Ss, complex(1.0), w(K, N), Rs0, Rs1),
+		list__append(Rs0, Rs1, Outs)
+	).
+
+:- pred xform(list(complex), list(complex), complex, complex,
+		list(complex), list(complex)).
+:- mode xform(in, in, in, in, out, out) is det.
+
+xform([], [_|_], _WJN, _WKN, _, _) :-
+	error("ran out of Di's").
+xform([_|_], [], _WJN, _WKN, _, _) :-
+	error("ran out of Dj's").
+xform([], [], _WJN, _WKN, [], []).
+xform([Di0|Dis0], [Dj0|Djs0], WJN, WKN, [Di|Dis], [Dj|Djs]) :-
+	Tmp = WJN * Dj0,
+	Di = Di0 + Tmp,
+	Dj = Di0 - Tmp,
+	xform(Dis0, Djs0, WJN*WKN, WKN, Dis, Djs).
+
+:- pred divide(int, list(T), list(T), list(T)).
+:- mode divide(in, in, out, out) is det.
+
+divide(R, Ins, Fs, Ss) :-
+	(
+		R =< 0
+	->
+		Fs = [],
+		Ss = Ins
+	;
+		(
+			Ins = [],
+			error("divide error!")
+		;
+			Ins = [F|Ins1],
+			divide(R-1, Ins1, Fs1, Ss),
+			Fs = [F|Fs1]
+		)
+	).
+
+:- func w(float, float) = complex.
+
+w(J, N) = cis(Theta) :-
+	Theta = J * 2.0 * pi / N.
+
+%------------------------------------------------------------------------------%
-- 
Fergus Henderson <fjh at cs.mu.oz.au>   |  "I have always known that the pursuit
WWW: <http://www.cs.mu.oz.au/~fjh>   |  of excellence is a lethal habit"
PGP: finger fjh at 128.250.37.3         |     -- the last words of T. S. Garp.
Previous message: diff: changes to CVSROOT admin files
Next message: diff: type_name for higher order, round 2
Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]
More information about the developers mailing list