Pythagorean-Triple-Based Table Generation for Trigonometric Function - - PowerPoint PPT Presentation

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation

MetaLibm Project Annual Meeting’2015 Hugues de Lassus Saint-Geniès Advisors: David Defour and Guillaume Revy

DALI project-team, University of Perpignan Via Domitia LIRMM, University of Montpellier 2 CNRS UMR 5506

18th February 2015

Who I am in 30”

Paris

→ 2008: bac S 2008 → 2011: CPGE

Grenoble

2 1 1

2011 → 2013: Grenoble INP S2 2013: univ. autonome de Basse-Californie

Ensenada, Mexico

August 2013

February 2014: internships (Vérimag, Argosim)

February 2014

October 2014: DALI/UPVD/LIRMM

Perpignan

October 2014

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 1 / 24

1

DALI in MetaLibm

2

Background on Function Evaluation

3

Primitive Pythagorean Triples

4

Experimental Results

5

Conclusion and Perspectives

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 2 / 24

DALI in the MetaLibm Project

Architecture-aware code generation:

◮ Parallelism, vector operations? ◮ FPU? ◮ Memory management, caches? ◮ Features like FMA, if-conversion?

Focus on elementary functions Contributions: high quality code generation chain MetaLibm function interval architecture accuracy

• ptimized

code

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 2 / 24

Basic Scheme for Function Evaluation

x no yes Special input Floating-point number unpacking Normalization Range redution Sign and exponent computation Significand approximation Rounding decision Correct rounding computation Special output selection Result construction r {±0, ±∞, NaN, . . .}

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 3 / 24

Basic Scheme for Function Evaluation

x no yes Special input Floating-point number unpacking Normalization Range redution Sign and exponent computation Significand approximation Rounding decision Correct rounding computation Special output selection Result construction r {±0, ±∞, NaN, . . .}

We focus on: sine and cosine range reduction based on tabulated values combined with polynomial approximations

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 3 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

2 Keeping only values in [0, π/4]:

x ← s · x with s = sign(x)

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

2 Keeping only values in [0, π/4]:

x ← s · x with s = sign(x)

3 Splitting x into high and low parts xh and xℓ

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

2 Keeping only values in [0, π/4]:

x ← s · x with s = sign(x)

3 Splitting x into high and low parts xh and xℓ 4 Table lookup of precomputed sine and cosine of xh

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

2 Keeping only values in [0, π/4]:

x ← s · x with s = sign(x)

3 Splitting x into high and low parts xh and xℓ 4 Table lookup of precomputed sine and cosine of xh 5 Polynomial approximations of sine and cosine of xℓ

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

2 Keeping only values in [0, π/4]:

x ← s · x with s = sign(x)

3 Splitting x into high and low parts xh and xℓ 4 Table lookup of precomputed sine and cosine of xh 5 Polynomial approximations of sine and cosine of xℓ 6 Reconstruction using basic trigonometric relations:

sin(x) = sin(xh) · cos(xℓ) + cos(xh) · sin(xℓ) Then sin(y) depending on value of k ≡ 4 and s.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Basic Scheme for Sine and Cosine Evaluation

Let y be a floating-point number, y ∈ F.

1 Range reduction

x = y − k π 2 , such that x ∈ (−π 4 , π 4 ], k ∈ Z

2 Keeping only values in [0, π/4]:

x ← s · x with s = sign(x)

3 Splitting x into high and low parts xh and xℓ 4 Table lookup of precomputed sine and cosine of xh 5 Polynomial approximations of sine and cosine of xℓ 6 Reconstruction using basic trigonometric relations:

sin(x) = sin(xh) · cos(xℓ) + cos(xh) · sin(xℓ) Then sin(y) depending on value of k ≡ 4 and s.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 4 / 24

Sources of Error in the Classical Method

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ) ⊖ sinh ⊗ sin(xℓ) sin(x) = sinh ⊗ cos(xℓ) ⊕ cosh ⊗ sin(xℓ) cosh sinh cos(xℓ + corr) sin(xℓ + corr) cos(xℓ) sin(xℓ) xh xℓ

Errors appear in: x, inherently by the range reduction

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 5 / 24

Sources of Error in the Classical Method

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ) ⊖ sinh ⊗ sin(xℓ) sin(x) = sinh ⊗ cos(xℓ) ⊕ cosh ⊗ sin(xℓ) cosh sinh cos(xℓ + corr) sin(xℓ + corr) cos(xℓ) sin(xℓ) xh xℓ

Errors appear in: x, inherently by the range reduction sinh and cosh, which are rounded approximations

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 5 / 24

Sources of Error in the Classical Method

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ) ⊖ sinh ⊗ sin(xℓ) sin(x) = sinh ⊗ cos(xℓ) ⊕ cosh ⊗ sin(xℓ) cosh sinh cos(xℓ + corr) sin(xℓ + corr) cos(xℓ) sin(xℓ) xh xℓ

Errors appear in: x, inherently by the range reduction sinh and cosh, which are rounded approximations sin(xℓ) and cos(xℓ), which are computed by polynomial evaluation

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 5 / 24

Sources of Error in the Classical Method

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ) ⊖ sinh ⊗ sin(xℓ) sin(x) = sinh ⊗ cos(xℓ) ⊕ cosh ⊗ sin(xℓ) cosh sinh cos(xℓ + corr) sin(xℓ + corr) cos(xℓ) sin(xℓ) xh xℓ

Errors appear in: x, inherently by the range reduction sinh and cosh, which are rounded approximations sin(xℓ) and cos(xℓ), which are computed by polynomial evaluation in the reconstruction process

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 5 / 24

Sources of Error in the Classical Method

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ) ⊖ sinh ⊗ sin(xℓ) sin(x) = sinh ⊗ cos(xℓ) ⊕ cosh ⊗ sin(xℓ) cosh sinh cos(xℓ + corr) sin(xℓ + corr) cos(xℓ) sin(xℓ) xh xℓ

Errors appear in: x, inherently by the range reduction sinh and cosh, which are rounded approximations sin(xℓ) and cos(xℓ), which are computed by polynomial evaluation in the reconstruction process Drawback: additional bits in computation steps for accuracy purpose

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 5 / 24

Reducing the Error on Tabulated Values1

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ + corr) ⊖ sinh ⊗ sin(xℓ + corr) sin(x) = sinh ⊗ cos(xℓ + corr) ⊕ cosh ⊗ sin(xℓ + corr) cosh sinh cos(xℓ + corr) sin(xℓ + corr) xh xℓ corr

Gal’s accurate table method 1 more source of error: correction term corr reduces the error on sinh, cosh, corr by making them very close to floating-point numbers

1Gal and Bachelis, “An Accurate Elementary Mathematical Library for the IEEE

Floating Point Standard”.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 6 / 24

Our Method for Reducing the Error on Tabulated Values

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ + corr) ⊖ sinh ⊗ sin(xℓ + corr) sin(x) = sinh ⊗ cos(xℓ + corr) ⊕ cosh ⊗ sin(xℓ + corr) cosh sinh cos(xℓ + corr) sin(xℓ + corr) xh xℓ corr

Our method: still 1 more source of error: correction term corr but sinh and cosh stored exactly, that is, without any rounding error

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 7 / 24

Our Method for Reducing the Error on Tabulated Values

x ∈ [0, π

4 ]

Splitting x into (xh, xℓ) Table lookup Polynomial approximations Reconstruction cos(x) = cosh ⊗ cos(xℓ + corr) ⊖ sinh ⊗ sin(xℓ + corr) sin(x) = sinh ⊗ cos(xℓ + corr) ⊕ cosh ⊗ sin(xℓ + corr) cosh sinh cos(xℓ + corr) sin(xℓ + corr) xh xℓ corr

Our method: still 1 more source of error: correction term corr but sinh and cosh stored exactly, that is, without any rounding error How? With Pythagorean triples.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 7 / 24

Pythagorean Triples

Let TP and TPP be the sets of pythagorean triples and prime pythagorean triples, respectively. TP = {(a, b, c) ∈ N∗3 | a2 + b2 = c2} TPP = {(a, b, c) ∈ TP | a, b, c coprime}

α b a c

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 8 / 24

Pythagorean Triples

Let TP and TPP be the sets of pythagorean triples and prime pythagorean triples, respectively. TP = {(a, b, c) ∈ N∗3 | a2 + b2 = c2} TPP = {(a, b, c) ∈ TP | a, b, c coprime}

α b a c

For each pythagorean triple (a, b, c): sin(α) = a/c cos(α) = b/c

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 8 / 24

Pythagorean Triples

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 π/4

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 9 / 24

Pythagorean Triples

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 π/4

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 9 / 24

Pythagorean Triples

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 π/4

corr

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 9 / 24

Pythagorean Triples

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 π/4

corr b a c

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 9 / 24

Primitive Pythagorean Triples Generation

One easy way: ternary-trees with root (3, 4, 5)2 3 linear relationships 3 children/node, e.g:3

  1 −2 2 2 −1 2 2 −2 3   ,   −1 2 2 2 −1 2 −2 2 3   ,   1 2 2 2 1 2 2 2 3  

2Price, “The Pythagorean Tree: A New Species”. 3Barning, “On Pythagorean and quasi-Pythagorean triangles and a generation

process with the help of unimodular matrices”.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 10 / 24

Primitive Pythagorean Triples Generation

One easy way: ternary-trees with root (3, 4, 5)2 3 linear relationships 3 children/node, e.g:3

  1 −2 2 2 −1 2 2 −2 3   ,   −1 2 2 2 −1 2 −2 2 3   ,   1 2 2 2 1 2 2 2 3  

3, 4, 5 5, 12, 13 15, 8, 17 21, 20, 29

2Price, “The Pythagorean Tree: A New Species”. 3Barning, “On Pythagorean and quasi-Pythagorean triangles and a generation

process with the help of unimodular matrices”.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 10 / 24

Primitive Pythagorean Triples with c ≤ 212

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 11 / 24

Primitive Pythagorean Triples with c ≤ 212

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 π/4

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 11 / 24

Primitive Pythagorean Triples with c ≤ 212

500 1000 1500 2000 2500 3000 3500 4000 500 1000 1500 2000 2500 3000 3500 4000 π/4

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 11 / 24

Utility of Primitive Pythagorean Triples: Solution 1

For each table entry: store a, b, c as float or double exactly compute and store corr = xh − arcsin(a/c)

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 12 / 24

Utility of Primitive Pythagorean Triples: Solution 1

For each table entry: store a, b, c as float or double exactly compute and store corr = xh − arcsin(a/c) Reconstruction steps:

1 Table lookup of a, b, c, corr 2 Polynomial approximations

Cℓ ≈ cos(xℓ + corr) and Sℓ ≈ sin(xℓ + corr)

3 Result reconstruction:

sin(x) = (a ⊗ Cℓ ⊕ b ⊗ Sℓ) ⊘ c

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 12 / 24

Utility of Primitive Pythagorean Triples: Solution 1

For each table entry: store a, b, c as float or double exactly compute and store corr = xh − arcsin(a/c) Reconstruction steps:

1 Table lookup of a, b, c, corr 2 Polynomial approximations

Cℓ ≈ cos(xℓ + corr) and Sℓ ≈ sin(xℓ + corr)

3 Result reconstruction:

sin(x) = (a ⊗ Cℓ ⊕ b ⊗ Sℓ) ⊘ c ⊖ Additional cost of 1 division by c slower, more rounding.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 12 / 24

Utility of Primitive Pythagorean Triples: Solution 2

For each table entry: store exactly A = a

c · k and B = b c · k such that A, B ∈ N

include 1

k in polynomial coefficients

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 13 / 24

Utility of Primitive Pythagorean Triples: Solution 2

For each table entry: store exactly A = a

c · k and B = b c · k such that A, B ∈ N

include 1

k in polynomial coefficients

⇒ k should be a small lowest common multiple (lcm) of all possible combinations of one triple hypothenuse (c) per table entry

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 13 / 24

Utility of Primitive Pythagorean Triples: Solution 2

For each table entry: store exactly A = a

c · k and B = b c · k such that A, B ∈ N

include 1

k in polynomial coefficients

⇒ k should be a small lowest common multiple (lcm) of all possible combinations of one triple hypothenuse (c) per table entry Reconstruction steps:

1 Table-lookup of A, B, corr 2 Polynomial approximations

Sℓ ≈ sin(xℓ + corr) k and Cℓ ≈ cos(xℓ + corr) k

3 Result reconstruction:

sin(x) = A ⊗ Cℓ ⊕ B ⊗ Sℓ

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 13 / 24

Table Generation Algorithm

Inputs

◮ p: table index size ◮ n: log2(·) of the bound on considered hypothenuse triples

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 14 / 24

Table Generation Algorithm

Inputs

◮ p: table index size ◮ n: log2(·) of the bound on considered hypothenuse triples

Algorithm steps

1

Compute triples with c ≤ 2n

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 14 / 24

Table Generation Algorithm

Inputs

◮ p: table index size ◮ n: log2(·) of the bound on considered hypothenuse triples

Algorithm steps

1

Compute triples with c ≤ 2n

2

Determine all the possible lcm among computed hypotenuses

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 14 / 24

Table Generation Algorithm

Inputs

◮ p: table index size ◮ n: log2(·) of the bound on considered hypothenuse triples

Algorithm steps

1

Compute triples with c ≤ 2n

2

Determine all the possible lcm among computed hypotenuses

3

Search for a lcm k common to all table entries

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 14 / 24

Table Generation Algorithm

Inputs

◮ p: table index size ◮ n: log2(·) of the bound on considered hypothenuse triples

Algorithm steps

1

Compute triples with c ≤ 2n

2

Determine all the possible lcm among computed hypotenuses

3

Search for a lcm k common to all table entries

4

For each entry, search for a triple with hypothenuse k and minimizing correction term

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 14 / 24

Table Generation Algorithm

Inputs

◮ p: table index size ◮ n: log2(·) of the bound on considered hypothenuse triples

Algorithm steps

1

Compute triples with c ≤ 2n

2

Determine all the possible lcm among computed hypotenuses

3

Search for a lcm k common to all table entries

4

For each entry, search for a triple with hypothenuse k and minimizing correction term

Outputs

◮ table with entries of the form

table[i] = {Ai, Bi, corri}

◮ polynomial approximants of sin(xℓ+corr)

k

and cos(xℓ+corr)

k

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 14 / 24

Experimental Protocol

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first n p

means:

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 15 / 24

Experimental Protocol

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first n p

means:

1 ⌈ π

4 · 2p⌉ table entries

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 15 / 24

Experimental Protocol

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first n p

means:

1 ⌈ π

4 · 2p⌉ table entries

2 search k among generated hypotenuses less than 2n

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 15 / 24

Experimental Protocol

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first n p

means:

1 ⌈ π

4 · 2p⌉ table entries

2 search k among generated hypotenuses less than 2n 3 discard first entry: triple (0, 1, 1) is fine

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 15 / 24

Example for (n, p) = (5, 2)

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first 5 2
• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 16 / 24

Example for (n, p) = (5, 2)

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first 5 2

Output table:

static const double lcm = 25; static const double table [4][3] = { { /* A = */ 0 , /* B = */ 25, /* corr = */ 0 x0000000000000000 }, /* corr ~ 0 */ /* angle ~ 0 */ /* ppt = (0,1,1) */ { /* A = */ 7 , /* B = */ 24, /* corr = */ 0 xbfa14d7623b1c6e8 }, /* corr ~

• 0.0337941 */

/* angle ~ 0.283794 */ /* ppt = (7 ,24 ,25) */ { /* A = */ 15, /* B = */ 20, /* corr = */ 0 xbfc25e3e8c9a7b84 }, /* corr ~

• 0.143501 */

/* angle ~ 0.643501 */ /* ppt = (3,4,5) */ { /* A = */ 20, /* B = */ 15, /* corr = */ 0 xbfc6b19c1586ed40 } /* corr ~

• 0.177295 */

/* angle ~ 0.927295 */ /* ppt = (4,3,5) */ };

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 16 / 24

Example for (n, p) = (5, 2)

./search_lcm --address=1 --algorithm=exhaustive \

• -omit-first 5 2

Output table:

static const double lcm = 25; static const double table [4][3] = { { /* A = */ 0 , /* B = */ 25, /* corr = */ 0 x0000000000000000 }, /* corr ~ 0 */ /* angle ~ 0 */ /* ppt = (0,1,1) */ { /* A = */ 7 , /* B = */ 24, /* corr = */ 0 xbfa14d7623b1c6e8 }, /* corr ~

• 0.0337941 */

/* angle ~ 0.283794 */ /* ppt = (7 ,24 ,25) */ { /* A = */ 15, /* B = */ 20, /* corr = */ 0 xbfc25e3e8c9a7b84 }, /* corr ~

• 0.143501 */

/* angle ~ 0.643501 */ /* ppt = (3,4,5) */ { /* A = */ 20, /* B = */ 15, /* corr = */ 0 xbfc6b19c1586ed40 } /* corr ~

• 0.177295 */

/* angle ~ 0.927295 */ /* ppt = (4,3,5) */ }; 5 10 15 20 25 5 10 15 20 25 Stored PT Selected PPT π/4

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 16 / 24

Exhaustive Search

Intel(R) Xeon(R) CPU E5-2650 v2 @ 2.60 GHz (32 cores) 125 GB RAM p n kmin time (s) kmax t (s) hypotenuses 3 10 725 ≪ 1 725 ≪ 1 125 4 14 10, 625 0.01 13, 325 0.02 1, 712 5 17 130, 645 0.14 130, 645 0.08 12, 463 6 21 1, 676, 285 6 1, 698, 385 3.8 179, 632 7 25 32, 846, 125 1, 000 32, 846, 125 65 2, 635, 323

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 17 / 24

Exhaustive Search

Intel(R) Xeon(R) CPU E5-2650 v2 @ 2.60 GHz (32 cores) 125 GB RAM p n kmin time (s) kmax t (s) hypotenuses 3 10 725 ≪ 1 725 ≪ 1 125 4 14 10, 625 0.01 13, 325 0.02 1, 712 5 17 130, 645 0.14 130, 645 0.08 12, 463 6 21 1, 676, 285 6 1, 698, 385 3.8 179, 632 7 25 32, 846, 125 1, 000 32, 846, 125 65 2, 635, 323 Current bottlenecks: searching the lcm and memory Finding k for a 10-bit addressed table is desperate with this exhaustive technique

◮ our estimations show that n ≈ 37

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 17 / 24

Focus on Exhaustive Search for (n, p) = (25, 7)

[Sun Feb 15 11:04:05 2015] Start Pythagorean -triple -based table generation tool *** table address size (# bits): 7 *** maximum hypothenuse size (# bits): 25 ***

• mit

first table cell: yes *** algorithm: exhaustive search *** addressing mode: with

• ffset

*** max angle: 0.789062 rad *** entry width: 0.0078125 rad ***

• utput: cout

[Sun Feb 15 11:04:06 2015] Triple generation

• n 25 bits ...

done. [Sun Feb 15 11:04:12 2015] Triple sort by angle division ... done. *** 101 angle divisions are used [Sun Feb 15 11:04:22 2015] Lowest common multiple computation ... done. *** computed lcm: 32846125 [Sun Feb 15 11:05:09 2015] Triple table building ... done. [Sun Feb 15 11:05:10 2015] Table printing

• ut ...

done [Sun Feb 15 11:05:10 2015] End Pythagorean -triple -based table generation tool

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 18 / 24

Optimization of LCM Searching

For kmin: search lcm in [2n−1, 2n] algorithm called big_c

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 19 / 24

Optimization of LCM Searching

For kmin: search lcm in [2n−1, 2n] algorithm called big_c p n kmin time (s) hypotenuses 5 17 130, 645 0.05 6, 044 6 21 1, 676, 285 4.5 87, 633 7 25 32, 846, 125 670 1, 290, 663 ⊕ ≈ twice less hypotenuses ⊖ still a bottleneck: too many triples/hypotenuses!

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 19 / 24

First Heuristic

p n kmin prime factorization triples 5 17 130, 645 5 · 17 · 29 · 53 21, 588 6 21 1, 676, 285 5 · 13 · 17 · 37 · 41 338, 660 7 25 32, 846, 125 53 · 13 · 17 · 29 · 41 5, 365, 290 For hypotenuses in [2n−1, 2n], only keep those ending by 5

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 20 / 24

First Heuristic

p n kmin prime factorization triples 5 17 130, 645 5 · 17 · 29 · 53 21, 588 6 21 1, 676, 285 5 · 13 · 17 · 37 · 41 338, 660 7 25 32, 846, 125 53 · 13 · 17 · 29 · 41 5, 365, 290 For hypotenuses in [2n−1, 2n], only keep those ending by 5 p n kmin time (s) triples hypotenuses 6 21 1, 676, 285 1.6 225, 788 18, 505 7 25 32, 846, 125 213 3, 576, 819 270, 360 Roughly: ×3 speedup compared to big_c

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 20 / 24

Second heuristic

Observation: k is always a product of small prime factors of the form 4n + 1 Pythagorean primes4 Pythagorean primes less than 70: p = {5, 13, 17, 29, 37, 41, 53, 61}

4Fässler, “Multiple Pythagorean Number Triples”.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 21 / 24

Second heuristic

Observation: k is always a product of small prime factors of the form 4n + 1 Pythagorean primes4 Pythagorean primes less than 70: p = {5, 13, 17, 29, 37, 41, 53, 61} During the generation, only keep triples with hypotenuse c =

• p[i]ri · Q

with Q < 5 and ri = 0

• r

1 if p[i] = 5

4Fässler, “Multiple Pythagorean Number Triples”.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 21 / 24

Second heuristic

Observation: k is always a product of small prime factors of the form 4n + 1 Pythagorean primes4 Pythagorean primes less than 70: p = {5, 13, 17, 29, 37, 41, 53, 61} During the generation, only keep triples with hypotenuse c =

• p[i]ri · Q

with Q < 5 and ri = 0

• r

1 if p[i] = 5 n p kmin time (s) triples hypotenuses 21 6 1,676,285 0.09 1,566 43 25 7 32,846,125 0.8 3,308 48 28 8 243,061,325 4.3 5,217 49

4Fässler, “Multiple Pythagorean Number Triples”.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 21 / 24

Limit of our Implementation

With such results. . . why not try with p = 9?

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 22 / 24

Limit of our Implementation

With such results. . . why not try with p = 9? For p = 9, n ≈ 33 from our estimations.

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 22 / 24

Limit of our Implementation

With such results. . . why not try with p = 9? For p = 9, n ≈ 33 from our estimations. ./search_lcm --address=1 --omit-first \

• -algorithm=big_c 33 9

n p kmin time (s) triples hypotenuses 33 9 SIGSEGV – – –

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 22 / 24

Limit of our Implementation

With such results. . . why not try with p = 9? For p = 9, n ≈ 33 from our estimations. ./search_lcm --address=1 --omit-first \

• -algorithm=big_c 33 9

n p kmin time (s) triples hypotenuses 33 9 SIGSEGV – – – Limit: recursion depth in generation function for p = 9, it was 40,310 when program crashed

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 22 / 24

Conclusion

A new approach using primitive pythagorean triples for sine and cosine evaluation:

⊲ removing the error on two tabulated terms ⊲ concentrating the error on half the terms in the reconstruction

A prototype able to compute:

⊲ size-variable tables up to 8 indexing bits ⊲ size-matching approximation polynomials with Sollya

• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 23 / 24

Perspectives

Generation of tables addressed by more than 8 bits Evaluation of the actual accuracy and performance of sine and cosine implementations using this approach Integration in a full code-generation chain Number theory might help us: ignore bad triples directly build hypotenuses composed of Pythagorean primes select the best triple for each table entry, reducing the lcm-search to

• ne computation
• H. de Lassus Saint-Geniès

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation 24 / 24

Pythagorean-Triple-Based Table Generation for Trigonometric Function Evaluation

MetaLibm Project Annual Meeting’2015 Hugues de Lassus Saint-Geniès Advisors: David Defour and Guillaume Revy

DALI project-team, University of Perpignan Via Domitia LIRMM, University of Montpellier 2 CNRS UMR 5506

18th February 2015

Perspectives

Experiment shows that for p = 10, n needs only be 20 to get at least

• ne triple/entry (forgetting first entry).

50 100 150 200 250 300 350 400 100 200 300 400 500 600 700 800 Number of PPT wrt entry number

⊖ Still, computing ≈ 200800 combinations is not yet foreseen. . .