The Accuracy of l3fp: a Series of Test Cases

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:

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:

  1. test_fp_sum: In each run, the sum of 300 random floating point numbers is computed.
  2. test_fp_mult: In each run, the product of 300 random floating point numbers is computed.
  3. test_fp_div: In each run, the quotient of a number divided by 300 random floating point numbers is computed.
  4. 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

Details

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')