ch24safe.c

稀疏矩阵、链表、图、队列、二叉树、多叉树、排序、遗传算法等的实现

 * double operands, which may raise error conditions, and indicates to *
 * the caller the precision of the result.                             *
 ***********************************************************************/

/***********************************************************************
 * Name:          flpSafeMultiplySensitive
 * 
 * Description:   This function performs a safe floating-point multiply
 *                operation, with detection and signalling of errors
 *                
 * Parameters: 
 *    flpFirstMultiplicand    the first Multiplicand
 *    flpFirstMultiplicandPrecision precision of the first multiplicand
 *    flpSecondMultiplicand   the second Multiplcand
 *    flpSecondMultiplicandPrecision precision of the second multiplicand
 *    pflpResultPrecision     pointer to a double which will contain
 *                            the resultant precision after the multiply.
 *                            If this pointer is NULL it is not used
 *    RaiseError              an int containing a logical value. 
 *                            If TRUE then errors will be Raised, if
 *                            FALSE they will not, even if detected
 *    pErrorRaised            an int * pointing to a word containing
 *                            error indication. If this pointer is NULL
 *                            then it is not used
 * 
 * External/Global variables:
 *    N/A
 *
 * Returns:
 *    the product, as a double
 *    the result precision, if requested
 *    the error condition, if requested
 *    Raises a signal error, if requested, and required
 * 
 **********************************************************************/

double flpSafeMultiplySensitive ( double flpFirstMultiplicand,
                                  double flpFirstMultiplicandPrecision,
                                  double flpSecondMultiplicand,
                                  double flpSecondMultiplicandPrecision,
                                  double * pflpResultPrecision,
                                  int RaiseError, int * pErrorRaised )
{
   double flpReturnValue           = 0.0;
   double flpZero                  = 0.0;
   double flpLocalResultPrecision  = 0.0;
   double flpLocalFirstPrecision   = 0.0;
   double flpLocalSecondPrecision  = 0.0;
   long double flpLocalReturnValue = 0.0;
   
   /* Set the outgoing precision initially to its minimum:            */
   if ( pflpResultPrecision != NULL )
   {
      flpLocalFirstPrecision = DBL_EPSILON * fabs( flpFirstMultiplicand );
      if ( DBL_EPSILON > flpLocalFirstPrecision )
         flpLocalFirstPrecision = DBL_EPSILON;
      if ( flpFirstMultiplicandPrecision > flpLocalFirstPrecision )
         flpLocalFirstPrecision = flpFirstMultiplicandPrecision;
      flpLocalSecondPrecision = DBL_EPSILON * fabs( flpSecondMultiplicand );
      if ( DBL_EPSILON > flpLocalSecondPrecision )
         flpLocalSecondPrecision = DBL_EPSILON;
      if ( flpSecondMultiplicandPrecision > flpLocalSecondPrecision )
         flpLocalSecondPrecision = flpSecondMultiplicandPrecision;
      flpLocalResultPrecision = ( flpLocalFirstPrecision   * 
                                  fabs ( flpSecondMultiplicand ) )
                              + ( flpLocalSecondPrecision * 
                                  fabs ( flpFirstMultiplicand ) );
      if ( flpLocalResultPrecision < DBL_EPSILON )
           flpLocalResultPrecision = DBL_EPSILON;
      *pflpResultPrecision = flpLocalResultPrecision;
   }
   
   /* Set the outgoing error condition to no error:                   */
   if ( pErrorRaised != NULL )
       *pErrorRaised = 0;

   /* If either of the incoming operands is close to zero, set the    */
   /* result to zero, and compute the resultant precision as being the*/
   /* largest of the incoming precisions:                             */
   if (    flpcmp ( flpFirstMultiplicand,  "==", 0.0 ) 
        || flpcmp ( flpSecondMultiplicand, "==", 0.0 ) )
   {
      if ( pflpResultPrecision != NULL )
      {
        if ( flpFirstMultiplicandPrecision > flpSecondMultiplicandPrecision ) 
              *pflpResultPrecision = flpFirstMultiplicandPrecision;
         else *pflpResultPrecision = flpSecondMultiplicandPrecision;
      }
      /* Return to caller with the computed value:                    */
      return ( flpReturnValue );
   }
   
   /* Neither of the incoming multiplicands is (close to) zero. So    */
   /* compute the resultant precision of the product:                 */
   if ( pflpResultPrecision != NULL )
   {
      *pflpResultPrecision = (  ( flpFirstMultiplicandPrecision  * 
                                  fabs ( flpSecondMultiplicand ) )
                              + ( flpSecondMultiplicandPrecision * 
                                  fabs ( flpFirstMultiplicand  ) ) );
   }
   
   /* Now compute the product, but in long double format. Look at the */
   /* result to determine whether we have overshot the allowable      */
   /* size for "double". The same comments apply here as are stated in*/
   /* the routine flpSafeDivideSensitive:                             */
   /* In the following expression the casts are made explicitly to    */
   /* ensure that the conversions to (long double) are performed early*/
   /* enough:                                                         */
   flpLocalReturnValue = (long double) flpFirstMultiplicand * 
                         (long double) flpSecondMultiplicand;
   /* Now we know the magnitude of the result, adjust the outgoing    */
   /* precision, to account for the minimum error for numbers of this */
   /* value-range:                                                    */
   if ( flpLocalResultPrecision < DBL_EPSILON * fabs ( flpLocalReturnValue ) )
        flpLocalResultPrecision = DBL_EPSILON * fabs ( flpLocalReturnValue ) ;

   if ( fabs ( flpLocalReturnValue ) > DBL_MAX )
   {
      /* Indicate that the precision is extremely low:                */
      if ( pflpResultPrecision != NULL )
          *pflpResultPrecision = DBL_MAX;
      
      /* Indicate the error we have detected:                         */
      if ( pErrorRaised != NULL )
          *pErrorRaised = ARITHMETIC_PRECISION_BAD;
      
      /* Set the return value to the maximum possible, with the       */
      /* the sign of the real result:                                 */
      if ( flpLocalReturnValue < 0.0L )
         flpReturnValue = - DBL_MAX;
      else flpReturnValue = DBL_MAX;
      
      /* If we are raising errors, then raise the error, otherwise    */
      /* return to the caller with a zero value:                      */
      if ( RaiseError )
      {
         raise ( ARITHMETIC_PRECISION_BAD );
      } 

      /* return to caller with computed value:                        */
      return ( flpReturnValue );

   } else {

      /* All is well: set the return value to the computed value:     */
      flpReturnValue = flpLocalReturnValue;

      /* Look at the resultant precision and compare it in magnitude  */
      /* to the result. If better than ARITHMETIC_ACCEPT_PRECISION    */
      /* then accept it, otherwise raise an error:                    */
      /* In any case, return it to the caller:                        */
      if ( pflpResultPrecision != NULL )
          *pflpResultPrecision = flpLocalResultPrecision;

      /* Note that the following calculation is done as a multiplica- */
      /* tion to avoid having to do all the division tests. The cross-*/
      /* multiplication of the inequality                             */
      /*      (error / value) > acceptable_precision                  */
      /* works for all values of "value"                              */
      if ( flpLocalResultPrecision  >
           ( fabs ( flpReturnValue ) > DBL_EPSILON ?
               ( fabs ( flpReturnValue ) * ARITHMETIC_ACCEPT_PRECISION ) : 
               DBL_EPSILON ) )
      {

         /* Indicate the error we have detected:                      */
         if ( pErrorRaised != NULL )
             *pErrorRaised = ARITHMETIC_BEYOND_RANGE;
      
         /* If we are raising errors, then raise the error:           */
         if ( RaiseError )
         {
            raise ( ARITHMETIC_BEYOND_RANGE );
         }  /* end of "if/then" raising ARITHMETIC_BEYOND_RANGE error */ 
      }     /* end of "if/then" precision is bad                      */
   }        /* end of "if/else" result lies with acceptable range     */
   
   /* Return to caller with the computed value:                       */
   return ( flpReturnValue );

}

/***********************************************************************
 * This is the first ("simple") safe multiply procedure for floating   *
 * (double) values. It calls the "complex" safe multiply procedure     *
 * (flpSafeMultiplySensitive), indicating that no errors are to be     *
 * raised, to perform the body of the calculation.                     *
 ***********************************************************************/

/***********************************************************************
 * Name:          flpSafeMultiply
 * 
 * Description:   This function performs a safe floating-point multiply
 *                operation, without signalling any error
 *                
 * Parameters: 
 *    flpFirstMultiplicand    the first Multiplicand
 *    flpSecondMultiplicand   the second Multiplcand
 * 
 * External/Global variables:
 *    N/A
 *
 * Returns:
 *    the product, as a double
 * 
 **********************************************************************/

double flpSafeMultiply (double flpFirstMultiplicand,
                        double flpSecondMultiplicand)
{
  double flpReturnValue = 0.0;
   
  flpReturnValue =  flpSafeMultiplySensitive ( flpFirstMultiplicand, 0.0,
                                               flpSecondMultiplicand, 0.0,
                                               NULL, FALSE, NULL);
  
  return (flpReturnValue);
 
} 

/***********************************************************************/

/***********************************************************************
 * This is the second ("complex") level Safe Addition function for     *
 * double operands, which may raise error conditions, and indicates to *
 * the caller the precision of the result.                             *
 ***********************************************************************/

/***********************************************************************
 * Name:          flpSafeAddSensitive
 * 
 * Description:   This function performs a safe floating-point add
 *                operation, with detection and signalling of errors
 *                
 * Parameters: 
 *    flpFirstAddend          the first addend
 *    flpFirstAddendPrecision precision of the first addend
 *    flpSecondAddend         the second addend
 *    flpSecondAddendPrecision precision of the second addend
 *    pflpResultPrecision     pointer to a double which will contain
 *                            the resultant precision after the add.
 *                            If this pointer is NULL it is not used
 *    RaiseError              An int containing a logical value. 
 *                            If TRUE then errors will be Raised, if
 *                            FALSE they will not, even if detected
 *    pErrorRaised            An int * pointing to a word containing
 *                            error indication. If this pointer is NULL
 *                            then it is not used
 * 
 * External/Global variables:
 *    N/A
 *
 * Returns:
 *    the sum, as a double
 *    the result precision, if requested
 *    the error condition, if requested
 *    Raises a signal error, if requested, and required
 * 
 **********************************************************************/

double flpSafeAddSensitive ( double flpFirstAddend,
                             double flpFirstAddendPrecision,
                             double flpSecondAddend,
                             double flpSecondAddendPrecision,
                             double * pflpResultPrecision,
                             int RaiseError, int * pErrorRaised )
{
   double flpReturnValue           = 0.0;
   double flpZero                  = 0.0;
   double flpLocalResultPrecision  = 0.0;
   double flpLocalFirstPrecision   = 0.0;
   double flpLocalSecondPrecision  = 0.0;
   long double flpLocalReturnValue = 0.0;
   
   /* Set the outgoing precision initially to its minimum:            */
   if ( pflpResultPrecision != NULL )
   {
      flpLocalFirstPrecision = DBL_EPSILON * fabs( flpFirstAddend );
      if ( DBL_EPSILON > flpLocalFirstPrecision )
         flpLocalFirstPrecision = DBL_EPSILON;
      if ( flpFirstAddendPrecision > flpLocalFirstPrecision )
         flpLocalFirstPrecision = flpFirstAddendPrecision;
      flpLocalSecondPrecision = DBL_EPSILON * fabs( flpSecondAddend );
      if ( DBL_EPSILON > flpLocalSecondPrecision )
         flpLocalSecondPrecision = DBL_EPSILON;
      if ( flpSecondAddendPrecision > flpLocalSecondPrecision )
         flpLocalSecondPrecision = flpSecondAddendPrecision;
      flpLocalResultPrecision = flpLocalFirstPrecision    
                              + flpLocalSecondPrecision  ;
      if ( flpLocalResultPrecision < DBL_EPSILON )
         flpLocalResultPrecision = DBL_EPSILON;
      *pflpResultPrecision = flpLocalResultPrecision;
   }
   
   /* Set the outgoing error condition to no error:                   */
   if ( pErrorRaised != NULL )
       *pErrorRaised = 0;
   
   /* Now compute the sum, but in long double format. Look at the     */
   /* result to determine whether we have overshot the allowable      */
   /* size for "double". The same comments apply here as are stated in*/
   /* the routine flpSafeDivideSensitive:                             */
   /* In the following expression the casts are made explicitly to    */
   /* ensure that the conversions to (long double) are performed early*/
   /* enough:                                                         */
   flpLocalReturnValue = (long double) flpFirstAddend +
                         (long double) flpSecondAddend;
   /* Now we know the magnitude of the result, adjust the outgoing    */
   /* precision, to account for the minimum error for numbers of this */
   /* value-range:                                                    */
   if ( flpLocalResultPrecision < DBL_EPSILON * fabs ( flpLocalReturnValue ) )
        flpLocalResultPrecision = DBL_EPSILON * fabs ( flpLocalReturnValue ) ;

   if ( fabs ( flpLocalReturnValue ) > DBL_MAX )
   {
      /* Indicate that the precision is extremely low:                */
      if ( pflpResultPrecision != NULL )
          *pflpResultPrecision = DBL_MAX;
      
      /* Indicate the error we have detected:                         */
      if ( pErrorRaised != NULL )
          *pErrorRaised = ARITHMETIC_PRECISION_BAD;
      
      /* Set the return value to the maximum possible, with the       */
      /* the sign of the real result:                                 */
      if ( flpLocalReturnValue < 0.0L )
         flpReturnValue = - DBL_MAX;
      else flpReturnValue = DBL_MAX;