Logo Search packages:      
Sourcecode: calcoo version File versions

aux.c

/*
 * Calcoo: aux.c
 *
 * Copyright (C) 2001, 2002 Alexei Kaminski
 *
 * auxuliary general-purpose functions
 * calcoo code translators
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "basic.h"
#include "codes.h"
#include "defaults.h"
#include "displays.h"
#include "const.h"
#include "aux_headers.h"

/*******************************************************************
 *
 * Calcoo-specific auxiliary functions
 *
 *******************************************************************/

int digit_to_code(int d)
{
      switch (d) {
      case 0: return CODE_0; break;
      case 1: return CODE_1; break;
      case 2: return CODE_2; break;
      case 3: return CODE_3; break;
      case 4: return CODE_4; break;
      case 5: return CODE_5; break;
      case 6: return CODE_6; break;
      case 7: return CODE_7; break;
      case 8: return CODE_8; break;
      case 9: return CODE_9; break;
      }
      error_occured("digit_to_code(): request to transate a non-digit",
                  FALSE);
      return CODE_0;
}

int code_to_digit(int c)
{
      switch (c) {
      case CODE_0: return 0; break;
      case CODE_1: return 1; break;
      case CODE_2: return 2; break;
      case CODE_3: return 3; break;
      case CODE_4: return 4; break;
      case CODE_5: return 5; break;
      case CODE_6: return 6; break;
      case CODE_7: return 7; break;
      case CODE_8: return 8; break;
      case CODE_9: return 9; break;
      }
      error_occured("code_to_digit(): request to transate a non-digit code",
                  FALSE);
      return 0;
}

int binop_to_od(int c)
{
      switch (c) {
      case CODE_ADD: return OD_ADD; break;
      case CODE_SUB: return OD_SUB; break;
      case CODE_MUL: return OD_MUL; break;
      case CODE_DIV: return OD_DIV; break;
      case CODE_POW: return OD_POW; break;
      }
      error_occured("binop_to_od(): request to transate a non-binop code",
                  FALSE);
      return OD_ADD;
      
}

int inverse_sign(int s)
{
      switch (s) {
      case SIGN_PLUS:  return SIGN_MINUS; break;
      case SIGN_MINUS: return SIGN_PLUS;  break;
      }
      error_occured("inverse_sign(): request to invert a non-sign",
                  FALSE);
      return SIGN_PLUS;
}

int priority(int op)
{
      switch (op) {
      case CODE_ADD:
      case CODE_SUB:
            return PRIORITY_CODE_ADD; break;
      case CODE_MUL:
      case CODE_DIV:
            return PRIORITY_CODE_MUL; break;
      case CODE_POW:
            return PRIORITY_CODE_POW; break;
      }
      error_occured("priority(): request of priority of a non-binop code",
                  FALSE);
      return PRIORITY_CODE_ADD;
}

/*******************************************************************
 *
 * General-purpose auxiliary and math functions
 *
 *******************************************************************/

int last_digit(double a)
{
      return round_double_to_int(a - 10.0 * floor(a / 10.0));
}

int almost_integer(double a, double precision)
{
      if (fabs(a - round_double_to_int(a)) < precision * fabs(a))
            return TRUE;
      else 
            return FALSE;
}

/* the name is weird on purpose; in order to prevent possible
 * name conflicts */
double sign_of_double(double a)
{
      if (a > 0.0)
            return 1.0;
      else if (a < 0.0)
            return -1.0;
      else 
            return 0.0;
}

double max_fabs(double a, double b) 
{
      double fa, fb;
      fa = fabs(a);
      fb = fabs(b);
      if (fa > fb)
            return fa;
      else 
            return fb;
}


double smart_sum(double a, double b, double precision)
{
      double sumabs, sumsign, cutoff;

      if ( fabs(a + b) < fabs(a) * precision ){
            /* hack to have 100.1 - 100 - 0.1 == 0 */
            return 0.0;
      } else {
            sumabs = fabs(a + b);
            sumsign = sign_of_double(a + b); 
            if (floor(log10(max_fabs(a, b))) <= ceil(log10(sumabs)))
                  return sumsign * sumabs;
            cutoff = pow(10, ceil(log10(max_fabs(a, b)))) * precision;
            if (cutoff < sumabs) {
                  return sumsign * cutoff * rint(sumabs / cutoff);
            } else {
                  return 0.0;
            }
      }

}


/* Yes, this is one weird name. It used to be just "round", 
 * but it conflicted with a function with the same name
 * defined in math.h on Mac OS X. In order to prevent this and 
 * possible future name conflicts, the present name was chosen */
int round_double_to_int(double a)
{
      if (a >= 0.0)
            return (int)(a + 0.5);
      else
            return (int)(a - 0.5);
}

/* the formulas used in the functions below are taken from 
 * "Handbook on mathematic functions", ed. Abramovitz and Stegun.
 * The relative precision may be as low as 10^(-7) for smaller numbers, 
 * but it is FAST */

int fact_too_large(double x)
/* determines if the factorial is too large to be shown by calcoo
 * this is only an estimate, it still can turn out to be larger than
 * the calcoo's max number, but if it returns FALSE it is at least guaranteed
 * that the factorial will fit into double */
{
      double log_10_x_fact;
      if (x <= 0.0 )
            return TRUE;
      log_10_x_fact = (x + 0.5) * log10(x) - x * log10(exp(1.0));
      return (log_10_x_fact > pow(10, EXP_INPUT_LENGTH) );
}

double fact_function(double x)
/* if calcoo is unable to show all the meaningful digits of the result,
 * there is no sense in the exact calculation, so we can use the
 * Stirling formula for the approximate calculation. Otherwise, we
 * call fact_function_jr(x), which calculates the factorial of integer
 * numbers exactly, but is slower */
{
      double log_10_x_fact, a;
      log_10_x_fact = (x + 0.5) * log10(x) - x * log10(exp(1.0));
      if (log_10_x_fact > INPUT_LENGTH ) {
            a = 1.0/(x+1.0);
            return exp(-x-1.0) * sqrt(2.0*PI) * pow(x+1.0, (x+0.5))
                  * ( 1.0 
                      + a * ((1.0/12.0) 
                      + a * ((1.0/288.0)
                      + a * ((-139.0/51840.0)
                      + a * ((-571.0/2488320.0)
                        )))));
      } else
            return fact_function_jr(x);
}

double fact_function_jr(double x)
{
      if (x >= 1.0)
            return x * fact_function_jr(x-1.0);
      else
            return ( 1
                   + x * ((-0.577191652) 
                       + x * ((0.988205891)
                           + x * ((-0.897056937) 
                               + x * ((0.918206857)
                                   + x * ((-0.756704078)
                                       + x * ((0.482199394)
                                           + x * ((-0.193527818) 
                                               + x * (0.035868343)
                                           ))))))));
}


/*******************************************************************
 *
 * Functions used in debugging
 *
 *******************************************************************/

/* this is a kind of assert() */
/* REPORT_NONCRITICAL_ERRORS is defined in defaults.h */
void error_occured(char *message, int is_critical)
{
      if (is_critical) {
            fprintf(stderr, "Calcoo critical error: %s, exiting\n", 
                  message);
            exit(1);
      } else 
            if (REPORT_NONCRITICAL_ERRORS)
                  fprintf(stderr, "Calcoo non-critical error: %s\n", 
                        message);
}

/* DEBUG_MESSAGES_ON is defined in defaults.h */
void mess(char *message, int number)
{
      if(DEBUG_MESSAGES_ON)
            fprintf(stderr,   "%s %d\n", message, number);
}

void messd(char *message, double number)
{
      if(DEBUG_MESSAGES_ON)
            fprintf(stderr,   "%s %f\n", message, number);
}

void verify_malloc(void *p)
{
      if (p == NULL){
            fprintf(stderr, "malloc() failed, exiting\n");
            exit(1);
      }
}

Generated by  Doxygen 1.6.0   Back to index