[gambit-list] magnitude vs hypot

Bradley Lucier lucier at math.purdue.edu
Sat Mar 30 13:27:48 EDT 2019


On 3/29/19 11:58 PM, Marc Feeley wrote:
> Were you thinking of using glibc for better performance?

Generally, I have tried not to rewrite code from glibc in the Gambit 
runtime, but the complex absolute value code for ##cabs was in the same 
Kahan paper as the complex transcendental functions so I didn't notice 
that hypot() could do the same job.

When I looked into it, both codes achieve < 1ulp error bounds, but I 
wanted to test speed and accuracy before replacing the code.

So the accuracy testing showed that hypot() gives an answer that has > 
1/2 ulp error about 12.9% of the time, while Gambit's algorithm gives an 
error > 1/2 ulp about 4.5% of the time:

(n: 100000 hypot-errors: 12906 magnitude-errors: 4557)

<Somewhat later ...>

Gambit's algorithm for ##cabs is from the paper:

https://people.freebsd.org/~das/kahan86branch.pdf

I translated the Gambit algorithm into C, and tested it on

Ubuntu 18.04
model name	: Intel(R) Core(TM) i7-4770 CPU @ 3.40GHz
gcc version 7.3.0 (Ubuntu 7.3.0-27ubuntu1~18.04)
gcc -O3 -fno-trapping-math -fno-math-errno -march=native -W -Wall 
test-hypot.c -lm

for which I got the timings (first is hypot, second is Gambit's 
algorithm translated to C:

firefly:~/text/courses/computation/computational-reals/src/FFT> ./a.out
41209517.673652 0.164801
41209517.673652 0.255371

so 0.164801 seconds for hypot and 0.255371 seconds for Gambit's 
algorithm.  (The first number is just something so that gcc wouldn't 
eliminate the loops.)

I also tested it on my Raspberry Pi:

model name	: ARMv7 Processor rev 4 (v7l)
Pi 3 Model B
gcc version 4.9.2 (Raspbian 4.9.2-10+deb8u2)
gcc -O3 -fno-trapping-math -fno-math-errno -mcpu=cortex-a53 
-mfpu=neon-vfpv4 -mfloat-abi=hard -std=c99 -W -Wall test-hypot.c -lm

which gave the timings

sweety-pi:~/text/courses/computation/computational-reals/src/FFT> ./a.out
41209517.673652 5.460862
41209517.673652 2.310178

So now Gambit's algorithm is faster.

Of course, with Gambit's boxing of floating-point variables and the 
trampoline, etc., Gambit's actual code will be slower, indeed we find:

(declare (standard-bindings)
          (extended-bindings)
          (block)
          (not safe))

(define arg1 (make-rectangular 1.2345678 (fl- 1.2345678 0.1)))
(define arg2 (make-rectangular 1.2345678 0.1))

(time
  (do ((i 0 (fx+ i 1)))
      ((fx= i 10000000))
    (##cabs arg1)
    (##cabs arg2)))

which after compilation with Gambit configured with

firefly:~/text/courses/computation/computational-reals/src/FFT> gsi -v
v4.9.1 20180930122740 x86_64-unknown-linux-gnu "./configure 'CC=gcc 
-march=native -D___CAN_IMPORT_CLIB_DYNAMICALLY' '--enable-single-host' 
'--enable-shared' '--enable-multiple-versions'"

gives on my x86-64 box

firefly:~/text/courses/computation/computational-reals/src/FFT> gsi 
time-cabs
(time (do ((i 0 (fx+ i 1))) ((fx= i 10000000)) (##cabs arg1) (##cabs arg2)))
     0.752108 secs real time
     0.752072 secs cpu time (0.751388 user, 0.000684 system)
     528 collections accounting for 0.052368 secs real time (0.052385 
user, 0.000095 system)
     2560000000 bytes allocated
     1154 minor faults
     no major faults

and after compilation with Gambit configured with

v4.9.3 20180930122740 armv7l-unknown-linux-gnueabihf "./configure 
'CC=gcc -mcpu=cortex-a53 -mfpu=neon-vfpv4 -mfloat-abi=hard 
-D___CAN_IMPORT_CLIB_DYNAMICALLY' '--enable-single-host' 
'--enable-shared' '--enable-multiple-versions'"

gives on my Raspberry Pi:

sweety-pi:~/text/courses/computation/computational-reals/src/FFT> gsi 
-:m100000 time-cabs
(time (do ((i 0 (fx+ i 1))) ((fx= i 10000000)) (##cabs arg1) (##cabs arg2)))
     7.193316 secs real time
     7.190000 secs cpu time (7.070000 user, 0.120000 system)
     24 collections accounting for 0.024148 secs real time (0.020000 
user, 0.000000 system)
     2559955416 bytes allocated
     25212 minor faults
     no major faults

Here's the code for test-hypot.c:


#include <math.h>
#include <stdio.h>
#include <time.h>

double myhypot(double x, double y)
{
   double
     r2   = 1.4142135623730951,
     r2p1 = 2.414213562373095,
     t2p1 = 1.2537167179050217e-16;
   double t;

   x = fabs(x);
   y = fabs(y);

   if (x < y) {
     double temp = x;
     x = y;
     y = temp;
   }
   if (__builtin_expect(isinf(y),0))
     x = y;
   t = x - y;
   if (__builtin_expect(!(x == INFINITY) && !(t == x), 1)) {
     if (t > y) {
       double s = x / y;
       s += sqrt(1.0 + s*s);
       return x + y/s;
     } else {
       double s = t/y;
       t = s * (2.0 + s);
       s = r2p1 + (s + (t2p1 + t/(r2 + sqrt(2.0 + t))));
       return x + y/s;
     }
   }
   else
     return x;
}

int main()
{
   double sum;
   double x1;
   clock_t start, end;
   long N = 10000000;

   x1 = 1.2345678;
   sum = 0.;
   start = clock();
   for (int i = 0; i < N; i++) {
     /* arguments exercise both main paths through code */
     sum += hypot(x1, x1-0.1) + hypot(x1, 0.1);
     x1 += 0.0000001;
   }
   end = clock();
   printf("%lf %lf \n", sum, ((double) (end - start)) / CLOCKS_PER_SEC);

   x1 = 1.2345678;
   sum = 0.;
   start = clock();
   for (int i = 0; i < N; i++) {
     /* arguments exercise both main paths through code */
     sum += myhypot(x1, x1-0.1) + myhypot(x1, 0.1);
     x1 += 0.0000001;
   }
   end = clock();
   printf("%lf %lf \n", sum, ((double) (end - start)) / CLOCKS_PER_SEC);


   return 1;
}



More information about the Gambit-list mailing list