Computing correctly rounded logarithms with fixed-point operations - - PowerPoint PPT Presentation

▶
computing correctly rounded logarithms with fixed point
SMART_READER_LITE
LIVE PREVIEW

Computing correctly rounded logarithms with fixed-point operations - - PowerPoint PPT Presentation

May 14, 2023 195 likes •1.08k views

Computing correctly rounded logarithms with fixed-point operations Julien Le Maire, Florent de Dinechin, Jean-Michel Muller and Nicolas Brunie Outline Introduction and context Algorithm Results and comparisons for libm log Bonus: a

slide-1
SLIDE 1

Computing correctly rounded logarithms with fixed-point operations

Julien Le Maire, Florent de Dinechin, Jean-Michel Muller and Nicolas Brunie

slide-2
SLIDE 2

Outline

Introduction and context Algorithm Results and comparisons for libm log Bonus: a floating-point in, fixed-point out variant Conclusions

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 2

slide-3
SLIDE 3

Software evaluating elementary function

This work is about libm functions: prototype: floating-point in, floating-point out, e.g. double log(double x);

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 3

slide-4
SLIDE 4

Software evaluating elementary function

This work is about libm functions: prototype: floating-point in, floating-point out, e.g. double log(double x); implementation: Nearly 100% of the implementations and literature use floating-point

Integer-based implementations: only on processors without FPU (StrongArm, ST200)

1960 1980 2000 IEEE-754 (64 bits) mainstream floating-point 32-bits mainstream integer

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 3

slide-5
SLIDE 5

The times they are a-changing

1960 1980 2000 IEEE-754 (64 bits) mainstream floating-point 32-bits mainstream integer

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 4

slide-6
SLIDE 6

The times they are a-changing

1960 1980 2000 IEEE-754 (64 bits) mainstream floating-point 32-bits 64-bits mainstream integer

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 4

slide-7
SLIDE 7

The times they are a-changing

1960 1980 2000 IEEE-754 (64 bits) mainstream floating-point 32-bits 64-bits mainstream integer

At the same time, architectural changes such as unified integer/floating-point registers.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 4

slide-8
SLIDE 8

The times they are a-changing

1960 1980 2000 IEEE-754 (64 bits) mainstream floating-point 32-bits 64-bits mainstream integer

At the same time, architectural changes such as unified integer/floating-point registers.

This work:

Re-evaluate the idea of implementing floating-point functions using integer arithmetic.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 4

slide-9
SLIDE 9

Integer now seems better than floating-point

most operations are faster on integers, especially addition

3-5 cycles in floating point 1 cycle in integer (more or less defines the processor cycle time)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 5

slide-10
SLIDE 10

Integer now seems better than floating-point

most operations are faster on integers, especially addition

3-5 cycles in floating point 1 cycle in integer (more or less defines the processor cycle time)

64 bits of significand is better than 52

if you can predict the value of the exponent, exponent bits are wasted bits convert precision to speed?

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 5

slide-11
SLIDE 11

Integer now seems better than floating-point

most operations are faster on integers, especially addition

3-5 cycles in floating point 1 cycle in integer (more or less defines the processor cycle time)

64 bits of significand is better than 52

if you can predict the value of the exponent, exponent bits are wasted bits convert precision to speed?

modern 64-bit machines offer all the integer instructions we need

addition multiplication 64x64 → 128 (mulq) count leading zeroes, shifts (lzcnt, bsr)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 5

slide-12
SLIDE 12

Integer now seems better than floating-point

most operations are faster on integers, especially addition

3-5 cycles in floating point 1 cycle in integer (more or less defines the processor cycle time)

64 bits of significand is better than 52

if you can predict the value of the exponent, exponent bits are wasted bits convert precision to speed?

modern 64-bit machines offer all the integer instructions we need

addition multiplication 64x64 → 128 (mulq) count leading zeroes, shifts (lzcnt, bsr)

Fast small multiprecision out of the box:

mainstream compilers (gcc, clang, icc) support __int_128 addition 128x128 → 128 (add, adc) shift on two registers (shld, shrd) multiplication, etc...

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 5

slide-13
SLIDE 13

No vectorization yet

Integer SIMD/vector support still lagging behind FP until recently, no vector multiplication AVX512: 52-bit vector multiplication (recycling mantissa multiplier) So all hope is not left.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 6

slide-14
SLIDE 14

This work: an experiment

Implementing the floating-point logarithm function using only integer arithmetic for performance (previous work motivated by lack of FP hardware) with state of the art accuracy: correct rounding

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 7

slide-15
SLIDE 15

This work: an experiment

Implementing the floating-point logarithm function using only integer arithmetic for performance (previous work motivated by lack of FP hardware) with state of the art accuracy: correct rounding Why the log? Because it seemed the easiest function for this.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 7

slide-16
SLIDE 16

Main results

worst-case execution time now in a factor 3 of the best faithful implementations

An improvement of a factor 5 over previous state of the art

average time almost twice better than current glibc. proposal of a floating-point in, fix-point out variant of the log function

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 8

slide-17
SLIDE 17

Outline

Introduction and context Algorithm Results and comparisons for libm log Bonus: a floating-point in, fixed-point out variant Conclusions

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 9

slide-18
SLIDE 18

Logarithm, the mathematical version

1 2 3 4 5 6 7 −2 −1 1 2 y x y = ln(x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 10

slide-19
SLIDE 19

Logarithm, the mathematical version

ln (a × b) = ln (a) + ln (b) 1 2 3 4 5 6 7 −2 −1 1 2 y x y = ln(x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 10

slide-20
SLIDE 20

Logarithm, the mathematical version

ln (a × b) = ln (a) + ln (b) ln (ba) = a × ln(b) 1 2 3 4 5 6 7 −2 −1 1 2 y x y = ln(x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 10

slide-21
SLIDE 21

Logarithm, the mathematical version

ln (a × b) = ln (a) + ln (b) ln (ba) = a × ln(b) Taylor: for x small, ln(1 + x) ≈ x − x2/2 + x3/3... 1 2 3 4 5 6 7 −2 −1 1 2 y x y = ln(x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 10

slide-22
SLIDE 22

Logarithm, the floating-point version

The floating point version of the natural logarithm is called log (you will also find log2 and log10 and a few others) ∀x ∈ F64 log(x) = ◦ (ln(x))

1 2 3 4 5 6 7 −2 −1 1 2 y x y = ln(x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 11

slide-23
SLIDE 23

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-24
SLIDE 24

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-25
SLIDE 25

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-26
SLIDE 26

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin CRLibm refinement of Ziv’s technique: First step: quick-and-dirty evaluation of ln(x) (just accurate enough to ensure correct rounding in most cases)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-27
SLIDE 27

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin CRLibm refinement of Ziv’s technique: First step: quick-and-dirty evaluation of ln(x) (just accurate enough to ensure correct rounding in most cases) test if rounding can be decided

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-28
SLIDE 28

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin CRLibm refinement of Ziv’s technique: First step: quick-and-dirty evaluation of ln(x) (just accurate enough to ensure correct rounding in most cases) test if rounding can be decided if not (rarely), recompute ln(x) with the worst-case accuracy

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-29
SLIDE 29

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin CRLibm refinement of Ziv’s technique: First step: quick-and-dirty evaluation of ln(x) (just accurate enough to ensure correct rounding in most cases) test if rounding can be decided if not (rarely), recompute ln(x) with the worst-case accuracy Trade-off between first and second steps: MeanTime = Time(1st step) + Pr[need 2nd step] · Time(2nd step)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-30
SLIDE 30

On-demand accuracy

Muller and Lefèvre solved the table maker dilema for log

Computing the log with an error ≤ 2−113 enables correct rounding real numbers two consecutive floating-point numbers computed logarithm, with error margin CRLibm refinement of Ziv’s technique: First step: quick-and-dirty evaluation of ln(x) (just accurate enough to ensure correct rounding in most cases) test if rounding can be decided if not (rarely), recompute ln(x) with the worst-case accuracy Trade-off between first and second steps: MeanTime = Time(1st step) + Pr[need 2nd step] · Time(2nd step) Best so far: Time(2nd step) ≈ 10 × Time(1st step) This work: Time(2nd step) ≈ 2 × Time(1st step)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 12

slide-31
SLIDE 31

The big picture

  • 1. Filter special cases (negative numbers, ∞, ...)
  • 2. Argument range reduction
  • 3. Polynomial approximation
  • 4. Solution reconstruction
  • 5. Error evaluation and rounding test
  • 6. If more accuracy needed:

Rerun the steps 3 and 4 with the worst-case accuracy.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 13

slide-32
SLIDE 32

First argument range reduction

input = 2E · (1 + x) ln(input) = E · ln (2) + ln (1 + x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 14

slide-33
SLIDE 33

First argument range reduction

input = 2E · (1 + x) ln(input) = E · ln (2) + ln (1 + x) Evaluation algorithm: approximate ln (1 + x) with a polynomial p(x) degree needed: at least 26 evaluate E · ln (2) add both terms

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 14

slide-34
SLIDE 34

Tang’s range reduction

1 + x:

1.

fractional part on 52 bits

1 + x1: 1.

fractional part on 7 bits

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 15

slide-35
SLIDE 35

Tang’s range reduction

1 + x:

1.

fractional part on 52 bits

1 + x1: 1.

fractional part on 7 bits

rx:

0.

fractional part on 18 bits

A table, addressed by x1, the 7 most significand bits of x, stores rx ≈ 1 1 + x1 ≈ 1 1 + x and ln 1 rx

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 15

slide-36
SLIDE 36

Tang’s range reduction

1 + x:

1.

fractional part on 52 bits

1 + x1: 1.

fractional part on 7 bits

rx:

0.

fractional part on 18 bits

1 + y:

1000000

.

fractional part on 64 bits plus 6 implicit zeros

  • 64
  • 70
  • 6

A table, addressed by x1, the 7 most significand bits of x, stores rx ≈ 1 1 + x1 ≈ 1 1 + x and ln 1 rx

  • Define

1 + y = rx · (1 + x) ≈ 1 Then ln(1 + x) = ln(1 + y) + ln( 1

rx )

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 15

slide-37
SLIDE 37

Tang’s range reduction algorithm

1 + x:

1.

fractional part on 52 bits

1 + x1: 1.

fractional part on 7 bits

rx:

0.

fractional part on 18 bits

1 + y:

1000000

.

fractional part on 64 bits plus 6 implicit zeros

  • 64
  • 70
  • 6

Extract the index x1 Read, from a table addressed by x1, both rx and ln( 1

rx )

compute y = rx · (1 + x) − 1 (exactly) approximate ln (1 + y) with a polynomial p(y) Degree needed: 8 add it all: ln(input) ≈ E · ln (2) + p (y) + ln( 1 rx )

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 16

slide-38
SLIDE 38

Tang’s range reduction

1 + x:

1.

fractional part on 52 bits

1 + x1: 1.

fractional part on 7 bits

rx:

0.

fractional part on 18 bits

1 + y:

1000000

.

fractional part on 64 bits plus 6 implicit zeros

  • 64
  • 70
  • 6

y = rx · (1 + x) − 1 With 1 + x on 53 bits we can tabulate rx on 18 bits: the exact product would need 71 bits but we can predict the 7 leading bits ... so we can let them overflow quietly and use a 64 × 64 → 64 multiplication.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 17

slide-39
SLIDE 39

Two levels of Tang reduction

  • 64
  • 76
  • 12

1 + x:

1.

fractional part on 52 bits

1 + x1: 0.

fractional part on 7 bits

rx:

0.

fractional part on 9 bits

1 + y:

100000

.

fractional part on 60 bits including 5 zeros

1 + y1: 100000

.

fractional part on 13 bits

ry:

0.

fractional part on 15 bits

1 + z:

1000000000000

.

fractional part on 64 bits plus 12 implicit zeros

x ∈ [0, 1) y ∈

  • 0, 2−6.41504

z ∈

  • 0, 2−12.6747

x1 takes 64 different values y1 takes 96 different values the whole reduction of x to z is computed exactly in 64-bit int.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 18

slide-40
SLIDE 40

A few Pareto points in the design space

Table size (bytes) degree 1st degree 2nd 39,936 3 5 12,288 3 6 4,032 4 7 2,240 4 8 2,016 4 9 900 5 10 594 6 12 298 7 14

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 19

slide-41
SLIDE 41

Why stop at two levels of reduction?

Answer is: diminushing return. For a target accuracy of 2−60: interval of x degree needed No reduction [−1/2, 1/2] 29 1 level [−2−7, 2−7] 8 2 levels [−2−12, 2−12] 4 3 levels [−2−18, 2−18] 3 Adding more levels will cost more operations than it saves...

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 20

slide-42
SLIDE 42

Polynomial approximation (advertisement)

We want to approximate log(1 + z) on an interval around 0. Use the (now standard) tool set to obtain it. Sollya:

finds a machine-efficient polynomial P(z) computes a safe bound on the approximation error P(z) − ln(1 + z)

Gappa: bounds the accumulation of rounding errors when evaluating P(z) in C We obtain a Coq proof of the error: real numbers computed approximation of ln(1 + z), with relative error margin

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 21

slide-43
SLIDE 43

Reconstructing the solution

input = 2e · (1 + x)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-44
SLIDE 44

Reconstructing the solution

input = 2e · 1 rx · (1 + y)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-45
SLIDE 45

Reconstructing the solution

input = 2e · 1 rx · 1 ry · (1 + z)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-46
SLIDE 46

Reconstructing the solution

input = 2e · 1 rx · 1 ry · (1 + z) ln(input) = e · ln(2) + ln(rx

−1)

+ ln(ry

−1)

+ ln(1 + z)

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-47
SLIDE 47

Reconstructing the solution

input = 2e · 1 rx · 1 ry · (1 + z) ln(input) = e · ln(2) + ln(rx

−1)

+ ln(ry

−1)

+ ln(1 + z)

11

  • 53
  • 117

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

“If we can predict the exponents, exponent bits are wasted bits”

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-48
SLIDE 48

Reconstructing the solution

input = 2e · 1 rx · 1 ry · (1 + z) ln(input) = e · ln(2) + ln(rx

−1)

+ ln(ry

−1)

+ ln(1 + z)

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

“If we can predict the exponents, exponent bits are wasted bits”

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-49
SLIDE 49

Reconstructing the solution

input = 2e · 1 rx · 1 ry · (1 + z) ln(input) = e · ln(2) + ln(rx

−1)

+ ln(ry

−1)

+ ln(1 + z)

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

“If we can predict the exponents, exponent bits are wasted bits”

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-50
SLIDE 50

Reconstructing the solution

input = 2e · 1 rx · 1 ry · (1 + z) ln(input) = e · ln(2) + ln(rx

−1)

+ ln(ry

−1)

+ ln(1 + z)

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

“If we can predict the exponents, exponent bits are wasted bits”

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 22

slide-51
SLIDE 51

Error evaluation

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 23

slide-52
SLIDE 52

Error evaluation

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

ǫ < (|e|) · 2−117

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 23

slide-53
SLIDE 53

Error evaluation

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

ǫ < (|e| + 1 + 1) · 2−117

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 23

slide-54
SLIDE 54

Error evaluation

11

  • 53
  • 117
  • 12

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

ǫ <

  • |e| + 1 + 1 + P(z) · 2−59

· 2−117

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 23

slide-55
SLIDE 55

Rounding test

Simple technique: compute the two bounds of the interval, and see if they round to the same mantissa (two additions, a xor and a shift) real numbers

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 24

slide-56
SLIDE 56

Rounding test

Simple technique: compute the two bounds of the interval, and see if they round to the same mantissa (two additions, a xor and a shift) real numbers

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 24

slide-57
SLIDE 57

Rounding test

Simple technique: compute the two bounds of the interval, and see if they round to the same mantissa (two additions, a xor and a shift) real numbers For comparison, the proof of the floating-point-based rounding test (invented by Ziv and used in CRLibm) is an 18-page paper that took 20 years to publish...

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 24

slide-58
SLIDE 58

Second step

Use 3 words instead of 2 for the precomputed ln(2) Use a much more accurate polynomial:

with coefficients on 128 bits instead of 64 (but z is still only a 64-bit number) and using a higher degree (7) polynomial

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 25

slide-59
SLIDE 59

Outline

Introduction and context Algorithm Results and comparisons for libm log Bonus: a floating-point in, fixed-point out variant Conclusions

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 26

slide-60
SLIDE 60

Implementation parameters

  • f correctly rounded implementations

glibc crlibm-td crlibm-de cr-FixP degree pol. 1 3/8 6 7 4 degree pol. 2 20 12 14 7 tables size 13 Kb 8192 bytes 6144 bytes 4032 bytes % accurate phase N/A 1.5 0.4 4.4

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 27

slide-61
SLIDE 61

Average and max runing time (in processor cycles)

Pentium timing (lower is better) cycles MKL glibc crlibm cr-de cr-FixP avg time 25 90 69 46 49 max time 25 11,554 642 410 79 Timing breakdown on two processors cycles Core i5 Bostan System glibc newlib 90 105 quick phase alone 42 94 accurate phase alone 74 181 both phases (avg time) 49 121 both phases (max time) 79 225 Slanted means: no correct rounding

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 28

slide-62
SLIDE 62

Average and max runing time (in processor cycles)

Pentium timing (lower is better) cycles MKL glibc crlibm cr-de cr-FixP avg time 25 90 69 46 49 max time 25 11,554 642 410 79 Timing breakdown on two processors cycles Core i5 Bostan System glibc newlib 90 105 quick phase alone 42 94 accurate phase alone 74 181 both phases (avg time) 49 121 both phases (max time) 79 225 Slanted means: no correct rounding

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 28

slide-63
SLIDE 63

Average and max runing time (in processor cycles)

Pentium timing (lower is better) cycles MKL glibc crlibm cr-de cr-FixP avg time 25 90 69 46 49 max time 25 11,554 642 410 79 Timing breakdown on two processors cycles Core i5 Bostan System glibc newlib 90 105 quick phase alone 42 94 accurate phase alone 74 181 both phases (avg time) 49 121 both phases (max time) 79 225 Slanted means: no correct rounding

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 28

slide-64
SLIDE 64

Average and max runing time (in processor cycles)

Pentium timing (lower is better) cycles MKL glibc crlibm cr-de cr-FixP avg time 25 90 69 46 49 max time 25 11,554 642 410 79 Timing breakdown on two processors cycles Core i5 Bostan System glibc newlib 90 105 quick phase alone 42 94 accurate phase alone 74 181 both phases (avg time) 49 121 both phases (max time) 79 225 Slanted means: no correct rounding

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 28

slide-65
SLIDE 65

Outline

Introduction and context Algorithm Results and comparisons for libm log Bonus: a floating-point in, fixed-point out variant Conclusions

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 29

slide-66
SLIDE 66

Floating-point in, fixed-point out

  • utput: fixed-point, 11 bits integer part, 53 bit fractional part

integer part fraction

  • target absolute accuracy 2−52
  • utput

absolute table Core i5 Bostan format accuracy size cycles cycles Fix64 2−52 2304 24 66 Fix128 2−116 4032 60 179 double (libm) 2−42 90 105 Fix64 is the code of the first step only, without the conversion to float.

tweak: poly degree 3 only for abs. accuracy 2−59

Fix128 is the code of the second step only, without the conversion to float.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 30

slide-67
SLIDE 67

Motivation

TKF91 : DNA sequence alignment algorithm dynamic programming algorithm: alignment as a path within a 2D array. borders of an array initialized with log-likelihoods then array filled using recurrence formulae that involve only max and + operations. All current implementations of this algorithm use a floating-point array, but int64 + and max are 1-cycle, vectorizable, and exact operations; absolute accuracy of initialization logs: up to 2−42 with FP log, 2−52 with FixP log.

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 31

slide-68
SLIDE 68

Only partial experiments

Improvement in accuracy measured No noticeable improvement in performance

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 32

slide-69
SLIDE 69

Outline

Introduction and context Algorithm Results and comparisons for libm log Bonus: a floating-point in, fixed-point out variant Conclusions

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 33

slide-70
SLIDE 70

Conclusions

Better range reduction thanks to a wider format ... leading to improvements in polynomial degree and table size Faster multiprecision due to higher precision Able to reuse some computations of the fast step in the accurate step Alternative rounding test for acurate step Probability to launch 2nd step is high, but this is acceptable since 2nd step is so cheap

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 34

slide-71
SLIDE 71

Conclusions

Better range reduction thanks to a wider format ... leading to improvements in polynomial degree and table size Faster multiprecision due to higher precision Able to reuse some computations of the fast step in the accurate step Alternative rounding test for acurate step Probability to launch 2nd step is high, but this is acceptable since 2nd step is so cheap Competitive against state-of-the-art Worst case improved compared to others implementations Second step-only version is a viable alternative

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 34

slide-72
SLIDE 72

Conclusions

Better range reduction thanks to a wider format ... leading to improvements in polynomial degree and table size Faster multiprecision due to higher precision Able to reuse some computations of the fast step in the accurate step Alternative rounding test for acurate step Probability to launch 2nd step is high, but this is acceptable since 2nd step is so cheap Competitive against state-of-the-art Worst case improved compared to others implementations Second step-only version is a viable alternative

Limitations:

Require 64 bits integer support No support for vectorization

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 34

slide-73
SLIDE 73

Thanks for your attention

Any question ?

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 35

slide-74
SLIDE 74

Outline

Some code details on the solution reconstruction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 36

slide-75
SLIDE 75

Special cases: businesss as usual

/∗ r e i n t e r p r e t x to manipulate i t s b i t s more e a s i l y ∗/ uint64_t i n p u t b i t s = (( union { double d ; uint64_t u ; }){ i n p u t } ) . u ; i n t xe = i n p u t b i t s > > 52; /∗ f i l t e r the s p e c i a l c a s e s : ! ( x i s normalized and 0 < x < +I n f ) ∗/ i f (0 x7FEu <= ( unsigned ) xe − 1u ) { /∗ x = + − 0 : r a i s e a DivideByzero , r e t u r n −I n f ∗/ i f (( x b i t s & ~(1 u l l < < 63)) == 0) return −1.0/0.0; /∗ x < 0 . 0 : r a i s e a I n v a l i d O p e r a t i o n , r e t u r n a qNaN ∗/ i f (( x b i t s & (1 u l l < < 63)) != 0) return ( x−x ) / 0 ; /∗ x = qNaN : r e t u r n a qNaN x = sNaN : r a i s e a I n v a l i d O p e r a t i o n , r e t u r n a qNaN x = +I n f : r e t u r n +I n f ∗/ i f ( xe != 0) return x+x ; /∗ x subnormal : change x to a normalized number ∗/ e l s e { i n t u = c l z 6 4 ( x b i t s ) − 12; x b i t s < <= u + 1; xe −= u ; } } xe −= 1023;

Only interesting line: the subnormal management

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 37

slide-76
SLIDE 76

Argument reduction

/∗ i n p u t = 2^ xe ∗ (1 + X) ∗/ uint64_t x = i n p u t b i t s & 0xFFFFFFFFFFFFFull ; /∗ 1 + X = (1/ Ri ) ∗ Y ∗/ uint8_t x1 = x > > (52 − ARG_REDUC_1_PREC) ; uint64_t y = argReduc1_inv [ x1 ] ∗ ( x + (1 u l l < < 5 2 ) ) ; _ _ b u i l t i n _ p r e f e t c h (& argReduc1_log [ x1 ] , 0 , 0 ) ; /∗ Y = (1/S) ∗ (1 + Z) ∗/ /∗ with dZ = dz /2^(52 + ARG_REDUC_1_SIZE + ARG_REDUC_2_SIZE) and 1/S = argReduc2 [ s i ] . v a l /2^ARG_REDUC_2_SIZE ∗/ uint8_t y1 = ( y > > (52 + ARG_REDUC_1_SIZE − ARG_REDUC_2_PREC)) − (1 u < < ARG_REDUC_2_PREC) ; uint64_t z = argReduc2_inv [ y1 ] ∗ y ; // +1 part removed by

  • v e r f l o w

_ _ b u i l t i n _ p r e f e t c h (& argReduc2_log [ y1 ] , 0 , 0 ) ;

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 38

slide-77
SLIDE 77

Polynomial evaluation

/∗ Polynomial approximation

  • f

log (1+Z)/Z ~= P(Z) and Z∗P(Z) ∗/ s t a t i c const uint64_t a4 = UINT64_C(0 x3ffc147cb4539237 ) ; s t a t i c const uint64_t a3 = UINT64_C(0 x5555553dc70f0dfd ) ; s t a t i c const uint64_t a2 = UINT64_C(0 x 7 f f f f f f f f f d 5 7 4 f d ) ; s t a t i c const uint64_t a1 = UINT64_C(0 x f f f f f f f f f f f f f f f a ) ; uint64_t pz = a1 − ( highmul ( z , a2 − ( highmul ( z , a3 − ( highmul ( z , a4 ) > > 12) ) > > 12) ) > > 1 2 ) ; uint128_t zpz = f u l l m u l ( z , pz ) ; /∗ Polynomial approximation without s h i f t s : r e p l a c e highmul ( z , a ) > > 12 with highmul ( z , a > > 12) ∗/ s t a t i c const uint64_t a4 = UINT64_C(0 x0000000003ffc147 ) ; s t a t i c const uint64_t a3 = UINT64_C(0 x0000005555553dc6 ) ; s t a t i c const uint64_t a2 = UINT64_C(0 x 0 0 0 7 f f f f f f f f f d 5 7 ) ; s t a t i c const uint64_t a1 = UINT64_C(0 x f f f f f f f f f f f f f f f a ) ; uint64_t pz = a1 − highmul ( z , a2 − highmul ( z , a3 − highmul ( z , a4 ) ) ) ; uint128_t zpz = f u l l m u l ( z , pz ) ;

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 39

slide-78
SLIDE 78

Reconstructing the solution

/∗ Compute part

  • f

the r e s u l t that don ’ t depend on Z ( xe ∗ log (2) + lo g (1/ Ri ) + log (1/ Si )) ∗/ uint128_t c s t p a r t = f u l l i m u l ( xe , log2fw_mid ) + UINT128 (( int64_t ) xe ∗ log2fw_high , 0) // f u l l m u l not needed + UINT128 ( argReduc1 [ r i ] . log_hi , argReduc1 [ r i ] . log_lo ) + UINT128 ( argReduc2 [ s i ] . log_hi , argReduc2 [ s i ] . log_lo ) ; /∗ Polynomial approximation

  • f

log (1+Z)/Z ~= P(Z) and Z∗P(Z) ∗/ . . . /∗ Assemble the two parts , compute sign , mantissa and exponent ∗/ uint128_t l o n g r e s = c s t p a r t + ( zpz > > (11 + IMPLICIT_ZEROS ) ) ; uint64_t s i g n = − ( HI ( l o n g r e s ) > > 6 3 ) ; // 0

  • r

~0 // i f s i g n != 0 , t h i s i s l o n g r e s = ~ l o n g r e s (−a = ~a + 1) // to avoid the +1 approx , do : // l o n g r e s = (( int64_t ) s i g n + l o n g r e s ) ^ UINT128 ( sign , s i g n ) ; l o n g r e s ^= UINT128 ( sign , s i g n ) ; i n t u = c l z 6 4 ( HI ( l o n g r e s )) + 1 ; i n t exponent = 11 − u ; uint64_t mantissa = ( HI ( l o n g r e s ) < < u ) | (LO( l o n g r e s ) > > (64 − u ) ) ;

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 40

slide-79
SLIDE 79

Rounding test and conversion

/∗ Compute the maximal a b s o l u t e e r r o r ( a l i g n e d with l o n g r e s ) : + 2 + abs ( xe ) f o r xe ∗ log (2 ) , log (1/ Ri ) and log (1/ Si ) + 1 + zpz>>(POLYNOMIAL_PREC+IMPLICIT_ZEROS+11) f o r the polynom I f r e s u l t ∗(1 + − maxRelErr ) are not rounded to the same number , we uint64_t maxAbsErr = 3 + abs ( xe ) + ( HI ( zpz ) > > (POLYNOMIAL_PREC + I uint64_t maxRelErr = ( maxAbsErr > > (64 − u )) + 1; i f ( ( ( mantissa + maxRelErr ) ^ ( mantissa − maxRelErr )) > > 11) { return log_rn_accurate ( c s t p a r t , z , xe ) ; } /∗ Assemble the computed r e s u l t ∗/ uint64_t r e s u l t b i t s = (( uint64_t ) s i g n < < 63) + (( uint64_t )( exponent +1023) < < 52) + ( mantissa > > 12) + (( mantissa > > 11) & 1 ) ; /∗ round to n e a r e s t ∗/ return ( union { uint64_t u ; double d ; }){ r e s u l t b i t s } . d ;

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 41

slide-80
SLIDE 80

Outline

Some code details on the solution reconstruction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 42

slide-81
SLIDE 81

Reconstructing the solution

11

  • 53
  • 117
  • 11

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 43

slide-82
SLIDE 82

Reconstructing the solution

11

  • 53
  • 117
  • 11

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

1 floating-point fraction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 43

slide-83
SLIDE 83

Reconstructing the solution

11

  • 53
  • 117
  • 11

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

1 floating-point fraction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 43

slide-84
SLIDE 84

Reconstructing the solution

11

  • 53
  • 117
  • 11

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

1 floating-point fraction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 43

slide-85
SLIDE 85

Reconstructing the solution

11

  • 53
  • 117
  • 11

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

−1 floating-point fraction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 43

slide-86
SLIDE 86

Reconstructing the solution

11

  • 53
  • 117
  • 11

e · ln(2): ln(rx

−1):

ln(ry

−1):

P(z) ≈ ln(1 + z): sum:

1 floating-point fraction

  • J. Le Maire, F. de Dinechin, J.-M. Muller and N. Brunie

Computing correctly rounded logarithm with fixed-point operations 43