A C library for fast multiprecision evaluation of transcendental functions, using fixed-point arithmetic directly on top of GMP/MPIR. The goal is to provide optimized base implementations of transcendental functions, with minimal overhead at low precision as well as asymptotic speed. It should be usable directly for computations that work well in fixed point (i.e. involving both inputs and outputs of order of magnitude close to unity), or as a base library for implementation of floating-point arithmetic with bells and whistles.

Current status

The current target is to test various algorithms, mainly for elementary functions. There is no code with a usable library interface yet. The svn trunk contains test implementations for exp (the same code may also compute cos or sin), log, and gamma.

Exp

The test program computes exp(0.37) and compares both accuracy and speed to MPFR-2.4.1. This implementation is uniformly faster than MPFR up to at least several hundred bits of precision on one AMD64 system. Accuracy is within 1-2 bits of the correctly rounded value (as computed by MPFR). Times (5th and 6th columns) are in nanoseconds.

fredrik@airy:~/src/fastfunlib/fastfunlib/exptest$ ./exptest
 prec   acc   J   r     mpfr     this   faster
   53    52   1   2     5070     1240   4.089
   66    65   1   6     5200     1390   3.741
   82    81   1   8     5590     1460   3.829
  102   101   1   6     5810     1670   3.479
  127   126   2   8     7170     2280   3.145
  158   157   2   7     8390     2500   3.356
  197   196   1  11     8610     3030   2.842
  246   246   1   7     9870     3900   2.531
  307   307   2  12    12260     4440   2.761
  383   381   3   9    12860     5220   2.464
  478   477   2  11    16260     6580   2.471
  597   596   2  14    19220     8039   2.391
  746   746   3  13    23400    10500   2.229
  932   931   3   8    27900    14100   1.979
 1165  1164   3  14    35400    18700   1.893
 1456  1455   3  14    48500    27000   1.796
 1820  1819   3  12    67000    38000   1.763
 2275  2274   4  11    93000    56500   1.646
 2843  2842   4  17   130500    84000   1.554
 3553  3553   5  14   190500   124000   1.536
 4441  4440   6  14   278000   187000   1.487
 5551  5550   6  17   400500   279500   1.433
 6938  6937   6  19   599000   429500   1.395
 8672  8670   7  18   891000   645000   1.381
10840 10839   7  25  1310000   973000   1.346
13550 13549   8  23  2024500  1472000   1.375
16937 16935   9  25  3121000  2299000   1.358
21171 21170   9  27  4647000  3446000   1.349
26463 26461   9  36  5335000  5209500   1.024

The J and r parameters are algorithmic tuning parameters. The test program performs an exhaustive search and reports the best combination. Eventually a program for pre-tuning will have to be written.

This is still a naive implementation and optimistically the speed at low precision could be improved by perhaps an additional factor two.

Log

Logarithm using Taylor series and argument reduction is twice as fast as MPFR up to a few hundred bits:

fredrik@airy:~/src/fastfunlib/fastfunlib/logtest$ ./logtest
 prec   acc   J   r     mpfr     this   faster
   53    53   1   1     4760     1750   2.720
   66    66   2   1     4800     1960   2.449
   82    82   1   3     5370     2390   2.247
  102   102   2   3     6630     2560   2.590
  127   127   3   2     8150     3310   2.462
  158   158   2   2     8160     3960   2.061
  197   197   2   2     9530     4760   2.002
  246   246   4   2    13430     6150   2.184
  307   307   4   2    15180     7300   2.079
  383   383   4   3    18120     9180   1.974
  478   478   4   3    23040    11340   2.032
  597   597   5   3    27760    14820   1.873
  746   746   6   4    30900    18900   1.635
  932   932   5   4    39800    25100   1.586
 1165  1165   5   5    49400    34600   1.428
 1456  1456   5   5    59500    46000   1.293
 1820  1820   5   7    79500    65000   1.223
 2275  2275   6   7   108000    92000   1.174
 2843  2843   7   9   133000   131500   1.011
 3553  3552   6  11   189500   188000   1.008
 4441  4441   7  11   253500   278500   0.910
 5551  5551   7  12   350500   409000   0.857

Logarithm using a 256KB (2^9 values x 4096 bits) precomputed lookup table buys another 2x speedup at low precision:

fredrik@airy:~/src/fastfunlib/fastfunlib/logtest2$ ./logtest2
 prec   acc   J   r     mpfr     this   faster
   53    53   1   0     4710      900   5.233
   66    66   1   0     4770      930   5.129
   82    82   1   0     5700     1040   5.481
  102   102   1   0     6660     1200   5.550
  127   127   1   0     7800     1430   5.455
  158   158   2   0     8660     1750   4.949
  197   197   1   0     9280     2170   4.276
  246   246   2   0    13480     2670   5.049
  307   307   2   0    15020     3120   4.814
  383   383   3   0    18100     4080   4.436
  478   478   3   0    23640     5120   4.617
  597   597   3   0    28620     6880   4.160
  746   746   4   0    31000     9400   3.298
  932   932   4   0    39800    12200   3.262
 1165  1165   4   0    50200    17600   2.852
 1456  1456   4   0    59000    25500   2.314
 1820  1820   4   0    80000    38500   2.078
 2275  2275   6   0   109000    56500   1.929
 2843  2843   6   0   134500    87500   1.537
 3553  3553   6   0   191000   133000   1.436
 4441  4441   7  11   254500   278500   0.914
 5551  5550   7  12   352000   410500   0.857

Gamma function

The gammatest program implements the gamma function using Taylor series, designed for small arguments, e.g. x < 100 (for very large x, a second implementation based on Stirling's series should be used). This approach pays off greatly, giving 20-100x speedups over MPFR.

Times are for computing gamma(5.7) (Columns: prec, accurate, MPFR time, this time, speedup.)

fredrik@airy:~/src/fastfunlib/fastfunlib/gammatest$ ./gammatest
   53    53      76350       3150   24.238
   66    66      83210       3530   23.572
   82    82      90030       3770   23.881
  102   102     108540       4220   25.720
  127   127     121530       5850   20.774
  158   158     140270       6920   20.270
  197   197     174500       9090   19.197
  246   246     238250      12100   19.690
  307   307     309100      15800   19.563
  383   383     409400      19720   20.761
  478   478     604360      25260   23.926
  597   597     847160      33280   25.456
  746   746    1292500      48000   26.927
  932   932    1872700      66200   28.289
 1165  1165    3168700     111200   28.496
 1456  1456    5332000     171000   31.181
 1820  1820    9475000     288000   32.899
 2275  2275   16993000     474000   35.850
 2843  2843   30980500     809000   38.295
 3553  3553   56078000    1362000   41.173
 4441  4441  109397000    2362500   46.306
 5551  5551  219725000    3975500   55.270
 6938  6938  416183500    6976500   59.655
 8672  8672  835161500   11983500   69.693
10840 10840 1652142000   20552500   80.386
13550 13550 3450440500   34704000   99.425

Generating the coefficients for 1000-digit precision from scratch should take about a second (a small fraction of a second for lower precisions). Currently, the test code must read pre-generated coefficients from a file. The svn repository includes a data file generated with mpmath that works up to about 15000 bits.

There should be an interface for loading and saving coefficient data to files between sessions. The fastfunlib project could also provide a downloadable repository of expensive data, so that not every user who needs many extremely high-precision function evaluations has to generate the same data from scratch (which could take hours).

Other functions

Initially exp, cosh/sinh, cos/sin, log, atan, zeta constants/Bernoulli numbers, gamma and digamma functions should be added; more functions may be added later on.