psmile_ddadd.c

Go to the documentation of this file.

/*-----------------------------------------------------------------------
* Copyright 2008-2010, DKRZ, Hamburg, Germany.
* All rights reserved. Use is subject to OASIS4 license terms.
*-----------------------------------------------------------------------*/

#include <stddef.h>
#include "PSMILe_f2c.h"

/*This function is based on the code from the MPFUN90 library
 *http://crd.lbl.gov/~dhbailey/mpdist/
 *
 *It was necessary to rewrite this function in C, because F90
 *does not support the volatile keyword. This keyword is
 *necessary to avoid compiler optimisations which use more than
 *64bit for intermediate floating point results.
 *However, this has some performance drawbacks. Another possible
 *solution is to force the compiler to use 64bit arithmetics for
 *this function. Unfortunately, this would require compiler
 *specific flags. Therefore, I will use this general solution for
 *now and may add optimisations for different compilers later.
 */

/* This function add two double-double precision numbers.*/

#ifdef ANSI_C

void psmile_ddadd (double * dda, double * ddb)

#else

void psmile_ddadd (dda, ddb)

double * dda, * ddb;

#endif

{
   volatile double e, t1, t2, t2_1, t2_2, t2_3;
/* Compute dda + ddb using Knuth's trick.*/
   t1 = dda[0] + ddb[0];
   e = t1 - dda[0];
   t2_1 = t1 - e;
   t2_2 = dda[0] - t2_1;
   t2_3 = ddb[0] - e;
   t2   = t2_3 + t2_2;
   t2  += dda[1];
   t2  += ddb[1];

/* The result is t1 + t2, after normalization.*/
   ddb[0] = t1 + t2;
   e = t1 + t2;
   e -= t1;
   ddb[1] = t2 - e;
   
   return;
}