CPSC 355: Tutorial 17
Floating Point Operations
PhD Student
Fall 2017
Floating Point Arithmetic
ARMv8 has 32 floating point registers (sort of, they're actually 128bit, but ignore that for now)
- d0-d31 are 64 bit double precision registers
- s0-s31 are the lower 32 bits of these registers and are used for single precision FP operations
- d8-d15 are callee saved registers
- d0-d7 and d16-d31 may be overwritten by subroutines.
- d0-d7 are used to pass floating point arguments into a function
Floating Point Arithmetic
An example function...
double approximate(double x, double y, int iterations) {
// d0 contains x
// d1 contains y
// w0 contains iterations
// do some calculations
// return value goes in d0
}
double approximate(int iterations, double x, double y) {
// d0 contains x
// d1 contains y
// w0 contains iterations
// do some calculations
// return value goes in d0
}
Floating Point Arithmetic
An example function...
int approximate(double x, double y, int iterations) {
// d0 contains x
// d1 contains y
// w0 contains iterations
// do some calculations
// return value goes in w0
}
int approximate(int iterations, double x, double y) {
// d0 contains x
// d1 contains y
// w0 contains iterations
// do some calculations
// return value goes in w0
}
Floating Point Arithmetic
- fmov S/Dd, S/Dn: move S/Dn into S/Dd
- fmov S/Dd, #fpimm: move immediate value into S/Dd (pretty limited must be of the form \(\pm n / 16 \times 2^r\) where \(n \in [16, 31]\) and \(r \in [-3, 4]\))
- fcvt Sd, Dn: Convert Dn to single prec. Sd
- fcvt Dd, Sn: Convert Sn to double prec. Dd
- fadd, fsub S/Dd, S/Dn, S/Dm: Add or subtract S/Dn, S/Dm, place result in S/Dd
- fmul S/Dd, S/Dn, S/Dm: S/Dd = S/Dn * S/Dm
- fdiv S/Dd, S/Dn, S/Dm: S/Dd = S/Dn / S/Dm
- fmadd S/Dd, S/Dn, S/Dm, S/Da: S/Dd = S/Da + S/Dn * S/Dm
Floating Point Arithmetic
- fabs S/Dd, S/Dn: S/Dd = abs(S/Dn)
- fneg S/Dd, S/Dn: S/Dd = -(S/Dn)
- fcmp S/Dn, S/Dm: Compares two registers, stores result in condition flags
- fcmp S/Dn, #fpimm: Compares register and immediate value, stores result in condition flags
An example...
Taylor Series
Recall the taylor series of an infinitely differentiable function (i.e. holomorphic or analytic function)
f(x) = \displaystyle\sum^\infty_{n=0} \frac{f^{(n)}(a)}{n!}(x-a)^n
Or its Maclaurin expansion
f(x) = \displaystyle\sum^\infty_{n=0} \frac{f^{(n)}(0)}{n!}x^n
For \(e^x\) we have
e^x = \displaystyle\sum^\infty_{n=0} \frac{x^n }{n!}
Taylor Series
We can stop after some amount of terms in the sum
e^x \approx \displaystyle\sum^j_{n=0} \frac{x^n }{n!}
In a machine, eventually we start to lose precision, so we stop after the terms stop decreasing the approximation beyond a threshold
Taylor Series
#include <math.h>
#include <stdio.h>
double taylor_exp(double x) {
double n = 0;
double nfact = 1;
double exp, exp_old;
double xx = 1;
exp = 1.0;
do {
n += 1.0;
nfact *= n;
xx *= x;
exp_old = exp;
exp += xx/nfact;
} while(fabs(exp_old - exp) > 1e-10);
return exp;
}
int main(int argc, char *argv[]) {
printf("%f\n", taylor_exp(0.0));
printf("%f\n", taylor_exp(1.0));
printf("%f\n", taylor_exp(2.0));
}
Taylor Series
.data
fmt1: .string "%f\n"
threshold: .double 0r1e-10 // A double prec. constant
.text
.balign 4
.global main
main:
stp x29, x30, [sp, -16]!
mov x29, sp
fmov d0, xzr
bl taylor_exp
adrp x0, fmt1
add x0, x0, :lo12:fmt1
bl printf
fmov d0, 1.0
bl taylor_exp
adrp x0, fmt1
add x0, x0, :lo12:fmt1
bl printf
fmov d0,2.0
bl taylor_exp
adrp x0, fmt1
add x0, x0, :lo12:fmt1
bl printf
ldp x29, x30, [sp], 16
ret
Taylor Series
taylor_exp: // d0 is the first argument
n_s = 16
nfact_s = 24
exp_s = 32
exp_old_s = 40
xx_s = 48
stp x29, x30, [sp, -64]! // Standard stack setup/teardown
mov x29, sp
fmov d1, xzr // Convert 0 to a double
str d1, [x29, n_s] // n = 0.0
fmov d1, 1.0
str d1, [x29, nfact_s] // nfact = 1.0
str d1, [x29, xx_s] // xx = 1.0
str d1, [x29, exp_s] // exp = 1.0
exp_do_top:
ldr d1, [x29, n_s] // loading and storing to addrs works just as it does for x
fmov d2, 1.0
fadd d1, d1, d2 // d1 = d1 + 1
str d1, [x29, n_s] // n = d1
// Build up the factorial over each iteration of the loop
ldr d2, [x29, nfact_s]
fmul d2, d2, d1
str d2, [x29, nfact_s] // nfact = nfact * n = n * (n-1)!
ldr d3, [x29, xx_s]
fmul d3, d0, d3
str d3, [x29, xx_s] // xx = x*xx = x^n
Taylor Series
ldr d4, [x29, exp_s]
str d4, [x29, exp_old_s] // exp_old = exp
fdiv d2, d3, d2 // d2 = x^n / nfact
fadd d5, d4, d2 // d5 = exp + x^n/nfact
str d5, [x29, exp_s] // exp = d5
ldr d1, [x29, exp_s]
ldr d2, [x29, exp_old_s]
fsub d1, d1, d2 // d1 = exp - exp_old
fabs d1, d1 // d1 = fabs(exp - exp_old)
adrp x0, threshold // Load the addr of the threshold
add x0, x0, :lo12: threshold
ldr d2, [x0] // load threshold value into d2
fcmp d1, d2 // fabs(exp - exp_old) cmp threshold
b.gt exp_do_top // if it still contributes to the upper 10 digits, keep going
ldr d0, [x29, exp_s] // load ~exp(x) into d0 (return value)
ldp x29, x30, [sp], 64 // stack teardown
ret
Next Day
Work Day
CPSC 355: Tutorial 17
By Joshua Horacsek
CPSC 355: Tutorial 17
- 1,571