more stable implementation of the Law of Cosines.

lawofcosine.py

import numpy as np

method_list = [1, 2, 3]

for method in method_list:
    for i in range(2000):
        length = 10000
        a = np.random.randn(length).astype(np.float32)
        b = np.random.randn(length).astype(np.float32)
        # non-stable when theta is small
        theta = (np.random.rand(length)*2*np.pi - np.pi).astype(np.float32) * 0.0001

        if method == 1:
            # non-stable using float32
            # you may get runtime warnning
            c = np.sqrt(a**2 + b**2 - 2*a*b*np.cos(theta))
        elif method == 2:
            # more stable1: use float64
            a, b, theta = a.astype(np.float64), b.astype(np.float64), theta.astype(np.float64)
            c = np.sqrt(a**2 + b**2 - 2*a*b*np.cos(theta))
        elif method == 3:
            # more stable2: limit a and b between 0 and 1
            max_val = np.maximum(a, b)
            a /= max_val
            b /= max_val
            c = np.sqrt(a**2 + b**2 - 2*a*b*np.cos(theta))
            c *= max_val