Luxophia

LUX.Curve.Poly.pas
553 строки · 13.7 Кб
Перенос по словам
1
unit LUX.Curve.Poly;
2

3
interface //#################################################################### ■
4

5
uses LUX,
6
     LUX.D1,
7
     LUX.D2,
8
     LUX.D2x2,
9
     LUX.D3,
10
     LUX.D3x3,
11
     LUX.D4,
12
     LUX.M4,
13
     LUX.D5,
14
     LUX.DN;
15

16
//type //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【型】
17

18
     //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【レコード】
19

20
     //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【クラス】
21

22
//const //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【定数】
23

24
//var //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【変数】
25

26
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【ルーチン】
27

28
function Legendre( const X_:Single; const N_:Cardinal ) :Single; overload;
29
function Legendre( const X_:Double; const N_:Cardinal ) :Double; overload;
30

31
function Chebyshev1( const X_:Single; const N_:Cardinal ) :Single; overload;
32
function Chebyshev1( const X_:Double; const N_:Cardinal ) :Double; overload;
33

34
function SumLegendre( const X_:Single; const Ws_:array of Single ) :Single; overload;
35
function SumLegendre( const X_:Double; const Ws_:array of Double ) :Double; overload;
36

37
function SumChebyshev1( const X_:Single; const Ws_:array of Single ) :Single; overload;
38
function SumChebyshev1( const X_:Double; const Ws_:array of Double ) :Double; overload;
39

40
function Poly( const X_:Single; const Ks_:TSingle2D ) :Single; overload;
41
function Poly( const X_:Double; const Ks_:TDouble2D ) :Double; overload;
42

43
function Poly( const X_:Single; const Ks_:TSingle3D ) :Single; overload;
44
function Poly( const X_:Double; const Ks_:TDouble3D ) :Double; overload;
45

46
function Poly( const X_:Single; const Ks_:TSingle4D ) :Single; overload;
47
function Poly( const X_:Double; const Ks_:TDouble4D ) :Double; overload;
48

49
function Poly( const X_:Single; const Ks_:TSingle5D ) :Single; overload;
50
function Poly( const X_:Double; const Ks_:TDouble5D ) :Double; overload;
51

52
function Poly( const X_:Single; const Ks_:TSingleND ) :Single; overload;
53
function Poly( const X_:Double; const Ks_:TDoubleND ) :Double; overload;
54

55
function FitPoly( const P_:TSingle4D ) :TSingle4D; overload;
56
function FitPoly( const P_:TDouble4D ) :TDouble4D; overload;
57

58
function FitPoly( const P_:TSingle5D ) :TSingle5D; overload;
59
function FitPoly( const P_:TDouble5D ) :TDouble5D; overload;
60

61
procedure RandPoly( out Ks_:TSingle4D ); overload;
62
procedure RandPoly( out Ks_:TDouble4D ); overload;
63

64
procedure RandPoly( out Ks_:TSingle5D ); overload;
65
procedure RandPoly( out Ks_:TDouble5D ); overload;
66

67
implementation //############################################################### ■
68

69
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【レコード】
70

71
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【クラス】
72

73
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$【ルーチン】
74

75
function Legendre( const X_:Single; const N_:Cardinal ) :Single;
76
var
77
   I :Cardinal;
78
   P0, P1, P2 :Single;
79
begin
80
     if N_ = 0 then Result := 1
81
     else
82
     begin
83
          P1 := 1 ;
84
          P2 := X_;
85
          for I := 2 to N_ do
86
          begin
87
               P0 := P1;  P1 := P2;
88

89
               P2 := ( ( 2 * I - 1 ) * X_ * P1 - ( I - 1 ) * P0 ) / I;
90
          end;
91

92
          Result := P2;
93
     end;
94
end;
95

96
function Legendre( const X_:Double; const N_:Cardinal ) :Double;
97
var
98
   I :Cardinal;
99
   P0, P1, P2 :Double;
100
begin
101
     if N_ = 0 then Result := 1
102
     else
103
     begin
104
          P1 := 1 ;
105
          P2 := X_;
106
          for I := 2 to N_ do
107
          begin
108
               P0 := P1;  P1 := P2;
109

110
               P2 := ( ( 2 * I - 1 ) * X_ * P1 - ( I - 1 ) * P0 ) / I;
111
          end;
112

113
          Result := P2;
114
     end;
115
end;
116

117
//------------------------------------------------------------------------------
118

119
function Chebyshev1( const X_:Single; const N_:Cardinal ) :Single;
120
var
121
   I :Cardinal;
122
   T0, T1, T2 :Single;
123
begin
124
     if N_ = 0 then Result := 1
125
     else
126
     begin
127
          T1 := 1 ;
128
          T2 := X_;
129
          for I := 2 to N_ do
130
          begin
131
               T0 := T1;  T1 := T2;
132

133
               T2 := 2 * X_ * T1 - T0;
134
          end;
135

136
          Result := T2;
137
     end;
138
end;
139

140
function Chebyshev1( const X_:Double; const N_:Cardinal ) :Double;
141
var
142
   I :Cardinal;
143
   T0, T1, T2 :Double;
144
begin
145
     if N_ = 0 then Result := 1
146
     else
147
     begin
148
          T1 := 1 ;
149
          T2 := X_;
150
          for I := 2 to N_ do
151
          begin
152
               T0 := T1;  T1 := T2;
153

154
               T2 := 2 * X_ * T1 - T0;
155
          end;
156

157
          Result := T2;
158
     end;
159
end;
160

161
//------------------------------------------------------------------------------
162

163
function SumLegendre( const X_:Single; const Ws_:array of Single ) :Single;
164
var
165
   H, I :Integer;
166
   P0, P1, P2 :Single;
167
begin
168
     Result := Ws_[ 0 ];
169

170
     H := High( Ws_ );
171

172
     if H > 0 then
173
     begin
174
          P1 := 1 ;
175
          P2 := X_;
176

177
          Result := Result + Ws_[ 1 ] * P2;
178

179
          for I := 2 to H do
180
          begin
181
               P0 := P1;  P1 := P2;
182

183
               P2 := ( ( 2 * I - 1 ) * X_ * P1 - ( I - 1 ) * P0 ) / I;
184

185
               Result := Result + Ws_[ I ] * P2;
186
          end;
187
     end;
188
end;
189

190
function SumLegendre( const X_:Double; const Ws_:array of Double ) :Double;
191
var
192
   H, I :Integer;
193
   P0, P1, P2 :Double;
194
begin
195
     Result := Ws_[ 0 ];
196

197
     H := High( Ws_ );
198

199
     if H > 0 then
200
     begin
201
          P1 := 1 ;
202
          P2 := X_;
203

204
          Result := Result + Ws_[ 1 ] * P2;
205

206
          for I := 2 to H do
207
          begin
208
               P0 := P1;  P1 := P2;
209

210
               P2 := ( ( 2 * I - 1 ) * X_ * P1 - ( I - 1 ) * P0 ) / I;
211

212
               Result := Result + Ws_[ I ] * P2;
213
          end;
214
     end;
215
end;
216

217
//------------------------------------------------------------------------------
218

219
function SumChebyshev1( const X_:Single; const Ws_:array of Single ) :Single;
220
var
221
   H, I :Integer;
222
   P0, P1, P2 :Single;
223
begin
224
     Result := Ws_[ 0 ];
225

226
     H := High( Ws_ );
227

228
     if H > 0 then
229
     begin
230
          P1 := 1 ;
231
          P2 := X_;
232

233
          Result := Result + Ws_[ 1 ] * P2;
234

235
          for I := 2 to H do
236
          begin
237
               P0 := P1;  P1 := P2;
238

239
               P2 := 2 * X_ * P1 - P0;
240

241
               Result := Result + Ws_[ I ] * P2;
242
          end;
243
     end;
244
end;
245

246
function SumChebyshev1( const X_:Double; const Ws_:array of Double ) :Double;
247
var
248
   H, I :Integer;
249
   P0, P1, P2 :Double;
250
begin
251
     Result := Ws_[ 0 ];
252

253
     H := High( Ws_ );
254

255
     if H > 0 then
256
     begin
257
          P1 := 1 ;
258
          P2 := X_;
259

260
          Result := Result + Ws_[ 1 ] * P2;
261

262
          for I := 2 to H do
263
          begin
264
               P0 := P1;  P1 := P2;
265

266
               P2 := 2 * X_ * P1 - P0;
267

268
               Result := Result + Ws_[ I ] * P2;
269
          end;
270
     end;
271
end;
272

273
//------------------------------------------------------------------------------
274

275
function Poly( const X_:Single; const Ks_:TSingle2D ) :Single;
276
begin
277
     with Ks_ do Result := _2 * X_ + _1;
278
end;
279

280
function Poly( const X_:Double; const Ks_:TDouble2D ) :Double;
281
begin
282
     with Ks_ do Result := _2 * X_ + _1;
283
end;
284

285
//------------------------------------------------------------------------------
286

287
function Poly( const X_:Single; const Ks_:TSingle3D ) :Single;
288
begin
289
     with Ks_ do Result := ( _3 * X_ + _2 ) * X_ + _1;
290
end;
291

292
function Poly( const X_:Double; const Ks_:TDouble3D ) :Double;
293
begin
294
     with Ks_ do Result := ( _3 * X_ + _2 ) * X_ + _1;
295
end;
296

297
//------------------------------------------------------------------------------
298

299
function Poly( const X_:Single; const Ks_:TSingle4D ) :Single;
300
begin
301
     with Ks_ do Result := ( ( _4 * X_ + _3 ) * X_ + _2 ) * X_ + _1;
302
end;
303

304
function Poly( const X_:Double; const Ks_:TDouble4D ) :Double;
305
begin
306
     with Ks_ do Result := ( ( _4 * X_ + _3 ) * X_ + _2 ) * X_ + _1;
307
end;
308

309
//------------------------------------------------------------------------------
310

311
function Poly( const X_:Single; const Ks_:TSingle5D ) :Single;
312
begin
313
     with Ks_ do Result := ( ( ( _5 * X_ + _4 ) * X_ + _3 ) * X_ + _2 ) * X_ + _1;
314
end;
315

316
function Poly( const X_:Double; const Ks_:TDouble5D ) :Double;
317
begin
318
     with Ks_ do Result := ( ( ( _5 * X_ + _4 ) * X_ + _3 ) * X_ + _2 ) * X_ + _1;
319
end;
320

321
//------------------------------------------------------------------------------
322

323
function Poly( const X_:Single; const Ks_:TSingleND ) :Single;
324
var
325
   I :Integer;
326
begin
327
     Result := Ks_[ Ks_.DimN-1 ];
328

329
     for I := Ks_.DimN-2 downto 0 do Result := Result * X_ + Ks_[ I ];
330
end;
331

332
function Poly( const X_:Double; const Ks_:TDoubleND ) :Double;
333
var
334
   I :Integer;
335
begin
336
     Result := Ks_[ Ks_.DimN-1 ];
337

338
     for I := Ks_.DimN-2 downto 0 do Result := Result * X_ + Ks_[ I ];
339
end;
340

341
//------------------------------------------------------------------------------
342

343
function FitPoly( const P_:TSingle4D ) :TSingle4D;
344
const
345
     M :TSingleM3 = ( _11:+09  ;  _12:-09/2;  _13:+01  ;
346
                      _21:-45/2;  _22:+18  ;  _23:-09/2;
347
                      _31:+27/2;  _32:-27/2;  _33:+09/2; );
348
var
349
   P, K :TSingle3D;
350
begin
351
     with P_ do
352
     begin
353
          P._1 := _2 - _1;
354
          P._2 := _3 - _1;
355
          P._3 := _4 - _1;
356
     end;
357

358
     K := M * P;
359

360
     with Result do
361
     begin
362
          _1 := P_._1;
363
          _2 := K ._1;
364
          _3 := K ._2;
365
          _4 := K ._3;
366
     end;
367
end;
368

369
function FitPoly( const P_:TDouble4D ) :TDouble4D;
370
const
371
     M :TDoubleM3 = ( _11:+09  ;  _12:-09/2;  _13:+01  ;
372
                      _21:-45/2;  _22:+18  ;  _23:-09/2;
373
                      _31:+27/2;  _32:-27/2;  _33:+09/2; );
374
var
375
   P, K :TDouble3D;
376
begin
377
     with P_ do
378
     begin
379
          P._1 := _2 - _1;
380
          P._2 := _3 - _1;
381
          P._3 := _4 - _1;
382
     end;
383

384
     K := M * P;
385

386
     with Result do
387
     begin
388
          _1 := P_._1;
389
          _2 := K ._1;
390
          _3 := K ._2;
391
          _4 := K ._3;
392
     end;
393
end;
394

395
//------------------------------------------------------------------------------
396

397
function FitPoly( const P_:TSingle5D ) :TSingle5D;
398
const
399
     M :TSingleM4 = ( _11:+016  ;  _12:-012  ;  _13:+016/3;  _14:-001  ;
400
                      _21:-208/3;  _22:+076  ;  _23:-112/3;  _24:+022/3;
401
                      _31:+096  ;  _32:-128  ;  _33:+224/3;  _34:-016  ;
402
                      _41:-128/3;  _42:+064  ;  _43:-128/3;  _44:+032/3; );
403
var
404
   P, K :TSingle4D;
405
begin
406
     with P_ do
407
     begin
408
          P._1 := _2 - _1;
409
          P._2 := _3 - _1;
410
          P._3 := _4 - _1;
411
          P._4 := _5 - _1;
412
     end;
413

414
     K := M * P;
415

416
     with Result do
417
     begin
418
          _1 := P_._1;
419
          _2 := K ._1;
420
          _3 := K ._2;
421
          _4 := K ._3;
422
          _5 := K ._4;
423
     end;
424
end;
425

426
function FitPoly( const P_:TDouble5D ) :TDouble5D;
427
const
428
     M :TDoubleM4 = ( _11:+016  ;  _12:-012  ;  _13:+016/3;  _14:-001  ;
429
                      _21:-208/3;  _22:+076  ;  _23:-112/3;  _24:+022/3;
430
                      _31:+096  ;  _32:-128  ;  _33:+224/3;  _34:-016  ;
431
                      _41:-128/3;  _42:+064  ;  _43:-128/3;  _44:+032/3; );
432
var
433
   P, K :TDouble4D;
434
begin
435
     with P_ do
436
     begin
437
          P._1 := _2 - _1;
438
          P._2 := _3 - _1;
439
          P._3 := _4 - _1;
440
          P._4 := _5 - _1;
441
     end;
442

443
     K := M * P;
444

445
     with Result do
446
     begin
447
          _1 := P_._1;
448
          _2 := K ._1;
449
          _3 := K ._2;
450
          _4 := K ._3;
451
          _5 := K ._4;
452
     end;
453
end;
454

455
//------------------------------------------------------------------------------
456

457
procedure RandPoly( out Ks_:TSingle4D );
458
var
459
   Ws :array [ 0..3 ] of Single;
460
   Ps :TSingle4D;
461
begin
462
     Ws[ 0 ] := 2 * Random - 1;
463
     Ws[ 1 ] := 2 * Random - 1;
464
     Ws[ 2 ] := 2 * Random - 1;
465
     Ws[ 3 ] := 2 * Random - 1;
466

467
     with Ps do
468
     begin
469
          _1 := SumChebyshev1( 0/3, Ws );
470
          _2 := SumChebyshev1( 1/3, Ws );
471
          _3 := SumChebyshev1( 2/3, Ws );
472
          _4 := SumChebyshev1( 3/3, Ws );
473
     end;
474

475
     Ks_ := FitPoly( Ps / 2 );
476
end;
477

478
procedure RandPoly( out Ks_:TDouble4D );
479
var
480
   Ws :array [ 0..3 ] of Double;
481
   Ps :TDouble4D;
482
begin
483
     Ws[ 0 ] := 2 * Random - 1;
484
     Ws[ 1 ] := 2 * Random - 1;
485
     Ws[ 2 ] := 2 * Random - 1;
486
     Ws[ 3 ] := 2 * Random - 1;
487

488
     with Ps do
489
     begin
490
          _1 := SumChebyshev1( 0/3, Ws );
491
          _2 := SumChebyshev1( 1/3, Ws );
492
          _3 := SumChebyshev1( 2/3, Ws );
493
          _4 := SumChebyshev1( 3/3, Ws );
494
     end;
495

496
     Ks_ := FitPoly( Ps / 2 );
497
end;
498

499
//------------------------------------------------------------------------------
500

501
procedure RandPoly( out Ks_:TSingle5D );
502
var
503
   Ws :array [ 0..4 ] of Single;
504
   Ps :TSingle5D;
505
begin
506
     Ws[ 0 ] := 2 * Random - 1;
507
     Ws[ 1 ] := 2 * Random - 1;
508
     Ws[ 2 ] := 2 * Random - 1;
509
     Ws[ 3 ] := 2 * Random - 1;
510
     Ws[ 4 ] := 2 * Random - 1;
511

512
     with Ps do
513
     begin
514
          _1 := SumChebyshev1( 0/4, Ws );
515
          _2 := SumChebyshev1( 1/4, Ws );
516
          _3 := SumChebyshev1( 2/4, Ws );
517
          _4 := SumChebyshev1( 3/4, Ws );
518
          _5 := SumChebyshev1( 4/4, Ws );
519
     end;
520

521
     Ks_ := FitPoly( Ps / Roo2(5) );
522
end;
523

524
procedure RandPoly( out Ks_:TDouble5D );
525
var
526
   Ws :array [ 0..4 ] of Double;
527
   Ps :TDouble5D;
528
begin
529
     Ws[ 0 ] := 2 * Random - 1;
530
     Ws[ 1 ] := 2 * Random - 1;
531
     Ws[ 2 ] := 2 * Random - 1;
532
     Ws[ 3 ] := 2 * Random - 1;
533
     Ws[ 4 ] := 2 * Random - 1;
534

535
     with Ps do
536
     begin
537
          _1 := SumChebyshev1( 0/4, Ws );
538
          _2 := SumChebyshev1( 1/4, Ws );
539
          _3 := SumChebyshev1( 2/4, Ws );
540
          _4 := SumChebyshev1( 3/4, Ws );
541
          _5 := SumChebyshev1( 4/4, Ws );
542
     end;
543

544
     Ks_ := FitPoly( Ps / Roo2(5) );
545
end;
546

547
//############################################################################## □
548

549
initialization //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 初期化
550

551
finalization //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 最終化
552

553
end. //######################################################################### ■
554
Luxophia

Использование cookies
