Slow? Ada??

[Bill Findlay <yaldnif.w@blueyonder.co.uk>]

I was amused to read this ...

On 11/07/2013 00:21, in article krkq9u$fo9$1@loke.gir.dk, "Randy Brukardt"
<randy@rrsoftware.com> wrote:
...
> For things like Sine, the best one can do is copy an implementation from a
> cookbook -- it makes no sense to even imagine trying to create a new one.

...for reasons nothing to do with that controversy.

When debugging Whetstone Algol under my KDF9 emulator I looked at the KDF9
assembly code programming of the arctan function and was completely baffled
by it, so I asked my former colleague, Michael Jamieson to investigate.
He successfully reconstructed the mathemetics behind it, which are very
non-obvious (to put it mildly).

A couple of days ago I idly wondered how well this 50 years old algorithm
would compare with modern implementations, so I wrote a test program to race
it against GNAT GPL 2013's arctan.  I expected to find that 2013 would spank
1963's botty.  To my astonishment, the old algorithm was faster.

Digging down, I found that Ada.Numerics.Long_Elementary_Functions.arctan
does quite a bit of argument range reduction and result error checking, but
finally invokes the "fpatan" opcode of the x86_64 CPU.

The 1963 algorithm, expressed in Ada 2012 but compiled with aggressive
optimization, is only 17% slower than that single CPU opcode!

Note also that, although it was designed for a machine with 48-bit floating
point, it gets the same accuracy as the hardware method on a machine with
64-bit floating point.

Here is the code, with the observed results included as commentary:

with Ada.Numerics.Long_Elementary_Functions;
with Ada.Text_IO;
with CPU_Timing;
with System.Machine_Code;

use  Ada.Numerics.Long_Elementary_Functions;
use  Ada.Text_IO;
use  CPU_Timing;
use  System.Machine_Code;

procedure arctan_test64 is

   -- typical output:
   -- R = 100000000 evaluations
   -- checksum = 5.00000005000000E+07, loop   time per repetition = 6 ns
   -- checksum = 4.38824577044418E+07, P51V15 time per evaluation = 49 ns
   -- checksum = 4.38824577044418E+07, arctan time per evaluation = 53 ns
   -- checksum = 4.38824577044418E+07, fpatan time per evaluation = 41 ns

   -- P51V15 is the KDF9 algorithm used in Whetstone Algol, ca. 1963
   -- See http://www.findlayw.plus.com/KDF9/Arctan%20Paper.pdf

   type vector is array (0..5) of Long_Float;

   V : constant vector := (  28165298.0 / 1479104550.0,
                             28165300.0 / 1479104550.0,
                             56327872.0 / 1479104550.0,
                            113397760.0 / 1479104550.0,
                            179306496.0 / 1479104550.0,
                           1073741824.0 / 1479104550.0);

   function P51V15 (x : Long_Float) return Long_Float with Inline;

   function P51V15 (x : Long_Float)
   return Long_Float is
      A : Long_Float := 1.0;
      S : Long_Float := V(0);
      B : vector;
   begin
       B(0) := sqrt(x*x + 1.0);
       -- 4 AGM (Arithmetic-Geometric Mean) cycles give a set of values ...
       for i in 0..3 loop
          A := (A+B(i)) * 0.5;
          B(i+1) := sqrt(A*B(i));
       end loop;
       -- ... that is subjected to a convergence acceleration process:
       for i in 1..5 loop
          S := S + V(i)*B(i-1);
       end loop;
       return x / S;
   end P51V15;

   -- this is the hardware arctan function of the x86_64 Core i7 CPU
   function fpatan (X : Long_Float) return Long_Float with Inline;

   function fpatan (X : Long_Float)
   return Long_Float is
      Result  : Long_Float;
   begin
      Asm(Template  => "fld1"
                     & Character'Val(10) -- LF
                     & Character'Val(9)  -- HT
                     & "fpatan",
          Outputs   => Long_Float'Asm_Output ("=t", Result),
          Inputs    => Long_Float'Asm_Input  ("0", X));
      return Result;
   end fpatan;

   R : constant := 1e8; -- number of loop repetitions

   function ns_per_rep (c : CPU_Usage_in_Microseconds)
   return Natural is
   begin
      return Natural(c * 1e3 / R);
   end ns_per_rep;

   x : Long_Float;
   c : CPU_Timer;
   l : CPU_Usage_in_Microseconds;
   t : CPU_Usage_in_Microseconds;

begin
   Put_Line("R =" & Integer'Image(R) & " evaluations");

   -- determine the fixed overhead time
   x := 0.0;
   Reset_Timer(c);
   for i in 1..R loop
      x := x + Long_Float(i)/Long_Float(R);
   end loop;
   l := User_CPU_Time_Since(c);
   Put_Line("checksum =" & x'Img & ", "
          & "loop   time per repetition ="
          & Natural'Image(ns_per_rep(l)) & " ns");

   x := 0.0;
   Reset_Timer(c);
   for i in 1..R loop
      x := x + P51V15(Long_Float(i)/Long_Float(R));
   end loop;
   t := User_CPU_Time_Since(c) - l;
   Put_Line("checksum =" & x'Img & ", "
          & "P51V15 time per evaluation ="
          & Natural'Image(ns_per_rep(t)) & " ns");

   x := 0.0;
   Reset_Timer(c);
   for i in 1..R loop
      x := x + arctan(Long_Float(i)/Long_Float(R));
   end loop;
   t := User_CPU_Time_Since(c) - l;
   Put_Line("checksum =" & x'Img & ", "
          & "arctan time per evaluation ="
          & Natural'Image(ns_per_rep(t)) & " ns");

   x := 0.0;
   Reset_Timer(c);
   for i in 1..R loop
      x := x + fpatan(Long_Float(i)/Long_Float(R));
   end loop;
   t := User_CPU_Time_Since(c) - l;
   Put_Line("checksum =" & x'Img & ", "
          & "fpatan time per evaluation ="
          & Natural'Image(ns_per_rep(t)) & " ns");

end arctan_test64;

The difference in time between Ada.Numerics.Long_Elementary_Functions.arctan
and  the fpatan function above is due to the afore-mentioned range reduction
and checking.  The KDF9 programmers in 1963 were less punctillious about
such matters than we rightly expect Ada to be nowadays.

-- 
Bill Findlay
with blueyonder.co.uk;
use  surname & forename;