<commit_before>---
title: 'Recursive Doubling'
date: 2019-10-03
permalink: /posts/2019/10/03/recursive-doubling
categories:
  - math
tags:
  - computer science
  - math
  - parallel computing
---


Consider a problem of the form 

$$ x_i = a_i x_{i-1} + b_i $$

where $x_1=b_1$ and $a_i,b_i,x_i \in \mathbb{R}$. When $a$ and $b$ are known computing $x$ appears to be an intrinsically serial problem. To compute $x_i$ you must know $x_{i-1}$. The code would look something like below.

```c++
x[0] = b[0];
for (size_t i = 1; i < n; i++) {
  x[i] = a[i] * x[i-1] + b[i];
}
```

However, we can approach in a parallel manner when we consider some observations.

Observation #1: $x_i$ can be rewritten in terms of $x_{i-2}$. 

$$x_i = a_i \left( x_{i-1} \right) + b_i = a_i \left( a_{i-1} x_{i-2} + b_{i-1} \right) + b_i $$

$$ = a_i a_{i-1} x_{i-2} + a_i b_{i-1} + b_i  $$


Just using Observation #1 we can already see room for parallelization. The odd and even indices of $x$ can be computed in parallel. Using OpenMP here is the code using Observation #1 alone.

``` c++
#pragma omp parallel {
  int thread_idx = omp_get_thread_num();
  if (thread_idx == 0) {
    for (size_t i = 1; i < n; i += 2) {
      x[i] = a[i] * a[i-1] * x[i-2] + a[i] * b[i-1] + b[i];
    }
  } else if (thread_idx == 1) {<commit_msg>fixed error<commit_after>---
title: 'Recursive Doubling'
date: 2019-10-03
permalink: /posts/2019/10/03/recursive-doubling
categories:
  - math
tags:
  - computer science
  - math
  - parallel computing
---


Consider a problem of the form 

$$ x_i = a_i x_{i-1} + b_i $$

where $x_1=b_1$ and $a_i,b_i,x_i \in \mathbb{R}$. When $a$ and $b$ are known computing $x$ appears to be an intrinsically serial problem. To compute $x_i$ you must know $x_{i-1}$. The code would look something like below.

```c++
x[0] = b[0];
for (size_t i = 1; i < n; i++) {
  x[i] = a[i] * x[i-1] + b[i];
}
```

However, we can approach in a parallel manner when we consider some observations.

Observation #1: $x_i$ can be rewritten in terms of $x_{i-2}$. 

$$x_i = a_i \left( x_{i-1} \right) + b_i = a_i \left( a_{i-1} x_{i-2} + b_{i-1} \right) + b_i $$

$$ = a_i a_{i-1} x_{i-2} + a_i b_{i-1} + b_i  $$


Just using Observation #1 we can already see room for parallelization. The odd and even indices of $x$ can be computed in parallel. Using OpenMP here is the code using Observation #1 alone.

``` c++
x[0] = b[0];
x[1] = a[1] * x[0] + b[1];
#pragma omp parallel {
  int thread_idx = omp_get_thread_num();
  if (thread_idx == 0) {
    for (size_t i = 3; i < n; i += 2) {
      x[i] = a[i] * a[i-1] * x[i-2] + a[i] * b[i-1] + b[i];
    }
  } else if (thread_idx == 1) {
{'fraction_non_alphanumeric': 0.17394578313253012, 'fraction_numerical': 0.0338855421686747, 'mean_word_length': 3.1003086419753085, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3385075', 'n_tokens_mistral': 1168, 'n_tokens_neox': 1109, 'n_words': 365}