The Accuracy of l3fp: a Series of Test Cases

03 Aug 2020

Representation of floating point numbers in l3fp
Test cases & results
Conclusion
Details

I was reading LaTeX3’s documentation for \dim_to_fp:n, where I encountered the following description:

Expands to an internal floating point number equal to the value of the <dimexpr> in pt. Since dimension expressions are evaluated much faster than their floating point equivalent, \dim_to_fp:n can be used to speed up parts of a computation where a low precision and a smaller range are acceptable.

At first, I was confused by the subject of this paragraph and thought the precision of l3fp is not ideal. Because many data processing packages (e.g. datatool) rely on l3fp to do floating point arithmetic, I had this idea of testing l3fp’s accuracy against IEEE 754 double precision floating point arithmetic. The result shows that the error of l3fp is very small compared to IEEE 754 and is negligible in everyday applications. However, the trigonometry functions of l3fp seems to be significantly less precise compared to other operations.

Representation of floating point numbers in l3fp

As this post points out, the secret of l3fp’s floating point representation can be exposed by simply showing the meaning of floating point macros. For example,

\fp_set:Nn \l_tmpa_fp {pi}
\fp_use:N \l_tmpa_fp
\par\cs_meaning:N \l_tmpa_fp

gives:

3.141592653589793
macro:->\s__fp \__fp_chk:w 10{1}{3141}{5926}{5358}{9793};

It can be seen that:

Different from IEEE 754 (binary), the number is represented in decimal format.
It has 16 significant digits. It is widely known that IEEE 754 double precision gives 15-17 significant digits. (source)

Test cases & results

I have created a series of test cases which are inspired by daily application scenarios. They do not evaluate the performance of l3fp rigorously. My intentions are to compare their performances in everyday scenarios and prove their reliability. There are 4 types of floating point arithmetic tests:

test_fp_sum: In each run, the sum of 300 random floating point numbers is computed.
test_fp_mult: In each run, the product of 300 random floating point numbers is computed.
test_fp_div: In each run, the quotient of a number divided by 300 random floating point numbers is computed.
test_fp_(func)_sum: In each run, apply (func) to 300 random floating point numbers and compute their sum.

Each test is repeated 200 times, then the relative differece w.r.t. to IEEE 754 double precision baseline is computed. The result is shown in the following box plot:

result

Conclusion

Compared to IEEE 754, l3fp is actually very precise: the relative difference never exceeds \(10^{-13}\) across all test cases. This can be considered insignificant for most users. For most functions except for trigonometry, the relative difference is actually much smaller than this value. It can be concluded that we can comfortably use l3fp in most everyday scenarios.
The paragraph at the opening of this post actually means that dimension expressions in LaTeX are faster but less accurate, which inspires the design idea of l3fparray.

Details

TeX distribution: TexLive 2020 on Windows
Python version: Python 3.7 on Windows

Python main test utility

import numpy as np
import subprocess
import os
import pickle
import string
import random
from multiprocessing.pool import ThreadPool

with open('tex_template.tex', 'r') as infile:
    tex_template = infile.read()

def get_random_name():
    k = 15
    return ''.join(random.choices(string.ascii_letters + string.digits, k=k))

class TeXExecutor:

    def __init__(self):
        self.tex_infile = get_random_name() + '.txt'
        self.tex_outfile = get_random_name() + '.txt'
        self.tex_source_prefix = get_random_name()
        self.tex_source = self.tex_source_prefix + '.tex'

    def delete_files(self):
        self._del_if_exists(self.tex_infile)
        self._del_if_exists(self.tex_outfile)
        self._del_if_exists(self.tex_source)
        self._del_if_exists(self.tex_source_prefix + '.aux')
        self._del_if_exists(self.tex_source_prefix + '.log')

    def _del_if_exists(self, fn):
        if os.path.exists(fn):
            os.remove(fn)

    def execute_tex_code(self, tex_src):
        tex_src = tex_src.replace('%%%%infile', self.tex_infile).replace('%%%%outfile', self.tex_outfile)
        with open(self.tex_source, 'w') as outfile:
            outfile.write(tex_src)
        subprocess.run(['pdflatex', '-interaction=nonstopmode', self.tex_source], stdout=subprocess.DEVNULL)

    def load_tex_result(self):
        assert os.path.exists(self.tex_outfile)
        with open(self.tex_outfile, 'r') as infile:
            num = float(infile.read().strip())
            return num

    def save_numbers_to_txt(self, arr):
        with open(self.tex_infile, 'w') as outfile:
            for num in arr:
                outfile.write('{:.18e}\n'.format(num))

def test_fp_sum(seed=0):
    np.random.seed(seed)
    numarr = np.concatenate(
        [np.random.uniform(-10000.0, 10000.0, 300)]
    )
    np.random.shuffle(numarr)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_sum = np.sum(numarr)

    initcode = r'''
    \fp_zero:N \l_result_fp
    '''
    loopcode = r'''
    \fp_add:Nn \l_result_fp {\l_tmpa_str}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_sum = tex_ex.load_tex_result()
    #print(actual_sum, tex_sum)
    tex_ex.delete_files()
    return (actual_sum, tex_sum)

def test_fp_mult(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(1.01, 1.1, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_prod = np.prod(numarr)

    initcode = r'''
    \fp_set:Nn \l_result_fp {1.0}
    '''
    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_set:Nn \l_result_fp {\fp_eval:n {\l_result_fp * \l_tmpa_fp}}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_prod = tex_ex.load_tex_result()
    #print(actual_prod, tex_prod)
    tex_ex.delete_files()
    return (actual_prod, tex_prod)


def test_fp_div(seed=0):
    np.random.seed(seed)
    start_num = 123456789.0
    numarr = np.random.uniform(1.01, 1.1, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_quot = start_num
    for num in numarr:
        actual_quot /= num

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % start_num

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_set:Nn \l_result_fp {\fp_eval:n {\l_result_fp / \l_tmpa_fp}}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_quot = tex_ex.load_tex_result()
    tex_ex.delete_files()
    #print(actual_quot, tex_quot)
    return (actual_quot, tex_quot)


def test_fp_sqrt_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(0.0, 10000.0, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.sqrt(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { sqrt(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    #print(actual_val, tex_val)
    return (actual_val, tex_val)

def test_fp_sin_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(-180.0, 180.0, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.sin(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { sin(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    print(actual_val, tex_val)
    return (actual_val, tex_val)

def test_fp_log_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(1.0, 1.0e9, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.log(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { ln(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    print(actual_val, tex_val)
    return (actual_val, tex_val)

def test_fp_tan_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(-np.pi / 2.0 + 0.01, np.pi / 2.0 - 0.01, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.tan(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { tan(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    print(actual_val, tex_val)
    return (actual_val, tex_val)

def test_fp_asin_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(-1.0, 1.0, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.arcsin(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { asin(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    print(actual_val, tex_val)
    return (actual_val, tex_val)

def test_fp_atan_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(-1.0e4, 1.0e4, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.arctan(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { atan(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    print(actual_val, tex_val)
    return (actual_val, tex_val)

def test_fp_exp_sum(seed=0):
    np.random.seed(seed)
    numarr = np.random.uniform(0.1, 4.0, 300)
    tex_ex = TeXExecutor()
    tex_ex.save_numbers_to_txt(numarr)
    actual_val = np.sum(np.exp(numarr))

    initcode = r'''
    \fp_set:Nn \l_result_fp {%.18e}
    ''' % (0, )

    loopcode = r'''
    \fp_set:Nn \l_tmpa_fp {\l_tmpa_str}
    \fp_add:Nn \l_result_fp {\fp_eval:n { exp(\l_tmpa_fp) }}
    '''

    # calculate tex sum
    tex_code = tex_template.replace('%%%%init', initcode).replace('%%%%loopbody', loopcode)
    tex_ex.execute_tex_code(tex_code)

    tex_val = tex_ex.load_tex_result()
    tex_ex.delete_files()
    print(actual_val, tex_val)
    return (actual_val, tex_val)

def run_test_case(test_func):
    print('testing {}...'.format(test_func.__name__))

    n_repeat = 200
    seeds = list(range(n_repeat))

    pool = ThreadPool(15)
    result = pool.map(test_func, seeds, chunksize=5)

    if os.path.exists('test_log.bin'):
        # retrieve result dict
        with open('test_log.bin', 'rb') as infile:
            log = pickle.load(infile)
    else:
        log = dict()

    log[test_func.__name__] = result

    with open('test_log.bin', 'wb') as outfile:
        pickle.dump(log, outfile)
        outfile.flush()


#run_test_case(test_func=test_fp_sum)
#run_test_case(test_func=test_fp_mult)
#run_test_case(test_func=test_fp_div)
#run_test_case(test_func=test_fp_sqrt_sum)
#run_test_case(test_func=test_fp_sin_sum)
#run_test_case(test_func=test_fp_log_sum)
#run_test_case(test_func=test_fp_tan_sum)
#run_test_case(test_func=test_fp_exp_sum)
#run_test_case(test_func=test_fp_asin_sum)
#run_test_case(test_func=test_fp_atan_sum)

LaTeX template file (template.tex)

\documentclass{article}
\usepackage{expl3}

\begin{document}

\ExplSyntaxOn

\fp_new:N \l_result_fp

\ior_open:Nn \g_tmpa_ior {
%%%%infile
}

%%%%init

\ior_str_map_variable:NNn \g_tmpa_ior \l_tmpa_str {
%%%%loopbody
}

\iow_open:Nn \g_tmpa_iow {
%%%%outfile
}
\iow_now:Nx \g_tmpa_iow {\fp_use:N \l_result_fp}

\end{document}

Python plotting

import numpy as np
import matplotlib.pyplot as plt
import pickle

with open('test_log.bin', 'rb') as infile:
    log = pickle.load(infile)

label_conv = {
    'test_fp_sum' : 'Sum',
    'test_fp_mult' : 'Mult',
    'test_fp_div' : 'Div',
    'test_fp_sqrt_sum' : 'Sqrt-Sum',
    'test_fp_sin_sum' : 'Sin-Sum',
    'test_fp_log_sum' : 'Log-Sum',
    'test_fp_tan_sum' : 'Tan-Sum',
    'test_fp_exp_sum' : 'Exp-Sum',
    'test_fp_asin_sum' : 'ArcSin-Sum',
    'test_fp_atan_sum' : 'ArcTan-Sum'
}
labels = []
data = []

fig = plt.figure(figsize=(6, 10))
plt.suptitle('Relative Difference of $l3fp$')

for key, val in log.items():
    val = np.asarray(val)
    labels.append(label_conv[key])
    if val.shape[1] == 2:
        rel_diff = np.abs(val[:, 0] - val[:, 1]) / val[:, 0]
        data.append(rel_diff)

plt.boxplot(data, sym='')
plt.ylabel('Relative difference w.r.t. to IEEE 754 double precision baseline')
plt.xlabel('Tests')
plt.xticks(np.arange(1, 1 + len(labels)), labels, rotation=45)
plt.savefig('result.svg')

Alan Xiang's Blog