mathlib.cpp

Go to the documentation of this file.

#include "mathlib.h"

#include "gelib.h"

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


namespace midaspace {

void shaker (long &n,long *a)
{
  long i,j,k;
  
  j=0;
  for (i=0;i<n;i++) {
    if (a[i]==-1) continue;
    a[j]=a[i];
    if (i!=j) a[i]=-1;
    for (k=i+1;k<n;k++)
      if (a[j]==a[k]) a[k]=-1;
    
    j++;
  }
  n=j;
}


void sort_4 (long *x)
{
  long a;
  
  if (x[0] < x[1])  if (x[1] > x[2])  if (x[1] > x[3]) {  a=x[1];  x[1]=x[3];  x[3]=a;  sort_3(x);  return; }
                                        else           {  a=x[1];  x[1]=x[2];  x[2]=a;  }
                    else              if (x[2] > x[3]) {  a=x[2];  x[2]=x[3];  x[3]=a;  sort_3(x);  return; }
                                      else             {  return; }
  else              if (x[0] > x[2])  if (x[0] > x[3]) {  a=x[0];  x[0]=x[3];  x[3]=a;  sort_3(x);  return; }
                                      else             {  a=x[0];  x[0]=x[2];  x[2]=a;  }
                    else              if (x[2] > x[3]) {  a=x[2];  x[2]=x[3];  x[3]=a;  sort_3(x);  return; }
  
  if (x[0]>x[1]) { a=x[0]; x[0]=x[1]; x[1]=a; }
}


void sort_3 (long *x)
{
  long a;
  
  if (x[0] > x[1]) { if (x[0] > x[2])  if (x[1] > x[2]) {  a=x[0];  x[0]=x[2];              x[2]=a; }
                                       else             {  a=x[0];  x[0]=x[1];  x[1]=x[2];  x[2]=a; }
                     else                               {  a=x[0];  x[0]=x[1];              x[1]=a; } }
  else             { if (x[1] > x[2]) { if (x[0] > x[2]) {  a=x[0];  x[0]=x[2];  x[2]=x[1];  x[1]=a; }
                                        else             {  a=x[1];  x[1]=x[2];              x[2]=a; } } }
}


void sort_2 (long *x)
{
  if (x[0] < x[1])
    return;
  long a = x[0];  x[0] = x[1];  x[1] = a;
}


void solve_3 (const double *m,const double *r,double *l)
{
  double detall,sd[6];
  
  // subdeterminants
  sd[0] = m[4]*m[8] - m[7]*m[5];
  sd[1] = m[3]*m[8] - m[6]*m[5];
  sd[2] = m[3]*m[7] - m[6]*m[4];
  
  sd[3] = m[4]*r[2] - m[7]*r[1];
  sd[4] = r[1]*m[8] - r[2]*m[5];
  sd[5] = m[3]*r[2] - m[6]*r[1];
  
  // total determinant
  detall = m[0] * sd[0] - m[1] * sd[1] + m[2] * sd[2] ;
  
  l[0] = ( r[0] * sd[0] - m[1] * sd[4] - m[2] * sd[3] ) / detall ;
  l[1] = ( m[0] * sd[4] - r[0] * sd[1] + m[2] * sd[5] ) / detall ;
  l[2] = ( m[0] * sd[3] - m[1] * sd[5] + r[0] * sd[2] ) / detall ;
}


void solve_3 (const double **m,const double *r,double *l)
{
  double detall,sd[6];
  
  // subdeterminants
  sd[0] = m[1][1]*m[2][2] - m[2][1]*m[1][2];
  sd[1] = m[1][0]*m[2][2] - m[2][0]*m[1][2];
  sd[2] = m[1][0]*m[2][1] - m[2][0]*m[1][1];
  
  sd[3] = m[1][1]*r[2]    - m[2][1]*r[1]   ;
  sd[4] = r[1]   *m[2][2] - r[2]   *m[1][2];
  sd[5] = m[1][0]*r[2]    - m[2][0]*r[1]   ;
  
  // total determinant
  detall = m[0][0] * sd[0] - m[0][1] * sd[1] + m[0][2] * sd[2] ;
  
  l[0] = ( r[0]    * sd[0] - m[0][1] * sd[4] - m[0][2] * sd[3] ) / detall ;
  l[1] = ( m[0][0] * sd[4] - r[0]    * sd[1] + m[0][2] * sd[5] ) / detall ;
  l[2] = ( m[0][0] * sd[3] - m[0][1] * sd[5] + r[0]    * sd[2] ) / detall ;
}


long solv_2nle (double zero,const double *a,const double *b,const double *c,const double *d,double *x,double *y)
{
  long ii,jj,nr;
  
  if ( nillret(zero,a[0])==0.0 && nillret(zero,a[1])==0.0 )
    return (solv_2le (zero,b,c,d,x[0],y[0]));
  
  if ( nillret(zero,a[0])==0.0 || nillret(zero,a[1])==0.0 ){
    if (nillret(zero,a[1])) { ii=0; jj=1; }
    else                    { ii=1; jj=0; }
    
    if (!nillret(zero,b[ii]))
      if (!nillret(zero,c[ii]))
    if (!nillret(zero,d[ii]))
      return -1;
    else
      return  0;
      else{
    y[0] = -d[ii]/c[ii];
    return solv_1le (zero, a[jj]*y[0]+b[jj] , c[jj]*y[0]+d[jj] ,x[0]);
      }
    
    else
      if (!nillret(zero,c[ii])){
    x[0] = -d[ii]/b[ii];
    return solv_1le (zero, a[jj]*x[0]+c[jj] , b[jj]*x[0]+d[jj] ,y[0]);
      }
      else{
    nr = solv_polynom_2 ( zero*zero , -a[jj]*b[ii] , -a[jj]*d[ii] + c[ii]*b[jj]-c[jj]*b[ii] , c[ii]*d[jj]-c[jj]*d[ii] ,x[0],x[1]);
    if (nr>0)
      y[0] = (-b[ii]*x[0]-d[ii]) / c[ii];
    if (nr>1)
      y[1] = (-b[ii]*x[1]-d[ii]) / c[ii];
    return nr;
      }
    
  }
  
  nr = solv_polynom_2 ( zero*zero , a[0]*b[1]-a[1]*b[0] , a[0]*d[1]-a[1]*d[0] + c[0]*b[1]-c[1]*b[0] , c[0]*d[1]-c[1]*d[0] ,x[0],x[1]);
  
  for (long i=0;i<nr;i++)
    if ((y[i] = nillret(zero,a[0]*x[i]+c[0])))
      y[i] = (-b[0]*x[i]-d[0]) / y[i] ;
    else
      if (nillret(zero,b[0]*x[i]+d[0]))
    //errol;
    return -2;
      else
    return -1;
  
  return nr;
}

long solv_2le (double zero,const double *a,const double *b,const double *c,double &x,double &y)
{
  if (nillret(zero,a[0])){
    // test of linear dependence
    if (!nillret ( zero*zero , a[1]*(a[0]+b[0])-a[0]*(a[1]+b[1]) )) {
      if (nillret ( zero*zero , a[1]*c[0]-a[0]*c[1] ))
    return  0;
      else
    return -1;
    }
  }
  else{
    // test of linear dependence
    if (!nillret ( zero*zero , b[1]*(a[0]+b[0])-b[0]*(a[1]+b[1]) )) {
      if (nillret ( zero*zero , b[1]*c[0]-b[0]*c[1] ) || ( nillret ( zero*zero , b[1]*c[0]-b[0]*c[1] )                                      ))
    return  0;
      else
    if (nillret(zero,b[0])==0.0 && nillret(zero,b[1])==0.0 && nillret(zero,c[0])!=0.0)
      return 0;
    else
      return -1;
    }
  }
  
  if (nillret(zero,a[0])){
    y = (c[0]*a[1]-c[1]*a[0]) / (a[0]*b[1]-a[1]*b[0]);
    x = (-c[0]-b[0]*y) / a[0];
  }
  else{
    y = -c[0]/b[0];
    x = (-c[1]-b[1]*y) / a[1];
  }
  
  return 1;
}

long solv_1le (double zero,double a,double b,double &x)
{
  if (!nillret(zero,a)) {
    if (!nillret(zero,b))
      return -1;
    else
      return 0;
  }
  x = -b/a;
  
  return 1;
}

void nilling (double zero, double &a)
{
  if (-zero<a && a<zero) a = 0.0;
}

double nillret (double zero, double a)
{
  if (-zero<a && a<zero) return 0.0;
  else return a;
}

bool isZero (double zero, double a)
{
  if (a == 0.0)          return true;
  if (-zero<a && a<zero) return true;
  else                   return false;
}
bool isZero (double abszero, double relzero, double a)
{
  return isZero (relzero > abszero ? relzero : abszero, a);
}

long solv_polynom_2 (double zero,double a,double b,double c,double &r1,double &r2)
{
  if (!nillret(zero,a)) {
    if (!nillret(zero,b))
      if (!nillret(zero,c))
    return -1;
      else
    return 0;
    else {
      r1 = -c/b;
      return 1;
    }
  }
  
  double D = nillret ( zero*zero , b*b - 4*a*c );
  
  if (D < 0.0)
    return 0;
  
  if (D == 0.0){
    r1 = -b/2.0/a;
    return 1;
  }
  
  r1 = (-b + sqrt(D))/2.0/a;
  r2 = (-b - sqrt(D))/2.0/a;
  return 2;
}

long div_dd (double d1, double d2, const char* s1, const char* s2, const long line)
{
  long answer = lround (d1/d2);
  if ( fabs(answer-d1/d2) > 0.00001 ) { fprintf (stderr,"%s = %f is not dividable by %s = %f (%f), line %ld\n",s1,d1,s2,d2,d1/d2,line); abort (); }
  return answer;
}

long div_dd (long &answer, double d1, double d2)
{
  answer = lround (d1/d2);
  if ( fabs(answer-d1/d2) > 0.00001 ) return 1;
  return 0;
}

long div_dd (int &answer, double d1, double d2)
{
  long aux;
  int ret = div_dd (aux, d1, d2);
  answer = (int) aux;
  
  return ret;
}


int       decomp_int (int answer[], int n, int       l, int       rad)
{
  for (int i=n; i>0; i--) { answer[i-1] = l%rad;  l = l/rad; }
  return l;
}
long      decomp_int (int answer[], int n, long      l, long      rad)
{
  for (int i=n; i>0; i--) { answer[i-1] = l%rad;  l = l/rad; }
  return l;
}
long long decomp_int (int answer[], int n, long long l, long long rad)
{
  for (int i=n; i>0; i--) { answer[i-1] = l%rad;  l = l/rad; }
  return l;
}


bool isNonZero (double x, double tolerance)
{
  if ( fabs(x) > tolerance )  return 1;
  else                        return 0;
}


} // namespace midaspace