First, if anyone knows a clean way to add quadruple precision to numerical recipes in C++ please tell me (or I suppose I could even learn Fortran if that's easier).

I'm trying to change some calculations to quaduple precision. I expected it to be straightforward to redefine DP to be quadruple as opposed to double. So I downloaded a package from http://www.nersc.gov/~dhb/mpdist/mpdist.html and have successfully implemented simple calculations (e.g., 4*4 etc).

Here is a sample program, so you can see how the code works:

CODE:
**************************************
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include "qd.h"
#include <cmath>
#include "x86.h"

using std::cout;
using std::cin;
using std::endl;

int main() {
unsigned short old_cw;
x86_fix_start(&old_cw);

qd_real a, b, c;

cout << "Enter a:";
cin >> a;
cout << "Enter b:";
cin >> b;

c = a + b;

qd_real asquared = a*a;
qd_real bsquared = sqr(b);

cout << "a = " << a << endl;
cout << "b = " << b << endl;
cout << "2ab = " << 2*a*b << endl;
cout << "a*a = " << asquared << endl;
cout << "b*b = " << bsquared << endl;
cout << "a*a-b*b = " << asquared -bsquared << endl;
cout << "a + b = " << c << endl;
cout << "sin(c) = " << sin(c) << endl;

x86_fix_end(&old_cw);

return 0;
}

********************************************
OUTPUT:

Enter a:5
Enter b:3
a = 5.000000000000000000000000000000000000000000000000 000000000000000E0
b = 3.000000000000000000000000000000000000000000000000 000000000000000E0
2ab = 3.000000000000000000000000000000000000000000000000 000000000000000E1
a*a = 2.500000000000000000000000000000000000000000000000 000000000000000E1
b*b = 9.000000000000000000000000000000000000000000000000 000000000000000E0
a*a-b*b = 1.600000000000000000000000000000000000000000000000 000000000000000E1
a + b = 8.000000000000000000000000000000000000000000000000 000000000000000E0
sin(c) = 9.893582466233817778081235982452886721164190808857 612628177152199E-1

*********************************************

I'm now trying to redefine 'DP' to be quadruple precision by changing the lines in nrtypes.h and nrutil.h mentioned on page 17 of the text to 'typedef qd_real DP;'

However, it is getting errors in compilation:
**********************************************
In file included from nrutil.h:1,
from nr.h:5,
from bsstep.cpp:2:
nrutil_nr.h:10: syntax error before `;'
nrutil_nr.h:416: `DP' was not declared in this scope
nrutil_nr.h:416: template argument 1 is invalid
nrutil_nr.h:416: ISO C++ forbids declaration of `cc_p' with no type
nrutil_nr.h:416: ISO C++ forbids declaration of `cr_p' with no type

etc.
********************************************
Line 10 of nrutil_nr.h was the line that was changed to 'typedef qd_real DP;'

I suspect that it has something to do with the fact that x86_fix_start and _end are in main. Somehow it doesn't know about qd_real when it goes through nr.h

Any ideas/suggestions would be greatly appreciated.

Have you tried

#include "qd.h"

in the 2 nr files where the typedefs are?

Well, I thought I had...

That seems to have cleared up most of the problems. Still getting a few errors, but I think I can track them down.

thanks,

Joel

Code changes:
Replace double by long double when declaring data type.
All mathematical functions, such as pow, fabs etc., take a double argument and return a double. The quadruple precision equivalents have an added l at the end of the function name, eg. powl, fabsl.
When printing or scanning a long double, you must use %Le or %Lf.

Now when compiling:
Use cc128 instead of cc. If using a Makefile, cc128 is located at /usr/bin/cc128 on the IBM SP2.
Make sure you are not linking in the library libc.a (denoted by -lc when linking). cc128 uses libc128.a automatically (needed for printf to understand %Le etc.) unless told to do otherwise.

Note:
Optimiser's may and can cause a problem when using quadruple precision. Check results from programs run with and without optimisation.
When compiling with cc128, (http://www.eatsum.com) you don't need to use the -qlongdouble flag.