Prosthaphaeresis - an ancient algorithm for approximating products

The year is 1601. You and your buddies have invested heavily and bought a new sleek ship which you all agreed deserved the name Helena. You are all pumped to sail the world, to meet new people and to see new places.

Having prepared everything, you are on the final step of doing calculations to guarantee the safe navigation of Helena and you come down to the expression $ 197598759 \times 84738549 $ that needs to be evaluated. Unfortunately, you all really suck at multiplication, logarithms are yet to be invented and you’re stuck. Oh, no! What do you do?

The Solution

Now, if you may forgive the rather fictatious depiction of you and your buddies, this was interestingly a problem that sailors and mathematicians faced pre-logarithms invention – that of estimating products. To solve their problems, they came up with an ingenious way of relating products to trigonometric identities, known as prosthaphaeresis.

We recently covered an article on how to derive trigonometric identities and in it, we came across the following identities,

$$ \tag{1} \sin(a)\cos(b) = \frac{\sin(a+b) + \sin(a-b)}{2} $$ $$ \tag{2} \cos(a)\cos(b) = \frac{\cos(a+b) + \cos(a-b)}{2} $$ $$ \tag{3} \sin(a)\sin(b) = \frac{\cos(a-b) - \cos(a+b)}{2} $$

Clearly, this may serve as a way to reduce difficult multiplication problems to calculating sums, differences and the easy division by two if we can tame any values $ x, y $ to be constrained in the domain $ -1 \leq x, y \leq 1 $ (the domain of cosine and sine trigonometric functions).

Luckily, we can always rewrite any numbers in the scientific notation but this time with a twist – we force the coefficient to be within the range of $ -1 $ and $ 1 $. Using the numbers we came across earlier, we have,

$$ \begin{aligned} 197598759 \times 84738549 \; & = \; (1.97598759 \times 10^{8}) \times (8.4738549 \times 10^{7}) \\ \; & = \; 1.97598759 \times 8.4738549 \times 10^{15} \\ \; & = \; 0.197598759 \times 0.84738549 \times 10^{17} \end{aligned} $$

Now the problem comes down to calculating $ 0.197598759 \times 0.84738549 $ and without loss of generality, we can relate the problem to the identity $ (2) $ above. Therefore, we let,

$$ \begin{aligned} \cos(a) & = 0.197598759 \\ \; \Rightarrow a & = \cos^{-1}(0.197598759) = 1.371888552 \\ \\ \; \cos(b) & = 0.84738549 \\ \; \Rightarrow b & = \cos^{-1}(0.84738549) = 0.559754503 \end{aligned} $$

With $ a $ and $ b $ in place, the problem is now reduced to evaluating,

$$ \begin{aligned} \frac{\cos(a+b) + \cos(a-b)}{2} & = \frac{\cos(1.931643055) + \cos(0.812134049)}{2} \\ \; & = 0.16744232131184622 \end{aligned} $$

Which on multiplying by $ 10^{17} $ we get $ 16744232131184622 $ which is a rather close approximation of the actual product!

The Implementation

To those of us who would rather test the algorithm out on their own, here’s an implementation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34


import math

# Convert a value x to the notation a * 10^b
# where  a <= 1. Note that we require x > -1;
def scale_coefficient(x: float) -> (float, float):
    assert x > -1
    exponent = 0
    while x > 1:
        x /= 10
        exponent += 1
    return (x, exponent)

def prosthaphaeresis(a: float, b: float) -> float:
    sign = 1
    if a < 0:
        a = math.abs(a)
        sign = -1
    if b < 0:
        b = math.abs(b)
        sign *= -1
    a_scaled, a_exp = scale_coefficient(a)
    b_scaled, b_exp = scale_coefficient(b)
    a_cos_inverse = math.acos(a_scaled)
    b_cos_inverse = math.acos(b_scaled)
    return (
        sign * 
        (math.cos(a_cos_inverse + b_cos_inverse) + 
         math.cos(a_cos_inverse - b_cos_inverse)) / 
        2.0 *
        math.pow(10, a_exp + b_exp)
    )

if __name__ == '__main__':
    print(prosthaphaeresis(197598759, 84738549))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43


#include <cassert>
#include <cmath>
#include <iostream>

// Convert a value x to the notation a * 10^b
// where  a <= 1. Note that we require x > -1;
std::pair<double, unsigned int> scaleCoefficient(double x) {
  assert(x > -1);
  int exponent = 0;
  while (x > 1) {
    x /= 10;
    exponent += 1;
  }
  return std::make_pair(x, exponent);
}

double prosthaphaeresis(double a, double b) {
  int sign = 1;
  if (a < 0) {
    a = std::abs(a);
    sign = -1;
  }
  if (b < 0) {
    b = std::abs(b);
    sign *= -1;
  }
  auto aScaled = scaleCoefficient(a);
  auto bScaled = scaleCoefficient(b);
  auto aCosInverse = std::acos(aScaled.first);
  auto bCosInverse = std::acos(bScaled.first);
  return (
    sign *
    (std::cos(aCosInverse + bCosInverse) +
     std::cos(aCosInverse - bCosInverse)) / 
    2.0 * 
    std::pow(10, aScaled.second + bScaled.second)
  );
}

int main() {
  std::cout << prosthaphaeresis(197598759, 84738549);
  return 0;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59


public class Main {
  private class Pair<K, V> {
    private K key;
    private V value;

    Pair(K key, V value) {
      this.key = key;
      this.value = value;
    }

    public K getKey() {
      return this.key;
    }

    public V getValue() {
      return this.value;
    }
  }
  
  // Convert a value x to the notation a * 10^b
  // where  a <= 1. Note that we require x > -1;
  private Pair<Double, Integer> scaleCoefficient(Double x) {
    assert x > -1;
    int exponent = 0;
    while (x > 1) {
      x /= 10;
      ++exponent;
    }
    return new Pair<>(x, exponent);
  }

  private double prosthaphaeresis(double a, double b) {
    int sign = 1;
    if (a < 0) {
      a = Math.abs(a);
      sign = -1;
    }
    if (b < 0) {
      b = Math.abs(b);
      sign *= -1;
    }
    Pair<Double, Integer> aScaled = scaleCoefficient(a);
    Pair<Double, Integer> bScaled = scaleCoefficient(b);
    double aCosInverse = Math.acos(aScaled.getKey());
    double bCosInverse = Math.acos(bScaled.getKey());
    return (
      sign *
      (Math.cos(aCosInverse + bCosInverse) + 
       Math.cos(aCosInverse - bCosInverse)) /
      2.0 *
      Math.pow(10, aScaled.getValue() + bScaled.getValue())
    );
  }

  public static void main (String[] args) {
    Main runner = new Main();
    System.out.println(runner.prosthaphaeresis(197598759, 84738549));
  }
}