{************************************************} { } { Turbo Pascal } { Programme IVEC } { Unit_ CMath } { Sources r_cup_r_s dans une librairie } { math_mathique en libre service aux _tats-unis} { Copyright (c) 1992 Odent Jean Philippe } { } {************************************************} Unit CMath; { The setting of the N compiler directive within this include file will } { determine whether your toolbox programs will use the standard Turbo } { Pascal 6 byte real number or the Turbo Pascal double precision real } { number. In order to use the double precision real number type, you } { must have an 8087 math coprocessor installed in your computer. } {$N-} { Change to $N+ if you want to use the 8 byte double precision real } { Change to $N- if you want to use the 6 byte non-8087 real number } interface uses Wintypes, WinProcs; { TPoint et MessageBeep } const PiRad = 180; PI = 3.14159265358979323846; { pi } PI2 = PI*2; {$IFOPT N+} type Float = Double; { 8 byte real, requires 8087 math chip } const Epsilon2 = 1e-6; Tolerance = 1E-4; Epsilon = 1E-015; Mant = 25; Deci = 25; {$ELSE} type Float = Real; { 6 byte real, no math chip required } const Tolerance = 1E-6; Epsilon2 = 1e-9; Epsilon = 1E-07; Mant = 15; Deci = 14; {$ENDIF} const DimenMat = 3; { Calcul sur des matrices 3*3 } TNArraySize = 4; { Size of the matrix } LimitTabVal = 4100; LimitCalculTabVal = 4003; MaxLongInt = 655350000; PartPi = 90; { Pas de rotation } Plus = 1; Moins = -1; Modul3 = 1.73205080757; { racine de 3 } DeciSurface = 3; CoefSurface = 5*5*0.000001; PasScrollHorizontal = 8; PasScrollVertical = 16; { 15 laisse une moirure } FlechePleine = 1; FlecheVide = 2; RAnatomique = 0; RAxesPropres = 1; RObservateur = 2; type TNvector = array[1..TNArraySize] of Float; TNmatrix = array[1..TNArraySize] of TNvector; TParamV = array[1..6] of real; P_TabVal = ^_TabVal; _TabVal = array[0..LimitTabVal] of Integer; P_TabFloat = ^_TabFloat; _TabFloat = array[0..LimitTabVal] of Float; P_TabPoints = ^_TabPoints; _TabPoints = array[0..LimitTabVal] of TPoint; P_TabPVal = ^_TabPVal; _TabPVal = array[1..3] of P_TabVal; TRectReel = record Left, Top, Bottom, Right : Float; end; RepereVector = record V : TNVector; Indice : Integer; Nom : pChar; end; TEchelle = record X, Y, Z : Real; end; function MathError:Integer; Function HexaSt(W:LongInt;n:byte):String; Function StHexa(St:String):LongInt; function ASin(x:Float):Float; function ACos(x:Float):Float; function Rad(x:Float):Float; function Deg(x:Float):Float; function Sgn(x:Float):Integer; function Sign(var N:Real):Integer; function DivNonZ(a,b : Float):Float; function Limite(N:Real):Real; function MaxI(A,B : LongInt):LongInt; function MinI(A,B : LongInt):LongInt; function MaxR(A,B : Float):Float; function MinR(A,B : Float):Float; function Atan2( y, x: Float):Float; function fabs(x: Float):Float; function floor(x:Float):Float; function powi(x : Float; nn:Integer):Float; function pow(x, y : Float):Float; function RoundR(x:Float):Float; function Sin(x:Float):Float; function Cos(x:Float):Float; function Exp(x:Float):Float; function Log(x:Float):Float; function Log10(x:Float):Float; function Exp10(x:Float):Float; function Cosh(x:Float):Float; function ACosh(x:Float):Float; function Sinh(x:Float):Float; function ASinh(x:Float):Float; function Cot(x:Float):Float; function Tan(x:Float):Float; function Atan(x:Float):Float; function Tanh(x:Float):Float; function ATanh(x:Float):Float; function sqrt(x:Float):Float; function cbrt(x:Float):Float; function NonZ(x:Float):Float; function Encadre(a, b, c: LongInt):LongInt; function Sqr3(a,b,c:Float):Float; function Module(var V:TNVector):Float; function Longitude(var V:TNVector):Float; function Latitude(var V:TNVector):Float; function Site(var V:TNVector):Float; function Elevation(var V:TNVector):Float; function Azimut(var V:TNVector):Float; procedure MesureVecteur(var Vecteur:TNVector; var Param:TParamV); procedure CoordonneeMatrice(var Mat:TNMatrix; var V:TNVector; IsRepereObj:Boolean); procedure OptimiseRotation(var Vn:TNVector); function VerifAngle(var V:TNVector; var Mat:TNMatrix):Boolean; procedure Ajuste(var V:TNVector; var Mat:TNMatrix); procedure Initial(Dimen : Integer; var Data: TNmatrix); procedure DivMatFloat(Dimen : integer; var Data : TNmatrix; Divisor : Float); procedure ClearVect(var Data : TNVector); procedure AddVecteur(Dimen:Integer; var A, B, C : TNVector); procedure SubVecteur(Dimen:Integer; var A, B, C : TNVector); procedure ClearMat(var Data : TNmatrix); procedure MultMat(Dimen : Integer; MatA, MatB:TNMatrix; var MatC:TNMatrix); procedure MultMatVect(Dimen : Integer; var Mat:TNMatrix; var V:TNVector); procedure TransposeMat(Dimen : Integer; var Data:TNMatrix); procedure InverseMatRot(DimenMat: Integer; var MatI:TNMatrix); procedure AffecteMat(p:Integer; Phi:Float; var Mat:TNMatrix); procedure InitialiseXmat; procedure RotationMatrice(var PhiRot:TNvector; var MatRot:TNMatrix; Repere:Byte); procedure IncrementeMatrice(NoAxe, Sens : Integer; var MatRot, MatVP:TNMatrix; Rotation:Boolean); procedure CSpline(var AA : _TabVal; N : integer; X1, Xm, Echelle : Float; var BB : _TabPoints; M, Offsx, Offsy : integer); procedure Bezier(var A : _TabVal; MaxContrPoints : integer; var B : _TabPoints; MaxIntPoints : integer); procedure CalculVecteurPropre(var EigenVectors : TNMatrix; var PVecteur : _TabPVal; Taille : Word); procedure Jacobi(Dimen : integer; Mat : TNmatrix; MaxIter : integer; Tolerance : Float; var Eigenvalues : TNvector; var Eigenvectors : TNmatrix; var Iter : integer; var Error : byte); function SurfaceOmbre(var Mat:TNMatrix; var Points:_TabPVal; Modulo, OffsetS, Count:Integer):float; { Les proc_dures de validation des r_sultats } type StErr20 = String[20]; var MathStrErreur : StErr20; XMat : array[1..6] of TNMatrix; RotationPropre : Boolean; RepereActuel : Byte; implementation uses Verif; procedure Jacobi; external 'IVECdll' index 3; procedure CalculVecteurPropre; external 'IVECdll' index 4; const CarHexa : set of char =['0'..'9','a'..'f','A'..'F']; BaseHexa = 16; AlphabHexa : array[0..15] of char = '0123456789ABCDEF'; const PIO2 = 1.57079632679489661923; { pi/2 } PIO4 = 7.85398163397448309616E-1; { pi/4 } SQRT2 = 1.41421356237309504880; { sqrt(2) } SQRTH = 7.07106781186547524401E-1; { sqrt(2)/2 } LOG2E = 1.4426950408889634073599; { 1/log(2) } SQ2OPI = 7.9788456080286535587989E-1; { sqrt( 2/pi )} LOGE2 = 6.93147180559945309417E-1; { log(2) } LOGSQ2 = 3.46573590279972654709E-1; { log(2)/2 } THPIO4 = 2.35619449019234492885; { 3*pi/4 } TWOOPI = 6.36619772367581343075535E-1; { 2/pi } lossth = 1.073741824e9; { This routine may be called to report one of the following error conditions (in the include file mconf.h). Mnemonic Value Significance } DOMAIN = 1; { argument domain error } SING = 2; { function singularity } OVERFLOW = 3; { overflow range error } UNDERFLOW = 4; { underflow range error } TLOSS = 5; { total loss of precision } PLOSS = 6; { partial loss of precision } EDOM = 33; { Unix domain error code } ERANGE = 34; { Unix range error code } MACHEP = 1.38777878078144567553E-17; { 2**-56 } MAXLOG = 8.8029691931113054295988E1; { log(2**127) } MINLOG = -8.872283911167299960540E1; { log(2**-128) } MAXNUM = 1.7e38; { MAXNUM = 1.701411834604692317316873e38; { 2**127 } type TabCoef = array[0..6] of Float; var MathErreur : Integer; function MathError:Integer; begin MathError := MathErreur; MathErreur := 0; end; procedure mtherr(St: StErr20; Error:Integer); begin if MathErreur = 0 then begin MathErreur := Error; case Error of DOMAIN : MathStrErreur := 'argument domain error'; SING : MathStrErreur := 'function singularity'; OVERFLOW : MathStrErreur := 'overflow range error'; UNDERFLOW : MathStrErreur := 'underflow range error'; TLOSS : MathStrErreur := 'total loss of precision'; PLOSS : MathStrErreur := 'partial loss of precision'; EDOM : MathStrErreur := 'Unix domain error code'; ERANGE : MathStrErreur := 'Unix range error code'; end; MathStrErreur := ST + MathStrErreur; end; end; Function HexaSt(W:LongInt;n:byte):String; var St : String; Tc : Byte; begin St := ''; Tc:=N; Repeat St := AlphabHexa[W mod BaseHexa] + St; W := W div BaseHexa; Dec(Tc); until Tc=0; St:='0000'+St; HexaSt := Copy(St, length(St)-n+1, n); end; Function StHexa(St:String):LongInt; var Tot : LongInt; i : byte; begin Tot := 0; for i:=1 to length(St) do Tot := Tot*BaseHexa + Pos(St[i],AlphabHexa)-1; StHexa := Tot; end; function MinI(A,B : LongInt):LongInt; begin if A>B then MinI:=B else MinI:=A; end; function MaxI(A,B : LongInt):LongInt; begin if A>B then MaxI:=A else MaxI:=B; end; function MaxR(A,B : Float):Float; begin if A>B then MaxR:=A else MaxR:=B; end; function MinR(A,B : Float):Float; begin if A>B then MinR:=B else MinR:=A; end; function fabs(x:Float):Float; begin if( x < 0.0 ) then fabs := -x else fabs := x ; end; function floor(x:Float):Float; begin Floor:=Trunc(x); end; function ldexp( x: Float; N: Integer):Float; { ldexp() multiplies x by 2**n. } var r : Float; begin R := 1; if N>0 then while N>0 do begin R:=R*2; Dec(N); end else while N<0 do begin R:=R/2; Inc(N); end; ldexp := x * R; end; function frexp(x:Float; var e:Integer ):Float; { frexp() extracts the exponent from x. It returns an integer power of two to expnt and the significand between 0.5 and 1 to y. Thus x = y * 2**expn. } begin e :=0; if (fabs(x)<0.5) then While (fabs(x)<0.5) do begin x := x*2; Dec(e); end else While (fabs(x)>1) do begin x := x/2; Inc(e); end; frexp := x; end; function polevl(var x:Float; var Coef:TabCoef; N:Integer):Float; { Evaluate polynomial SYNOPSIS: int N; double x, y, coef[N+1], polevl[]; y = polevl( x, coef, N ); DESCRIPTION: Evaluates polynomial of degree N: 2 N y = C + C x + C x +...+ C x 0 1 2 N Coefficients are stored in reverse order: coef[0] = C , ..., coef[N] = C . N 0 The function p1evl() assumes that coef[N] = 1.0 and is omitted from the array. Its calling arguments are otherwise the same as polevl(). SPEED: In the interest of speed, there are no checks for out of bounds arithmetic. This routine is used by most of the functions in the library. Depending on available equipment features, the user may wish to rewrite the program in microcode or assembly language. } var ans : Float; i : Integer; begin ans := Coef[0]; for i:=1 to N do ans := ans * x + Coef[i]; polevl:=ans; end; function p1evl(var x:Float; var Coef:TabCoef; N:Integer):Float; { Evaluate polynomial when coefficient of x is 1.0. Otherwise same as polevl. } var ans : Float; i : Integer; begin ans := x + Coef[0]; for i:=1 to N-1 do ans := ans * x + Coef[i]; p1evl := ans; end; function powi(x : Float; nn:Integer):Float; { Real raised to integer power SYNOPSIS: double x, y, powi(); int n; y = powi( x, n ); DESCRIPTION: Returns argument x raised to the nth power. The routine efficiently decomposes n as a sum of powers of two. The desired power is a product of two-to-the-kth powers of x. Thus to compute the 32767 power of x requires 28 multiplications instead of 32767 multiplications. } var n, e, sign, asign, lx : Integer; p : byte; w, s, y : Real; Label done; begin if( x = 0.0 ) then if( nn = 0 ) then begin powi := 1.0 ; exit; end else if( nn < 0 ) then begin powi := MAXNUM ; exit; end else begin powi := 0.0 ; exit; end; if( nn = 0 ) then begin powi := 1.0 ; exit; end; if( x < 0.0 ) then begin asign := -1; x := -x; end else asign := 0; if( nn < 0 ) then begin sign := -1; n := -nn; end else begin sign := 1; n := nn; end; { Overflow detection } { Calculate approximate logarithm of answer } s := frexp( x, lx ); e := (lx - 1)*n; if( (e = 0) or (e > 64) or (e < -64) ) then begin s := (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1); s := (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2; end else s := LOGE2 * e; if( s > MAXLOG ) then begin mtherr( 'powi', OVERFLOW ); y := MAXNUM; goto done; end; if( s < MINLOG ) then begin powi := 0.0; exit; end; { Handle tiny denormal answer, but with less accuracy since roundoff error in 1.0/x will be amplified. The precise demarcation should be the gradual underflow threshold. } if( (s < (-MAXLOG+2.0)) and (sign < 0) ) then begin x := 1.0/x; sign := -sign; end; { First bit of the power } if odd( n ) then y := x else begin y := 1.0; asign := 0; end; w := x; n := n shr 1; while( n<>0 ) do begin w := w * w; { arg to the 2-to-the-kth power } if odd( n ) then { if that bit is set, then include in product } y := y*w; n := n shr 1; end; done: if( asign<>0 ) then y := -y; { odd power of negative number } if( sign < 0 ) then y := 1.0/y; powi := y; end; function reduc(x: Float):Float; { Find a multiple of 1/16 that is within 1/16 of x. } var t:Float; begin t := ldexp( x, 4 ); t := floor( t ); t := ldexp( t, -4 ); reduc := t; end; function pow(x, y : Float):Float; { Power function SYNOPSIS: double x, y, z, pow(); z = pow( x, y ); DESCRIPTION: Computes x raised to the yth power. Analytically, x**y = exp( y log(x) ). Following Cody and Waite, this program uses a lookup table of 2**-i/16 and pseudo extended precision arithmetic to obtain an extra three bits of accuracy in both the logarithm and the exponential. } type Tab17Coef = array[0..16] of Float; const P : TabCoef = ( 4.97778295871696322025E-1, 3.73336776063286838734E0, 7.69994162726912503298E0, 4.66651806774358464979E0, 0, 0, 0); Q : TabCoef = ( 9.33340916416696166113E0, 2.79999886606328401649E1, 3.35994905342304405431E1, 1.39995542032307539578E1, 0, 0, 0); A : Tab17Coef = ( 1.00000000000000000000E0, 9.57603280698573700036E-1, 9.17004043204671215328E-1, 8.78126080186649726755E-1, 8.40896415253714502036E-1, 8.05245165974627141736E-1, 7.71105412703970372057E-1, 7.38413072969749673113E-1, 7.07106781186547572737E-1, 6.77127773468446325644E-1, 6.48419777325504820276E-1, 6.20928906036742001007E-1, 5.94603557501360513449E-1, 5.69394317378345782288E-1, 5.45253866332628844837E-1, 5.22136891213706877402E-1, 5.00000000000000000000E-1); B : Tab17Coef = ( 0.00000000000000000000E0, 1.64155361212281360176E-17, 4.09950501029074826006E-17, 3.97491740484881042808E-17, -4.83364665672645672553E-17, 1.26912513974441574796E-17, 1.99100761573282305549E-17, -1.52339103990623557348E-17, 0.00000000000000000000E0, 0, 0, 0, 0, 0, 0, 0, 0); R : TabCoef = ( 1.49664108433729301083E-5, 1.54010762792771901396E-4, 1.33335476964097721140E-3, 9.61812908476554225149E-3, 5.55041086645832347466E-2, 2.40226506959099779976E-1, 6.93147180559945308821E-1); { log2(e) - 1 } LOG2EA = 0.44269504088896340736; MEXP = 16383.0; MNEXP = -16383.0; { -#define MEXP 16383.0 -#if DENORMAL -#define MNEXP -17183.0 -#else -#define MNEXP -16383.0 -#endif -#endif } var w, z, Wa, Wb, ya, yb, u : Float; F, Fa, Fb, G, Ga, Gb, H, Ha, Hb : Float; e, i : Integer; nflg : Boolean; begin nflg := False; { flag = 1 if x<0 raised to integer power } w := floor(y); if( (w = y) and (fabs(w) < 32768.0) ) then begin i := Trunc(w); w := powi( x, i ); pow := w ; Exit; end; if( x <= 0.0 ) then if( x = 0.0 ) then if( y = 0.0 ) then begin pow := 1.0; { 0**0 } Exit; end else begin pow := 0.0; { 0**y } Exit; end else { x<0.0 } begin if( w <> y ) then begin { noninteger power of negative number } mtherr( 'pow', DOMAIN ); pow := 0.0; Exit; end; nflg := True; x := fabs(x); end; { separate significand from exponent } x := frexp( x, e ); { find significand in antilog table A[] } i := 1; if( x <= A[9] ) then i := 9; if( x <= A[i+4] ) then inc(i, 4); if( x <= A[i+2] ) then inc(i, 2); if( x >= A[1] ) then i := -1; inc(i); { Find (x - A[i])/A[i] in order to compute log(x/A[i]): log(x) = log( a x/a ) = log(a) + log(x/a) log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a } x := x - A[i]; x := x - B[i div 2]; x := x / A[i]; { rational approximation for log(1+v): log(1+v) = v - v**2/2 + v**3 P(v) / Q(v) } z := x*x; w := x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) ); w := w - ldexp( z, -1 ); { w - 0.5 * z } { Convert to base 2 logarithm: multiply by log2(e) } w := w + LOG2EA * w; { Note x was not yet added in to above rational approximation, so do it now, while multiplying by log2(e). } z := w + LOG2EA * x; z := z + x; { Compute exponent term of the base 2 logarithm. } w := -i; w := ldexp( w, -4 ); { divide by 16 } w := w + e; { Now base 2 log of x is w + z. } { Multiply base 2 log by y, in extended precision. separate y into large part ya and small part yb less than 1/16 } ya := reduc(y); yb := y - ya; F := z * y + w * yb; Fa := reduc(F); Fb := F - Fa; G := Fa + w * ya; Ga := reduc(G); Gb := G - Ga; H := Fb + Gb; Ha := reduc(H); w := ldexp( Ga+Ha, 4 ); { Test the power of 2 for overflow } if( w > MEXP ) then begin mtherr( 'pow', OVERFLOW ); Pow := MAXNUM ; Exit; end; if( w < MNEXP ) then begin mtherr( 'pow', UNDERFLOW ); pow := 0.0; Exit; end; e := trunc(w); Hb := H - Ha; if( Hb > 0.0 ) then begin inc(e); Hb := Hb - 0.0625; end; { Now the product y * log2(x) = Hb + e/16.0. Compute base 2 exponential of Hb, where -0.0625 <= Hb <= 0. } z := Hb * polevl( Hb, R, 6 ); { z = 2**Hb - 1 } { Express e/16 as an integer plus a negative number of 16ths. Find lookup table entry for the fractional power of 2. } if( e < 0 ) then i := 0 else i := 1; inc(i, e div 16); e := 16*i - e; w := A[e]; z := w + w * z; { 2**-e * ( 1 + (2**Hb-1) ) } z := ldexp( z, i ); { multiply by integer power of 2 } if( nflg ) then begin { For negative x, find out if the integer exponent is odd or even. } w := ldexp( y, -1 ); w := floor(w); w := ldexp( w, 1 ); if( w <> y ) then z := -z; { odd exponent } end; pow := z ; end; function RoundR(x:Float):Float; { Round double to nearest or even integer valued double SYNOPSIS: double x, y, round(); y = round(x); DESCRIPTION: Returns the nearest integer to x as a double precision floating point result. If x ends in 0.5 exactly, the nearest even integer is chosen. ACCURACY: If x is greater than 1/(2*MACHEP), its closest machine representation is already an integer, so rounding does not change it. } var y, r : Float; label rndup ; begin { Largest integer <= x } y := floor(x); { Fractional part } r := x - y; { Round up to nearest. } if( r > 0.5 ) then goto rndup; { Round to even } if( r = 0.5 ) then begin r := y - 2.0 * floor( 0.5 * y ); if( r = 1.0 ) then rndup: y := y + 1.0; end; { Else round down. } roundR := y; end; function Rad(x:Float):Float; begin Rad:= (x*Pi)/PiRad; end; function Sgn(x:Float):Integer; begin if x<0 then Sgn:=-1 else Sgn:=1; end; function Deg(x:Float):Float; begin x := (x*PiRad)/Pi; While fabs(x)>PiRad do x := x - Sgn(x)*PiRad; Deg := x; end; function DivNonZ(a,b : Float):Float; var Test : Byte; begin Test := 0; if abs(a)1 then Limite := Sgn(N) else Limite := N; end; function Sign(var N:Real):Integer; begin if Abs(N)0 then Sign:=1 else Sign:=-1; end; { PremiŠre m‚thode de calcul des arcsinus } function ASin1(x:Float):Float; begin { ASin:=Pi/2-ACos(x); } if fabs(x)>=1 then x:=Sin(Pi/2)*Sgn(x); Asin1 := ArcTan(x/sqrt(1-x*x)); end; function ACos1(x:Float):Float; begin ACos1:=PIO2-ASin1(x); end; function ASin3(x:Float):Float; forward; function ACos3(x:Float):Float; begin if fabs(x)<0.5 then ACos3:=PIO2 - ASin3(x) else ACos3:=ArcTan(sqrt(1-x*x)/x); end; function ASin3(x:Float):Float; begin if fabs(x)>0.5 then ASin3 := PIO2 - ACos3(x) else ASin3:=ArcTan(x/sqrt(1-x*x)); end; function ASin4(x:Float):Float; { M‚thode de la suite de taylor } var A, C, K, R, R2 : Float; Count : LongInt; begin Count:=0; A:=1; C:=x; K:=1; R := x; Repeat K := K *(A/(A+1)); C := C*x*x; A:=A+2; R2 := R; R := R + (K*C)/A; Inc(Count); Until (Count=80000) or (fabs(R-R2) lossth ) then begin mtherr( 'sin', TLOSS ); sin := 0.0; exit; end; y := floor( x/PIO4 ); { integer part of x/PIO4 } { strip high bits of integer part to prevent integer overflow } z := ldexp( y, -4 ); z := floor(z); { integer part of y/8 } z := y - ldexp( z, 4 ); { y - 16 * (y/16) } j := Trunc(z); { convert to integer for tests on the phase angle } { map zeros to origin } if odd( j ) then begin inc(j); y := y + 1.0; end; j := j and 7; { octant modulo 360 degrees } { reflect in x axis } if( j > 3) then begin sign := -sign; dec(j, 4); end; { Extended precision modular arithmetic } z := ((x - y * DP1) - y * DP2) - y * DP3; zz := z * z; if( (j=1) or (j=2) ) then y := 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 ) else { y = z + z * (zz * polevl( zz, sincof, 5 )); } y := z + z * z * z * polevl( zz, sincof, 5 ); if(sign < 0) then y := -y; sin := y; end; function Cos(x:Float):Float; { Circular cosine SYNOPSIS: double x, y, cos(); y = cos( x ); DESCRIPTION: Range reduction is into intervals of pi/4. The reduction error is nearly eliminated by contriving an extended precision modular arithmetic. Two polynomial approximating functions are employed. Between 0 and pi/4 the cosine is approximated by 1 - x**2 Q(x**2). Between pi/4 and pi/2 the sine is represented as x + x**3 P(x**2). } var y, z, zz : Float; j, sign : Integer; i : LongInt; begin { make argument positive } sign := 1; if( x < 0 ) then x := -x; if( x > lossth ) then begin mtherr( 'cos', TLOSS ); cos := 0.0; exit; end; y := floor( x/PIO4 ); z := ldexp( y, -4 ); z := floor(z); { integer part of y/8 } z := y - ldexp( z, 4 ); { y - 16 * (y/16) } { integer and fractional part modulo one octant } i := Trunc(z); if odd( i ) then { map zeros to origin } begin inc(i); y := y + 1.0; end; j := i and 07; if( j > 3) then begin dec(j,4); sign := -sign; end; if( j > 1 ) then sign := -sign; { Extended precision modular arithmetic } z := ((x - y * DP1) - y * DP2) - y * DP3; zz := z * z; if( (j=1) or (j=2) ) then { y = z + z * (zz * polevl( zz, sincof, 5 )); } y := z + z * z * z * polevl( zz, sincof, 5 ) else y := 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 ); if(sign < 0) then y := -y; cos := y ; end; function ASin(x:Float):Float; const P : TabCoef = ( -6.96822599948686174217E-1, 1.01598286089872099722E1, -3.97340771391578757294E1, 5.72912144709846496134E1, -2.74148200465925708020E1, 0,0); Q : TabCoef = ( -2.38368245005177488242E1, 1.51095072703128995631E2, -3.82340216045978957023E2, 4.17767300951716199422E2, -1.64488920279555473283E2,0,0); var Sign : Integer; a,z,zz,pp : Float; Flag : byte; begin if x>0 then begin Sign:=1; a:=x; end else begin Sign:=-1; a:=-x; end; if (a>1.0) then begin mtherr( 'ASin', OVERFLOW ); ASin :=0.0; end else begin if (a0.5) then begin zz:=(0.5-a+0.5)/2.0; z:=sqrt(zz); flag:=1; end else begin z:=a; zz:=z*z; flag:=0; end; z := z * zz * polevl( zz, P, 4)/p1evl( zz, Q, 5) + z; if( flag<>0 ) then begin z := z + z; z := PIO2 - z; end; end; if (sign < 0 ) then z := -z; ASin:=z; end; end; function ACos(x:Float):Float; { Inverse circular cosine SYNOPSIS: double x, y, acos(); y = acos( x ); DESCRIPTION: Returns radian angle between -pi/2 and +pi/2 whose cosine is x. Analytically, acos(x) = pi/2 - asin(x). However if |x| is near 1, there is cancellation error in subtracting asin(x) from pi/2. Hence if x < -0.5, acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) ); or if x > +0.5, acos(x) = 2.0 * asin( sqrt((1-x)/2) ). } label domerr; begin if fabs(x)>1 then begin mtherr( 'Acos', OVERFLOW ); ACos := x; end else if (x<-0.5) then ACos := pi - 2.0 * Asin( sqrt((1+x)*0.5) ) else if (x > 0.5) then ACos := 2.0 * Asin( sqrt((1-x)*0.5) ) else ACos := PIO2 - Asin(x) end; function Exp(x:Float):Float; { Exponential function SYNOPSIS: double x, y, exp(); y = exp( x ); DESCRIPTION: Returns e (2.71828...) raised to the x power. Range reduction is accomplished by separating the argument into an integer k and fraction f such that x k f e = 2 e. A Pade' form of degree 2/3 is used to approximate exp(f) - 1 in the basic range [-0.5 ln 2, 0.5 ln 2]. } const P : TabCoef = ( 1.26183092834458542160E-4, 3.02996887658430129200E-2, 1.00000000000000000000E0, 0, 0, 0, 0); Q : TabCoef = ( 3.00227947279887615146E-6, 2.52453653553222894311E-3, 2.27266044198352679519E-1, 2.00000000000000000005E0, 0 ,0 ,0); C1 = 6.9335937500000000000E-1; C2 = 2.1219444005469058277E-4; var n : Integer; px, qx, xx : Float; begin if( x > MAXLOG) then begin mtherr( 'Exp', OVERFLOW ); Exp := MAXNUM ; end else if( x < MINLOG ) then begin mtherr('Exp', UNDERFLOW ); Exp := 0.0; end else begin { Express e**x = e**g 2**n } { = e**g e**( n loge(2) ) } { = e**( g + n loge(2) ) } px := x * LOG2E; qx := floor( px + 0.5 ); { floor() truncates toward -infinity. } n := Trunc(qx); x := x - qx * C1; x := x + qx * C2; { rational approximation for exponential } { of the fractional part: } { e**x - 1 = 2x P(x**2)/( Q(x**2) - P(x**2) ) } xx := x * x; px := x * polevl( xx, P, 2 ); x := px/( polevl( xx, Q, 3 ) - px ); x := ldexp( x, 1 ); x := x + 1.0; x := ldexp( x, n ); Exp := x; end; end; function Log(x:Float):Float; { Natural logarithm SYNOPSIS: double x, y, log(); y = log( x ); DESCRIPTION: Returns the base e (2.718...) logarithm of x. The argument is separated into its exponent and fractional parts. If the exponent is between -1 and +1, the logarithm of the fraction is approximated by log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x). Otherwise, setting z = 2(x-1)/x+1), log(x) = z + z**3 P(z)/Q(z). } const P : TabCoef = ( { Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) 1/sqrt(2) <= x < sqrt(2) } 4.58482948458143443514E-5, 4.98531067254050724270E-1, 6.56312093769992875930E0, 2.97877425097986925891E1, 6.06127134467767258030E1, 5.67349287391754285487E1, 1.98892446572874072159E1); Q : TabCoef = ( 1.50314182634250003249E1, 8.27410449222435217021E1, 2.20664384982121929218E2, 3.07254189979530058263E2, 2.14955586696422947765E2, 5.96677339718622216300E1, 0); { Coefficients for log(x) = z + z**3 P(z)/Q(z), where z = 2(x-1)/(x+1) 1/sqrt(2) <= x < sqrt(2) } R : TabCoef = ( -7.89580278884799154124E-1, 1.63866645699558079767E1, -6.41409952958715622951E1, 0, 0, 0, 0); S : TabCoef = ( -3.56722798256324312549E1, 3.12093766372244180303E2, -7.69691943550460008604E2, 0, 0, 0, 0); var e : Integer; z, y : Float; Label Ldone; { All four routines return a double precision floating point result. floor() returns the largest integer less than or equal to x. It truncates toward minus infinity. ceil() returns the smallest integer greater than or equal to x. It truncates toward plus infinity. frexp() extracts the exponent from x. It returns an integer power of two to expnt and the significand between 0.5 and 1 to y. Thus x = y * 2**expn. ldexp() multiplies x by 2**n. These functions are part of the standard C run time library for many but not all C compilers. The ones supplied are written in C for either DEC or IEEE arithmetic. They should be used only if your compiler library does not already have them. The IEEE versions assume that denormal numbers are implemented in the arithmetic. Some modifications will be required if the arithmetic has abrupt rather than gradual underflow. } begin if( x <= 0.0 ) then begin if( x = 0.0 ) then mtherr( 'Log', SING ) else mtherr( 'Log', DOMAIN ); Log := MINLOG; end; x := frexp( x, e ); { logarithm using log(x) = z + z**3 P(z)/Q(z), where z = 2(x-1)/x+1) } if( (e > 2) or (e < -2) ) then begin if( x < SQRTH ) then begin { 2( 2x-1 )/( 2x+1 ) } Dec(e, 1); z := x - 0.5; y := 0.5 * z + 0.5; end else begin { 2 (x-1)/(x+1) } z := x - 0.5; z := z - 0.5; y := 0.5 * x + 0.5; end; x := z / y; { /* rational form */ } z := x*x; z := x + x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) ); goto ldone; end; { logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) } if( x < SQRTH ) then begin Dec(e, 1); x := ldexp( x, 1 ) - 1.0; { 2x - 1 } end else x := x - 1.0; { rational form } z := x*x; y := x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) ); y := y - ldexp( z, -1 ); { y - 0.5 * z } z := x + y; ldone: { recombine with exponent term } if( e <> 0 ) then begin y := e; z := z - y * 2.121944400546905827679e-4; z := z + y * 0.693359375; end; Log:= z; end; function Log10(x:Float):Float; { Common logarithm SYNOPSIS: double x, y, log10(); y = log10( x ); DESCRIPTION: Returns logarithm to the base 10 of x. The argument is separated into its exponent and fractional parts. The logarithm of the fraction is approximated by log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x). } { Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x) } { 1/sqrt(2) <= x < sqrt(2) } const P : TabCoef = ( 4.58482948458143443514E-5, 4.98531067254050724270E-1, 6.56312093769992875930E0, 2.97877425097986925891E1, 6.06127134467767258030E1, 5.67349287391754285487E1, 1.98892446572874072159E1); Q : TabCoef = ( 1.50314182634250003249E1, 8.27410449222435217021E1, 2.20664384982121929218E2, 3.07254189979530058263E2, 2.14955586696422947765E2, 5.96677339718622216300E1, 0); SQRTH = 0.70710678118654752440; L102A = 3.0078125E-1; L102B = 2.48745663981195213739E-4; L10EA = 4.3359375E-1; L10EB = 7.00731903251827651129E-4; var e : Integer; w, y, z : Float; begin if( x <= 0.0 ) then begin if( x = 0.0 ) then mtherr( 'Log10', SING ) else mtherr( 'Log10', DOMAIN ); Log10 := MINLOG; end else begin { separate mantissa from exponent } x := frexp( x, e ); { logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) } if( x < SQRTH ) then begin dec(e); x := ldexp( x, 1 ) - 1.0; { 2x - 1 } end else x := x - 1.0; { rational form } z := x*x; y := x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) ); y := y - ldexp( z, -1 ); { y - 0.5 * x**2 } { multiply log of fraction by log10(e) } { and base 2 exponent by log10(2) } z := (x + y) * L10EB; { accumulate terms in order of size } z := z + y * L10EA; z := z + x * L10EA; w := e; z := z + w * L102B; z := z + w * L102A; Log10 := z ; end; end; function Exp10(x:Float):Float; { Base 10 exponential function (Common antilogarithm) SYNOPSIS: double x, y, exp10(); y = exp10( x ); DESCRIPTION: Returns 10 raised to the x power. Range reduction is accomplished by expressing the argument as 10**x = 2**n 10**f, with |f| < 0.5 log10(2). The Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) is used to approximate 10**f. } const P : TabCoef = ( 4.09962519798587023075E-2, 1.17452732554344059015E1, 4.06717289936872725516E2, 2.39423741207388267439E3, 0, 0 ,0); Q : TabCoef = ( 8.50936160849306532625E1, 1.27209271178345121210E3, 2.07960819286001865907E3, 0, 0, 0, 0); LOG210 = 3.32192809488736234787e0; LG102A = 3.01025390625000000000E-1; LG102B = 4.60503898119521373889E-6; MAXL10 = 38.230809449325611792; var px, xx : Float; n : Integer; begin if( x > MAXL10 ) then begin mtherr( 'exp10', OVERFLOW ); Exp10 := MAXNUM; end else if( x < -MAXL10 ) then { Would like to use MINLOG but can't } begin mtherr( 'exp10', UNDERFLOW ); Exp10 := 0.0; end else begin { Express 10**x = 10**g 2**n } { = 10**g 10**( n log10(2) ) } { = 10**( g + n log10(2) ) } px := floor( LOG210 * x + 0.5 ); n := Trunc(px); x := x - px * LG102A; x := x - px * LG102B; { rational approximation for exponential } { of the fractional part: } { 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) ) } xx := x * x; px := x * polevl( xx, P, 3 ); x := px/( p1evl( xx, Q, 3 ) - px ); x := 1.0 + ldexp( x, 1 ); { multiply by power of 2 } x := ldexp( x, n ); Exp10 := x; end end; function Cosh(x:Float):Float; { Hyperbolic cosine SYNOPSIS: double x, y, cosh(); y = cosh( x ); DESCRIPTION: Returns hyperbolic cosine of argument in the range MINLOG to MAXLOG. cosh(x) = ( exp(x) + exp(-x) )/2. } var y : Float; begin if( x < 0 ) then x := -x; if( x > MAXLOG ) then begin mtherr( 'cosh', OVERFLOW ); Cosh := MAXNUM ; end else begin y := exp(x); y := y + 1.0/y; y := ldexp( y, -1 ); Cosh := y ; end; end; function ACosh(x:Float):Float; { Inverse hyperbolic cosine SYNOPSIS: double x, y, acosh(); y = acosh( x ); Returns inverse hyperbolic cosine of argument. If 1 <= x < 1.5, a rational approximation sqrt(z) * P(z)/Q(z) where z = x-1, is used. Otherwise, acosh(x) = log( x + sqrt( (x-1)(x+1) ). } const P : TabCoef = ( 1.18801130533544501356E2, 3.94726656571334401102E3, 3.43989375926195455866E4, 1.08102874834699867335E5, 1.10855947270161294369E5, 0, 0); Q : TabCoef = ( 1.86145380837903397292E2, 4.15352677227719831579E3, 2.97683430363289370382E4, 8.29725251988426222434E4, 7.83869920495893927727E4, 0, 0); var z,a : Float; begin { acosh(z) = sqrt(x) * R(x), z = x + 1, interval 0 < x < 0.5 } if( x < 1.0 ) then begin mtherr( 'acosh', DOMAIN ); ACosh :=0.0; end else if( x > 1.0e8 ) then ACosh :=( log(x) + LOGE2 ) else begin z := x - 1.0; if( z < 0.5 ) then begin a := sqrt(z) * (polevl(z, P, 4) / p1evl(z, Q, 5) ); ACosh := a; end else begin a := sqrt( z*(x+1.0) ); ACosh :=( log(x + a) ); end; end; end; function Sinh(x:Float):Float; { Hyperbolic sine SYNOPSIS: double x, y, sinh(); y = sinh( x ); DESCRIPTION: Returns hyperbolic sine of argument in the range MINLOG to MAXLOG. The range is partitioned into two segments. If |x| <= 1, a rational function of the form x + x**3 P(x)/Q(x) is employed. Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2. } const P : TabCoef = ( -7.89474443963537015605E-1, -1.63725857525983828727E2, -1.15614435765005216044E4, -3.51754964808151394800E5, 0, 0, 0); Q : TabCoef = ( -2.77711081420602794433E2, 3.61578279834431989373E4, -2.11052978884890840399E6, 0, 0, 0, 0); var a : Float; begin a := fabs(x); if( (a > MAXLOG) or (a > -MINLOG) ) then begin mtherr( 'sinh', DOMAIN ); if( x > 0 ) then Sinh := MAXNUM else Sinh := -MAXNUM; end else if( a > 1.0 ) then begin a := exp(a); a := 0.5*a - (0.5/a); if( x < 0 ) then a := -a; Sinh := a; end else begin a := a*a; Sinh := x + x * a * (polevl(a,P,3)/p1evl(a,Q,3)) ; end; end; function ASinh(x:Float):Float; { Inverse hyperbolic sine SYNOPSIS: double x, y, asinh(); y = asinh( x ); DESCRIPTION: Returns inverse hyperbolic sine of argument. If |x| < 0.5, the function is approximated by a rational form x + x**3 P(x)/Q(x). Otherwise, asinh(x) = log( x + sqrt(1 + x*x) ). } const P : TabCoef = ( -4.33231683752342103572E-3, -5.91750212056387121207E-1, -4.37390226194356683570E0, -9.09030533308377316566E0, -5.56682227230859640450E0, 0, 0); Q : TabCoef = ( 1.28757002067426453537E1, 4.86042483805291788324E1, 6.95722521337257608734E1, 3.34009336338516356383E1, 0, 0, 0); var Sign : Integer; a,z,zz,pp : Float; Flag : byte; begin if( x < 0.0 ) then begin sign := -1; x := -x end else sign := 1; if( x > 1.0e8 ) then ASinh := sign * (log(x) + LOGE2) else begin z := x * x; if( x < 0.5 ) then begin a := ( polevl(z, P, 4)/p1evl(z, Q, 4) ) * z; a := a * x + x; if( sign < 0 ) then a := -a; ASinh := a; end else begin a := sqrt( z + 1.0 ); ASinh :=( sign * log(x + a) ); end; end; end; function TanCot(xx:Float; cotflg:Boolean):Float; const P : TabCoef = ( -1.30936939181383777646E4, 1.15351664838587416140E6, -1.79565251976484877988E7, 0, 0, 0, 0); Q : TabCoef = ( 1.36812963470692954678E4, -1.32089234440210967447E6, 2.50083801823357915839E7, -5.38695755929454629881E7, 0, 0, 0); DP1 = 7.853981554508209228515625E-1; DP2 = 7.94662735614792836714E-9; DP3 = 3.06161699786838294307E-17; var x,y,z,zz : Float; Sign,j : Integer; begin { make argument positive but save the sign } if( xx < 0 ) then begin x := -xx; sign := -1; end else begin x := xx; sign := 1; end; if( x > lossth ) then begin if( cotflg ) then mtherr( 'cot', TLOSS ) else mtherr( 'tan', TLOSS ); TanCot := 0.0; end; { compute x mod PIO4 } y := floor( x/PIO4 ); { strip high bits of integer part } z := ldexp( y, -3 ); z := floor(z); { integer part of y/8 } z := y - ldexp( z, 3 ); { y - 16 * (y/16) } { integer and fractional part modulo one octant } j := Trunc(z); { map zeros and singularities to origin } if odd(j) then begin inc(j); y := y+1.0; end; z := ((x - y * DP1) - y * DP2) - y * DP3; zz := z * z; if( zz > 1.0e-14 ) then y := z + z * (zz * polevl( zz, P, 2 )/p1evl(zz, Q, 4)) else y := z; if( j and 2 )>0 then begin if( cotflg ) then y := -y else y := -1.0/y; end else if( cotflg ) then y := 1.0/y; if( sign < 0 ) then y := -y; Tancot := y ; end; function Cot(x:Float):Float; { Circular cotangent SYNOPSIS: double x, y, cot(); y = cot( x ); DESCRIPTION: Returns the circular cotangent of the radian argument x. Range reduction is modulo pi/4. A rational function x + x**3 P(x**2)/Q(x**2) is employed in the basic interval [0, pi/4]. } begin if( x = 0.0 ) then begin mtherr( 'cot', SING ); Cot := MAXNUM; end else Cot := tancot(x, True); end; function Tan(x:Float):Float; { Circular tangent SYNOPSIS: double x, y, tan(); y = tan( x ); DESCRIPTION: Returns the circular tangent of the radian argument x. Range reduction is modulo pi/4. A rational function x + x**3 P(x**2)/Q(x**2) is employed in the basic interval [0, pi/4]. } begin Tan := tancot(x, False); end; function Atan(x:Float):Float; { Inverse circular tangent (arctangent) SYNOPSIS: double x, y, atan(); y = atan( x ); DESCRIPTION: Returns radian angle between -pi/2 and +pi/2 whose tangent is x. Range reduction is from four intervals into the interval from zero to tan( pi/8 ). The approximant uses a rational function of degree 3/4 of the form x + x**3 P(x)/Q(x). } const P : TabCoef = ( -8.40980878064499716001E-1, -8.83860837023772394279E0, -2.18476213081316705724E1, -1.48307050340438946993E1, 0, 0, 0); Q : TabCoef = ( 1.54974124675307267552E1, 6.27906555762653017263E1, 9.22381329856214406485E1, 4.44921151021319438465E1, 0, 0, 0); { tan( 3*pi/8 ) } T3P8 = 2.41421356237309504880; { tan( pi/8 ) } TP8 = 0.41421356237309504880; var y,z : Float; Sign : Integer; begin { make argument positive and save the sign } sign := 1; if( x < 0.0 ) then begin sign := -1; x := -x; end; { range reduction } if( x > T3P8 ) then begin y := PIO2; x := -( 1.0/x ); end else if( x > TP8 ) then begin y := PIO4; x := (x-1.0)/(x+1.0); end else y := 0.0; { rational form in x**2 } z := x * x; y := y + ( polevl( z, P, 3 ) / p1evl( z, Q, 4 ) ) * z * x + x; if( sign < 0 ) then y := -y; Atan := y; end; function Atan2( y, x: Float):Float; { Quadrant correct inverse circular tangent SYNOPSIS: double x, y, z, atan2(); z = atan2( y, x ); DESCRIPTION: Returns radian angle whose tangent is y/x. Define compile time symbol ANSIC = 1 for ANSI standard, range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range 0 to 2PI, args (x,y). } var z, w : Float; code : Integer; begin code := 0; if ( x < 0.0 ) then code := 2; if ( y < 0.0 ) then code := code or 1; if ( x = 0.0 ) then begin if ( code and 1 )>0 then Atan2 := -PIO2 else if ( y = 0.0 ) then Atan2 := 0.0 else Atan2 := PIO2; Exit; end; if ( y = 0.0 ) then begin if ( code and 2 )>0 then Atan2 := PI else Atan2 := 0.0; Exit; end; case code of 0, 1: w := 0.0; 2: w := PI; 3: w := -PI; end; Atan2 := w + Arctan( y/x ) ; { Utilise le Arctan r_sident, prend moins de place } end; function Tanh(x:Float):Float; { Hyperbolic tangent SYNOPSIS: double x, y, tanh(); y = tanh( x ); DESCRIPTION: Returns hyperbolic tangent of argument in the range MINLOG to MAXLOG. A rational function is used for |x| < 0.625. The form x + x**3 P(x)/Q(x) of Cody _& Waite is employed. Otherwise, tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1). } const P : TabCoef = ( -9.64399179425052238628E-1, -9.92877231001918586564E1, -1.61468768441708447952E3, 0 ,0, 0, 0); Q : TabCoef = ( 1.12811678491632931402E2, 2.23548839060100448583E3, 4.84406305325125486048E3, 0 ,0, 0, 0); var z, s : Float; begin z := fabs(x); if( z > 0.5 * MAXLOG ) then begin if( x > 0 ) then Tanh := 1.0 else Tanh := -1.0 ; end else if( z >= 0.625 ) then begin s := exp(2.0*z); z := 1.0 - 2.0/(s + 1.0); if( x < 0 ) then z := -z; end else begin s := x * x; z := polevl( s, P, 2 )/p1evl(s, Q, 3); z := x * s * z; z := x + z; end; Tanh := z ; end; function ATanh(x:Float):Float; { Inverse hyperbolic tangent SYNOPSIS: double x, y, atanh(); y = atanh( x ); DESCRIPTION: Returns inverse hyperbolic tangent of argument in the range MINLOG to MAXLOG. If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is employed. Otherwise, atanh(x) = 0.5 * log( (1+x)/(1-x) ). } const P : TabCoef = ( -8.54074331929669305196E-1, 1.20426861384072379242E1, -4.61252884198732692637E1, 6.54566728676544377376E1, -3.09092539379866942570E1, 0, 0); Q : TabCoef = ( -1.95638849376911654834E1, 1.08938092147140262656E2, -2.49839401325893582852E2, 2.52006675691344555838E2, -9.27277618139601130017E1, 0, 0); var z, s : Float; begin z := fabs(x); if( z >= 1.0 ) then begin if( x = 1.0 ) then Atanh := MAXNUM else if( x = -1.0 ) then Atanh := -MAXNUM else { |x| > 1.0 } begin mtherr( 'atanh', DOMAIN ); Atanh := MAXNUM ; end end else if( z < 1.0e-7 ) then Atanh := x else if( z < 0.5 ) then begin z := x * x; s := x + x * z * (polevl(z, P, 4) / p1evl(z, Q, 5)); Atanh := s; end else Atanh := 0.5 * log((1.0+x)/(1.0-x)) ; end; function sqrt(x:Float):Float; { Square root SYNOPSIS: double x, y, sqrt(); y = sqrt( x ); DESCRIPTION: Returns the square root of x. Range reduction involves isolating the power of two of the argument and using a polynomial approximation to obtain a rough value for the square root. Then Heron's iteration is used three times to converge to an accurate value. } var e : Integer; w,z : Float; begin if( x <= 0.0 ) then begin if( x < 0.0 ) then mtherr( 'sqrt', DOMAIN ); sqrt := 0.0; end else begin w := x; { separate exponent and significand } z := frexp( x, e ); { approximate square root of number between 0.5 and 1 } { relative error of approximation = 7.47e-3 } x := 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z; { adjust for odd powers of 2 } if odd(e) then x := x*SQRT2; { re-insert exponent } x := ldexp( x, (e div 2) ); { Newton iterations: } x := 0.5*(x + w/x); x := 0.5*(x + w/x); x := 0.5*(x + w/x); x := 0.5*(x + w/x); x := 0.5*(x + w/x); x := 0.5*(x + w/x); sqrt := x; end; end; function cbrt(x:Float):Float; { Cube root SYNOPSIS: double x, y, cbrt(); y = cbrt( x ); DESCRIPTION: Returns the cube root of the argument, which may be negative. Range reduction involves determining the power of 2 of the argument. A polynomial of degree 2 applied to the mantissa, and multiplication by the cube root of 1, 2, or 4 approximates the root to within about 0.1%. Then Newton's iteration is used three times to converge to an accurate result. } const CBRT2 = 1.25992104989487316477; CBRT4 = 1.58740105196819947475; var e, rem, Sign : Integer; z : Float; begin if( x = 0 ) then cbrt := 0.0 else if( x > 0 ) then sign := 1 else begin sign := -1; x := -x; end; z := x; { extract power of 2, leaving } { mantissa between 0.5 and 1 } x := frexp( x, e ); { Approximate cube root of number between .5 and 1, } { peak relative error = 6.36e-4 } x := (-1.9150215751434832257e-1 * x + 6.9757045195246484393e-1) * x + 4.9329566506409572486e-1; { exponent divided by 3 } if( e >= 0 ) then begin rem := e; e := e div 3; rem := rem - 3*e; if( rem = 1 ) then x := x*CBRT2 else if( rem = 2 ) then x := x*CBRT4; end { argument less than 1 } else begin e := -e; rem := e; e := e div 3; rem := rem - 3*e; if( rem = 1 ) then x := x/CBRT2 else if( rem = 2 ) then x := x/CBRT4; e := -e; end; { multiply by power of 2 } x := ldexp( x, e ); { Newton iteration } x := x - ( x - (z/(x*x)) )/3.0; x := x - ( x - (z/(x*x)) )/3.0; x := x - ( x - (z/(x*x)) )/3.0; if( sign < 0 ) then x := -x; cbrt := x; end; { Imp_ratif pour _liminer les divisions par 0 ! } function NonZ(x:Float):Float; begin if Abs(x)0.01 then begin V[2] := PI-V[2]; V[1] := PI-V[1]; end; V[3] := Acos(Limite(DivNonZ(Mat[1,1], Cos(V[2])))); if Abs(Cos(V[1])*Sin(V[3])+Sin(V[1])*Mat[1,3]*Cos(V[3])-Mat[2,1])>0.01 then begin V[1] := V[1]-PI; V[2] := PI-V[2]; V[3] := PI-V[3]; end; { !!! Attention ceci marche tr_s bien pour Z !!! } if (Abs(Mat[2,3])<0.001) then { Z est sur le plan horizontal } if (Abs(Mat[3,3])<0.001) then { Z est sur l'axe horizontal } begin V[3] := Atan2(Mat[2,1], Mat[2,2]); { Bon pour horizontal } V[1] := 0; end; OptimiseRotation(V); { VerifAngle(V, Mat); } (* if MathError>0 then Averti(MathStrErreur); *) end else { Rep_re observateur = Latitude, Longitude, Rotation propre } begin { premi_re rotation autour de X0 } Mat1:=Mat; TransposeMat(DimenMat, Mat1); { pour avoir les axes de l'observateur } V[1] := Latitude(Mat1[3]); V[2] := Longitude(Mat1[3]); V[3] := 0; Initial(DimenMat, Mat1); RotationMatrice(V, Mat1, RObservateur); { Construction d'une base observateur } MultMat(DimenMat, Mat, Mat1, Mat1); V[3] := ACos(Limite(Mat1[2,2]))*Sgn(Mat1[1,2]); end; end; function Egal(A, B:Float):boolean; begin if abs(A-B)<=Epsilon then Egal := True else Egal := False; end; procedure OptimiseRotation(var Vn:TNVector); const Epsi = PI-100000*Epsilon2; var No, i, j, k:Byte; Max : Float; begin for No:=1 to DimenMat do if abs(Vn[No])>=PI2 then Vn[No]:=Vn[No]-PI2*Sgn(Vn[No]); for No:=1 to DimenMat do begin for i:=1 to DimenMat do if i<>No then begin k:=j; j:=i; end; if (abs(Vn[k])>=Epsi) and (Abs(Vn[j])>=Epsi) then for i:=1 to DimenMat do Vn[i] := Vn[i]-PI*Sgn(Vn[i]); end; end; function VerifAngle(var V:TNVector; var Mat:TNMatrix):Boolean; var C,S : TNVector; i : Byte; A, B, D: Real; Erreur : Boolean; begin for i:=1 to 3 do begin C[i] := Cos(V[i]); S[i] := Sin(V[i]); end; { Etape 4 } { x""= x(cos(2)cos(3)) + y(cos(1)sin(3)-sin(1)sin(2)cos(3)) + z(cos(1)sin(2)cos(3) + sin(1)sin(3)) y""= x(-cos(2)sin(3))+ y(cos(1)cos(3)-sin(1)sin(2)sin(3)) + z(cos(1)sin(2)sin(3) + sin(1)cos(3)) z""= x(-sin(2)) + y(-sin(1)cos(2)) + zcos(1)cos(2) } A := C[2]*C[3] - Mat[1,1] + C[1]*S[3]-S[1]*S[2]*C[3] - Mat[2,1] - S[1]*C[2] - Mat[2,3] { + C[1]*S[2]*C[3]-S[1]*S[3] - Mat[3,1] } + C[1]*S[2]*S[3]-S[1]*C[3] - Mat[3,2] { - S[2] - Mat[1,3] } { + C[1]*C[3]-S[1]*S[2]*S[3] - Mat[2,2] } + C[1]*C[2] - Mat[3,3]; Erreur := (abs(A)<0.01) ; VerifAngle := Erreur; end; procedure Ajuste(var V:TNVector; var Mat:TNMatrix); var W : TNVector; i : Byte; StTest : array[0..80] of char; begin i := 0; Repeat Case (i and 3) of 0 : W[1] := V[1]; 1 : W[1] := -V[1]; 2 : W[1] := PI-V[1]; 3 : W[1] := -PI-V[1]; end; Case ((i and 12) shr 2) of 0 : W[2] := V[2]; 1 : W[2] := -V[2]; 2 : W[2] := PI-V[2]; 3 : W[2] := -PI-V[2]; end; Case ((i and 48) shr 4) of 0 : W[3] := V[3]; 1 : W[3] := -V[3]; 2 : W[3] := PI-V[3]; 3 : W[3] := -PI-V[3]; end; Inc(i); Until (i=48) or VerifAngle(W, Mat); if i=48 then begin MessageBeep(0); { Erreur } { W := V; } end; { Str(i, StTest); TestTexte(StTest); } { V := W; } end; {----------------------------------------------------------} { } { Routines de calcul matriciel } { } {----------------------------------------------------------} procedure ClearVect(var Data : TNVector); begin FillChar(Data, SizeOf(Data), 0); end; procedure AddVecteur(Dimen:Integer; var A, B, C : TNVector); var i:Integer; begin for i:=1 to Dimen do C[i] := A[i]+B[i]; end; procedure SubVecteur(Dimen:Integer; var A, B, C : TNVector); var i:Integer; begin for i:=1 to Dimen do C[i] := A[i]-B[i]; end; procedure ClearMat(var Data : TNmatrix); begin FillChar(Data, SizeOf(Data), 0); end; procedure Initial(Dimen : integer; var Data : TNmatrix); {----------------------------------------------------------} {- Output: Data -} {- -} {- This procedure intializes the above variables to zero. -} {----------------------------------------------------------} var i:Byte; begin ClearMat(Data); for i:=1 to Dimen do Data[i,i]:=1; end; { procedure Initial } procedure MultMat(Dimen : Integer; MatA, MatB:TNMatrix; var MatC:TNMatrix); var i,j: Integer; begin for i:=1 to Dimen do for j:=1 to Dimen do MatC[i,j] := MatA[i,1]*MatB[1,j] + MatA[i,2]*MatB[2,j] + MatA[i,3]*MatB[3,j]; end; procedure MultMatVect(Dimen : Integer; var Mat:TNMatrix; var V:TNVector); var i,j: Integer; V2 : TNVector; begin ClearVect(V2); for i:=1 to Dimen do for j:=1 to Dimen do V2[i] := V2[i] + Mat[i,j]*V[j]; V:=V2; end; procedure DivMatFloat(Dimen : integer; var Data : TNmatrix; Divisor : Float); var i,j: Integer; begin if abs(Divisor)P then begin K:=J; J:=i; end; Mat[J,J]:= CosT; Mat[j,K]:= SinT; Mat[K,J]:= -SinT; Mat[K,K]:= CosT; end; procedure InitialiseXmat; var i:Integer; begin for i:=1 to 3 do AffecteMat(i, Pi/PartPi, XMat[i]); for i:=4 to 6 do AffecteMat(i-3, -Pi/PartPi, XMat[i]); end; procedure RotationMatrice(var PhiRot:TNvector; var MatRot:TNMatrix; Repere:Byte); { Angles en radian } var MatT : TNMatrix; V1,V2,V3,V4,V5,V6,V7,V8, V0, Phi, S0, C0, Rx, Ry, Rz : Real; i : Integer; begin Case Repere of RAnatomique : for i:=1 to DimenMat do { Rotations X, Y, Z } if PhiRot[i]<>0.0 then begin {!!! Rep_re Objet = Post-Multiplier !!!} AffecteMat(i, PhiRot[i], MatT); MultMat(DimenMat, MatRot, MatT, MatRot); end; Robservateur : for i:=DimenMat downto 1 do { Rotations X, Y, Z } if PhiRot[i]<>0.0 then { Evite les effacements de matrice } begin {!!! Rep_re Observateur = Pr_-Multiplier !!!} AffecteMat(i, PhiRot[i], MatT); MultMat(DimenMat, MatT, MatRot, MatRot); end; RAxesPropres : begin { Pour la description de cette m_thode, consulter : "Rotation autour d'un vecteur arbitraire" ROBOTICS, Control, Sensing, Vision and Intelligence Edition McGraw-Hill Book Company (1987) page 20-22 } Rx:= PhiRot[1]; Ry:= PhiRot[2]; Rz:= PhiRot[3]; Phi := PhiRot[4]; C0:=Cos(Phi); S0:=Sin(Phi); V0:=1-C0; { M_thode l_g_re, peu rapide : 24 multiplications} { MatT[1,1] := Rx*Rx*V0+C0; MatT[1,2] := Rx*Ry*V0-Rz*S0; MatT[1,3] := Rx*Rz*V0+Ry*S0; MatT[2,1] := Rx*Ry*V0+Rz*S0; MatT[2,2] := Ry*Ry*V0+C0; MatT[2,3] := Ry*Rz*V0-Rx*S0; MatT[3,1] := Rx*Rz*V0-Ry*S0; MatT[3,2] := Ry*Rz*V0+Rx*S0; MatT[3,3] := Rz*Rz*V0+C0; } { M_thode plus rapide, moins l_g_re : 12 multiplications} V1 := Rx*V0; V2 := Ry*V1; V3 := Rz*V1; V4 := Rx*S0; V5 := Ry*S0; V6 := Rz*S0; V7 := Ry*V0; V8 := Rz*V7; MatT[1,1] := Rx*V1+C0; MatT[1,2] := V2-V6; MatT[1,3] := V3+V5; MatT[2,1] := V2+V6; MatT[2,2] := Ry*V7+C0; MatT[2,3] := V8-V4; MatT[3,1] := V3-V5; MatT[3,2] := V8+V4; MatT[3,3] := Rz*Rz*V0+C0; MultMat(DimenMat, MatRot, MatT, MatRot); { Dans ses axes propres } end; end; end; procedure IncrementeMatrice(NoAxe, Sens : Integer; var MatRot, MatVP:TNMatrix; Rotation:Boolean); begin if Sens=Moins then { Change de matrice } Inc(NoAxe,3); if RotationPropre then case RepereActuel of RAnatomique : {!!! Rep_re Objet = Post-Multiplier !!!} MultMat(DimenMat, MatRot, XMat[NoAxe], MatRot); RAxesPropres : begin if NoAxe>3 then Dec(NoAxe,3); if NoAxe=2 then Sens:=-Sens; MatVP[NoAxe][4]:=Pi/PartPi*Sens; RotationMatrice(MatVP[NoAxe], MatRot, RepereActuel); end; Robservateur : {!!! Rep_re Observateur = Pr_-Multiplier !!!} MultMat(DimenMat, XMat[NoAxe], MatRot, MatRot); end else {!!! Rep_re Observateur = Pr_-Multiplier !!!} MultMat(DimenMat, XMat[NoAxe], MatRot, MatRot); end; { Entr_e : AA : Tableau de valeurs enti_res N : Taille du tableau AA Offsx : Offset x _ rajouter _ toutes les valeurs de BB.x Offsy : Offset y _ rajouter _ toutes les valeurs de BB.y X1, Xm : Indice de d_but et de fin du tableau _ agrandir Sortie : BB : Tableau de points sur l'_cran M : Taille tu tableau BB _ remplir } procedure CSpline(var AA:_TabVal; N : integer; X1, Xm, Echelle : Float; var BB : _TabPoints; M, Offsx, Offsy : integer); var I, K : integer; Dx, T, P: Float; B, C, D : P_TabFloat; { AA : _TabVal absolute AAA; } label 1,2,3; function SplineEval( T : Float; var I : integer) : Float; var J, K : integer; Dx : Float; begin if I >= N then I := 1; if (T < I) or (T > (I+1)) then begin I := 1; J := N + 1; repeat K := (I + J) div 2; if T < K then J := K; if T >= K then I := K; until J <= (I + 1); end; Dx := T - I; SplineEval := AA[I] + Dx * (B^[I] + Dx * (C^[I] + Dx * D^[I])); end; { SplineEval } begin { CSpline } if ((N+1)*SizeOf(B^[1]))>MaxAvail then Goto 3; GetMem(B, (N+1)*SizeOf(B^[1])); if ((N+1)*SizeOf(C^[1]))>MaxAvail then Goto 2; GetMem(C, (N+1)*SizeOf(C^[1])); if ((N+1)*SizeOf(D^[1]))>MaxAvail then Goto 1; GetMem(D, (N+1)*SizeOf(D^[1])); if N >= 3 then begin D^[1] := 2 - 1; C^[2] := (AA[2] - AA[1]) / D^[1]; for I := 2 to N-1 do begin D^[I] := I+1 - I; B^[I] := 2.0 * (D^[I-1] + D^[I]); C^[I+1] := (AA[I+1] - AA[I]) / D^[I]; C^[I] := C^[I+1] - C^[I]; end; B^[1] := -D^[1]; B^[N] := -D^[N-1]; C^[1] := 0.0; C^[N] := 0.0; if N > 3 then begin C^[1] := C^[3] / (4 - 2) - C^[2] / (3 - 1); C^[N] := C^[N-1] / (N - (N-2)) - C^[N-2] / ((N-1) - (N-3)); C^[1] := C^[1] * Sqr(D^[1]) / (4 - 1); C^[N] := -C^[N] * Sqr(D^[N-1]) / (N- (N-3)); end; for I := 2 to N do begin T := D^[I-1] / B^[I-1]; B^[I] := B^[I] - T * D^[I-1]; C^[I] := C^[I] - T * C^[I-1]; end; C^[N] := C^[N] / B^[N]; for I := N-1 downto 1 do C^[I] := (C^[I] - D^[I] * C^[I+1]) / B^[I]; B^[N] := (AA[N] - AA[N-1]) / D^[N-1] + D^[N-1] * (C^[N-1] + 2.0 * C^[N]); for I := 1 to N-1 do begin B^[I] := (AA[I+1] - AA[I]) / D^[I] - D^[I] * (C^[I+1] + 2.0 * C^[I]); D^[I] := (C^[I+1] - C^[I]) / D^[I]; C^[I] := 3.0 * C^[I]; end; C^[N] := 3.0 * C^[N]; D^[N] := D^[N-1]; end else if N = 2 then begin B^[1] := (AA[2] - AA[1]) / (2 - 1); C^[1] := 0.0; D^[1] := 0.0; B^[2] := B^[1]; C^[2] := 0.0; D^[2] := 0.0; end; if (N >= 2) and (M >= 2) then if (X1 >= 1) and (Xm <= N) then begin Dx := (Xm - X1) / (M - 1); K := 1; for I := 1 to M do begin P := X1 + (I - 1) * Dx; BB[I].y := -Round(SplineEval(P, K)*Echelle)+Offsy; BB[I].x := I+Offsx; end; end else { MtError(20, 7) } else { MtError(20, 4)}; FreeMem(D, (N+1)*SizeOf(D^[1])); 1: FreeMem(C, (N+1)*SizeOf(C^[1])); 2: FreeMem(B, (N+1)*SizeOf(B^[1])); 3: end; { CSpline } procedure Bezier(var A : _TabVal; MaxContrPoints : integer; var B : _TabPoints; MaxIntPoints : integer); const MaxControlPoints = 25; type CombiArray = array[0..MaxControlPoints] of Float; var N : integer; ContrPoint, IntPoint : integer; T, SumX, SumY, Prod, DeltaT, Quot : Float; Combi : CombiArray; begin MaxContrPoints := MaxContrPoints - 1; DeltaT := 1.0 / (MaxIntPoints - 1); Combi[0] := 1; Combi[MaxContrPoints] := 1; for N := 0 to MaxContrPoints - 2 do Combi[N + 1] := Combi[N] * (MaxContrPoints - N) / (N + 1); for IntPoint := 1 to MaxIntPoints do begin T := (IntPoint - 1) * DeltaT; if T <= 0.5 then begin Prod := 1.0 - T; Quot := Prod; for N := 1 to MaxContrPoints - 1 do Prod := Prod * Quot; Quot := T / Quot; SumX := MaxContrPoints + 1; SumY := A[MaxContrPoints + 1]; for N := MaxContrPoints downto 1 do begin SumX := Combi[N - 1] * N + Quot * SumX; SumY := Combi[N - 1] * A[N] + Quot * SumY; end; end else begin Prod := T; Quot := Prod; for N := 1 to MaxContrPoints - 1 do Prod := Prod * Quot; Quot := (1 - T) / Quot; SumX := 1; SumY := A[1]; for N := 1 to MaxContrPoints do begin SumX := Combi[N] * (N + 1) + Quot * SumX; SumY := Combi[N] * A[N + 1] + Quot * SumY; end; end; B[IntPoint].x := Round(SumX * Prod); B[IntPoint].y := Round(SumY * Prod); end; end; { Bezier } {===================================================================} { Procedure jacobi + ordonne d_j_ _crite le 31/8/82 1982 (voir documents) } {===================================================================} procedure CalculVecteurPropreB(var EigenVectors : TNmatrix; var PVecteur : _TabPVal; Taille : Word); const Precis=1.6E-9; type VCTI = array[1..DimenMat] of Integer; var i,j : Integer; SinT, CosT, Controle, Nu : Float; VCI,VCIA,WCIR:TNvector; SI, SCI, SCRI, TSCRI, AMat, SMat : TNMatrix; Indicat, QX,P,Q : Integer; AXIND, ColInd : VCTI; procedure CMATA0D; var i,j :Integer; Reduc:Float; begin ClearMat(AMat); for i:=0 to Taille-1 do begin Amat[1,1]:=AMat[1,1]+Int(PVecteur[1]^[i])*Int(PVecteur[1]^[i]); Amat[1,2]:=AMat[1,2]+Int(PVecteur[1]^[i])*Int(PVecteur[2]^[i]); Amat[1,3]:=AMat[1,3]-Int(PVecteur[1]^[i])*Int(PVecteur[3]^[i]); Amat[2,1]:=AMat[2,1]+Int(PVecteur[2]^[i])*Int(PVecteur[1]^[i]); Amat[2,2]:=AMat[2,2]+Int(PVecteur[2]^[i])*Int(PVecteur[2]^[i]); Amat[2,3]:=AMat[2,3]-Int(PVecteur[2]^[i])*Int(PVecteur[3]^[i]); Amat[3,1]:=AMat[3,1]-Int(PVecteur[3]^[i])*Int(PVecteur[1]^[i]); Amat[3,2]:=AMat[3,2]-Int(PVecteur[3]^[i])*Int(PVecteur[2]^[i]); Amat[3,3]:=AMat[3,3]+Int(PVecteur[3]^[i])*Int(PVecteur[3]^[i]); end; {$IFNDEF OPE } WriteMat('MatA=', AMat); {$ENDIF} Reduc:=0; for i:= 1 to DimenMat do for j:= 1 to DimenMat do begin if AMat[i,j]>Reduc then Reduc:=AMat[i,j]; end; DivMatFloat(DimenMat, AMat, Reduc); {$IFNDEF OPE } WriteMat('MatA normalis_e =', AMat); {$ENDIF} end; procedure Trigo; var Omega, OmegaP, Lambda, Mu : Float; begin Lambda:=-AMat[P,Q]; MU:=(AMat[P,P]-AMat[Q,Q])/2; OmegaP:=Lambda/Sqrt(sqr(Lambda)+sqr(MU)); if MU<0 then Omega:=-OmegaP else Omega:=OmegaP; SinT:=Omega/Sqrt(2*(1+sqrt(1-sqr(Omega)))); CosT:=sqrt(1-sqr(SinT)); Controle:=sqr(SinT)+sqr(CosT); end; procedure Clcas; var BMat, SP : TNMatrix; CC,SS,SC : Float; I : Integer; begin for i:=1 to DimenMat do begin BMat[i,P]:=AMat[i,P]*CosT-AMat[i,Q]*SinT; BMat[i,Q]:=AMat[i,P]*SinT+AMat[i,Q]*CosT; SP[i,P]:=SMat[i,p]*CosT-SMat[i,Q]*SinT; SP[i,Q]:=SMat[i,p]*SinT+SMat[i,Q]*CosT; end; CC:=sqr(CosT); SS:=sqr(SinT); SC:=SinT*CosT; BMat[P,P]:=AMat[P,P]*CC+AMat[Q,Q]*SS-2*AMAt[P,Q]*SC; BMat[Q,Q]:=AMat[P,P]*SS+AMat[Q,Q]*CC+2*AMAt[P,Q]*SC; BMat[P,Q]:=(AMat[P,P]-AMat[Q,Q])*SC+AMat[P,Q]*(CC-SS); BMat[Q,P]:=BMat[P,Q]; for i:=1 to DimenMat do begin AMat[i,P]:=BMat[i,P]; AMat[i,Q]:=BMat[i,Q]; SMat[i,P]:=SP[i,P]; SMat[i,Q]:=SP[i,Q]; end; for i:=1 to DimenMat do begin AMat[P,i]:=BMat[i,P]; AMat[Q,i]:=BMat[i,Q]; end; end; procedure ClcNU; { Racine de la somme des termes non diagonaux } var S2 : Float; i,j : Integer; begin S2:=0; for i:=1 to DimenMat-1 do for j:=i+1 to DimenMat do S2:=S2+AMat[i,j]*AMat[i,j]; NU:=sqrt(S2*2); Indicat:=0; end; procedure AXI; { On passe en coordonn_es inverses sur z } var i,j : Integer; begin {$IFNDEF OPE } WriteMat('SMat =', SMat); {$ENDIF} for i:=1 to DimenMat do for j:=1 to DimenMat do SI[i,j] := SMat[i,j]; for i:=1 to DimenMat do SI[DimenMat,i]:=-SMat[DimenMat,i]; {$IFNDEF OPE } WriteMat('SI =', SI); {$ENDIF} end; procedure CLASNBD(V:TNvector; var W:TNvector; var INDCL:VCTI); var Max:Float; i,j,iMax:Integer; begin for i:=1 to DimenMat do begin W[i]:=V[i]; INDCL[i]:=i; end; for i:=1 to DimenMat-1 do begin Max:=W[i]; for j:=i+1 to DimenMat do begin if W[j]>Max then begin Max:=W[j]; IMax:=INDCL[j]; W[j]:=W[i]; INDCL[j]:=INDCL[i]; W[i]:=MAx; INDCL[i]:=IMax; end; end; end; end; procedure CLASAX; { Classe les colonnes par ordre de grandeur de AMat[i,i] } var V,W:TNvector; i,j:Integer; begin for i:=1 to DimenMat do V[i]:=AMat[i,i]; {$IFNDEF OPE } WriteVect('EigenValues =', V); {$ENDIF} CLASNBD(V,W,Colind); for i:=1 to DimenMat do for j:=1 to DimenMat do SCI[i,j]:=SI[i, Colind[j]]; {$IFNDEF OPE } WriteMat('SCI =', SCI); {$ENDIF} end; procedure NMORAXE; var i,j,l,p:Integer; begin for j:=1 to DimenMat do begin for i:=1 to DimenMat do begin VCI[i]:=SCI[i,j]; VCIA[i]:=Abs(VCI[i]); end; CLASNBD(VCIA, WCIR, AXIND); l:=AXIND[1]; if SCI[l,j]>0 then for p:=1 to DimenMat do SCRI[p,j]:=SCI[p,j] else for p:=1 to DimenMat do SCRI[p,j]:=-SCI[p,j]; end; {$IFNDEF OPE } WriteMat('SCRI =', SCRI); {$ENDIF} end; procedure Creat; begin TSCRI:=SCRI; { TransposeMat(DimenMat, TSCRI); } EigenVectors:=TSCRI; {$IFNDEF OPE } WriteMat('EigenVector =', TSCRI); {$ENDIF} end; procedure NPoint; begin AXI; CLASAX; NMORAXE; CREAT; end; begin { Calcn } Initial(DimenMat, SMat); CMaTA0D; ClcNU; Nu:=Nu/DimenMat; QX:=2; P:=1; Q:=2; While NU>Precis do begin While P<=(DimenMat-1) do begin While QX<=DimenMat do begin Q:=QX; if Abs(AMat[P,Q])>=NU then begin Indicat:=1; Trigo; ClCas; end; QX:=QX+1; end; P:=P+1; QX:=P+1; end; NU:=NU/DimenMat; QX:=2; P:=1; end; NPoint; end; { Calcn } { Calcule la surface de l'ombre d'un objet projet_ sur l'_cran } function SurfaceOmbre(var Mat:TNMatrix; var Points:_TabPVal; Modulo, OffsetS, Count:Integer):float; const Org = 0; var i,index: integer; VA,VB,VC : P_TabVal; A11,A12,A13, A21,A22,A23, A31,A32,A33, x1, x2, y1, y2, A,B,C, Res : Float; { Surface du triangle } function S:Float; begin A := x2*y2; B := x1*y1; C := (y1-y2)*(x2-x1); S := (y1*x2 - (A + B + C) / 2); end; begin Res := 0; A11:=Mat[1,1]; A12:=Mat[1,2]; A13:=Mat[1,3]; A21:=Mat[2,1]; A22:=Mat[2,2]; A23:=Mat[2,3]; A31:=Mat[3,1]; A32:=Mat[3,2]; A33:=Mat[3,3]; { Rep_re les 3 vecteurs _ d_placer } VA:=Points[1]; { x } VB:=Points[2]; { y } VC:=Points[3]; { z } { Change les vecteurs } for index := OffsetS to (Count+OffsetS) do begin i := Index mod Modulo; A:= VA^[i]; B:= VB^[i]; C:=VC^[i]; x2:= x1; y2:= y1; x1:= (A*A11 + B*A12 + C*A13); y1:= (A*A21 + B*A22 + C*A23); if index>OffsetS then Res := Res + S ; end; { 5 ÝV/Bit } SurfaceOmbre := Abs(Res)*CoefSurface; end; begin MathError; end. ÿÿÿ