/ * Filename: MYLIB.H Author: Dcyu Date: 2002.9.9 * /
#include
#define infinity int_max # Define true 1 # define false 0
Typedef Double Graph; Typedef Int_Bool; TypeDef Struct Complex COMPLEX;
Void Error (Char * Message); void _pause (void); double _time (void); double _dif (double start, double end); double ** _ alloc2 (int R, int C); int ** _ alloc2_int (int R, int C); void _alloc2free (double ** x); void _dtoa (double value, char * string); void _putstring (int x, int y, char * msg, int color); long _fac (int N); int _fac2 (int N, int * a); Double _random (void); int _rand (int seed); double _avg (double A, double b); double _sta (double mu, double sigma); double; DOUBLE _STA2 (Double Mu, Double Sigma); Double _kai (int N); double _tdis (int N); double _f (int N1, int n2); double _poisson (int K, double lam); double _mean (double * a, INT Start, int end); Double _variance (Double * a, int start, int end); double _linear (double * x, double * y, double * a, double * b, int N); double _floyd (int i, INT J, INT K, DOUBLE ** D); Double _Dijkstra (graph ** g, int N, int * count); double _0618 (double x) , Double Start, Double End, Double EPS); Double_INTEGRAL (Double X) (Double X), Doubl e start, double end, int n; double _integral2 (double x) (double x), double start, double end, int div); long _comb (int N, int m); long _rank (int N, int m); void _allcomb_dg (int * ps, int * PE, int Elems, int * buf, int buffs, int **); void _allcomb (int * list, int n, int elems, int ** Comb); void _allrank_dg (int M, INT K, INT S, INT * A, INT * FLAG, INT * ICOUNT, INT ** RANK); void _allrank (int M, int ** rank); int _bolziman (Double DE , Double T); INT _ISPRIME (Long P); long _prime (int N); void _Swapi (int * a, int * b); void _swapl (long * a, long * b); void _swapd (double * a, Double * b); void _inv (int * x, int N); void _sort (double * x, int N); long _max2 (long n, long m); long _min2 (long n, long m);
Double _max (double * array, int N, int * id); double _min (double * array, int N, int * id); long _gcd (long M, long n); long _lcm (long M, long n); void _initgraph (void); void _SetPlotDefault (void); void _plot (double x) (double x) (double x); double _line_equation (double x, double a, double b); void _fit (double * X, Double * Y, INT N; Complex _Cplus (Complex Z1, Complex Z2); Complex_Cminus (Complex Z1, Complex Z2); Complex_cmul (Complex Z1, Complex Z2); Complex_CDIV (Complex Z1, Complex Z2); Double _CABS (Complex Z); Void _Croot (Complex Z1, INT N, Complex * Z); Complex_Cpow (Complex Z1, Double W); Double_Lag (Double * X, Double * Y, Double T, INT N); Double _Newt (double * x, double * y, int n, double t); int _mgauss (double ** a, double * b, int n, double ep); int _mgauss2 (double ** a, double ** b, int n, int m, double ep); int _mdjn (double ** a, double ** C, int N, int m); double _nor (double x, int L); double _mdet (double ** a, int N) ; int _MINV (Double T0, Double * T, Double * TT, INT N, INT M, DOUB Le ** b); void error (char * message) {clrs CR (); fprintf (stderr, "error:% s / n", message); exit (1);
}
Void_pause (void) {while (Bioskey (1) == 0);
}
Void_delay (double delaytime) {long bios_time; double start, end;
BIOS_TIME = Biostime (0, 0L); start = bios_time / clk_tck;
While (1) {BIOS_TIME = Biostime (0, 0L); End = BIOS_TIME / CLK_TCK; if (end-start> = delaytime) return;
}
}
Double _time (void) {long bios_time; double cur;
BIOS_TIME = Biostime (0, 0L); CUR = BIOS_TIME / CLK_TCK; RETURN CUR;
}
Double_dif (Double Start, Double End) {Return End-Start;
}
Double ** _ alloc2 (int R, int C) {double * x, ** y; int N;
X = (double *) Calloc (r * c, sizeof (double)); y = (double **) Calloc (r, sizeof (double *)); for (n = 0; n <= r-1; n ) Y [n] = & x [c * n]; return y;
}
INT ** _ Alloc2_int (int R, INT C) {INT * x, ** y; int N;
x = (int *) Calloc (R * C, SIZEOF (INT)); y = (int **) Calloc (R, SIZEOF (INT *)); for (n = 0; n <= r-1; n ) Y [n] = & x [c * n];
Return Y;
}
Void_alloc2free (double ** x) {free (x [0]); free (x);
}
Void_dtoa (double value, char * string) {unsigned long wn, df; char * str1, * str2, * str3, * STR4;
STR1 = value <0.0? "-": "; value = value <0.0? (- value): value; wn = (unsigned long); DF = (unsigned long) (Florue" -wn) * 1000)); str2 = (char *) malloc (10 * sizeof (char)); str4 = (char *) malloc (10 * sizeof (char)); UltoA (Wn, Str2, 10); STR3 = "."; Ultoa (DF, STR4, 10); strcpy (string, str1); strcat (string, str2); strcat (string, str3); strcat (string, str4); STRCAT
}
Void _putString (int x, int y, char * msg, int color) {int Len; const Int sIDEX = 8, Sidey = 10;
SetColor; LEN = Strlen (MSG); x = x-len * sIDEX; Y = Y-Sidey; OutTextxy (x, y, msg);
}
Long _fac (int N) {if (n <0 || n> 12) error ("n is not valid in _fac!"); if (n == 0 || n == 1) Return 1; return_fac ( N-1) * n;
}
INT _FAC2 (int N, int * a) {INT M, I, J, C, T;
A [0] = 1; m = 1; for (i = 2; i <= n; i ) {for (c = 0, j = 0; j For (j = 0; j <= (m-1) / 2; j ) {t = a [j]; a [j] = a [m-1-j]; a [m-1-j] = t;} return m; } Double _random (void) {Int a; double r; A = rand ()% 32767; R = (A 0.00) / 32767.00; Return R; } INT _RAND (int seed) {return Rand ()% SEED; } Double _avg (double a, double b) {double _random (void); Return A _Random () * (b-a); } Double _sta (double mu, double sigma) {double _random (void); INT I; Double R, Sum = 0.0; IF (Sigma <= 0.0) Error ("SIGMA <= 0.0 in _sta!"); for (i = 1; i <= 12; i ) SUM = SUM _RANDOM (); r = (SUM-6.00) * Sigma MU; Return R; } Double _sta2 (Double Mu, Double Sigma) {Double R1, R2; R1 = _random (); r2 = _random (); RETURN SQRT (-2 * log (r1)) * COS (2 * m_pi * r2) * Sigma MU; } Double _kai (int N) {double _sta (double mu, double sigma); INT I; Double SUM = 0.0; for (i = 1; i <= n; i ) SUM = SUM POW (_sta (0, 1), 2); Return SUM; } Double _TDIS (int N); Double _STA (Double MU, Double Sigma); Double X, Y; X = _STA (0, 1); y = _kai (n); Return X / SQRT (Y / N); } Double _f (int N1, int N2) {double _kai (int N); Double U1, U2; U1 = _KAI (N1); u2 = _kai (n2); Return (u1 * n2) / (u2 * n1); } Double _poisson (int K, Double Lam) {IF (k> 12 || K <0) Error ("k is not valid in _poisson!"); Return Pow (LAM, K) * Exp (-LAM) / _ FAC (K); } Double _me (double * a, int start, int end) {INT i; double sum = 0.0; IF (start <0 || end-start <0) error ("START IS LARGER THAN END, OR Start Is Not Valid In _Mean!"); For (i = start; i <= end; i ) SUM = SUM A [I]; Return Sum / (end-start 1); } Double _variance (Double * a, int start, int end) {double _mean (double * a, int start, int end); INT I; Double Mean, Sum = 0.0; Mean = _me (a, start, end); for (i = start; i <= end; i ) SUM = SUM (a [i] -mean) * (a [i] -mean); Return Sum / (end-start); Double _Linear (double * x, double * y, double * a, double * b, int N) {Double SXX, SYY, SXY, Meanx, Meany, QE; INT I; IF (N <= 2) Error ("N TOO SMALL IN _LINEAR!"); meanx = _mean (x, 0, n-1); meany = _mean (y, 0, n-1); sxx = _variance (x , 0, N-1); Syy = _Variance (Y, 0, N-1) * (N-1); SXY = 0.0; for (i = 0; i } Double ** _ adj_form (Double (* a) [3], int N, int m) {double ** g; int i, j; G = _alloc2 (n, n); for (i = 0; i For (i = 0; i } Double _Floyd (INT I, INT J, INT K, Double ** D) {Double Temp_ij, Temp_ik, Temp_kj; if (k> 0) {TEMP_IJ = _FLOYD (i, J, K-1, D); Temp_ik = _FLOYD (I, K, K-1, D); TEMP_KJ = _FLOYD (K, J, K-1, D); RETURN TEMP_IJ <(TEMP_IK TEMP_KJ)? Temp_ij: (Temp_ik Temp_kj);} IF (k == 0) RETURN * (* (D I) J); Else {Error ("K <0 in _floyd!"); Return 1;} } Double _dijkstra (graph ** g, int N, int * count) {INT i, J, * Mark, * pathd, * pathbk; double w, minc, * d; int from , TO; D = (double *) malloc (n * sizeof (double)); mark = (int *) malloc (n * sizeof (int)); pathd = (int *) malloc (n * sizeof (int)); Pathbk = (int *) malloc (n * sizeof (int)); for (i = 0; i } } J = 0; for (i = N-1; I> = 0; I -) IF (Pathbk [i]! = - 1) {path [j] = pathbk [i]; j ;} path [j] = T; * count = j 1; Return D [T]; } Double _0618 (double x) (double x) (double x), double, double eps) {DOUBLE A, B, C, D; IF (EPS <0.0 || end-start ELSE IF ((* f) (b)> (* f) && (* f) (d)> (* f)) {b = D; D = C; c = a bd; } Else {= C; B = D; C = A (B-A) * 0.381966011; D = A (B-A) * 0.618033989;}}}}}}}} RETURN (A B) / 2; } INT _2DIV (Double X), Double A, Double B, Double H, Double EPS, DOUBLE * X, INT N, INT * M) {Double Z, Z0, Z1, Y, Y0, Y1; * m = 0; z = a; y = (* f) (z); while (1) {IF ((z> b h / 2) || (* m == n)) Return 1; IF ( FABS (Y) } Double _INTEGRAL (double x) (double x), double start, double end, int N) {Double H, Al, BE, SUM; INT I, M Download Adobe Reader IF (n <= 0 || End-start <= 0.0) error ("START IS LARGER THAN END, OR N <= 0 in _INTEGRAL!"); h = (end-start) / n; al = start h * 0.4226497308; BE = start h * 1.577350269; SUM = (* f) (al) (* f) (be); m = n / 2-1; for (i = 1; i <= m; i ) {Al = Al 2 * h; be = be 2 * h; SUM = SUM (* F) (Al) (* f) (BE); } Sum = sum * sum; returnium; } Double _INTEGRAL2 (Double X) (Double X) (Double X) (Double START, DOUBLE END, INT DIV) {Double S, H, Y; INT I; IF (DIV <= 0 || End-start <= 0.0) Error ("START IS LARGER THAN END, OR DIV <= 0 in _INTEGRAL2!"); s = ((* f) (start) (* f) (END) / 2.00; H = (end-start) / div; for (i = 1; i Long _comb (int N, int m) {long _fac (int N); INT I; long TEMP = 1; IF (n> 15 && m> 8 || n <= 0 || m <0 || n For (i = n; i> = N-m 1; I -) TEMP = TEMP * I; TEMP = TEMP / _FAC (M); Return Temp; } Long _rank (int N, int m) {INT i; long TEMP = 1; IF (N> 12 || n <= 0 || m <0 || n For (i = n; i> = N-m 1; I -) TEMP = TEMP * I; RETURN TEMP; } Void _allcomb_dg (int * ps, int * pe, int elems, int buf [], int buffs, int **, int * iCount) {INT * PI, I; IF (Elems == 1) {for (PI = PS; pi } Void_Allcomb (int * list, int n, int elems, int **) {Int * buf; int i, icount; ICOUNT = 0; buf = (int *) Malloc (Elems * sizeof (int)); _allcomb_dg (& List [0], & List [N], Elems, BUF, ELEMS, COMB, & ICOUNT); for (i = 0; i } Void_allrank_dg (int * a, int * flag, int * iCount, int ** rank) {INT i; if (s> = k) {for (i = 0; i * iCount = * iCount 1; else {for (i = 1; i <= m; i ) IF (Flag [i] == 0) {FLAG [i] = 1; a [s] = i; _allrank_dg (M, K, S 1, A, Flag, ICOUNT, RANK); A [S] = 0; FLAG [i] = 0;}} } Void_allrank (int M, int ** rank) {INT * flag, * a; int i, icount, n, k; n = m; k = 0; flag = (int *) malloc ((m 1) * sizeof (int)); a = (int *) malloc ((m 1) * sizeof (int)); for i = 1; i <= m; i ) FLAG [i] = 0; ICOUNT = 0; _allrank_dg (M, N, K, A, Flag, & ICOUNT, RANK); } INT _BOLZIMAN (Double de, Double T) {Double _random (void); Return de <0.0 || _random () } INT _ISPRIME (long p) {long i; For (i = 2; i <= SQRT (P); i ) IF (p% i == 0) Return 0; Return 1; } Long _prime (int N) {INT _ISPRIME (long P); Long i = 2; int Num = 1; While (NUM <= n) {if (_isprime (i )) Num ;} Return I-1; } Void _SWAPI (int * a, int * b) {int P; p = * a; * a = * b; * b = P; } Void _Swapl (long * a, long * b) {long p; p = * a; * a = * b; * b = P; } Void _Swapd (double * a, double * b) {double p; p = * a; * a = * b; * b = P; } Void _inv (int * x, int N) {INT * P, T, * I, * J, M M = (n-1) / 2; i = x; j = x N-1; p = x m; for (; i <= p; i , j -) {t = * i; * i = * j; * j = t;} } Void _sort (double * x, int N) {INT I, J, K; Double T; For (i = 0; i } Long _max2 (long n, long m) {return n> m? N: m; } Long _min2 (long n, long m) {RETURN N } Double _max (double * array, int n, int * id) {double * p, * array_end, dmax; Array_end = array n; dmax = * array; * id = 0; for (p = array 1; p } Double _min (Double * Array, Int n, int * id) {double * p, * array_end, dmin; Array_end = array n; Dmin = * array; * id = 0; for (p = array 1; p Return DMIN; } Long _GCD (long M, long n) {long _max2 (long n, long m); long _min2 (long n, long m); Long T; T = _MAX2 (N, M); m = _min2 (n, m); n = t; while (nm! = 0) {n = nm; t = _max2 (n, m); m = _min2 (n, m ); N = T; } Return t; } Long _lcm (long M, long n) {long _GCD (long n, long m); RETURN M * N / _GCD (M, N); } Void _initgraph (void) {INT GDRIVER = VGA, GMODE = VGAHI, ERRORCODE INITGRAPH (& gDriver, & gmode, ".//"); errorcode = graphresult (); if (ErrorCode! = GROK) {Printf ("/ NGRAPHICS ERROR:% S / N", GrapherrorMsg (ErrorCode)); Printf ("Press Any key to halt! "); getCH (); exit (1); } Void _SetPlotDefault (void) {const Int basex = 60, basey = 420; const Int basebkcolor = lightblue, basecolor = red; Setbkcolor (BaseBkcolor); SetColor (Basecolor); Line (0, Basey, Getmaxx (), Basey; Line (basex, 0, basex, getmaxy ()); SetFillStyle (Solid_Fill, Basecolor); Fillellipse (Basex, Basey, 3, 3); } Void _plot (double x) (double x) (double x) {const Int basex = 60, basey = 420, grid = 20; const double eps = 1e-5; double step, stepy, * point Double DMAX, DMIN, DI; INT I, X, Y, LX, LY, MAXID, MINID; IF (End-Start _SETPLOTDEFAULT (); Stepx = (end-start) / (getmaxx () - basex 0.0); Point = (Double *) Malloc (GetMaxx () - Basex) * Sizeof (Double *)); for (i = basex i <= getmaxx (); i ) Point [i-basex] = (* f) (Start (i-basex) * stepx); dmax = _max (Point, getmaxx () - basex, & maxid; dmin = _min (Point, getmaxx () - basex, & minid; if (DMAX-DMIN } SetColor (Blue); outtextxy (getMaxx () / 2 50, getmaxy () - 30, "press any key to start analyysis!"); If (getch () == 27) Return; Closegraph (); CLRSCR Printf ("curve analysis: / n"); Printf ("FROM% LF TO% LF / N", START, END); Printf ("The maximum of the curve IS% LF, WHEN X =% LF. / N ", DMAX, START (MAXID 0.0) / (639-Basex 0.0) * (end-start); Printf (" "The minimum of the curve is% lf, when x =% lf. / n", Dmin , start (minid 0.0) / (end-start); Stepx = (end-start) / (Grid 0.0); for (di = start; di <= end; di = Di Stepx) Printf ("F (% LF) =% LF / N", DI, (* f) (DI)); _Pause (); } Double _Line_Equation (double x, double a, double b) {return A b * x; } Void _fit (double * x, double * y, int N) {double _Linear (double * x, double * y, double * a, double * b, int N); const int basex = 60, basey = 420, grid = 20; const double eps = 1e-5; Double A, B; Double DMAXX, DMINX, DMAXY, DMINY, * PDX, * PDY; INT I, Maxidx, Minidx, Maxidy, MINIDY, * PX, * Py; Double Stepx, STEPY, * POINT; DOUBLE DMAX, DMIN, DI, START, END; INT CX, CY, LX, LY, MAXID, MINID; CHAR * STR _SETPLOTDEFAULT (); _LINEAR (X, Y, & A, & B, N); Dmaxx = _max (x, n, & maxidx); dminx = _min (x, n, & minidx); dmaxy = _max (y, n, & maxidy); dminy = _MIN (Y, n, & minidy); if (DMAXX-DMINX < EPS || DMAXY-DMINY START = DMINX; END = DMAXX; Stepx = (end-start) / (getmaxx () - Basex 0.0); Point = (Double *) Malloc (GetMaxX () - Basex) * sizeof (Double *); for (i = basex; i <= getmaxx (); i ) Point [i-basex] = a b * (start (i-basex) * stepx); dmax = _max (Point, getmaxx () - basex, & maxid) Dmin = _MIN (Point, getmaxx () - basex, & minid); if (DMAX-DMIN } Str = (char *) Malloc (10 * sizeof (char)); for (i = 0; i } } Complex _Cplus (Complex Z1, Complex Z2) {complex z; z.x = z1.x z2.x; z.y = z1.y z2.y; return z; } Complex _cminus (Complex Z1, Complex Z2) {complex z; z.x = z1.x-z2.x; z.y = z1.y-z2.y; returnz; } Complex _cmul (complex z1, complex z2) {complex z; zx = z1.x * z2.x-z1.y * z2.y; zy = z1.x * z2.y z1.y * z2.x; return z; } Complex _CDIV (Complex Z1, Complex Z2) {Double E, F; Complex Z; IF (FABS (Z2.x)> = Fabs (Z2.Y)) {E = Z2.Y / Z2.x; f = z2.x e * z2.y; zx = (z1.x z1.y * e) / f; zy = (z1.y-z1.x * e) / f;} else {e = z2.x / z2.y; f = z2.y E * z2.x; zx = Z1.x * E Z1.Y) / f; Zy = (z1.y * e-z1.x) / f;} returnz; } Double _CABS (Complex Z) {Return Hypot (z.x, z.y); } Void _croot (Complex Z1, INT N, Complex * Z) {Int J; Double R, C, C1, CK, CS, S, S1, SK, SC; IF (z1.y == 0.0) {r = EXP (log (FABS (Z1.x)) / N); CS = 0;} else if (z1.x == 0.0) {if (z1.y> 0 CS = M_PI_2; ELSE CS = -M_PI_2; R = Exp (log (FABS (Z1.Y)) / n);} else {r = exp (log (z1.x * z1.x z1.y * z1 .y) / n / 2); cs = atan (z1.y / z1.x);} if (z1.x <0) cs = m_pi; cs / = n; sc = 2 * m_pi / n; c = COS (CS); S = SiN (CS); C1 = COS (SC); S1 = SiN (SC); CK = 1; SK = 0; z [0] .x = c * r; z [0]. Y = S * r; for (j = 1; j Complex _cpow (complex z1, double w) {double r, t; complex z; IF (z1.x == 0 && z1.y == 0) RETURN Z1; if (z1.x == 0) {if (z1.y> 0) T = m_pi_2; ELSE T = -M_PI_2;} else {IF Z1.X> 0) T = Atan (z1.y / z1.x); Else {if (z1.y> = 0) T = atan (z1.y / z1.x) m_pi; else t = atan Z1.Y / z1.x) -m_pi;}} r = exp (w * log (sqrt (z1.x * z1.x z1.y * z1.y))); zx = r * cos (w * T); zy = r * sin (w * t); returnz; } Double _lag (double * x, double * y, double t, int N) {INT I, J; Double P, S = 0.0; For (i = 0; i } Return s; } Double _newt (double * x, double * y, int N, double t) {INT I, J; Double * S, P; S = (double *) Calloc (n 1, sizeof (double)); if (s == null) Error ("Calloc Error IN _NEWT"); for (i = 1; i <= n; i ) s I] = y [i]; for (j = 1; j <= n-1; j ) for (i = n; i> = j 1; I -) s [i] = (s [i] -S [i-1]) / (x [i] -X [ij]); p = y [n]; for (i = N-1; I> = 1; I -) p = p * ( TX [I]) S [i]; free (s); returnit p; } INT _MGAUSS (Double ** a, double * B, int N, double ep) {INT I, J, K, L; Double C, T; for (k = 1; k <= n; k ) {c = 0.0 ; For (i = k; i <= n; i ) IF (Fabs (a [i-1] [k-1])> Fabs (c)) {c = a [i-1] [k-1] ; L = i;} IF (FABS (C) <= EP) Return 0; if (l! = K) {for (j = k; j <= n; j ) {t = a [k-1] [ J-1]; a [k-1] [j-1] = a [l-1] [j-1]; A [L-1] [J-1] = T;} T = B [k- 1]; B [k-1] = b [l-1]; b [l-1] = T;} C = 1 / c; for (j = k 1; j <= n; j ) {a [K-1] = a [k-1] [j-1] * C; for (i = k 1; i <= n; i ) a [i-1] [J-1 ] - = a [i-1] [k-1] * a [k-1] [j-1];} B [K-1] * = C; for (i = k 1; i <= n ; i ) B [i-1] - = b [k-1] * a [i-1] [k-1];} for (i = n; i> = 1; I -) for (j = i 1; j <= n; j ) B [i-1] - = a [i-1] [j-1] * b [j-1]; Return 1;} INT _MGAUSS2 (Double ** a, double ** B, int N, int m, double ep) {INT I, J, K, L; Double C, T; For (k = 1; k <= n; k ) {c = 0; for (j = k; j <= n; j ) IF (Fabs (a [k-1] [j-1])> FABS ( c) {c = a [k-1] [j-1]; l = j; } IF (FABS (C) <= EP) RETURN 0; if (l! = K) {for (i = k; i <= n; i ) {t = a [i-1] [k-1]; A [i-1] [k-1] = a [i-1] [l-1]; A [i-1] [L-1] = T; } For (i = 1; i <= m; i ) {t = b [i-1] [k-1]; b [i-1] [K-1] = B [i-1] [L- 1]; B [i-1] [L-1] = T; } C = 1 / c; for (i = k 1; i <= n; i ) {a [i-1] [k-1] * = C; for (j = k 1; j <= n ; J ) a [i-1] [j-1] - = a [k-1] [j-1] * a [i-1] [k-1]; } For (i = 1; i <= m; i ) {b [i-1] [k-1] * = C; for (j = k 1; j <= n; j ) B [i-1 ] [J-1] - = a [k-1] [j-1] * b [i-1] [k-1]; } } For (i = n; i> = 1; I -) {for (j = i 1; j <= n; j ) for (k = 1; k <= m; k ) B [K-1 ] [i-1] - = a [j-1] [i-1] * b [k-1] [j-1]; } Return 1; } INT _MDJN (double ** a, double ** C, int N, int m) {INT I, J, K, K1, K2, K3; IF (Fabs (A [0] [0]) == 0) RETURN 0; for (i = 2; i <= n; i ) a [i-1] [0] / = a [0] [0] ; For (i = 2; i <= n-1; i ) {for (j = 2; j <= i; j ) a [i-1] [i-1] - = a [i-1] [J -2] * a [i-1] [j-2] * a [j-2] [j-2]; for (k = i 1; k <= n; k ) {for (j = 2; J <= i; j ) a [k-1] [i-1] - = a [k-1] [j-2] * a [i-1] [j-2] * a [j-2] [J-2]; IF (FABS (a [i-1] [i-1]) == 0) RETURN 0; A [K-1] [i-1] / = a [i-1] [i -1]; } } For (j = 2; j <= n; j ) a [n-1] [n-1] - = a [n-1] [j-2] * a [n-1] [j-2] * a [j-2] [j-2]; for (j = 1; j <= m; j ) for (i = 2; i <= n; i ) for (k = 2; k <= i; K ) C [i-1] [j-1] - = a [i-1] [K-2] * C [K-2] [J-1]; For (i = 2; i <= n; i ) for (j = i; j <= n; j ) a [i-2] [j-1] = a [i-2] [i-2] * a [j-1] [i-2]; IF (Fabs (a [N-1] [N-1]) == 0) RETURN 0; for (j = 1; j <= m; j ) {C [N-1] [J-1] / = A [N-1] [N-1]; for (k = 2; k <= n; k ) {k1 = N-K 2; for (K2 = K1; K2 <= N; K2 ) {K3 = N-K 1; C [K3-1] [J-1] - = a [K3-1] [K2-1] * C [K2-1] [J-1]; } C [K3-1] [J-1] / = a [k3-1] [k3-1]; } Return 1; } Double _NOR (Double X, INT L) {Double RM, RN, RP, RT, PQ, RRT, X1, X2, Y, S, H; Double P1, P2, Q1, Q2; IF (x == 0) RETURN 0.5; RP = 1; if (x * l <0) RP = -1; x1 = FABS (X); X2 = x1 * x1; y = 1 / sqrt (2 * m_pi) * Exp (-0.5 * x2); RN = Y / X1; IF (RP> = 0 && Fabs (RN) <1e-12) Return 1.0; if (rp <= 0 && rn == 0) Return 0.0; RRT = 3.5; if (rp == - 1) RRT = 2.32; if (x1 <= RRT) {x1 * = y; s = x1; rn = 1.0; rt = 0; DO {RN = 2; rt = s ; X1 * = x2 / rn; s = x1;} while (s! = RT); PQ = 0.5 S; if (rp == - 1) PQ = 0.5-S; return pq;} Q1 = x1; p2 = y * x1; rn = 1; p1 = y; Q2 = x2 1; if (rp! = - 1) {RM = 1-p1 / q1; s = rm; rt = 1-p2 /} else {RM = P1 / Q1; S = RM; RT = P2 / Q2;} while (1) {RN ; if (RM == RT || S == RT) RTURN RT; S = x * p2 rn * p1 P1 = p2; p2 = S; s = x * Q2 RN * Q1; Q1 = Q2; Q2 = S; S = Rm; Rm = r = r = 1-p2 / q2; IF (rp == - 1 Rt = p2 / q2;} RETURN RT; } Double _md (double ** a, int N) {INT I, J, K, I0, J0, SGN; Double Q, MP, PV, PR, T, DET; SGN = 1; Pr = 1.0; for (k = 1; k <= N-1; K ) {mp = 0.0; for (i = k; i <= n; i ) for (j = k; j <= n; j ) {PV = a [i-1] [j-1]; if (FABS (PV)> = Fabs (MP)) {mp = PV; I0 = i; j0 = j;}} IF (MP == 0.0) Error ("failed in _mdet!"); If (i0! = K) {SGN = -sgn; for (j = k; j <= n; j ) {t = a [i0-1] [ J-1]; a [i0-1] [j-1] = a [k-1] [j-1]; A [k-1] [j-1] = T;}} f (j0! = K) {SGN = -SGN; for (i = k; i <= n; i ) {t = a [i-1] [j0-1]; A [i-1] [j0-1] = a [ I-1] [K-1]; A [I-1] [K-1] = T;}}} PR = Pr * MP; MP = -1 / Mp; for (i = k 1; i <= n; i ) {q = a [i-1] [K-1] * mp; for (j = k 1; j <= n; j ) a [i-1] [j-1] = a [ I-1] [J-1] q * a [k-1] [j-1];}}}} DET = SGN * Pr * a [N-1] [N-1]; Return Det; } INT_MINV (Double T0, Double * T, Double * TT, INT N, INT M, DOUBLE ** B) {INT I, J, K; Double A, * C, * R, * P, S; C = Double *) Malloc (M * SizeOf (Double)); r = (double *) Malloc (M * Sizeof (double)); p = (double *) malloc (M * sizeof (double)); if (c == NULL || R == NULL || R == NULL) Error ("Calloc Error IN _MINV!"); If (FABS (T0) == 0) {Free (c); Free (r); Free (p) RETURN 0;} a = T0; C [0] = TT [0] / T0; R [0] = T [0] / t0; for (k = 0; k <= n-2; k ) {s = 0; for (j = 1; j <= k 1; j ) s = s c [k 1-j] * TT [J-1]; s = (S-TT [K 1]) / a; for (i = 1; i <= k 1; i ) p [i-1] = C [i-1] S * R [K 1-I]; C [k 1] = -s; s = 0; for (j = 1; j <= k 1; j ) s = s r [k 1-j] * t [j-1]; s = (ST [k 1 ]) / a; for (i = 1; i <= k 1; i ) {r [i-1] = r [i-1] S * C [K 1-I]; C [K 1-i] = p [k 1-i];} R [k 1] = - S; a = 0; for (j = 1; j <= k 2; j ) a = a t [ J-1] * C [J-1]; A = T0-a; IF (Fabs (a) == 0) Return 0; } B [0] [0] = 1 / a; for (i = 1; i <= n; i ) {b [0] [i] = - r [i-1] / a; b [i] [ 0] = - C [i-1] / a;} for (i = 1; i <= n; i ) for (j = 1; j <= n; j ) {b [i] [j] = b [i-1] -C [i-1] * b [0] [j]; b [i] [j] = b [i] [j] c [nj] * b [0 ] [n 1-j];} free (c); free (r); free (p); returnif 1; } Void _fft (double * fr, double * fi, int n, int flag) {Int MP, Arg, Q, CNTR, P1, P2; INT I, J, A, B, K, M; Double Sign, PR, Pi, Harm, X, Y, T; Double * Ca, * SA;
