Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/spline3.cs @ 3021

Visit:

Last change on this file since 3021 was 2806, checked in by gkronber, 15 years ago
Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 31.2 KB

Rev	Line
[2806]	1	/*************************************************************************
	2	Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
	3
	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
	9
	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
	14
	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
	17
	18	>>> END OF LICENSE >>>
	19	*************************************************************************/
	20
	21	using System;
	22
	23	namespace alglib
	24	{
	25	public class spline3
	26	{
	27	public static void buildlinearspline(double[] x,
	28	double[] y,
	29	int n,
	30	ref double[] c)
	31	{
	32	int i = 0;
	33	int tblsize = 0;
	34
	35	x = (double[])x.Clone();
	36	y = (double[])y.Clone();
	37
	38	System.Diagnostics.Debug.Assert(n>=2, "BuildLinearSpline: N<2!");
	39
	40	//
	41	// Sort points
	42	//
	43	heapsortpoints(ref x, ref y, n);
	44
	45	//
	46	// Fill C:
	47	// C[0] - length(C)
	48	// C[1] - type(C):
	49	// 3 - general cubic spline
	50	// C[2] - N
	51	// C[3]...C[3+N-1] - x[i], i = 0...N-1
	52	// C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
	53	//
	54	tblsize = 3+n+(n-1)*4;
	55	c = new double[tblsize-1+1];
	56	c[0] = tblsize;
	57	c[1] = 3;
	58	c[2] = n;
	59	for(i=0; i<=n-1; i++)
	60	{
	61	c[3+i] = x[i];
	62	}
	63	for(i=0; i<=n-2; i++)
	64	{
	65	c[3+n+4*i+0] = y[i];
	66	c[3+n+4*i+1] = (y[i+1]-y[i])/(x[i+1]-x[i]);
	67	c[3+n+4*i+2] = 0;
	68	c[3+n+4*i+3] = 0;
	69	}
	70	}
	71
	72
	73	public static void buildcubicspline(double[] x,
	74	double[] y,
	75	int n,
	76	int boundltype,
	77	double boundl,
	78	int boundrtype,
	79	double boundr,
	80	ref double[] c)
	81	{
	82	double[] a1 = new double[0];
	83	double[] a2 = new double[0];
	84	double[] a3 = new double[0];
	85	double[] b = new double[0];
	86	double[] d = new double[0];
	87	int i = 0;
	88	int tblsize = 0;
	89	double delta = 0;
	90	double delta2 = 0;
	91	double delta3 = 0;
	92
	93	x = (double[])x.Clone();
	94	y = (double[])y.Clone();
	95
	96	System.Diagnostics.Debug.Assert(n>=2, "BuildCubicSpline: N<2!");
	97	System.Diagnostics.Debug.Assert(boundltype==0 \| boundltype==1 \| boundltype==2, "BuildCubicSpline: incorrect BoundLType!");
	98	System.Diagnostics.Debug.Assert(boundrtype==0 \| boundrtype==1 \| boundrtype==2, "BuildCubicSpline: incorrect BoundRType!");
	99	a1 = new double[n-1+1];
	100	a2 = new double[n-1+1];
	101	a3 = new double[n-1+1];
	102	b = new double[n-1+1];
	103
	104	//
	105	// Special case:
	106	// * N=2
	107	// * parabolic terminated boundary condition on both ends
	108	//
	109	if( n==2 & boundltype==0 & boundrtype==0 )
	110	{
	111
	112	//
	113	// Change task type
	114	//
	115	boundltype = 2;
	116	boundl = 0;
	117	boundrtype = 2;
	118	boundr = 0;
	119	}
	120
	121	//
	122	//
	123	// Sort points
	124	//
	125	heapsortpoints(ref x, ref y, n);
	126
	127	//
	128	// Left boundary conditions
	129	//
	130	if( boundltype==0 )
	131	{
	132	a1[0] = 0;
	133	a2[0] = 1;
	134	a3[0] = 1;
	135	b[0] = 2*(y[1]-y[0])/(x[1]-x[0]);
	136	}
	137	if( boundltype==1 )
	138	{
	139	a1[0] = 0;
	140	a2[0] = 1;
	141	a3[0] = 0;
	142	b[0] = boundl;
	143	}
	144	if( boundltype==2 )
	145	{
	146	a1[0] = 0;
	147	a2[0] = 2;
	148	a3[0] = 1;
	149	b[0] = 3(y[1]-y[0])/(x[1]-x[0])-0.5boundl*(x[1]-x[0]);
	150	}
	151
	152	//
	153	// Central conditions
	154	//
	155	for(i=1; i<=n-2; i++)
	156	{
	157	a1[i] = x[i+1]-x[i];
	158	a2[i] = 2*(x[i+1]-x[i-1]);
	159	a3[i] = x[i]-x[i-1];
	160	b[i] = 3(y[i]-y[i-1])/(x[i]-x[i-1])(x[i+1]-x[i])+3(y[i+1]-y[i])/(x[i+1]-x[i])(x[i]-x[i-1]);
	161	}
	162
	163	//
	164	// Right boundary conditions
	165	//
	166	if( boundrtype==0 )
	167	{
	168	a1[n-1] = 1;
	169	a2[n-1] = 1;
	170	a3[n-1] = 0;
	171	b[n-1] = 2*(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
	172	}
	173	if( boundrtype==1 )
	174	{
	175	a1[n-1] = 0;
	176	a2[n-1] = 1;
	177	a3[n-1] = 0;
	178	b[n-1] = boundr;
	179	}
	180	if( boundrtype==2 )
	181	{
	182	a1[n-1] = 1;
	183	a2[n-1] = 2;
	184	a3[n-1] = 0;
	185	b[n-1] = 3(y[n-1]-y[n-2])/(x[n-1]-x[n-2])+0.5boundr*(x[n-1]-x[n-2]);
	186	}
	187
	188	//
	189	// Solve
	190	//
	191	solvetridiagonal(a1, a2, a3, b, n, ref d);
	192
	193	//
	194	// Now problem is reduced to the cubic Hermite spline
	195	//
	196	buildhermitespline(x, y, d, n, ref c);
	197	}
	198
	199
	200	public static void buildhermitespline(double[] x,
	201	double[] y,
	202	double[] d,
	203	int n,
	204	ref double[] c)
	205	{
	206	int i = 0;
	207	int tblsize = 0;
	208	double delta = 0;
	209	double delta2 = 0;
	210	double delta3 = 0;
	211
	212	x = (double[])x.Clone();
	213	y = (double[])y.Clone();
	214	d = (double[])d.Clone();
	215
	216	System.Diagnostics.Debug.Assert(n>=2, "BuildHermiteSpline: N<2!");
	217
	218	//
	219	// Sort points
	220	//
	221	heapsortdpoints(ref x, ref y, ref d, n);
	222
	223	//
	224	// Fill C:
	225	// C[0] - length(C)
	226	// C[1] - type(C):
	227	// 3 - general cubic spline
	228	// C[2] - N
	229	// C[3]...C[3+N-1] - x[i], i = 0...N-1
	230	// C[3+N]...C[3+N+(N-1)*4-1] - coefficients table
	231	//
	232	tblsize = 3+n+(n-1)*4;
	233	c = new double[tblsize-1+1];
	234	c[0] = tblsize;
	235	c[1] = 3;
	236	c[2] = n;
	237	for(i=0; i<=n-1; i++)
	238	{
	239	c[3+i] = x[i];
	240	}
	241	for(i=0; i<=n-2; i++)
	242	{
	243	delta = x[i+1]-x[i];
	244	delta2 = AP.Math.Sqr(delta);
	245	delta3 = delta*delta2;
	246	c[3+n+4*i+0] = y[i];
	247	c[3+n+4*i+1] = d[i];
	248	c[3+n+4i+2] = (3(y[i+1]-y[i])-2d[i]delta-d[i+1]*delta)/delta2;
	249	c[3+n+4i+3] = (2(y[i]-y[i+1])+d[i]delta+d[i+1]delta)/delta3;
	250	}
	251	}
	252
	253
	254	public static void buildakimaspline(double[] x,
	255	double[] y,
	256	int n,
	257	ref double[] c)
	258	{
	259	int i = 0;
	260	double[] d = new double[0];
	261	double[] w = new double[0];
	262	double[] diff = new double[0];
	263
	264	x = (double[])x.Clone();
	265	y = (double[])y.Clone();
	266
	267	System.Diagnostics.Debug.Assert(n>=5, "BuildAkimaSpline: N<5!");
	268
	269	//
	270	// Sort points
	271	//
	272	heapsortpoints(ref x, ref y, n);
	273
	274	//
	275	// Prepare W (weights), Diff (divided differences)
	276	//
	277	w = new double[n-2+1];
	278	diff = new double[n-2+1];
	279	for(i=0; i<=n-2; i++)
	280	{
	281	diff[i] = (y[i+1]-y[i])/(x[i+1]-x[i]);
	282	}
	283	for(i=1; i<=n-2; i++)
	284	{
	285	w[i] = Math.Abs(diff[i]-diff[i-1]);
	286	}
	287
	288	//
	289	// Prepare Hermite interpolation scheme
	290	//
	291	d = new double[n-1+1];
	292	for(i=2; i<=n-3; i++)
	293	{
	294	if( (double)(Math.Abs(w[i-1])+Math.Abs(w[i+1]))!=(double)(0) )
	295	{
	296	d[i] = (w[i+1]diff[i-1]+w[i-1]diff[i])/(w[i+1]+w[i-1]);
	297	}
	298	else
	299	{
	300	d[i] = ((x[i+1]-x[i])diff[i-1]+(x[i]-x[i-1])diff[i])/(x[i+1]-x[i-1]);
	301	}
	302	}
	303	d[0] = diffthreepoint(x[0], x[0], y[0], x[1], y[1], x[2], y[2]);
	304	d[1] = diffthreepoint(x[1], x[0], y[0], x[1], y[1], x[2], y[2]);
	305	d[n-2] = diffthreepoint(x[n-2], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
	306	d[n-1] = diffthreepoint(x[n-1], x[n-3], y[n-3], x[n-2], y[n-2], x[n-1], y[n-1]);
	307
	308	//
	309	// Build Akima spline using Hermite interpolation scheme
	310	//
	311	buildhermitespline(x, y, d, n, ref c);
	312	}
	313
	314
	315	public static double splineinterpolation(ref double[] c,
	316	double x)
	317	{
	318	double result = 0;
	319	int n = 0;
	320	int l = 0;
	321	int r = 0;
	322	int m = 0;
	323
	324	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
	325	n = (int)Math.Round(c[2]);
	326
	327	//
	328	// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
	329	//
	330	l = 3;
	331	r = 3+n-2+1;
	332	while( l!=r-1 )
	333	{
	334	m = (l+r)/2;
	335	if( (double)(c[m])>=(double)(x) )
	336	{
	337	r = m;
	338	}
	339	else
	340	{
	341	l = m;
	342	}
	343	}
	344
	345	//
	346	// Interpolation
	347	//
	348	x = x-c[l];
	349	m = 3+n+4*(l-3);
	350	result = c[m]+x(c[m+1]+x(c[m+2]+x*c[m+3]));
	351	return result;
	352	}
	353
	354
	355	public static void splinedifferentiation(ref double[] c,
	356	double x,
	357	ref double s,
	358	ref double ds,
	359	ref double d2s)
	360	{
	361	int n = 0;
	362	int l = 0;
	363	int r = 0;
	364	int m = 0;
	365
	366	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineInterpolation: incorrect C!");
	367	n = (int)Math.Round(c[2]);
	368
	369	//
	370	// Binary search
	371	//
	372	l = 3;
	373	r = 3+n-2+1;
	374	while( l!=r-1 )
	375	{
	376	m = (l+r)/2;
	377	if( (double)(c[m])>=(double)(x) )
	378	{
	379	r = m;
	380	}
	381	else
	382	{
	383	l = m;
	384	}
	385	}
	386
	387	//
	388	// Differentiation
	389	//
	390	x = x-c[l];
	391	m = 3+n+4*(l-3);
	392	s = c[m]+x(c[m+1]+x(c[m+2]+x*c[m+3]));
	393	ds = c[m+1]+2xc[m+2]+3AP.Math.Sqr(x)c[m+3];
	394	d2s = 2c[m+2]+6x*c[m+3];
	395	}
	396
	397
	398	public static void splinecopy(ref double[] c,
	399	ref double[] cc)
	400	{
	401	int s = 0;
	402	int i_ = 0;
	403
	404	s = (int)Math.Round(c[0]);
	405	cc = new double[s-1+1];
	406	for(i_=0; i_<=s-1;i_++)
	407	{
	408	cc[i_] = c[i_];
	409	}
	410	}
	411
	412
	413	public static void splineunpack(ref double[] c,
	414	ref int n,
	415	ref double[,] tbl)
	416	{
	417	int i = 0;
	418
	419	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineUnpack: incorrect C!");
	420	n = (int)Math.Round(c[2]);
	421	tbl = new double[n-2+1, 5+1];
	422
	423	//
	424	// Fill
	425	//
	426	for(i=0; i<=n-2; i++)
	427	{
	428	tbl[i,0] = c[3+i];
	429	tbl[i,1] = c[3+i+1];
	430	tbl[i,2] = c[3+n+4*i];
	431	tbl[i,3] = c[3+n+4*i+1];
	432	tbl[i,4] = c[3+n+4*i+2];
	433	tbl[i,5] = c[3+n+4*i+3];
	434	}
	435	}
	436
	437
	438	public static void splinelintransx(ref double[] c,
	439	double a,
	440	double b)
	441	{
	442	int i = 0;
	443	int n = 0;
	444	double v = 0;
	445	double dv = 0;
	446	double d2v = 0;
	447	double[] x = new double[0];
	448	double[] y = new double[0];
	449	double[] d = new double[0];
	450
	451	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
	452	n = (int)Math.Round(c[2]);
	453
	454	//
	455	// Special case: A=0
	456	//
	457	if( (double)(a)==(double)(0) )
	458	{
	459	v = splineinterpolation(ref c, b);
	460	for(i=0; i<=n-2; i++)
	461	{
	462	c[3+n+4*i] = v;
	463	c[3+n+4*i+1] = 0;
	464	c[3+n+4*i+2] = 0;
	465	c[3+n+4*i+3] = 0;
	466	}
	467	return;
	468	}
	469
	470	//
	471	// General case: A<>0.
	472	// Unpack, X, Y, dY/dX.
	473	// Scale and pack again.
	474	//
	475	x = new double[n-1+1];
	476	y = new double[n-1+1];
	477	d = new double[n-1+1];
	478	for(i=0; i<=n-1; i++)
	479	{
	480	x[i] = c[3+i];
	481	splinedifferentiation(ref c, x[i], ref v, ref dv, ref d2v);
	482	x[i] = (x[i]-b)/a;
	483	y[i] = v;
	484	d[i] = a*dv;
	485	}
	486	buildhermitespline(x, y, d, n, ref c);
	487	}
	488
	489
	490	public static void splinelintransy(ref double[] c,
	491	double a,
	492	double b)
	493	{
	494	int i = 0;
	495	int n = 0;
	496	double v = 0;
	497	double dv = 0;
	498	double d2v = 0;
	499	double[] x = new double[0];
	500	double[] y = new double[0];
	501	double[] d = new double[0];
	502
	503	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineLinTransX: incorrect C!");
	504	n = (int)Math.Round(c[2]);
	505
	506	//
	507	// Special case: A=0
	508	//
	509	for(i=0; i<=n-2; i++)
	510	{
	511	c[3+n+4i] = ac[3+n+4*i]+b;
	512	c[3+n+4i+1] = ac[3+n+4*i+1];
	513	c[3+n+4i+2] = ac[3+n+4*i+2];
	514	c[3+n+4i+3] = ac[3+n+4*i+3];
	515	}
	516	}
	517
	518
	519	public static double splineintegration(ref double[] c,
	520	double x)
	521	{
	522	double result = 0;
	523	int n = 0;
	524	int i = 0;
	525	int l = 0;
	526	int r = 0;
	527	int m = 0;
	528	double w = 0;
	529
	530	System.Diagnostics.Debug.Assert((int)Math.Round(c[1])==3, "SplineIntegration: incorrect C!");
	531	n = (int)Math.Round(c[2]);
	532
	533	//
	534	// Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
	535	//
	536	l = 3;
	537	r = 3+n-2+1;
	538	while( l!=r-1 )
	539	{
	540	m = (l+r)/2;
	541	if( (double)(c[m])>=(double)(x) )
	542	{
	543	r = m;
	544	}
	545	else
	546	{
	547	l = m;
	548	}
	549	}
	550
	551	//
	552	// Integration
	553	//
	554	result = 0;
	555	for(i=3; i<=l-1; i++)
	556	{
	557	w = c[i+1]-c[i];
	558	m = 3+n+4*(i-3);
	559	result = result+c[m]*w;
	560	result = result+c[m+1]*AP.Math.Sqr(w)/2;
	561	result = result+c[m+2]AP.Math.Sqr(w)w/3;
	562	result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
	563	}
	564	w = x-c[l];
	565	m = 3+n+4*(l-3);
	566	result = result+c[m]*w;
	567	result = result+c[m+1]*AP.Math.Sqr(w)/2;
	568	result = result+c[m+2]AP.Math.Sqr(w)w/3;
	569	result = result+c[m+3]*AP.Math.Sqr(AP.Math.Sqr(w))/4;
	570	return result;
	571	}
	572
	573
	574	public static void spline3buildtable(int n,
	575	int diffn,
	576	double[] x,
	577	double[] y,
	578	double boundl,
	579	double boundr,
	580	ref double[,] ctbl)
	581	{
	582	bool c = new bool();
	583	int e = 0;
	584	int g = 0;
	585	double tmp = 0;
	586	int nxm1 = 0;
	587	int i = 0;
	588	int j = 0;
	589	double dx = 0;
	590	double dxj = 0;
	591	double dyj = 0;
	592	double dxjp1 = 0;
	593	double dyjp1 = 0;
	594	double dxp = 0;
	595	double dyp = 0;
	596	double yppa = 0;
	597	double yppb = 0;
	598	double pj = 0;
	599	double b1 = 0;
	600	double b2 = 0;
	601	double b3 = 0;
	602	double b4 = 0;
	603
	604	x = (double[])x.Clone();
	605	y = (double[])y.Clone();
	606
	607	n = n-1;
	608	g = (n+1)/2;
	609	do
	610	{
	611	i = g;
	612	do
	613	{
	614	j = i-g;
	615	c = true;
	616	do
	617	{
	618	if( (double)(x[j])<=(double)(x[j+g]) )
	619	{
	620	c = false;
	621	}
	622	else
	623	{
	624	tmp = x[j];
	625	x[j] = x[j+g];
	626	x[j+g] = tmp;
	627	tmp = y[j];
	628	y[j] = y[j+g];
	629	y[j+g] = tmp;
	630	}
	631	j = j-1;
	632	}
	633	while( j>=0 & c );
	634	i = i+1;
	635	}
	636	while( i<=n );
	637	g = g/2;
	638	}
	639	while( g>0 );
	640	ctbl = new double[4+1, n+1];
	641	n = n+1;
	642	if( diffn==1 )
	643	{
	644	b1 = 1;
	645	b2 = 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-boundl);
	646	b3 = 1;
	647	b4 = 6/(x[n-1]-x[n-2])*(boundr-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
	648	}
	649	else
	650	{
	651	b1 = 0;
	652	b2 = 2*boundl;
	653	b3 = 0;
	654	b4 = 2*boundr;
	655	}
	656	nxm1 = n-1;
	657	if( n>=2 )
	658	{
	659	if( n>2 )
	660	{
	661	dxj = x[1]-x[0];
	662	dyj = y[1]-y[0];
	663	j = 2;
	664	while( j<=nxm1 )
	665	{
	666	dxjp1 = x[j]-x[j-1];
	667	dyjp1 = y[j]-y[j-1];
	668	dxp = dxj+dxjp1;
	669	ctbl[1,j-1] = dxjp1/dxp;
	670	ctbl[2,j-1] = 1-ctbl[1,j-1];
	671	ctbl[3,j-1] = 6*(dyjp1/dxjp1-dyj/dxj)/dxp;
	672	dxj = dxjp1;
	673	dyj = dyjp1;
	674	j = j+1;
	675	}
	676	}
	677	ctbl[1,0] = -(b1/2);
	678	ctbl[2,0] = b2/2;
	679	if( n!=2 )
	680	{
	681	j = 2;
	682	while( j<=nxm1 )
	683	{
	684	pj = ctbl[2,j-1]*ctbl[1,j-2]+2;
	685	ctbl[1,j-1] = -(ctbl[1,j-1]/pj);
	686	ctbl[2,j-1] = (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,j-2])/pj;
	687	j = j+1;
	688	}
	689	}
	690	yppb = (b4-b3ctbl[2,nxm1-1])/(b3ctbl[1,nxm1-1]+2);
	691	i = 1;
	692	while( i<=nxm1 )
	693	{
	694	j = n-i;
	695	yppa = ctbl[1,j-1]*yppb+ctbl[2,j-1];
	696	dx = x[j]-x[j-1];
	697	ctbl[3,j-1] = (yppb-yppa)/dx/6;
	698	ctbl[2,j-1] = yppa/2;
	699	ctbl[1,j-1] = (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]dx)dx;
	700	yppb = yppa;
	701	i = i+1;
	702	}
	703	for(i=1; i<=n; i++)
	704	{
	705	ctbl[0,i-1] = y[i-1];
	706	ctbl[4,i-1] = x[i-1];
	707	}
	708	}
	709	}
	710
	711
	712	public static double spline3interpolate(int n,
	713	ref double[,] c,
	714	double x)
	715	{
	716	double result = 0;
	717	int i = 0;
	718	int l = 0;
	719	int half = 0;
	720	int first = 0;
	721	int middle = 0;
	722
	723	n = n-1;
	724	l = n;
	725	first = 0;
	726	while( l>0 )
	727	{
	728	half = l/2;
	729	middle = first+half;
	730	if( (double)(c[4,middle])<(double)(x) )
	731	{
	732	first = middle+1;
	733	l = l-half-1;
	734	}
	735	else
	736	{
	737	l = half;
	738	}
	739	}
	740	i = first-1;
	741	if( i<0 )
	742	{
	743	i = 0;
	744	}
	745	result = c[0,i]+(x-c[4,i])(c[1,i]+(x-c[4,i])(c[2,i]+c[3,i]*(x-c[4,i])));
	746	return result;
	747	}
	748
	749
	750	private static void heapsortpoints(ref double[] x,
	751	ref double[] y,
	752	int n)
	753	{
	754	int i = 0;
	755	int j = 0;
	756	int k = 0;
	757	int t = 0;
	758	double tmp = 0;
	759	bool isascending = new bool();
	760	bool isdescending = new bool();
	761
	762
	763	//
	764	// Test for already sorted set
	765	//
	766	isascending = true;
	767	isdescending = true;
	768	for(i=1; i<=n-1; i++)
	769	{
	770	isascending = isascending & (double)(x[i])>(double)(x[i-1]);
	771	isdescending = isdescending & (double)(x[i])<(double)(x[i-1]);
	772	}
	773	if( isascending )
	774	{
	775	return;
	776	}
	777	if( isdescending )
	778	{
	779	for(i=0; i<=n-1; i++)
	780	{
	781	j = n-1-i;
	782	if( j<=i )
	783	{
	784	break;
	785	}
	786	tmp = x[i];
	787	x[i] = x[j];
	788	x[j] = tmp;
	789	tmp = y[i];
	790	y[i] = y[j];
	791	y[j] = tmp;
	792	}
	793	return;
	794	}
	795
	796	//
	797	// Special case: N=1
	798	//
	799	if( n==1 )
	800	{
	801	return;
	802	}
	803
	804	//
	805	// General case
	806	//
	807	i = 2;
	808	do
	809	{
	810	t = i;
	811	while( t!=1 )
	812	{
	813	k = t/2;
	814	if( (double)(x[k-1])>=(double)(x[t-1]) )
	815	{
	816	t = 1;
	817	}
	818	else
	819	{
	820	tmp = x[k-1];
	821	x[k-1] = x[t-1];
	822	x[t-1] = tmp;
	823	tmp = y[k-1];
	824	y[k-1] = y[t-1];
	825	y[t-1] = tmp;
	826	t = k;
	827	}
	828	}
	829	i = i+1;
	830	}
	831	while( i<=n );
	832	i = n-1;
	833	do
	834	{
	835	tmp = x[i];
	836	x[i] = x[0];
	837	x[0] = tmp;
	838	tmp = y[i];
	839	y[i] = y[0];
	840	y[0] = tmp;
	841	t = 1;
	842	while( t!=0 )
	843	{
	844	k = 2*t;
	845	if( k>i )
	846	{
	847	t = 0;
	848	}
	849	else
	850	{
	851	if( k<i )
	852	{
	853	if( (double)(x[k])>(double)(x[k-1]) )
	854	{
	855	k = k+1;
	856	}
	857	}
	858	if( (double)(x[t-1])>=(double)(x[k-1]) )
	859	{
	860	t = 0;
	861	}
	862	else
	863	{
	864	tmp = x[k-1];
	865	x[k-1] = x[t-1];
	866	x[t-1] = tmp;
	867	tmp = y[k-1];
	868	y[k-1] = y[t-1];
	869	y[t-1] = tmp;
	870	t = k;
	871	}
	872	}
	873	}
	874	i = i-1;
	875	}
	876	while( i>=1 );
	877	}
	878
	879
	880	private static void heapsortdpoints(ref double[] x,
	881	ref double[] y,
	882	ref double[] d,
	883	int n)
	884	{
	885	int i = 0;
	886	int j = 0;
	887	int k = 0;
	888	int t = 0;
	889	double tmp = 0;
	890	bool isascending = new bool();
	891	bool isdescending = new bool();
	892
	893
	894	//
	895	// Test for already sorted set
	896	//
	897	isascending = true;
	898	isdescending = true;
	899	for(i=1; i<=n-1; i++)
	900	{
	901	isascending = isascending & (double)(x[i])>(double)(x[i-1]);
	902	isdescending = isdescending & (double)(x[i])<(double)(x[i-1]);
	903	}
	904	if( isascending )
	905	{
	906	return;
	907	}
	908	if( isdescending )
	909	{
	910	for(i=0; i<=n-1; i++)
	911	{
	912	j = n-1-i;
	913	if( j<=i )
	914	{
	915	break;
	916	}
	917	tmp = x[i];
	918	x[i] = x[j];
	919	x[j] = tmp;
	920	tmp = y[i];
	921	y[i] = y[j];
	922	y[j] = tmp;
	923	tmp = d[i];
	924	d[i] = d[j];
	925	d[j] = tmp;
	926	}
	927	return;
	928	}
	929
	930	//
	931	// Special case: N=1
	932	//
	933	if( n==1 )
	934	{
	935	return;
	936	}
	937
	938	//
	939	// General case
	940	//
	941	i = 2;
	942	do
	943	{
	944	t = i;
	945	while( t!=1 )
	946	{
	947	k = t/2;
	948	if( (double)(x[k-1])>=(double)(x[t-1]) )
	949	{
	950	t = 1;
	951	}
	952	else
	953	{
	954	tmp = x[k-1];
	955	x[k-1] = x[t-1];
	956	x[t-1] = tmp;
	957	tmp = y[k-1];
	958	y[k-1] = y[t-1];
	959	y[t-1] = tmp;
	960	tmp = d[k-1];
	961	d[k-1] = d[t-1];
	962	d[t-1] = tmp;
	963	t = k;
	964	}
	965	}
	966	i = i+1;
	967	}
	968	while( i<=n );
	969	i = n-1;
	970	do
	971	{
	972	tmp = x[i];
	973	x[i] = x[0];
	974	x[0] = tmp;
	975	tmp = y[i];
	976	y[i] = y[0];
	977	y[0] = tmp;
	978	tmp = d[i];
	979	d[i] = d[0];
	980	d[0] = tmp;
	981	t = 1;
	982	while( t!=0 )
	983	{
	984	k = 2*t;
	985	if( k>i )
	986	{
	987	t = 0;
	988	}
	989	else
	990	{
	991	if( k<i )
	992	{
	993	if( (double)(x[k])>(double)(x[k-1]) )
	994	{
	995	k = k+1;
	996	}
	997	}
	998	if( (double)(x[t-1])>=(double)(x[k-1]) )
	999	{
	1000	t = 0;
	1001	}
	1002	else
	1003	{
	1004	tmp = x[k-1];
	1005	x[k-1] = x[t-1];
	1006	x[t-1] = tmp;
	1007	tmp = y[k-1];
	1008	y[k-1] = y[t-1];
	1009	y[t-1] = tmp;
	1010	tmp = d[k-1];
	1011	d[k-1] = d[t-1];
	1012	d[t-1] = tmp;
	1013	t = k;
	1014	}
	1015	}
	1016	}
	1017	i = i-1;
	1018	}
	1019	while( i>=1 );
	1020	}
	1021
	1022
	1023	private static void solvetridiagonal(double[] a,
	1024	double[] b,
	1025	double[] c,
	1026	double[] d,
	1027	int n,
	1028	ref double[] x)
	1029	{
	1030	int k = 0;
	1031	double t = 0;
	1032
	1033	a = (double[])a.Clone();
	1034	b = (double[])b.Clone();
	1035	c = (double[])c.Clone();
	1036	d = (double[])d.Clone();
	1037
	1038	x = new double[n-1+1];
	1039	a[0] = 0;
	1040	c[n-1] = 0;
	1041	for(k=1; k<=n-1; k++)
	1042	{
	1043	t = a[k]/b[k-1];
	1044	b[k] = b[k]-t*c[k-1];
	1045	d[k] = d[k]-t*d[k-1];
	1046	}
	1047	x[n-1] = d[n-1]/b[n-1];
	1048	for(k=n-2; k>=0; k--)
	1049	{
	1050	x[k] = (d[k]-c[k]*x[k+1])/b[k];
	1051	}
	1052	}
	1053
	1054
	1055	private static double diffthreepoint(double t,
	1056	double x0,
	1057	double f0,
	1058	double x1,
	1059	double f1,
	1060	double x2,
	1061	double f2)
	1062	{
	1063	double result = 0;
	1064	double a = 0;
	1065	double b = 0;
	1066
	1067	t = t-x0;
	1068	x1 = x1-x0;
	1069	x2 = x2-x0;
	1070	a = (f2-f0-x2/x1(f1-f0))/(AP.Math.Sqr(x2)-x1x2);
	1071	b = (f1-f0-a*AP.Math.Sqr(x1))/x1;
	1072	result = 2at+b;
	1073	return result;
	1074	}
	1075	}
	1076	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences