dhp - Dalsoft High Precision package allows to perform floating point calculations with the precision higher than that supported by most of the numeric hardware. It provides C++ wrapper which allows for easy incorporation with the existing code.

dhp may be particularly useful in number of ways:
  1. Technical specifications
  2. Installation
  3. High Precision values
  4. Compiling and linking with the dhp library
  5. Exceptions
  6. Programming with dhp
    1. Global Variables
    2. Using dhp in C++
      1. dhp provides the following C++ class functions
        1. get_double
        2. get_string
        3. set_string
        4. cast
        5. is0
        6. operator=
      2. dhp provides the following C++ external functions
        1. operator+
        2. operator-
        3. operator*
        4. operator/
        5. compare functions
        6. sqrt
        7. trunc
        8. round
        9. fabs
        10. cmp
    3. Using dhp in C
      1. Naming conventions
      2. Data types
      3. dhp C functions
        1. dhp_get_double, dhp1024_get_double
        2. dhp_get_string, dhp1024_get_string
        3. dhp_set_double, dhp1024_set_double
        4. dhp_set_string, dhp1024_set_string
        5. dhp_add, dhp1024_add
        6. dhp_sub, dhp1024_sub
        7. dhp_mul, dhp1024_mul
        8. dhp_div, dhp1024_div
        9. dhp_sqrt, dhp1024_sqrt
        10. dhp_trunc, dhp1024_trunc
        11. dhp_round, dhp1024_round
        12. dhp_fabs, dhp1024_fabs
        13. dhp_neg, dhp1024_neg
        14. dhp_is0, dhp1024_is0
        15. dhp_cmp, dhp1024_cmp
  7. Examples
    1. computing the sequence of Fibonacci numbers
      1. using dhp in C
      2. using dhp in C++
      3. Results of execution
    2. Rump's Example
    3. Using to test parallelization of a stencil

Technical specifications


Current version:
Maximum size of the fraction:
Source code in:
C++ wrapper:
Operating system supported:
Documentation:
Support:
 2.1
unlimited
C++
provided
64-bit Linux/Windows
manual
on-line
`unlimited' should be understood as `determined by
the amount of the available memory'


WE ASSUME NO LIABILITY WHATSOEVER, AND DISCLAIM ANY EXPRESS OR IMPLIED WARRANTY, RELATING TO SALE AND/OR USE OF OUR PRODUCTS INCLUDING LIABILITY OR WARRANTIES RELATING TO FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, OR INFRINGEMENT OF ANY PATENT, COPYRIGHT OR OTHER INTELLECTUAL PROPERTY RIGHT. Our products are not intended for use in medical, life saving, life sustaining, critical control or safety systems, or in nuclear facility applications The software described in this document may contain software defects which may cause it to deviate from expected behavior.

Installation

To install dhp on Linux, copy the provided library ( linux_libdhp.a and included file ( dhp.h ) into the appropriate working directories ( on Linux it may be /usr/local/lib64 for libdhp.a and /usr/local/include for dhp.h ) renaming the library as libdhp.a; e.g. on Linux do:
sudo cp linux_libdhp.a /usr/local/lib64/libdhp.a

To install dhp on Windows, rename the ( extracted ) file windows_libdhp.lib as necessary and move it into the appropriate working directory; move included file ( dhp.h ) into the appropriate working directory.

High Precision values

dhp uses high precision representation of the floating point numbers. This representation has the form ( for non-zero values ):
1.fraction × 2exponent

and provides 32 bit exponent ( compare with 11 bits for IEEE double precision ). For C programs, dhp supports high precision representations with the fraction size of 128 and 1024 bits ( compare with 52 bits for IEEE double precision ). For C++ programs, dhp supports high precision representations with any fraction size .

The fraction size if a high precision representation of the floating point number supported by dhp, is called it precision.

Compiling and linking with the dhp library

This section describes how to use the dhp with your programs.

To use dhp in C/C++ programs, you must #include the file dhp.h in every source file that calls dhp functions, accesses dhp global variables, or uses #defined constants provided by dhp.

The procedures for linking dhp with the rest of your program vary according to the compiler and method you are using. When using C compiler, you must include the standard C++ library. On a typical Linux system this can be done as following
gcc myprogram.c -ldhp -lstdc++

or
g++ myprogram.cpp -ldhp
 

Exceptions

dhp functions, in addition to producing result, may also generate error code ( exceptions ). The possible error codes defined in dhp.h file and are:
 In this document for every dhp function we list the possible exceptions it may generate. The global variable dhp_error contains an error code of the last dhp function executed. The global variable dhp_error_sum contains summary of all errors reported during the use of dhp. It is a good idea to periodically check this variables to ensure they are equal to DHP_OK. Note that in the case of exception, the results of dhp functions are undefined. Note also that dhp items are always valid high precision floating point values.

Programming with dhp

This section describes how to make dhp calls in your program.

Global Variables

dhp_error
unsigned int dhp_error;
dhp_error contains the error code from the last  executed dhp function. If the call succeeded, dhp_error will be 0 ( DHP_ERR_OK ). A list of error code values, their meanings, and #defined constants for them is in section Exceptions.

dhp_error_sum
unsigned int dhp_error_sum;
dhp_error_sum record bitwise OR of the error code from dhp calls. A list of error code values, their meanings, and #defined constants for them is in section Exceptions.

dhp_version
unsigned short dhp_version;
dhp_version contains the dhp's version in packed BCD format. For example, if the dhp's version is 2.1, the value of dhp_version will be 0x0201. The high byte of dhp_version represents two digits to the left of the decimal point (one digit per nibble) and the low byte represents two digits to the right of the decimal point (again, one digit per nibble).

Using dhp in C++

The C++ class supported by dhp is dhpreal. It can be defined in the following ways:

dhpreal data() - defines data to be high precision representation of the floating point numbers that provides provides 32-bit exponent and 128-bit fraction.
dhpreal data( unsigned int FractionSizeInBits ) - defines data to be high precision representation of the floating point numbers that provides provides 32 bits exponent and and at least  FractionSizeInBits bits fraction.
dhpreal data( const dhpreal& source ) - defines data to be high precision representation identical to source ( having the same fraction size and the same value ).

Note that some compilers require explicit cast for defining dhprel data items with a parameter, e.g.
dhpreal data( ( unsigned int )100 );

to define data to be high precision representation of the floating point numbers that provides 32 bits exponent and at least 100 bits fraction.


dhp supports all arithmetic operations on the variables of the type dhpreal and allows mixing them with argument of the type double. For example, the following:
dhpreal a, b(500);

a = 1.23;
b = a*10;

will define a to be a floating point high precision representation with the fraction size of 128 bits, b to be a floating point high precision representation with the fraction size of 500 bits, set a to 1.23 and b to 12.3

dhp provides the following C++ class functions

get_double

double get_double()

return the value of this as a double.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhpreal a;
double d;
     .
d = a.get_double();

get_string

int get_string( char *str, unsigned int strSz, unsigned int cntAfterDot )

set str to represent the value of this in the following format:
[-+][D1][.D2]

where D1 - is a decimal integer denoting the integer part of the floating point value D2 - is a decimal integer denoting the fractional part of the floating point value. At most strSz characters of str are affected ( including null character terminating the string ) and at most cntAfterDot characters are set after the decimal point '.'.

return not-zero value if str contains the integer part D1 ( and decimal point in the case this is not a whole number ), zero - otherwise

Possible Exceptions: none

Example:
dhpreal a;
char s[1024];
int stat;
    .
stat = a.get_string( s, 1024, 4 );

set_string

void set_string( const char *str )

set this to floating point number represented as a string str. Input string str has the following format:
[-+][D1][.D2]

where D1 - is a decimal integer denoting the integer part of the floating point value D2 - is a decimal integer denoting the fractional part of the floating point value Result of on inputs not confirming to the above description is undefined.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhpreal a;
     .
a.set_string( "-123.456" );

cast

void cast( unsigned int FractionSizeInBits )

set this to be high precision representation of the floating point numbers that provides provides 32 bit exponent and at least FractionSizeInBits bits fraction. The old value of this is attempted to be preserved.

Possible Exceptions: none

Example:

dhpreal a;
.
a.cast( 1024 );

is0

int is0()

return not-zero value if this is equal to 0, zero - if this has numeric representation of a not zero value.

Possible Exceptions: none

Example:
dhpreal a( 2000 );
     .
  if ( !a.is0() )
     .

operator=

dhpreal& operator( double src )

set this to the value of src.

Possible Exceptions: DHP_ERR_INV, DHP_ERR_OVF

Example:
dhpreal a( 250 );
double d;
     .
  a = d;


dhpreal& operator=( const dhpreal& src )

set this to the value of src; the precision of this ( the size of it fraction ) is not affected.

Possible Exceptions: none

Example:
dhpreal a( 250 ), b;
     .
  a = b;


dhp provides the following C++ external functions

operator+

const dhpreal operator+( const dhpreal& src1, const dhpreal& src2 )

return dhpreal value equal src1 + src2. The precision of the returned result is set to the maximum of precisions of src1 and src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

const dhpreal operator+( double src1, const dhpreal& src2 )
const dhpreal operator+( const dhpreal& src1, double src2 )

return dhpreal value equal src1 + src2. The precision of the returned result is set to the maximum of precisions of the dhpreal argument.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_INV

Example:
dhpreal a, b( 200 ), c;
     .
a = b + c;
c = a + 1.2;

operator-

const dhpreal operator-( const dhpreal& src1, const dhpreal& src2 )

return dhpreal value equal src1 - src2. The precision of the returned result is set to the maximum of precisions of src1 and src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

const dhpreal operator-( double src1, const dhpreal& src2 )
const dhpreal operator-( const dhpreal& src1, double src2 )

return dhpreal value equal src1 - src2. The precision of the returned result is set to the maximum of precisions of the dhpreal argument.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_INV

Example:
dhpreal a, b( 200 ), c;
     .
a = b - c;
c = a - 1.2;


operator*

const dhpreal operator*( const dhpreal& src1, const dhpreal& src2 )

return dhpreal value equal src1 * src2. The precision of the returned result is set to the maximum of precisions of src1 and src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

const dhpreal operator*( double src1, const dhpreal& src2 )
const dhpreal operator*( const dhpreal& src1, double src2 )

return dhpreal value equal src1 * src2. The precision of the returned result is set to the maximum of precisions of the dhpreal argument.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_INV

Example:
dhpreal a, b( 200 ), c;
     .
a = b * c;
b = c*3.4;

operator/

const dhpreal operator/( const dhpreal& src1, const dhpreal& src2 )

return dhpreal value equal src1 / src2. The precision of the returned result is set to the maximum of precisions of src1 and src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_DZE

const dhpreal operator/( double src1, const dhpreal& src2 )
const dhpreal operator/( const dhpreal& src1, double src2 )

return dhpreal value equal src1 / src2. The precision of the returned result is set to the maximum of precisions of the dhpreal argument.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_INV, DHP_ERR_DZE

Example:
dhpreal a, b( 200 ), c;
     .
a = b / c;
c = 5.6/b;

compare functions

const int operator<( const dhpreal& src1, const dhpreal& src2 )
const int operator<=( const dhpreal& src1, const dhpreal& src2 )
const int operator>( const dhpreal& src1, const dhpreal& src2 )
const int operator>=( const dhpreal& src1, const dhpreal& src2 )
const int operator==( const dhpreal& src1, const dhpreal& src2 )
const int operator!=( const dhpreal& src1, const dhpreal& src2 )

Possible Exceptions: none

const int operator<( const dhpreal& src1, double src2 )
const int operator<( double src1, const dhpreal& src2 )

const int operator<=( const dhpreal& src1, double src2 )
const int operator<=( double src1, const dhpreal& src2 )

 const int operator>( const dhpreal& src1, double src2 )
const int operator>( double src1, const dhpreal& src2 )

const int operator>=( const dhpreal& src1, double src2 )
const int operator>=( double src1, const dhpreal& src2 )

const int operator==( const dhpreal& src1, double src2 )
const int operator==( double src1, const dhpreal& src2 )


const int operator!=( const dhpreal& src1, double src2 )
const int operator!=( double src1, const dhpreal& src2 )

Possible Exceptions:  DHP_ERR_INV

return not-zero value if corresponding relation between src1 and src2 is true, zero - otherwise.

Example:
dhpreal a, b( 200 ), c;
     .
if ( ( a < b ) && ( c >= 2.3 ) )


sqrt

dhpreal sqrt( const dhpreal& src )

return dhpreal value equal SQRT( src ). The precision of the returned result is set to the precision of src.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_INV

Example:
dhpreal a, b( 200 );
     .
a = sqrt( b );


trunc

dhpreal trunc( const dhpreal& src )

return dhpreal value equal to the integer nearest to src and not larger in absolute value than src. The precision of the returned result is set to the precision of src.

Possible Exceptions: none

Example:
dhpreal a, b( 200 );
     .
a = trunc( b );


round

dhpreal round( const dhpreal& src )

return dhpreal value equal to the integer nearest to src. The precision of the returned result is set to the precision of src.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhpreal a, b( 200 );
     .
a = round( b );


fabs

dhpreal fabs( const dhpreal& src )

return dhpreal value equal to the absolute value of src.The precision of the returned result is set to the precision of src.

Possible Exceptions: none

Example:
dhpreal a, b( 200 );
     .
a = fabs( b );


cmp

int cmp( const dhpreal& src1, const dhpreal& src2 )

return and integer less than, equal to, or greater than zero  is src1 is , respectively, less than, equal, or greater than src2

Possible Exceptions: none

Example:
dhpreal a, b( 200 );
     .
if ( cmp( a, b ) < 0 )

Using dhp in C

Naming conventions


All functions, provided by dhp have name that begins with 'dhp_' or 'dhp1024_' ( depending on the type of it arguments/result ) and followed by the low-case mnemonic, e.g 'dhp_add'.

Data types


The C data types supported by dhp are dhp_t and dhp1024_t

dhp_t is a high precision representation of the floating point numbers that provides 32 bit exponent and 128 bits fraction.

dhp1024_t is a high precision representation of the floating point numbers that provides 32 bit exponent and 1024 bits fraction.

dhp C functions

dhp provides several sets of functions for each of the supported data types.
  1. Functions for copying data to and from formats supported by the dhp. These functions are:
  2. Functions which are intended to be orthogonal with the standard floating point operations ( +, -, *, / etc. ). These functions are:
  3. Comparison functions:

dhp_get_double, dhp1024_get_double

double dhp_get_double( dhp_t src )
double dhp1024_get_double( dhp1024_t src )

return the value of src as a double.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhp_t a;
double d;
     .
d = dhp_get_double( a );

dhp_get_string, dhp1024_get_string

int dhp_get_string( dhp_t, src, char *str, unsigned int strSz, unsigned int cntAfterDot )
int dhp1024_get_string( dhp1024_t, src, char *str, unsigned int strSz, unsigned int cntAfterDot )

set str to represent the value of src in the following format:
[-+][D1][.D2]

where D1 - is a decimal integer denoting the integer part of the floating point value D2 - is a decimal integer denoting the fractional part of the floating point value. At most strSz characters of str are affected ( including null character terminating the string ) and at most cntAfterDot characters are set after the decimal point '.'.

return not-zero value if str contains the integer part D1 ( and decimal point in the case this is not a whole number ), zero - otherwise

Possible Exceptions: none

Example:
dhp_t a;
char s[1024];
int stat;
    .
stat = dhp_get_string( a, s, 1024, 4 );

dhp_set_double, dhp1024_set_double

dhp_t dhp_set_double( double src )
dhp1024_t dhp1024_set_double( double src )

return the item of an appropriate data type that has high precision representation of src.

Possible Exceptions: DHP_ERR_INV, DHP_ERR_OVF

Example:
dhp1024_t a;
double d;
     .
  a = dhp1024_set_double( d );

dhp_set_string, dhp1024_set_string

dhp_t dhp_set_string( const char *str )
dhp1024_t dhp1024_set_string( const char *str )

return the item of an appropriate data type that has high precision floating point number with the value represented as a string str. Input string str has the following format:
[-+][D1][.D2]

where D1 - is a decimal integer denoting the integer part of the floating point value D2 - is a decimal integer denoting the fractional part of the floating point value Result of on inputs not confirming to the above description is undefined.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhp_t a;
     .
a = dhp_set_string( "-123.456" );

dhp_add, dhp1024_add

dhp_t dhp_add( dhp_t src1, dhp_t src2 )
dhp1024_t dhp1024_add( dhp1024_t src1, dhp1024_t src2 )
return the item of an appropriate data type that has high precision floating point number with the value equal to src1 + src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhp1024_t a, b, c;
     .
a = dhp1024_add( b, c );

dhp_sub, dhp1024_sub

dhp_t dhp_sub( dhp_t src1, dhp_t src2 )
dhp1024_t dhp1024_sub( dhp1024_t src1, dhp1024_t src2 )
return the item of an appropriate data type that has high precision floating point number with the value equal to src1 - src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhp_t a, b, c;
     .
a = dhp_sub( b, c );

dhp_mul, dhp1024_mul

dhp_t dhp_mul( dhp_t src1, dhp_t src2 )
dhp1024_t dhp1024_mul( dhp1024_t src1, dhp1024_t src2 )
return the item of an appropriate data type that has high precision floating point number with the value equal to src1 * src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF

Example:
dhp_t a, b, c;
     .
a = dhp_mul( b, c );

dhp_div, dhp1024_div

dhp_t dhp_div( dhp_t src1, dhp_t src2 )
dhp1024_t dhp1024_div( dhp1024_t src1, dhp1024_t src2 )
return the item of an appropriate data type that has high precision floating point number with the value equal to src1 / src2.

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_DZE
Example:
dhp_t a, b, c;
     .
a = dhp_div( b, c );

dhp_sqrt, dhp1024_sqrt

dhp_t dhp_sqrt( dhp_t src )
dhp1024_t dhp1024_sqrt( dhp1024_t src )

return the item of an appropriate data type that has high precision floating point number with the value equal to SQRT( src )

Possible Exceptions: DHP_ERR_OVF, DHP_ERR_UNF, DHP_ERR_INV

Example:
dhp_t a, b;
     .
a = dhp_sqrt( b );

dhp_trunc, dhp1024_trunc

dhp_t dhp_trunc( dhp_t src )
dhp_t dhp1024_trunc( dhp1024_t src )

return the item of an appropriate data type that has high precision floating point number with the value equal to the integer nearest to src and not larger in absolute value than src

Possible Exceptions: none

Example:
dhp_t a, b;
     .
a = dhp_trunc( b );

dhp_round, dhp1024_round

dhp_t dhp_round( dhp_t src )
dhp_t dhp1024_round( dhp1024_t src )

return the item of an appropriate data type that has high precision floating point number with the value equal to the integer nearest to src

Possible Exceptions: none

Example:
dhp_t a, b;
     .
a = dhp_trunc( b );

dhp_fabs, dhp1024_fabs

dhp_t dhp_fabs( dhp_t src )
dhp1024_t dhp1024_fabs( dhp1024_t src )

return the item of an appropriate data type that has high precision floating point number with the value equal the absolute value of src

Possible Exceptions: none

Example:
dhp_t a, b;
     .
a = dhp_fabs( b );

dhp_neg, dhp1024_neg

dhp_t dhp_neg( dhp_t src )
dhp1024_t dhp1024_neg( dhp1024_t src )

return the item of an appropriate data type that has high precision floating point number with the value equal to  -src

Possible Exceptions: none

Example:
dhp_t a, b;
     .
a = dhp_neg( b );

dhp_is0, dhp1024_is0

int dhp_is0( dhp_t src )
 int dhp1024_is0( dhp1024_t src )
return not-zero value if src is equal to 0, zero - if src has numeric representation of a not zero value.

Possible Exceptions: none

Example:
dhp_t a;
     .
  if ( !dhp_is0( a ) )
     .

dhp_cmp, dhp1024_cmp

int dhp_cmp( dhp_t src1, dhp_t src2 )
int dhp1024_cmp( dhp_t src1, dhp_t src2 )

return and integer less than, equal to, or greater than zero  is src1 is, respectively, less than, equal, or greater than src2

Possible Exceptions: none

Example:
dhp_t a, b;
     .
if ( dhp_cmp( a, b ) < 0 )


Examples

computing the sequence of Fibonacci numbers

Assume that the task is to compute the first 100 items from the sequence of Fibonacci numbers
F0 = 0, F1 = 1, Fn+2 = Fn+1 + Fn

using de Moivre's-Binet's formula
Fn = ( φn - φ'n )/ √5

where
φ = ( 1 + √5 )/2, φ' = ( 1 - √5 )/2


The following program uses standard double precision variables to achieve this goal:
double p, p1, pn, p1n, fn, sqrt5;
int i;

sqrt5 = sqrt( 5. );
p = ( 1. + sqrt5 )/2.;
p1 = ( 1. - sqrt5 )/2.;
pn = p;
p1n = p1;
for ( i = 0; i < 100; i++ ) {
fn = ( pn - p1n )/sqrt5;
printf( "F%d = %f\n", i, fn );
pn *= p;
p1n *= p1;
}

using dhp in C

In C the above program can be rewritten to utilize dhp as following:
#include <dhp.h>
    .
dhp_t p, p1, pn, p1n, fn, sqrt5;
int i;
char s[1024];

sqrt5 = dhp_sqrt( dhp_set_double( 5. ) );
p = dhp_div( dhp_add( dhp_set_double( 1. ), sqrt5 ), dhp_set_double( 2. ) );
p1 = dhp_div( dhp_add( dhp_set_double( 1. ), sqrt5 ), dhp_set_double( 2. ) );
pn = p;
p1n = p1;
for ( i = 0; i < 100; i++ ) {
fn = dhp_div( dhp_sub( pn, p1n ), sqrt5 );
dhp_get_string( fn, s, 1000, 6 );
printf( "F%d = %s\n", i, s );
pn = dhp_mul( pn, p );
p1n = dhp_mul( p1n, p1 );
}

This example shows that implementing dhp functionality in C requires significant modification of the original program.

using dhp in C++

The program that uses standard double precision variables ( listed above ) can be rewritten to utilize dhp as following:
#include <dhp.h>
     .
dhpreal p, p1, pn, p1n, fn, sqrt5;
int i;
char s[1000];

fn = 5.;
sqrt5 = sqrt( fn );
p = ( 1. + sqrt5 )/2.;
p1 = ( 1. - sqrt5 )/2.;
pn = p;
p1n = p1;
for ( i = 0; i < 100; i++ ) {
fn = ( pn - p1n )/sqrt5;
fn.get_string( s, 1000, 6 );
printf( "F%d = %s\n", i, s );
pn *= p;
p1n *= p1;
}

The changes from the original ( normal ) program are marked in this color. This example shows that implementing dhp functionality in C++ is straightforward and requires very minor modification of the original program.

Results of execution

The following table lists generated Fibonacci number for selected indexes:
  1. The column under # header presents the index of Fibonacci number in the following cells of the raw.
  2. The column under double header presents Fibonacci numbers generated by the program using standard double precision variables ( listed above ).
  3. The column under 64-bit uints header presents Fibonacci numbers generated by the program utilizing the original formula using 64-bit unsigned integer variables.
  4. The column under dhp header presents Fibonacci numbers generated by the program using dhp ( listed above, C and C++ versions generated the identical results ).
Results that are not accurate marked in this color.
# double 64-bit uints dhp
68 72723460248141.156250
 72723460248141 72723460248141.000000
69 117669030460994.250000 117669030460994 117669030460994.000000
70 190392490709135.406250 190392490709135 190392490709135.000000
71 308061521170129.625000 308061521170129 308061521170129.000000
72 498454011879265.125000 498454011879264 498454011879264.000000
73 806515533049394.750000 806515533049393 806515533049393.000000
74 1304969544928659.750000 1304969544928657 1304969544928657.000000
75 2111485077978054.750000 2111485077978050 2111485077978050.000000
90 2880067194370823680.000000 2880067194370816120 2880067194370816120.000000
91 4660046610375542784.000000 4660046610375530309 4660046610375530309.000000
92 7540113804746366976.000000
7540113804746346429 7540113804746346429.000000
93 12200160415121909760.000000 12200160415121876738 12200160415121876738.000000
94 19740274219868274688.000000 1293530146158671551 19740274219868223167.000000
95  31940434634990190592.000000 13493690561280548289 31940434634990099905.000000
96 51680708854858465280.000000 14787220707439219840 51680708854858323072.000000
97  83621143489848655872.000000 9834167195010216513 83621143489848422977.000000
98 135301852344707137536.000000 6174643828739884737 135301852344706746049.000000
99  218922995834555793408.000000 16008811023750101250 218922995834555169026.000000

Rump's Example

Rump's example is to compute the expression
f = ( 333.75 - a2 ) b6 + a2 (11. a2 b2 - 121. b4 - 2. ) + 5.5 b8 + a / ( 2. b )
with
a = 77617. and b = 33096.

Attempts to implement the above formula using the standard ( single or double or extended ) precision, available on the modern computers, produce wrong result for f.

The following program uses standard double precision variables to implement Rump's example:

double a, b, a2, b2, b4, b6, b8, f;

a = 77617.;
b = 33096.;
a2 = a*a;
b2 = b*b;
b4 = b2*b2;
b6 = b4*b2;
b8 = b4*b4;

f = ( 333.75 - a2 )*b6 + a2* ( 11.*a2*b2 - 121.*b4 -2. ) + 5.5*b8 + a/(2.*b);

printf( "IEEE\t%f\n", f );

Execution of this program generates ( wrong ) result:

IEEE 1.172604

The program that uses standard double precision variables ( listed above ) can be rewritten in C++ to utilize dhp as following:

#include <dhp.h>
.
dhpreal a, b, a2, b2, b4, b6, b8, f;

a = 77617.;
b = 33096.;
a2 = a*a;
b2 = b*b;
b4 = b2*b2;
b6 = b4*b2;
b8 = b4*b4;

f = ( 333.75 - a2 )*b6 + a2* ( 11.*a2*b2 - 121.*b4 -2. ) + 5.5*b8 + a/(2.*b);

char s[1024];
f.get_string( s, 1023, 64 );
printf( "dhp\t%s\n", s );

Execution of this program generates ( the correct ) result:

dhp -.8273960599468213681411650954798162919984358498609068290206066511

The changes from the original ( "normal" ) program are highlighted in this color. This example shows that implementing dhp functionality in C++ is straightforward and requires very minor modification of the original program. Note also how easy it is, using dhp, to print the data with "any" desirable precision.

Using to test parallelization of a stencil

See this on how dhp was used to test parallel code generated by the Dalsoft Code Optimizer ( dco ).