Path: utzoo!utgpu!news-server.csri.toronto.edu!rutgers!ucsd!sdd.hp.com!think.com!linus!linus!eachus From: eachus@linus.mitre.org (Robert I. Eachus) Newsgroups: comp.benchmarks Subject: Re: BYTE and benchmarks (anecdote) Message-ID: Date: 11 Jan 91 00:33:58 GMT References: <1991Jan8.101912.15157@odin.diku.dk> <1872@ruunsa.fys.ruu.nl> Sender: usenet@linus.mitre.org Distribution: comp Organization: The Mitre Corporation, Bedford, MA Lines: 301 Nntp-Posting-Host: aries.mitre.org In-reply-to: hooft@ruunsa.fys.ruu.nl's message of 10 Jan 91 08:47:19 GMT In article <1872@ruunsa.fys.ruu.nl> hooft@ruunsa.fys.ruu.nl (Rob Hooft) writes: program test do i=1,46 write(*,*) ifib(i) enddo end function ifib(i) c c The actual function is: c fib=(((1.+5.**.5)/2.)**i-((1.-5.**.5)/2.)**i)/5.**.5 c ifib=int(aint(((1.d0+5.d0**.5d0)/2.d0)**i)/5.d0**.5d0)+1 end Not bad, but did you check your output? A good compiler will figure out that the square root of 5 is a common subexpression, and good compilers will even evaluate it at compile time. But I'd have to run your code through a FORTRAN compiler to see if you got the truncation right--it looks wrong. (The second term gets small rapidly but it alternates in sign.) Also I don't think your function works correctly for the magic numbers of computer science -- zero and one. My program (in Ada for readability, it should easily translate to C or Pascal) is append below, together with output. Some explanations are in order. The 68020 is (relatively) slow on floating point operations compared to chips which have floating point on chip (I say the 68020, not the 68881, since most of the overhead is moving stuff on and off chip). Also the math package I used calls the Unix exp function, which is fairly slow. Doing the exp by calling the 68881 directly would probably speed up the result using Fib4 below immensely. If you are wondering about the complexity of the code for the "simple" solution, I had to make sure that it would not overflow for N=46... Solutions which happily return junk are not acceptable! These function all raise CONSTRAINT_ERROR for values out of range, and only for values out of range. (If you run the Ada on your machine, check that the accuracy of float is acceptable. I only used FLOAT here to allow me to use the Verdix Math package.) I also have another math package to try out. If you are wondering about the comparison to 42, that was just to make sure that the compiler didn't optimize stuff away... Robert I. Eachus package Fibonacci is function Fib1(N: in Natural) return Natural; function Fib2(N: in Natural) return Natural; function Fib3(N: in Natural) return Natural; function Fib4(N: in Natural) return Natural; function Fib5(N: in Natural) return Natural; function Fib6(N: in Natural) return Natural; pragma INLINE(Fib1,Fib2,Fib3,Fib4,Fib5,Fib6); -- Procedure call overhead will be larger than most of these -- functions. end Fibonacci; with Math; package body Fibonacci is SQRT5 : constant := 2.23606_79774_99789_69640_91736_68731_27624; GOLDEN_RATIO : constant := 1.61803_39887_49894_84820_45868_34365_63812; NAT_LOG_PHI : constant := 0.48121_18250_59603_44749_77589_13424_36842; -- borrowed from the Verdix supplied Math package, copied to make the -- code more portable. InvSqrt5: constant := 1.0/Sqrt5; Phi: constant Float := Golden_Ratio; Rho: constant Float := Sqrt5/2.0 - 0.5; -- also = 1.0/Phi. -- Phi is a common designator for the golden ratio, -- Rho is just my choice of name for the inverse. -- Phi and Rho are declared Float to determine the -- accuracy used in evaluating the expressions below. function Fib1(N: in Natural) return Natural is N1: Natural := 0; N2: Natural := 1; begin if N in 0..1 then return N; else if N mod 2 = 0 then for I in 1..N/2-1 loop N1 := N1 + N2; N2 := N2 + N1; end loop; return N1+N2; else for I in 1..N/2 loop N1 := N1 + N2; N2 := N2 + N1; end loop; return N2; end if; end if; end Fib1; function Fib2(N: in Natural) return Natural is begin return Natural((Phi**N-(-Rho)**N) / Sqrt5 ); end Fib2; -- Now to cheat. For large N, Rho**N < 0.5, so all we have to -- do is round..., How large must N be? Let's see, 1/Sqrt5 -- rounds ok, and so does Phi/Sqrt5, and Rho**2 < 0.5 so no -- check need be made... function Fib3(N: in Natural) return Natural is begin return Natural(Phi**N/Sqrt5); end Fib3; -- Now, since a float to an integer power is supposed to be -- calculated iteratively in Ada (or iteratively followed by an -- inverse, let's see what happens when we explicitly force the power -- operation to be done by exponentiation. function Fib4(N: in Natural) return Natural is begin return Natural(Math.Exp(Nat_Log_Phi*Float(N))/Sqrt5); end Fib4; -- See if explicit looping for exponentiation is faster or slower... function Fib5(N: in Natural) return Natural is Temp: Float := InvSqrt5; begin for I in 1..N loop Temp := Temp * Phi; end loop; return Natural(Temp); end Fib5; -- Last but not least, see if multiplication is faster... function Fib6(N: in Natural) return Natural is begin return Natural(InvSqrt5*Phi**N); end Fib6; end Fibonacci; with Text_IO; with Fibonacci; use Fibonacci; with Calendar; procedure Main is Start, Stop: Calendar.Time; function "-" (L,R: Calendar.Time) return Duration renames Calendar."-"; package Int_IO is new Text_IO.Integer_IO(Integer); package Duration_IO is new Text_IO.Fixed_IO(Duration); use Int_IO, Duration_IO; -- If you don't know Ada, ignore these lines. They are only here -- to make right aligned columnar output of integers and output -- of timestamps easier. begin -- First check that everything works: for I in 0..46 loop Put(I,5); Put(Fib1(I),12); Put(Fib2(I),12); Put(Fib3(I),12); Put(Fib4(I),12); Put(Fib5(I),12); Put(Fib6(I),12); Text_IO.New_Line; end loop; delay 0.5; -- wait for quiescence Start := Calendar.Clock; for I in 1..1_000_000 loop if Fib1(I mod 47) = 42 then Text_IO.Put_Line("Ooops. :-)"); end if; end loop; Stop := Calendar.Clock; Text_IO.Put(" Fib1 took, on average, "); Put(Stop-Start,4,2); Text_IO.Put_Line(" microseconds."); delay 0.5; -- wait for quiescence Start := Calendar.Clock; for I in 1..1_000_000 loop if Fib2(I mod 47) = 42 then Text_IO.Put_Line("Ooops. :-)"); end if; end loop; Stop := Calendar.Clock; Text_IO.Put(" Fib2 took, on average, "); Put(Stop-Start,4,2); Text_IO.Put_Line(" microseconds."); delay 0.5; -- wait for quiescence Start := Calendar.Clock; for I in 1..1_000_000 loop if Fib3(I mod 47) = 42 then Text_IO.Put_Line("Ooops. :-)"); end if; end loop; Stop := Calendar.Clock; Text_IO.Put(" Fib3 took, on average, "); Put(Stop-Start,4,2); Text_IO.Put_Line(" microseconds."); delay 0.5; -- wait for quiescence Start := Calendar.Clock; for I in 1..1_000_000 loop if Fib4(I mod 47) = 42 then Text_IO.Put_Line("Ooops. :-)"); end if; end loop; Stop := Calendar.Clock; Text_IO.Put(" Fib4 took, on average, "); Put(Stop-Start,4,2); Text_IO.Put_Line(" microseconds."); delay 0.5; -- wait for quiescence Start := Calendar.Clock; for I in 1..1_000_000 loop if Fib5(I mod 47) = 42 then Text_IO.Put_Line("Ooops. :-)"); end if; end loop; Stop := Calendar.Clock; Text_IO.Put(" Fib5 took, on average, "); Put(Stop-Start,4,2); Text_IO.Put_Line(" microseconds."); delay 0.5; -- wait for quiescence Start := Calendar.Clock; for I in 1..1_000_000 loop if Fib6(I mod 47) = 42 then Text_IO.Put_Line("Ooops. :-)"); end if; end loop; Stop := Calendar.Clock; Text_IO.Put(" Fib6 took, on average, "); Put(Stop-Start,4,2); Text_IO.Put_Line(" microseconds."); end Main; Output on a Sun 3/260, using VADS Ada 6.03: aries% fib 0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 1 1 1 1 1 1 3 2 2 2 2 2 2 4 3 3 3 3 3 3 5 5 5 5 5 5 5 6 8 8 8 8 8 8 7 13 13 13 13 13 13 8 21 21 21 21 21 21 9 34 34 34 34 34 34 10 55 55 55 55 55 55 11 89 89 89 89 89 89 12 144 144 144 144 144 144 13 233 233 233 233 233 233 14 377 377 377 377 377 377 15 610 610 610 610 610 610 16 987 987 987 987 987 987 17 1597 1597 1597 1597 1597 1597 18 2584 2584 2584 2584 2584 2584 19 4181 4181 4181 4181 4181 4181 20 6765 6765 6765 6765 6765 6765 21 10946 10946 10946 10946 10946 10946 22 17711 17711 17711 17711 17711 17711 23 28657 28657 28657 28657 28657 28657 24 46368 46368 46368 46368 46368 46368 25 75025 75025 75025 75025 75025 75025 26 121393 121393 121393 121393 121393 121393 27 196418 196418 196418 196418 196418 196418 28 317811 317811 317811 317811 317811 317811 29 514229 514229 514229 514229 514229 514229 30 832040 832040 832040 832040 832040 832040 31 1346269 1346269 1346269 1346269 1346269 1346269 32 2178309 2178309 2178309 2178309 2178309 2178309 33 3524578 3524578 3524578 3524578 3524578 3524578 34 5702887 5702887 5702887 5702887 5702887 5702887 35 9227465 9227465 9227465 9227465 9227465 9227465 36 14930352 14930352 14930352 14930352 14930352 14930352 37 24157817 24157817 24157817 24157817 24157817 24157817 38 39088169 39088169 39088169 39088169 39088169 39088169 39 63245986 63245986 63245986 63245986 63245986 63245986 40 102334155 102334155 102334155 102334155 102334155 102334155 41 165580141 165580141 165580141 165580141 165580141 165580141 42 267914296 267914296 267914296 267914296 267914296 267914296 43 433494437 433494437 433494437 433494437 433494437 433494437 44 701408733 701408733 701408733 701408733 701408733 701408733 45 1134903170 1134903170 1134903170 1134903170 1134903170 1134903170 46 1836311903 1836311903 1836311903 1836311903 1836311903 1836311903 Fib1 took, on average, 47.52 microseconds. Fib2 took, on average, 141.74 microseconds. Fib3 took, on average, 83.50 microseconds. Fib4 took, on average, 275.32 microseconds. Fib5 took, on average, 317.32 microseconds. Fib6 took, on average, 82.98 microseconds. aries% -- Robert I. Eachus When the dictators are ready to make war upon us, they will not wait for an act of war on our part." - Franklin D. Roosevelt