|
|
@@ -439,79 +439,11 @@ end;
|
|
|
{$endif not SYSTEM_HAS_LDEXP}
|
|
|
|
|
|
|
|
|
- function polevl(x:Real; Coef:PReal; N:sizeint):Real;
|
|
|
- {*****************************************************************}
|
|
|
- { 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 : Real;
|
|
|
- i : sizeint;
|
|
|
-
|
|
|
- begin
|
|
|
- ans := Coef[0];
|
|
|
- for i:=1 to N do
|
|
|
- ans := ans * x + Coef[i];
|
|
|
- polevl:=ans;
|
|
|
- end;
|
|
|
-
|
|
|
-
|
|
|
- function p1evl(x:Real; Coef:PReal; N:sizeint):Real;
|
|
|
- { }
|
|
|
- { Evaluate polynomial when coefficient of x is 1.0. }
|
|
|
- { Otherwise same as polevl. }
|
|
|
- { }
|
|
|
- var
|
|
|
- ans : Real;
|
|
|
- i : sizeint;
|
|
|
+ function floord(x: double): double; inline;
|
|
|
begin
|
|
|
- ans := x + Coef[0];
|
|
|
- for i:=1 to N-1 do
|
|
|
- ans := ans * x + Coef[i];
|
|
|
- p1evl := ans;
|
|
|
- end;
|
|
|
-
|
|
|
-
|
|
|
- function floord(x: double): double;
|
|
|
- var
|
|
|
- t: double;
|
|
|
- begin
|
|
|
- t := int(x);
|
|
|
- if (x>=0.0) or (x=t) then
|
|
|
- floord := t
|
|
|
- else
|
|
|
- floord := t - 1.0;
|
|
|
+ result := int(x);
|
|
|
+ if result>x then
|
|
|
+ result := result - 1.0;
|
|
|
end;
|
|
|
|
|
|
{*
|
|
|
@@ -1281,10 +1213,9 @@ end;
|
|
|
{ of the fractional part: }
|
|
|
{ e**x - 1 = 2x P(x**2)/( Q(x**2) - P(x**2) ) }
|
|
|
xx := d * d;
|
|
|
- px := d * polevl( xx, P, 2 );
|
|
|
- d := px/( polevl( xx, Q, 3 ) - px );
|
|
|
- d := ldexp( d, 1 );
|
|
|
- d := d + 1.0;
|
|
|
+ px := d * ((P[0] * xx + P[1]) * xx + P[2]);
|
|
|
+ d := px/( (((Q[0] * xx + Q[1]) * xx + Q[2]) * xx + Q[3]) - px );
|
|
|
+ d := 2 * d + 1.0;
|
|
|
d := ldexp( d, n );
|
|
|
result := d;
|
|
|
end;
|
|
|
@@ -1429,8 +1360,6 @@ end;
|
|
|
Lg6: double = 1.531383769920937332e-01; { 3FC39A09 D078C69F }
|
|
|
Lg7: double = 1.479819860511658591e-01; { 3FC2F112 DF3E5244 }
|
|
|
|
|
|
- zero: double = 0.0;
|
|
|
-
|
|
|
var
|
|
|
hfsq,f,s,z,R,w,t1,t2,dk: double;
|
|
|
k,hx,i,j: longint;
|
|
|
@@ -1538,27 +1467,23 @@ end;
|
|
|
{ Between pi/4 and pi/2 the cosine is represented as }
|
|
|
{ 1 - x**2 Q(x**2). }
|
|
|
{*****************************************************************}
|
|
|
- var y, z, zz : Real;
|
|
|
+ var z, zz : Real;
|
|
|
j : sizeint;
|
|
|
|
|
|
begin
|
|
|
{ This seemingly useless condition ensures that sin(-0.0)=-0.0 }
|
|
|
if (d=0.0) then
|
|
|
exit(d);
|
|
|
- j := rem_pio2(d,z) and 3;
|
|
|
+ j := rem_pio2(d,z);
|
|
|
|
|
|
zz := z * z;
|
|
|
-
|
|
|
- if( (j=1) or (j=3) ) then
|
|
|
- y := 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 )
|
|
|
+ if j and 1<>0 then { j and 3 = 1 or j and 3 = 3 }
|
|
|
+ result := 1.0 - zz * 0.5 + zz * zz * (((((coscof[0] * zz + coscof[1]) * zz + coscof[2]) * zz + coscof[3]) * zz + coscof[4]) * zz + coscof[5])
|
|
|
else
|
|
|
- { y = z + z * (zz * polevl( zz, sincof, 5 )); }
|
|
|
- y := z + z * z * z * polevl( zz, sincof, 5 );
|
|
|
+ result := z + zz * z * (((((sincof[0] * zz + sincof[1]) * zz + sincof[2]) * zz + sincof[3]) * zz + sincof[4]) * zz + sincof[5]);
|
|
|
|
|
|
- if (j > 1) then
|
|
|
- result := -y
|
|
|
- else
|
|
|
- result := y;
|
|
|
+ if j and 2<>0 then { j and 3 = 2 or j and 3 = 3 }
|
|
|
+ result := -result;
|
|
|
end;
|
|
|
{$endif}
|
|
|
|
|
|
@@ -1593,20 +1518,16 @@ end;
|
|
|
var y, z, zz : Real;
|
|
|
j : sizeint;
|
|
|
begin
|
|
|
- j := rem_pio2(d,z) and 3;
|
|
|
+ j := rem_pio2(d,z);
|
|
|
|
|
|
zz := z * z;
|
|
|
-
|
|
|
- if( (j=1) or (j=3) ) then
|
|
|
- { y = z + z * (zz * polevl( zz, sincof, 5 )); }
|
|
|
- y := z + z * z * z * polevl( zz, sincof, 5 )
|
|
|
+ if j and 1<>0 then { j and 3 = 1 or j and 3 = 3 }
|
|
|
+ result := z + zz * z * (((((sincof[0] * zz + sincof[1]) * zz + sincof[2]) * zz + sincof[3]) * zz + sincof[4]) * zz + sincof[5])
|
|
|
else
|
|
|
- y := 1.0 - ldexp(zz,-1) + zz * zz * polevl( zz, coscof, 5 );
|
|
|
+ result := 1.0 - zz * 0.5 + zz * zz * (((((coscof[0] * zz + coscof[1]) * zz + coscof[2]) * zz + coscof[3]) * zz + coscof[4]) * zz + coscof[5]);
|
|
|
|
|
|
- if (j = 1) or (j = 2) then
|
|
|
- result := -y
|
|
|
- else
|
|
|
- result := y ;
|
|
|
+ if (j+1) and 2<>0 then { j and 3 = 1 or j and 3 = 2 }
|
|
|
+ result := -result;
|
|
|
end;
|
|
|
{$endif}
|
|
|
|
|
|
@@ -1766,10 +1687,6 @@ end;
|
|
|
{$OPTIMIZATION NOFASTMATH}
|
|
|
{$endif VER3_2}
|
|
|
function fpc_frac_real(d : ValReal) : ValReal;compilerproc;
|
|
|
-{$ifndef VER3_2}
|
|
|
- const
|
|
|
- zero: ValReal = 0.0;
|
|
|
-{$endif VER3_2}
|
|
|
begin
|
|
|
{ Nan or +/-Inf }
|
|
|
if (float64high(d) and $7ff00000)=$7ff00000 then
|
Rika Ichinose