1  /*************************************************************************


2  AP library


3  Copyright (c) 20032009 Sergey Bochkanov (ALGLIB project).


4 


5  >>> SOURCE LICENSE >>>


6  This program is free software; you can redistribute it and/or modify


7  it under the terms of the GNU General Public License as published by


8  the Free Software Foundation (www.fsf.org); either version 2 of the


9  License, or (at your option) any later version.


10 


11  This program is distributed in the hope that it will be useful,


12  but WITHOUT ANY WARRANTY; without even the implied warranty of


13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


14  GNU General Public License for more details.


15 


16  A copy of the GNU General Public License is available at


17  http://www.fsf.org/licensing/licenses


18  >>> END OF LICENSE >>>


19  *************************************************************************/


20  using System;


21  public partial class alglib


22  {


23  /********************************************************************


24  Callback definitions for optimizers/fitters/solvers.


25 


26  Callbacks for unparameterized (general) functions:


27  * ndimensional_func calculates f(arg), stores result to func


28  * ndimensional_grad calculates func = f(arg),


29  grad[i] = df(arg)/d(arg[i])


30  * ndimensional_hess calculates func = f(arg),


31  grad[i] = df(arg)/d(arg[i]),


32  hess[i,j] = d2f(arg)/(d(arg[i])*d(arg[j]))


33 


34  Callbacks for systems of functions:


35  * ndimensional_fvec calculates vector function f(arg),


36  stores result to fi


37  * ndimensional_jac calculates f[i] = fi(arg)


38  jac[i,j] = df[i](arg)/d(arg[j])


39 


40  Callbacks for parameterized functions, i.e. for functions which


41  depend on two vectors: P and Q. Gradient and Hessian are calculated


42  with respect to P only.


43  * ndimensional_pfunc calculates f(p,q),


44  stores result to func


45  * ndimensional_pgrad calculates func = f(p,q),


46  grad[i] = df(p,q)/d(p[i])


47  * ndimensional_phess calculates func = f(p,q),


48  grad[i] = df(p,q)/d(p[i]),


49  hess[i,j] = d2f(p,q)/(d(p[i])*d(p[j]))


50 


51  Callbacks for progress reports:


52  * ndimensional_rep reports current position of optimization algo


53 


54  Callbacks for ODE solvers:


55  * ndimensional_ode_rp calculates dy/dx for given y[] and x


56 


57  Callbacks for integrators:


58  * integrator1_func calculates f(x) for given x


59  (additional parameters xminusa and bminusx


60  contain xa and bx)


61  ********************************************************************/


62  public delegate void ndimensional_func (double[] arg, ref double func, object obj);


63  public delegate void ndimensional_grad (double[] arg, ref double func, double[] grad, object obj);


64  public delegate void ndimensional_hess (double[] arg, ref double func, double[] grad, double[,] hess, object obj);


65 


66  public delegate void ndimensional_fvec (double[] arg, double[] fi, object obj);


67  public delegate void ndimensional_jac (double[] arg, double[] fi, double[,] jac, object obj);


68 


69  public delegate void ndimensional_pfunc(double[] p, double[] q, ref double func, object obj);


70  public delegate void ndimensional_pgrad(double[] p, double[] q, ref double func, double[] grad, object obj);


71  public delegate void ndimensional_phess(double[] p, double[] q, ref double func, double[] grad, double[,] hess, object obj);


72 


73  public delegate void ndimensional_rep(double[] arg, double func, object obj);


74 


75  public delegate void ndimensional_ode_rp (double[] y, double x, double[] dy, object obj);


76 


77  public delegate void integrator1_func (double x, double xminusa, double bminusx, ref double f, object obj);


78 


79  /********************************************************************


80  Class defining a complex number with double precision.


81  ********************************************************************/


82  public struct complex


83  {


84  public double x;


85  public double y;


86 


87  public complex(double _x)


88  {


89  x = _x;


90  y = 0;


91  }


92  public complex(double _x, double _y)


93  {


94  x = _x;


95  y = _y;


96  }


97  public static implicit operator complex(double _x)


98  {


99  return new complex(_x);


100  }


101  public static bool operator==(complex lhs, complex rhs)


102  {


103  return ((double)lhs.x==(double)rhs.x) & ((double)lhs.y==(double)rhs.y);


104  }


105  public static bool operator!=(complex lhs, complex rhs)


106  {


107  return ((double)lhs.x!=(double)rhs.x)  ((double)lhs.y!=(double)rhs.y);


108  }


109  public static complex operator+(complex lhs)


110  {


111  return lhs;


112  }


113  public static complex operator(complex lhs)


114  {


115  return new complex(lhs.x,lhs.y);


116  }


117  public static complex operator+(complex lhs, complex rhs)


118  {


119  return new complex(lhs.x+rhs.x,lhs.y+rhs.y);


120  }


121  public static complex operator(complex lhs, complex rhs)


122  {


123  return new complex(lhs.xrhs.x,lhs.yrhs.y);


124  }


125  public static complex operator*(complex lhs, complex rhs)


126  {


127  return new complex(lhs.x*rhs.xlhs.y*rhs.y, lhs.x*rhs.y+lhs.y*rhs.x);


128  }


129  public static complex operator/(complex lhs, complex rhs)


130  {


131  complex result;


132  double e;


133  double f;


134  if( System.Math.Abs(rhs.y)<System.Math.Abs(rhs.x) )


135  {


136  e = rhs.y/rhs.x;


137  f = rhs.x+rhs.y*e;


138  result.x = (lhs.x+lhs.y*e)/f;


139  result.y = (lhs.ylhs.x*e)/f;


140  }


141  else


142  {


143  e = rhs.x/rhs.y;


144  f = rhs.y+rhs.x*e;


145  result.x = (lhs.y+lhs.x*e)/f;


146  result.y = (lhs.x+lhs.y*e)/f;


147  }


148  return result;


149  }


150  public override int GetHashCode()


151  {


152  return x.GetHashCode() ^ y.GetHashCode();


153  }


154  public override bool Equals(object obj)


155  {


156  if( obj is byte)


157  return Equals(new complex((byte)obj));


158  if( obj is sbyte)


159  return Equals(new complex((sbyte)obj));


160  if( obj is short)


161  return Equals(new complex((short)obj));


162  if( obj is ushort)


163  return Equals(new complex((ushort)obj));


164  if( obj is int)


165  return Equals(new complex((int)obj));


166  if( obj is uint)


167  return Equals(new complex((uint)obj));


168  if( obj is long)


169  return Equals(new complex((long)obj));


170  if( obj is ulong)


171  return Equals(new complex((ulong)obj));


172  if( obj is float)


173  return Equals(new complex((float)obj));


174  if( obj is double)


175  return Equals(new complex((double)obj));


176  if( obj is decimal)


177  return Equals(new complex((double)(decimal)obj));


178  return base.Equals(obj);


179  }


180  }


181 


182  /********************************************************************


183  Class defining an ALGLIB exception


184  ********************************************************************/


185  public class alglibexception : System.Exception


186  {


187  public string msg;


188  public alglibexception(string s)


189  {


190  msg = s;


191  }


192 


193  }


194 


195  /********************************************************************


196  reverse communication structure


197  ********************************************************************/


198  public class rcommstate


199  {


200  public rcommstate()


201  {


202  stage = 1;


203  ia = new int[0];


204  ba = new bool[0];


205  ra = new double[0];


206  ca = new alglib.complex[0];


207  }


208  public int stage;


209  public int[] ia;


210  public bool[] ba;


211  public double[] ra;


212  public alglib.complex[] ca;


213  };


214 


215  /********************************************************************


216  internal functions


217  ********************************************************************/


218  public class ap


219  {


220  public static int len<T>(T[] a)


221  { return a.Length; }


222  public static int rows<T>(T[,] a)


223  { return a.GetLength(0); }


224  public static int cols<T>(T[,] a)


225  { return a.GetLength(1); }


226  public static void swap<T>(ref T a, ref T b)


227  {


228  T t = a;


229  a = b;


230  b = t;


231  }


232 


233  public static void assert(bool cond, string s)


234  {


235  if( !cond )


236  throw new alglibexception(s);


237  }


238 


239  public static void assert(bool cond)


240  {


241  assert(cond, "ALGLIB: assertion failed");


242  }


243 


244  /****************************************************************


245  returns dps (digitsofprecision) value corresponding to threshold.


246  dps(0.9) = dps(0.5) = dps(0.1) = 0


247  dps(0.09) = dps(0.05) = dps(0.01) = 1


248  and so on


249  ****************************************************************/


250  public static int threshold2dps(double threshold)


251  {


252  int result = 0;


253  double t;


254  for (result = 0, t = 1; t / 10 > threshold*(1+1E10); result++, t /= 10) ;


255  return result;


256  }


257 


258  /****************************************************************


259  prints formatted complex


260  ****************************************************************/


261  public static string format(complex a, int _dps)


262  {


263  int dps = Math.Abs(_dps);


264  string fmt = _dps>=0 ? "F" : "E";


265  string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps);


266  string fmty = String.Format("{{0:"+fmt+"{0}}}", dps);


267  string result = String.Format(fmtx, a.x) + (a.y >= 0 ? "+" : "") + String.Format(fmty, Math.Abs(a.y)) + "i";


268  result = result.Replace(',', '.');


269  return result;


270  }


271 


272  /****************************************************************


273  prints formatted array


274  ****************************************************************/


275  public static string format(bool[] a)


276  {


277  string[] result = new string[len(a)];


278  int i;


279  for(i=0; i<len(a); i++)


280  if( a[i] )


281  result[i] = "true";


282  else


283  result[i] = "false";


284  return "{"+String.Join(",",result)+"}";


285  }


286 


287  /****************************************************************


288  prints formatted array


289  ****************************************************************/


290  public static string format(int[] a)


291  {


292  string[] result = new string[len(a)];


293  int i;


294  for (i = 0; i < len(a); i++)


295  result[i] = a[i].ToString();


296  return "{" + String.Join(",", result) + "}";


297  }


298 


299  /****************************************************************


300  prints formatted array


301  ****************************************************************/


302  public static string format(double[] a, int _dps)


303  {


304  int dps = Math.Abs(_dps);


305  string sfmt = _dps >= 0 ? "F" : "E";


306  string fmt = String.Format("{{0:" + sfmt + "{0}}}", dps);


307  string[] result = new string[len(a)];


308  int i;


309  for (i = 0; i < len(a); i++)


310  {


311  result[i] = String.Format(fmt, a[i]);


312  result[i] = result[i].Replace(',', '.');


313  }


314  return "{" + String.Join(",", result) + "}";


315  }


316 


317  /****************************************************************


318  prints formatted array


319  ****************************************************************/


320  public static string format(complex[] a, int _dps)


321  {


322  int dps = Math.Abs(_dps);


323  string fmt = _dps >= 0 ? "F" : "E";


324  string fmtx = String.Format("{{0:"+fmt+"{0}}}", dps);


325  string fmty = String.Format("{{0:"+fmt+"{0}}}", dps);


326  string[] result = new string[len(a)];


327  int i;


328  for (i = 0; i < len(a); i++)


329  {


330  result[i] = String.Format(fmtx, a[i].x) + (a[i].y >= 0 ? "+" : "") + String.Format(fmty, Math.Abs(a[i].y)) + "i";


331  result[i] = result[i].Replace(',', '.');


332  }


333  return "{" + String.Join(",", result) + "}";


334  }


335 


336  /****************************************************************


337  prints formatted matrix


338  ****************************************************************/


339  public static string format(bool[,] a)


340  {


341  int i, j, m, n;


342  n = cols(a);


343  m = rows(a);


344  bool[] line = new bool[n];


345  string[] result = new string[m];


346  for (i = 0; i < m; i++)


347  {


348  for (j = 0; j < n; j++)


349  line[j] = a[i, j];


350  result[i] = format(line);


351  }


352  return "{" + String.Join(",", result) + "}";


353  }


354 


355  /****************************************************************


356  prints formatted matrix


357  ****************************************************************/


358  public static string format(int[,] a)


359  {


360  int i, j, m, n;


361  n = cols(a);


362  m = rows(a);


363  int[] line = new int[n];


364  string[] result = new string[m];


365  for (i = 0; i < m; i++)


366  {


367  for (j = 0; j < n; j++)


368  line[j] = a[i, j];


369  result[i] = format(line);


370  }


371  return "{" + String.Join(",", result) + "}";


372  }


373 


374  /****************************************************************


375  prints formatted matrix


376  ****************************************************************/


377  public static string format(double[,] a, int dps)


378  {


379  int i, j, m, n;


380  n = cols(a);


381  m = rows(a);


382  double[] line = new double[n];


383  string[] result = new string[m];


384  for (i = 0; i < m; i++)


385  {


386  for (j = 0; j < n; j++)


387  line[j] = a[i, j];


388  result[i] = format(line, dps);


389  }


390  return "{" + String.Join(",", result) + "}";


391  }


392 


393  /****************************************************************


394  prints formatted matrix


395  ****************************************************************/


396  public static string format(complex[,] a, int dps)


397  {


398  int i, j, m, n;


399  n = cols(a);


400  m = rows(a);


401  complex[] line = new complex[n];


402  string[] result = new string[m];


403  for (i = 0; i < m; i++)


404  {


405  for (j = 0; j < n; j++)


406  line[j] = a[i, j];


407  result[i] = format(line, dps);


408  }


409  return "{" + String.Join(",", result) + "}";


410  }


411 


412  /****************************************************************


413  checks that matrix is symmetric.


414  maxAA^T is calculated; if it is within 1.0E14 of maxA,


415  matrix is considered symmetric


416  ****************************************************************/


417  public static bool issymmetric(double[,] a)


418  {


419  int i, j, n;


420  double err, mx, v1, v2;


421  if( rows(a)!=cols(a) )


422  return false;


423  n = rows(a);


424  if( n==0 )


425  return true;


426  mx = 0;


427  err = 0;


428  for( i=0; i<n; i++)


429  {


430  for(j=i+1; j<n; j++)


431  {


432  v1 = a[i,j];


433  v2 = a[j,i];


434  if( !math.isfinite(v1) )


435  return false;


436  if( !math.isfinite(v2) )


437  return false;


438  err = Math.Max(err, Math.Abs(v1v2));


439  mx = Math.Max(mx, Math.Abs(v1));


440  mx = Math.Max(mx, Math.Abs(v2));


441  }


442  v1 = a[i,i];


443  if( !math.isfinite(v1) )


444  return false;


445  mx = Math.Max(mx, Math.Abs(v1));


446  }


447  if( mx==0 )


448  return true;


449  return err/mx<=1.0E14;


450  }


451 


452  /****************************************************************


453  checks that matrix is Hermitian.


454  maxAA^H is calculated; if it is within 1.0E14 of maxA,


455  matrix is considered Hermitian


456  ****************************************************************/


457  public static bool ishermitian(complex[,] a)


458  {


459  int i, j, n;


460  double err, mx;


461  complex v1, v2, vt;


462  if( rows(a)!=cols(a) )


463  return false;


464  n = rows(a);


465  if( n==0 )


466  return true;


467  mx = 0;


468  err = 0;


469  for( i=0; i<n; i++)


470  {


471  for(j=i+1; j<n; j++)


472  {


473  v1 = a[i,j];


474  v2 = a[j,i];


475  if( !math.isfinite(v1.x) )


476  return false;


477  if( !math.isfinite(v1.y) )


478  return false;


479  if( !math.isfinite(v2.x) )


480  return false;


481  if( !math.isfinite(v2.y) )


482  return false;


483  vt.x = v1.xv2.x;


484  vt.y = v1.y+v2.y;


485  err = Math.Max(err, math.abscomplex(vt));


486  mx = Math.Max(mx, math.abscomplex(v1));


487  mx = Math.Max(mx, math.abscomplex(v2));


488  }


489  v1 = a[i,i];


490  if( !math.isfinite(v1.x) )


491  return false;


492  if( !math.isfinite(v1.y) )


493  return false;


494  err = Math.Max(err, Math.Abs(v1.y));


495  mx = Math.Max(mx, math.abscomplex(v1));


496  }


497  if( mx==0 )


498  return true;


499  return err/mx<=1.0E14;


500  }


501 


502 


503  /****************************************************************


504  Forces symmetricity by copying upper half of A to the lower one


505  ****************************************************************/


506  public static bool forcesymmetric(double[,] a)


507  {


508  int i, j, n;


509  if( rows(a)!=cols(a) )


510  return false;


511  n = rows(a);


512  if( n==0 )


513  return true;


514  for( i=0; i<n; i++)


515  for(j=i+1; j<n; j++)


516  a[i,j] = a[j,i];


517  return true;


518  }


519 


520  /****************************************************************


521  Forces Hermiticity by copying upper half of A to the lower one


522  ****************************************************************/


523  public static bool forcehermitian(complex[,] a)


524  {


525  int i, j, n;


526  complex v;


527  if( rows(a)!=cols(a) )


528  return false;


529  n = rows(a);


530  if( n==0 )


531  return true;


532  for( i=0; i<n; i++)


533  for(j=i+1; j<n; j++)


534  {


535  v = a[j,i];


536  a[i,j].x = v.x;


537  a[i,j].y = v.y;


538  }


539  return true;


540  }


541  };


542 


543  /********************************************************************


544  math functions


545  ********************************************************************/


546  public class math


547  {


548  //public static System.Random RndObject = new System.Random(System.DateTime.Now.Millisecond);


549  public static System.Random rndobject = new System.Random(System.DateTime.Now.Millisecond + 1000*System.DateTime.Now.Second + 60*1000*System.DateTime.Now.Minute);


550 


551  public const double machineepsilon = 5E16;


552  public const double maxrealnumber = 1E300;


553  public const double minrealnumber = 1E300;


554 


555  public static bool isfinite(double d)


556  {


557  return !System.Double.IsNaN(d) && !System.Double.IsInfinity(d);


558  }


559 


560  public static double randomreal()


561  {


562  double r = 0;


563  lock(rndobject){ r = rndobject.NextDouble(); }


564  return r;


565  }


566  public static int randominteger(int N)


567  {


568  int r = 0;


569  lock(rndobject){ r = rndobject.Next(N); }


570  return r;


571  }


572  public static double sqr(double X)


573  {


574  return X*X;


575  }


576  public static double abscomplex(complex z)


577  {


578  double w;


579  double xabs;


580  double yabs;


581  double v;


582 


583  xabs = System.Math.Abs(z.x);


584  yabs = System.Math.Abs(z.y);


585  w = xabs>yabs ? xabs : yabs;


586  v = xabs<yabs ? xabs : yabs;


587  if( v==0 )


588  return w;


589  else


590  {


591  double t = v/w;


592  return w*System.Math.Sqrt(1+t*t);


593  }


594  }


595  public static complex conj(complex z)


596  {


597  return new complex(z.x, z.y);


598  }


599  public static complex csqr(complex z)


600  {


601  return new complex(z.x*z.xz.y*z.y, 2*z.x*z.y);


602  }


603 


604  }


605 


606  /********************************************************************


607  serializer object (should not be used directly)


608  ********************************************************************/


609  public class serializer


610  {


611  enum SMODE { DEFAULT, ALLOC, TO_STRING, FROM_STRING };


612  private const int SER_ENTRIES_PER_ROW = 5;


613  private const int SER_ENTRY_LENGTH = 11;


614 


615  private SMODE mode;


616  private int entries_needed;


617  private int entries_saved;


618  private int bytes_asked;


619  private int bytes_written;


620  private int bytes_read;


621  private char[] out_str;


622  private char[] in_str;


623 


624  public serializer()


625  {


626  mode = SMODE.DEFAULT;


627  entries_needed = 0;


628  bytes_asked = 0;


629  }


630 


631  public void alloc_start()


632  {


633  entries_needed = 0;


634  bytes_asked = 0;


635  mode = SMODE.ALLOC;


636  }


637 


638  public void alloc_entry()


639  {


640  if( mode!=SMODE.ALLOC )


641  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


642  entries_needed++;


643  }


644 


645  private int get_alloc_size()


646  {


647  int rows, lastrowsize, result;


648 


649  // check and change mode


650  if( mode!=SMODE.ALLOC )


651  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


652 


653  // if no entries needes (degenerate case)


654  if( entries_needed==0 )


655  {


656  bytes_asked = 1;


657  return bytes_asked;


658  }


659 


660  // nondegenerate case


661  rows = entries_needed/SER_ENTRIES_PER_ROW;


662  lastrowsize = SER_ENTRIES_PER_ROW;


663  if( entries_needed%SER_ENTRIES_PER_ROW!=0 )


664  {


665  lastrowsize = entries_needed%SER_ENTRIES_PER_ROW;


666  rows++;


667  }


668 


669  // calculate result size


670  result = ((rows1)*SER_ENTRIES_PER_ROW+lastrowsize)*SER_ENTRY_LENGTH;


671  result += (rows1)*(SER_ENTRIES_PER_ROW1)+(lastrowsize1);


672  result += rows*2;


673  bytes_asked = result;


674  return result;


675  }


676 


677  public void sstart_str()


678  {


679  int allocsize = get_alloc_size();


680 


681  // check and change mode


682  if( mode!=SMODE.ALLOC )


683  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


684  mode = SMODE.TO_STRING;


685 


686  // other preparations


687  out_str = new char[allocsize];


688  entries_saved = 0;


689  bytes_written = 0;


690  }


691 


692  public void ustart_str(string s)


693  {


694  // check and change mode


695  if( mode!=SMODE.DEFAULT )


696  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


697  mode = SMODE.FROM_STRING;


698 


699  in_str = s.ToCharArray();


700  bytes_read = 0;


701  }


702 


703  public void serialize_bool(bool v)


704  {


705  if( mode!=SMODE.TO_STRING )


706  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


707  bool2str(v, out_str, ref bytes_written);


708  entries_saved++;


709  if( entries_saved%SER_ENTRIES_PER_ROW!=0 )


710  {


711  out_str[bytes_written] = ' ';


712  bytes_written++;


713  }


714  else


715  {


716  out_str[bytes_written+0] = '\r';


717  out_str[bytes_written+1] = '\n';


718  bytes_written+=2;


719  }


720  }


721 


722  public void serialize_int(int v)


723  {


724  if( mode!=SMODE.TO_STRING )


725  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


726  int2str(v, out_str, ref bytes_written);


727  entries_saved++;


728  if( entries_saved%SER_ENTRIES_PER_ROW!=0 )


729  {


730  out_str[bytes_written] = ' ';


731  bytes_written++;


732  }


733  else


734  {


735  out_str[bytes_written+0] = '\r';


736  out_str[bytes_written+1] = '\n';


737  bytes_written+=2;


738  }


739  }


740 


741  public void serialize_double(double v)


742  {


743  if( mode!=SMODE.TO_STRING )


744  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


745  double2str(v, out_str, ref bytes_written);


746  entries_saved++;


747  if( entries_saved%SER_ENTRIES_PER_ROW!=0 )


748  {


749  out_str[bytes_written] = ' ';


750  bytes_written++;


751  }


752  else


753  {


754  out_str[bytes_written+0] = '\r';


755  out_str[bytes_written+1] = '\n';


756  bytes_written+=2;


757  }


758  }


759 


760  public bool unserialize_bool()


761  {


762  if( mode!=SMODE.FROM_STRING )


763  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


764  return str2bool(in_str, ref bytes_read);


765  }


766 


767  public int unserialize_int()


768  {


769  if( mode!=SMODE.FROM_STRING )


770  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


771  return str2int(in_str, ref bytes_read);


772  }


773 


774  public double unserialize_double()


775  {


776  if( mode!=SMODE.FROM_STRING )


777  throw new alglib.alglibexception("ALGLIB: internal error during (un)serialization");


778  return str2double(in_str, ref bytes_read);


779  }


780 


781  public void stop()


782  {


783  }


784 


785  public string get_string()


786  {


787  return new string(out_str, 0, bytes_written);


788  }


789 


790 


791  /************************************************************************


792  This function converts sixbit value (from 0 to 63) to character (only


793  digits, lowercase and uppercase letters, minus and underscore are used).


794 


795  If v is negative or greater than 63, this function returns '?'.


796  ************************************************************************/


797  private static char[] _sixbits2char_tbl = new char[64]{


798  '0', '1', '2', '3', '4', '5', '6', '7',


799  '8', '9', 'A', 'B', 'C', 'D', 'E', 'F',


800  'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N',


801  'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V',


802  'W', 'X', 'Y', 'Z', 'a', 'b', 'c', 'd',


803  'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l',


804  'm', 'n', 'o', 'p', 'q', 'r', 's', 't',


805  'u', 'v', 'w', 'x', 'y', 'z', '', '_' };


806  private static char sixbits2char(int v)


807  {


808  if( v<0  v>63 )


809  return '?';


810  return _sixbits2char_tbl[v];


811  }


812 


813  /************************************************************************


814  This function converts character to sixbit value (from 0 to 63).


815 


816  This function is inverse of ae_sixbits2char()


817  If c is not correct character, this function returns 1.


818  ************************************************************************/


819  private static int[] _char2sixbits_tbl = new int[128] {


820  1, 1, 1, 1, 1, 1, 1, 1,


821  1, 1, 1, 1, 1, 1, 1, 1,


822  1, 1, 1, 1, 1, 1, 1, 1,


823  1, 1, 1, 1, 1, 1, 1, 1,


824  1, 1, 1, 1, 1, 1, 1, 1,


825  1, 1, 1, 1, 1, 62, 1, 1,


826  0, 1, 2, 3, 4, 5, 6, 7,


827  8, 9, 1, 1, 1, 1, 1, 1,


828  1, 10, 11, 12, 13, 14, 15, 16,


829  17, 18, 19, 20, 21, 22, 23, 24,


830  25, 26, 27, 28, 29, 30, 31, 32,


831  33, 34, 35, 1, 1, 1, 1, 63,


832  1, 36, 37, 38, 39, 40, 41, 42,


833  43, 44, 45, 46, 47, 48, 49, 50,


834  51, 52, 53, 54, 55, 56, 57, 58,


835  59, 60, 61, 1, 1, 1, 1, 1 };


836  private static int char2sixbits(char c)


837  {


838  return (c>=0 && c<127) ? _char2sixbits_tbl[c] : 1;


839  }


840 


841  /************************************************************************


842  This function converts three bytes (24 bits) to four sixbit values


843  (24 bits again).


844 


845  src array


846  src_offs offset of threebytes chunk


847  dst array for ints


848  dst_offs offset of fourints chunk


849  ************************************************************************/


850  private static void threebytes2foursixbits(byte[] src, int src_offs, int[] dst, int dst_offs)


851  {


852  dst[dst_offs+0] = src[src_offs+0] & 0x3F;


853  dst[dst_offs+1] = (src[src_offs+0]>>6)  ((src[src_offs+1]&0x0F)<<2);


854  dst[dst_offs+2] = (src[src_offs+1]>>4)  ((src[src_offs+2]&0x03)<<4);


855  dst[dst_offs+3] = src[src_offs+2]>>2;


856  }


857 


858  /************************************************************************


859  This function converts four sixbit values (24 bits) to three bytes


860  (24 bits again).


861 


862  src pointer to four ints


863  src_offs offset of the chunk


864  dst pointer to three bytes


865  dst_offs offset of the chunk


866  ************************************************************************/


867  private static void foursixbits2threebytes(int[] src, int src_offs, byte[] dst, int dst_offs)


868  {


869  dst[dst_offs+0] = (byte)(src[src_offs+0]  ((src[src_offs+1]&0x03)<<6));


870  dst[dst_offs+1] = (byte)((src[src_offs+1]>>2)  ((src[src_offs+2]&0x0F)<<4));


871  dst[dst_offs+2] = (byte)((src[src_offs+2]>>4)  (src[src_offs+3]<<2));


872  }


873 


874  /************************************************************************


875  This function serializes boolean value into buffer


876 


877  v boolean value to be serialized


878  buf buffer, at least 11 characters wide


879  offs offset in the buffer


880 


881  after return from this function, offs points to the char's past the value


882  being read.


883  ************************************************************************/


884  private static void bool2str(bool v, char[] buf, ref int offs)


885  {


886  char c = v ? '1' : '0';


887  int i;


888  for(i=0; i<SER_ENTRY_LENGTH; i++)


889  buf[offs+i] = c;


890  offs += SER_ENTRY_LENGTH;


891  }


892 


893  /************************************************************************


894  This function unserializes boolean value from buffer


895 


896  buf buffer which contains value; leading spaces/tabs/newlines are


897  ignored, traling spaces/tabs/newlines are treated as end of


898  the boolean value.


899  offs offset in the buffer


900 


901  after return from this function, offs points to the char's past the value


902  being read.


903 


904  This function raises an error in case unexpected symbol is found


905  ************************************************************************/


906  private static bool str2bool(char[] buf, ref int offs)


907  {


908  bool was0, was1;


909  string emsg = "ALGLIB: unable to read boolean value from stream";


910 


911  was0 = false;


912  was1 = false;


913  while( buf[offs]==' '  buf[offs]=='\t'  buf[offs]=='\n'  buf[offs]=='\r' )


914  offs++;


915  while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )


916  {


917  if( buf[offs]=='0' )


918  {


919  was0 = true;


920  offs++;


921  continue;


922  }


923  if( buf[offs]=='1' )


924  {


925  was1 = true;


926  offs++;


927  continue;


928  }


929  throw new alglib.alglibexception(emsg);


930  }


931  if( (!was0) && (!was1) )


932  throw new alglib.alglibexception(emsg);


933  if( was0 && was1 )


934  throw new alglib.alglibexception(emsg);


935  return was1 ? true : false;


936  }


937 


938  /************************************************************************


939  This function serializes integer value into buffer


940 


941  v integer value to be serialized


942  buf buffer, at least 11 characters wide


943  offs offset in the buffer


944 


945  after return from this function, offs points to the char's past the value


946  being read.


947 


948  This function raises an error in case unexpected symbol is found


949  ************************************************************************/


950  private static void int2str(int v, char[] buf, ref int offs)


951  {


952  int i;


953  byte[] _bytes = System.BitConverter.GetBytes((int)v);


954  byte[] bytes = new byte[9];


955  int[] sixbits = new int[12];


956  byte c;


957 


958  //


959  // copy v to array of bytes, sign extending it and


960  // converting to little endian order. Additionally,


961  // we set 9th byte to zero in order to simplify


962  // conversion to sixbit representation


963  //


964  if( !System.BitConverter.IsLittleEndian )


965  System.Array.Reverse(_bytes);


966  c = v<0 ? (byte)0xFF : (byte)0x00;


967  for(i=0; i<sizeof(int); i++)


968  bytes[i] = _bytes[i];


969  for(i=sizeof(int); i<8; i++)


970  bytes[i] = c;


971  bytes[8] = 0;


972 


973  //


974  // convert to sixbit representation, output


975  //


976  // NOTE: last 12th element of sixbits is always zero, we do not output it


977  //


978  threebytes2foursixbits(bytes, 0, sixbits, 0);


979  threebytes2foursixbits(bytes, 3, sixbits, 4);


980  threebytes2foursixbits(bytes, 6, sixbits, 8);


981  for(i=0; i<SER_ENTRY_LENGTH; i++)


982  buf[offs+i] = sixbits2char(sixbits[i]);


983  offs += SER_ENTRY_LENGTH;


984  }


985 


986  /************************************************************************


987  This function unserializes integer value from string


988 


989  buf buffer which contains value; leading spaces/tabs/newlines are


990  ignored, traling spaces/tabs/newlines are treated as end of


991  the integer value.


992  offs offset in the buffer


993 


994  after return from this function, offs points to the char's past the value


995  being read.


996 


997  This function raises an error in case unexpected symbol is found


998  ************************************************************************/


999  private static int str2int(char[] buf, ref int offs)


1000  {


1001  string emsg = "ALGLIB: unable to read integer value from stream";


1002  string emsg3264 = "ALGLIB: unable to read integer value from stream (value does not fit into 32 bits)";


1003  int[] sixbits = new int[12];


1004  byte[] bytes = new byte[9];


1005  byte[] _bytes = new byte[sizeof(int)];


1006  int sixbitsread, i;


1007  byte c;


1008 


1009  //


1010  // 1. skip leading spaces


1011  // 2. read and decode sixbit digits


1012  // 3. set trailing digits to zeros


1013  // 4. convert to little endian 64bit integer representation


1014  // 5. check that we fit into int


1015  // 6. convert to big endian representation, if needed


1016  //


1017  sixbitsread = 0;


1018  while( buf[offs]==' '  buf[offs]=='\t'  buf[offs]=='\n'  buf[offs]=='\r' )


1019  offs++;


1020  while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )


1021  {


1022  int d;


1023  d = char2sixbits(buf[offs]);


1024  if( d<0  sixbitsread>=SER_ENTRY_LENGTH )


1025  throw new alglib.alglibexception(emsg);


1026  sixbits[sixbitsread] = d;


1027  sixbitsread++;


1028  offs++;


1029  }


1030  if( sixbitsread==0 )


1031  throw new alglib.alglibexception(emsg);


1032  for(i=sixbitsread; i<12; i++)


1033  sixbits[i] = 0;


1034  foursixbits2threebytes(sixbits, 0, bytes, 0);


1035  foursixbits2threebytes(sixbits, 4, bytes, 3);


1036  foursixbits2threebytes(sixbits, 8, bytes, 6);


1037  c = (bytes[sizeof(int)1] & 0x80)!=0 ? (byte)0xFF : (byte)0x00;


1038  for(i=sizeof(int); i<8; i++)


1039  if( bytes[i]!=c )


1040  throw new alglib.alglibexception(emsg3264);


1041  for(i=0; i<sizeof(int); i++)


1042  _bytes[i] = bytes[i];


1043  if( !System.BitConverter.IsLittleEndian )


1044  System.Array.Reverse(_bytes);


1045  return System.BitConverter.ToInt32(_bytes,0);


1046  }


1047 


1048 


1049  /************************************************************************


1050  This function serializes double value into buffer


1051 


1052  v double value to be serialized


1053  buf buffer, at least 11 characters wide


1054  offs offset in the buffer


1055 


1056  after return from this function, offs points to the char's past the value


1057  being read.


1058  ************************************************************************/


1059  private static void double2str(double v, char[] buf, ref int offs)


1060  {


1061  int i;


1062  int[] sixbits = new int[12];


1063  byte[] bytes = new byte[9];


1064 


1065  //


1066  // handle special quantities


1067  //


1068  if( System.Double.IsNaN(v) )


1069  {


1070  buf[offs+0] = '.';


1071  buf[offs+1] = 'n';


1072  buf[offs+2] = 'a';


1073  buf[offs+3] = 'n';


1074  buf[offs+4] = '_';


1075  buf[offs+5] = '_';


1076  buf[offs+6] = '_';


1077  buf[offs+7] = '_';


1078  buf[offs+8] = '_';


1079  buf[offs+9] = '_';


1080  buf[offs+10] = '_';


1081  offs += SER_ENTRY_LENGTH;


1082  return;


1083  }


1084  if( System.Double.IsPositiveInfinity(v) )


1085  {


1086  buf[offs+0] = '.';


1087  buf[offs+1] = 'p';


1088  buf[offs+2] = 'o';


1089  buf[offs+3] = 's';


1090  buf[offs+4] = 'i';


1091  buf[offs+5] = 'n';


1092  buf[offs+6] = 'f';


1093  buf[offs+7] = '_';


1094  buf[offs+8] = '_';


1095  buf[offs+9] = '_';


1096  buf[offs+10] = '_';


1097  offs += SER_ENTRY_LENGTH;


1098  return;


1099  }


1100  if( System.Double.IsNegativeInfinity(v) )


1101  {


1102  buf[offs+0] = '.';


1103  buf[offs+1] = 'n';


1104  buf[offs+2] = 'e';


1105  buf[offs+3] = 'g';


1106  buf[offs+4] = 'i';


1107  buf[offs+5] = 'n';


1108  buf[offs+6] = 'f';


1109  buf[offs+7] = '_';


1110  buf[offs+8] = '_';


1111  buf[offs+9] = '_';


1112  buf[offs+10] = '_';


1113  offs += SER_ENTRY_LENGTH;


1114  return;


1115  }


1116 


1117  //


1118  // process general case:


1119  // 1. copy v to array of chars


1120  // 2. set 9th byte to zero in order to simplify conversion to sixbit representation


1121  // 3. convert to little endian (if needed)


1122  // 4. convert to sixbit representation


1123  // (last 12th element of sixbits is always zero, we do not output it)


1124  //


1125  byte[] _bytes = System.BitConverter.GetBytes((double)v);


1126  if( !System.BitConverter.IsLittleEndian )


1127  System.Array.Reverse(_bytes);


1128  for(i=0; i<sizeof(double); i++)


1129  bytes[i] = _bytes[i];


1130  for(i=sizeof(double); i<9; i++)


1131  bytes[i] = 0;


1132  threebytes2foursixbits(bytes, 0, sixbits, 0);


1133  threebytes2foursixbits(bytes, 3, sixbits, 4);


1134  threebytes2foursixbits(bytes, 6, sixbits, 8);


1135  for(i=0; i<SER_ENTRY_LENGTH; i++)


1136  buf[offs+i] = sixbits2char(sixbits[i]);


1137  offs += SER_ENTRY_LENGTH;


1138  }


1139 


1140  /************************************************************************


1141  This function unserializes double value from string


1142 


1143  buf buffer which contains value; leading spaces/tabs/newlines are


1144  ignored, traling spaces/tabs/newlines are treated as end of


1145  the double value.


1146  offs offset in the buffer


1147 


1148  after return from this function, offs points to the char's past the value


1149  being read.


1150 


1151  This function raises an error in case unexpected symbol is found


1152  ************************************************************************/


1153  private static double str2double(char[] buf, ref int offs)


1154  {


1155  string emsg = "ALGLIB: unable to read double value from stream";


1156  int[] sixbits = new int[12];


1157  byte[] bytes = new byte[9];


1158  byte[] _bytes = new byte[sizeof(double)];


1159  int sixbitsread, i;


1160 


1161 


1162  //


1163  // skip leading spaces


1164  //


1165  while( buf[offs]==' '  buf[offs]=='\t'  buf[offs]=='\n'  buf[offs]=='\r' )


1166  offs++;


1167 


1168 


1169  //


1170  // Handle special cases


1171  //


1172  if( buf[offs]=='.' )


1173  {


1174  string s = new string(buf, offs, SER_ENTRY_LENGTH);


1175  if( s==".nan_______" )


1176  {


1177  offs += SER_ENTRY_LENGTH;


1178  return System.Double.NaN;


1179  }


1180  if( s==".posinf____" )


1181  {


1182  offs += SER_ENTRY_LENGTH;


1183  return System.Double.PositiveInfinity;


1184  }


1185  if( s==".neginf____" )


1186  {


1187  offs += SER_ENTRY_LENGTH;


1188  return System.Double.NegativeInfinity;


1189  }


1190  throw new alglib.alglibexception(emsg);


1191  }


1192 


1193  //


1194  // General case:


1195  // 1. read and decode sixbit digits


1196  // 2. check that all 11 digits were read


1197  // 3. set last 12th digit to zero (needed for simplicity of conversion)


1198  // 4. convert to 8 bytes


1199  // 5. convert to big endian representation, if needed


1200  //


1201  sixbitsread = 0;


1202  while( buf[offs]!=' ' && buf[offs]!='\t' && buf[offs]!='\n' && buf[offs]!='\r' && buf[offs]!=0 )


1203  {


1204  int d;


1205  d = char2sixbits(buf[offs]);


1206  if( d<0  sixbitsread>=SER_ENTRY_LENGTH )


1207  throw new alglib.alglibexception(emsg);


1208  sixbits[sixbitsread] = d;


1209  sixbitsread++;


1210  offs++;


1211  }


1212  if( sixbitsread!=SER_ENTRY_LENGTH )


1213  throw new alglib.alglibexception(emsg);


1214  sixbits[SER_ENTRY_LENGTH] = 0;


1215  foursixbits2threebytes(sixbits, 0, bytes, 0);


1216  foursixbits2threebytes(sixbits, 4, bytes, 3);


1217  foursixbits2threebytes(sixbits, 8, bytes, 6);


1218  for(i=0; i<sizeof(double); i++)


1219  _bytes[i] = bytes[i];


1220  if( !System.BitConverter.IsLittleEndian )


1221  System.Array.Reverse(_bytes);


1222  return System.BitConverter.ToDouble(_bytes,0);


1223  }


1224  }}

