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; }