Department of Scientific Computing    Institute for Numerical Simulation    University of Bonn

Documentation

Download

Programming References

Bug Reports / Suggestions

FAQ

Authors

           

---

Math.cc

#include "Math.hpp"

void swapbytes(void *p , int size, int nitem) {
 assert((size&1)==0) ;
 
 char *p1=(char *)p , a ;
  
 for (int j=0; j<nitem; j++) {
   for (int i=0; i<size/2; i++) {
     a           =p1[i] ;
     p1[i]       =p1[size-1-i] ;
     p1[size-1-i]=a ;
   }
   p1=&p1[size] ;
 } 
}

int  doublecmp(const void *pa , const void *pb) {
 double a=(*(double *)pa) , b=(*(double *)pb) ;
 if (a <b) return -1 ;
 if (a==b) return  0 ;
 if (a >b) return +1 ;
 return -MAXINT ;
}

// (l,s)->x
void code  (int l,unsigned long s , unsigned long *x , int jmx1, int L0) {
  if (l > L0) *x=(2*s+1)*(1<<(jmx1-l)) + 0 ;
  else        *x=     s *(1<<(jmx1-l)) + 0 ;
}

// x->(l,s) 
void decode(int *l,unsigned long *s , unsigned long x , int jmx1, int L0) {
  int q ;
  x -=0 ; // MatLab arrays start with 1 instead of 0 ....
  q=1<<(jmx1-L0) ; // initial divisor
 *l=L0           ; // inital level

  if ((x % q)==0) {
   *s=x/q ;
  } else {
   do {
     q/=2 ;
     (*l)++ ;
     (*s)=(x/q-1)/2 ;
    } while((x % q)!=0) ;
  };
}


#ifdef _CLN_
// Integral Sinus function
cl_F Si(cl_F x) {
 int l,L=8 ; // Abruch der Taylorreihe bei ... 
  
 cl_float_format_t prec = cl_float_format(_CLN_DIGITS_);
 cl_F z1=cl_float(0.0,prec) ;
 cl_F z3=cl_float(0.0,prec) ;
 cl_F y =x*x*x*x ;
 cl_RA r ;
 
 for (l=L; l>0; l--) {
   r=4*l+1 ; r=recip(r) ;
   z1 = z1 + cl_float(r,prec) ;

   r=(4*l+1)*4*l*(4*l-1)*(4*l-2) ; r=recip(r) ;
   z1 = z1*y*cl_float(r,prec) ; 

   r=4*l+3 ; r=recip(r) ;
   z3 = z3 + cl_float(r,prec) ;

   r=(4*l+3)*(4*l+2)*(4*l+1)*4*l ; r=recip(r) ;
   z3 = z3*y*cl_float(r,prec) ;
 }
 
 z1=(1+z1)*x ;
 
 r=3 ; r=recip(r) ;
 z3=(cl_float(r,prec) + z3)*x*x*x ;
 
 r=6 ; r=recip(r) ;
 z3=z3*cl_float(r,prec) ;

 return z1-z3 ;
}


// Integral Cosinus function
cl_F Ci(cl_F x) {
 cl_float_format_t prec = cl_float_format(_CLN_DIGITS_);
 return cl_eulerconst(prec) + ln(x) + CCi(x) ;
}

// Integral Cosinus function
cl_F CCi(cl_F x) {
 int l,L=8 ; // Abruch der Taylorreihe bei ... 

 cl_float_format_t prec = cl_float_format(_CLN_DIGITS_);
 cl_F z0=cl_float(0.0,prec) ;
 cl_F z2=cl_float(0.0,prec) ;
 cl_F y =x*x*x*x ;
 cl_RA r ;

 for (l=L; l>0; l--) {
   r=4*l ; r=recip(r) ;
   z0 = z0 + cl_float(r,prec) ; 

   r=(4*l)*(4*l-1)*(4*l-2)*(4*l-3) ; r=recip(r) ;
   z0 = z0*y*cl_float(r,prec) ;

   r=4*l+2 ; r=recip(r) ;
   z2 = z2 + cl_float(r,prec) ;

   r=(4*l+2)*(4*l+1)*(4*l)*(4*l-1) ; r=recip(r) ;
   z2 = z2*y*cl_float(r,prec) ;
 }

 r=2 ; r=recip(r) ;
 z2 = (cl_float(r,prec) + z2)*x*x*cl_float(r,prec) ;

 return z0-z2 ;
}

//
// interfaces fuer double 
//
double Si(double x) {
 cl_float_format_t prec = cl_float_format(_CLN_DIGITS_);
 cl_F cx=cl_float(x,prec) ;
 return cl_double_approx(Si(cx)) ;
}

double Ci(double x) {
 cl_float_format_t prec = cl_float_format(_CLN_DIGITS_);
 cl_F cx=cl_float(x,prec) ;
 return cl_double_approx(Ci(cx)) ;
}

double CCi(double x) {
 cl_float_format_t prec = cl_float_format(_CLN_DIGITS_);
 cl_F cx=cl_float(x,prec) ;
 return cl_double_approx(CCi(cx)) ;
}

#endif

---