Xref: utzoo sci.math:16705 comp.lang.c:38230
Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!usc!elroy.jpl.nasa.gov!decwrl!pa.dec.com!hollie.rdg.dec.com!ryn.mro4.dec.com!allvax.enet.dec.com!jroth
From: jroth@allvax.enet.dec.com (Jim Roth)
Newsgroups: sci.math,comp.lang.c
Subject: Re: fibonacci Numbers
Message-ID: <4369@ryn.mro4.dec.com>
Date: 13 Apr 91 22:29:26 GMT
Sender: guest@ryn.mro4.dec.com
Followup-To: sci.math
Organization: Digital Equipment Corporation
Lines: 90


In article <1991Apr12.051844.15063@milton.u.washington.edu>, amigo@milton.u.washington.edu (The Friend) writes...
> 
>     Can someone send me source code for a routine to calculate 
>Fibonacci numbers to _x_ (x being input from user)? Additionally it'd
>be great if this found the PRIME numbers in the set. It could print its
>results as it went through (if that makes it any easier).

Hmm... sounds like a tough homework assignment :-)

Anyhow, here's a FORTRAN program that computes a few Fibonacci numbers
using a Cauchy's residue theorem and the fast Fourier transform; you can
always translate it to c...

- Jim

	PROGRAM FIBONACCI
C
C Compute some Fibonacci numbers
C
	IMPLICIT NONE

	INTEGER*4 LOGN, NPTS, NFIB, IPASS, I, J, K, L
	COMPLEX*16 T, U, V, W, Z, Z0, A(64)
	REAL*8 PI, X, R

	DATA PI/3.14159265358979323846264338D0/
C
C Generating function for Fibonacci numbers
C
	COMPLEX*16 F
	F(Z) = 1.0/(1.0-Z-Z**2)

	LOGN = 6
	NPTS = 2**LOGN
	NFIB = 2**(LOGN-1)

	Z0 = 0.0
	R = 0.4
C
C Approximate contour integrals with trapezoidal rule using FFT
C
	W = DCMPLX(DCOS(2.0*PI/NPTS), DSIN(2.0*PI/NPTS))
	V = R
	A(1) = F(Z0+V)
	A(2) = F(Z0-V)
	J = 0
	DO I = 1,NPTS/2-1
	  V = V*W
	  K = NPTS
100	  CONTINUE
	  K = K/2
	  J = J.XOR.K
	  IF ((J.AND.K).EQ.0) GOTO 100
	  A(J+1) = F(Z0+V)
	  A(J+2) = F(Z0-V)
	ENDDO

	L = 1
	DO IPASS = 1,LOGN
	  U = 1.0
	  W = DCMPLX(DCOS(PI/L), -DSIN(PI/L))
	  DO J = 1,L
	    DO I = J,NPTS,2*L
	      K = I+L
	      T = A(K)*U
	      A(K) = A(I)-T
	      A(I) = A(I)+T
	    ENDDO
	    U = U*W
	  ENDDO
	  L = L+L
	ENDDO
C
C Denormalize series coefficients
C
	X = 1.0D0/NPTS
	DO I = 1,NPTS
	  A(I) = A(I)*X
	  X = X/R
	ENDDO
C
C Show the answer
C
	TYPE 101, (I, DREAL(A(I)), I=1,NFIB)
101	FORMAT (1X, I10, F28.1)

	STOP

	END