This is getting pretty deep into the internals, but ...
I was thinking of replacing my implementation of ##cabs (complex absolute value) in _num.scm with glibc's hypot. Turns out that ##cabs (which is based on one of Kahan's papers) returns the correctly-rounded answer more often than hypot (which Kahan also had a hand in): Here's a count of how many times hypot and magnitude returned answers that were off by between 1/2 and 1 ulp:
(time (test 100000))
(n: 100000 hypot-errors: 12906 magnitude-errors: 4557) (time (test 100000)) 3.600013 secs real time 3.599883 secs cpu time (3.592672 user, 0.007211 system) 401 collections accounting for 0.047208 secs real time (0.047259 user, 0.000060 system) 1942508792 bytes allocated no minor faults no major faults
These shouldn't really be called "errors" just because they have slightly bigger than optimal floating-point errors.
The correctly-rounded answers were computed using my computable reals package.
I haven't yet compared computation speed, but it appears I'll leave ##cabs as is.
This is on Ubuntu 18.04 on x86-64.
Brad
Afficher les réponses par date
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; }
On 3/30/19 1:27 PM, Bradley Lucier wrote:
After looking at this data, I recommend adding primitives
(flhypot x y) (##flhypot x y)
to _num.scm, adding
FLHYPOT
to gambit.h and changing the compiler to inline flhypot and ##flhypot in the right circumstances.
We could use the existing implementation of ##cabs as the implementation for backends without hypot.
Brad
On 3/30/19 4:06 PM, Bradley Lucier wrote:
I ran the simple code at the end of this message, and got the following times:
(time (do ((i 0 (fx+ i 1))) ((fx= i 10000000)) (##cabs arg1) (##cabs arg2))) 0.755936 secs real time 0.755946 secs cpu time (0.755289 user, 0.000657 system) 528 collections accounting for 0.028905 secs real time (0.028845 user, 0.000110 system) 2560000000 bytes allocated 1152 minor faults no major faults (time (do ((i 0 (fx+ i 1))) ((fx= i 10000000)) (hypot arg1-x arg1-y) (hypot arg2-x arg2-y))) 1.173139 secs real time 1.173110 secs cpu time (1.169110 user, 0.004000 system) 190 collections accounting for 0.372126 secs real time (0.372004 user, 0.000017 system) 1556944 bytes allocated 2053 minor faults no major faults
I don't generally use c-lambda, but this seems pretty heavy duty, when the same number C calls to hypot takes roughly 0.167994 seconds.
Brad
(declare (standard-bindings) (extended-bindings) (block) (not safe))
(define hypot (c-lambda (double double) double "hypot"))
(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)))
(define arg1-x (real-part arg1)) (define arg1-y (imag-part arg1)) (define arg2-x (real-part arg2)) (define arg2-y (imag-part arg2))
(time (do ((i 0 (fx+ i 1))) ((fx= i 10000000)) (hypot arg1-x arg1-y) (hypot arg2-x arg2-y) ))