diff --git a/rtl/template/math.inc b/rtl/template/math.inc deleted file mode 100644 index a4cfaafcc9..0000000000 --- a/rtl/template/math.inc +++ /dev/null @@ -1,885 +0,0 @@ -{ - $Id$ - This file is part of the Free Pascal run time library. - Copyright (c) 1999-2000 by xxxx - member of the Free Pascal development team - - See the file COPYING.FPC, included in this distribution, - for details about the copyright. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. - - **********************************************************************} -{*************************************************************************} -{ math.inc } -{*************************************************************************} -{ Copyright Abandoned, 1987, Fred Fish } -{ } -{ This previously copyrighted work has been placed into the } -{ public domain by the author (Fred Fish) and may be freely used } -{ for any purpose, private or commercial. I would appreciate } -{ it, as a courtesy, if this notice is left in all copies and } -{ derivative works. Thank you, and enjoy... } -{ } -{ The author makes no warranty of any kind with respect to this } -{ product and explicitly disclaims any implied warranties of } -{ merchantability or fitness for any particular purpose. } -{-------------------------------------------------------------------------} -{ Copyright (c) 1992 Odent Jean Philippe } -{ } -{ The source can be modified as long as my name appears and some } -{ notes explaining the modifications done are included in the file. } -{-------------------------------------------------------------------------} -{ Copyright (c) 1997 Carl Eric Codere } -{ } -{ This include implements a template for } -{ all real (whatever this type maps to) and fixed point standard } -{ rtl routines. } -{ NOTE: Trunc and Int must be implemented depending on the target } -{ for real values.Sqrt must also be implemented for fixed. } -{*************************************************************************} - - - -type - TabCoef = array[0..6] of Real; - - -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; - MAXLOG = 8.8029691931113054295988E1; { log(2**127) } - MINLOG = -8.872283911167299960540E1; { log(2**-128) } - - DP1 = 7.85398125648498535156E-1; - DP2 = 3.77489470793079817668E-8; - DP3 = 2.69515142907905952645E-15; - -const sincof : TabCoef = ( - 1.58962301576546568060E-10, - -2.50507477628578072866E-8, - 2.75573136213857245213E-6, - -1.98412698295895385996E-4, - 8.33333333332211858878E-3, - -1.66666666666666307295E-1, 0); - coscof : TabCoef = ( - -1.13585365213876817300E-11, - 2.08757008419747316778E-9, - -2.75573141792967388112E-7, - 2.48015872888517045348E-5, - -1.38888888888730564116E-3, - 4.16666666666665929218E-2, 0); - - - function int(d : real) : real; - { these routine should be implemented all depending on the } - { target processor/operating system. } - begin - end; - - function trunc(d : real) : longint; - { these routine should be implemented all depending on the } - { target processor/operating system. } - Begin - end; - - - function abs(d : Real) : Real; - begin - if( d < 0.0 ) then - abs := -d - else - abs := d ; - end; - - - function frexp(x:Real; var e:Integer ):Real; - {* 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 (abs(x)<0.5) then - While (abs(x)<0.5) do - begin - x := x*2; - Dec(e); - end - else - While (abs(x)>1) do - begin - x := x/2; - Inc(e); - end; - frexp := x; - end; - - - function ldexp( x: Real; N: Integer):Real; - {* ldexp() multiplies x by 2**n. *} - var r : Real; - 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 polevl(var x:Real; var Coef:TabCoef; N:Integer):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 : Integer; - - begin - ans := Coef[0]; - for i:=1 to N do - ans := ans * x + Coef[i]; - polevl:=ans; - end; - - - function p1evl(var x:Real; var Coef:TabCoef; N:Integer):Real; - { } - { Evaluate polynomial when coefficient of x is 1.0. } - { Otherwise same as polevl. } - { } - var - ans : Real; - i : Integer; - begin - ans := x + Coef[0]; - for i:=1 to N-1 do - ans := ans * x + Coef[i]; - p1evl := ans; - end; - - - - - - function sqr(d : Real) : Real; - begin - sqr := d*d; - end; - - - function pi : Real; - begin - pi := 3.1415926535897932385; - end; - - - function sqrt(x:Real):Real; - {*****************************************************************} - { 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 : Real; - begin - if( x <= 0.0 ) then - begin - if( x < 0.0 ) then - RunError(207); - 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 Exp(x:Real):Real; - {*****************************************************************} - { 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 : Real; - begin - if( x > MAXLOG) then - RunError(205) - else - if( x < MINLOG ) then - begin - Runerror(205); - 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 := Trunc( px + 0.5 ); { Trunc() 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 Round(x: Real): longint; - var - fr: Real; - tr: Real; - Begin - fr := Frac(x); - tr := Trunc(x); - if fr > 0.5 then - Round:=Trunc(x)+1 - else - if fr < 0.5 then - Round:=Trunc(x) - else { fr = 0.5 } - { check sign to decide ... } - { as in Turbo Pascal... } - if x >= 0.0 then - Round := Trunc(x)+1 - else - Round := Trunc(x); - end; - - - function Ln(x:Real):Real; - {*****************************************************************} - { Natural Logarithm } - {*****************************************************************} - { } - { SYNOPSIS: } - { } - { double x, y, log(); } - { } - { y = ln( 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 : Real; - - Label Ldone; - begin - if( x <= 0.0 ) then - RunError(207); - 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; - - Ln:= z; - end; - - - - function Sin(x:Real):Real; - {*****************************************************************} - { Circular Sine } - {*****************************************************************} - { } - { SYNOPSIS: } - { } - { double x, y, sin(); } - { } - { y = sin( 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 sine is approximated by } - { x + x**3 P(x**2). } - { Between pi/4 and pi/2 the cosine is represented as } - { 1 - x**2 Q(x**2). } - {*****************************************************************} - var y, z, zz : Real; - j, sign : Integer; - - begin - { make argument positive but save the sign } - sign := 1; - if( x < 0 ) then - begin - x := -x; - sign := -1; - end; - - { above this value, approximate towards 0 } - if( x > lossth ) then - begin - sin := 0.0; - exit; - end; - - y := Trunc( x/PIO4 ); { integer part of x/PIO4 } - - { strip high bits of integer part to prevent integer overflow } - z := ldexp( y, -4 ); - z := Trunc(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 frac(d : Real) : Real; - begin - frac := d - Int(d); - end; - - - function sqrt(d : fixed) : fixed; - begin - end; - - - function Cos(x:Real):Real; - {*****************************************************************} - { Circular cosine } - {*****************************************************************} - { } - { 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 : Real; - j, sign : Integer; - i : LongInt; - begin - { make argument positive } - sign := 1; - if( x < 0 ) then - x := -x; - - { above this value, round towards zero } - if( x > lossth ) then - begin - cos := 0.0; - exit; - end; - - y := Trunc( x/PIO4 ); - z := ldexp( y, -4 ); - z := Trunc(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 ArcTan(x:Real):Real; - {*****************************************************************} - { 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 : Real; - 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; - Arctan := y; - end; - - - function int(d : fixed) : fixed; - {*****************************************************************} - { Returns the integral part of d } - {*****************************************************************} - begin - int:=d and $ffff0000; { keep only upper bits } - end; - - - function trunc(d : fixed) : longint; - {*****************************************************************} - { Returns the Truncated integral part of d } - {*****************************************************************} - begin - trunc:=longint(integer(d shr 16)); { keep only upper 16 bits } - end; - - function frac(d : fixed) : fixed; - {*****************************************************************} - { Returns the Fractional part of d } - {*****************************************************************} - begin - frac:=d AND $ffff; { keep only decimal parts - lower 16 bits } - end; - - function abs(d : fixed) : fixed; - {*****************************************************************} - { Returns the Absolute value of d } - {*****************************************************************} - var - w: integer; - begin - w:=integer(d shr 16); - if w < 0 then - begin - w:=-w; { invert sign ... } - d:=d and $ffff; - d:=d or (fixed(w) shl 16); { add this to fixed number ... } - abs:=d; - end - else - abs:=d; { already positive... } - end; - - - function sqr(d : fixed) : fixed; - {*****************************************************************} - { Returns the Absolute squared value of d } - {*****************************************************************} - begin - {16-bit precision needed, not 32 =)} - sqr := d*d; -{ sqr := (d SHR 8 * d) SHR 8; } - end; - - - function Round(x: fixed): longint; - {*****************************************************************} - { Returns the Rounded value of d as a longint } - {*****************************************************************} - var - lowf:integer; - highf:integer; - begin - lowf:=x and $ffff; { keep decimal part ... } - highf :=integer(x shr 16); - if lowf > 5 then - highf:=highf+1 - else - if lowf = 5 then - begin - { here we must check the sign ... } - { if greater or equal to zero, then } - { greater value will be found by adding } - { one... } - if highf >= 0 then - Highf:=Highf+1; - end; - Round:= longint(highf); - end; - - function power(bas,expo : real) : real; - begin - if bas=0 then - begin - if expo<>0 then - power:=0.0 - else - HandleError(207); - end - else if expo=0 then - power:=1 - else - { bas < 0 is not allowed } - if bas<0 then - handleerror(207) - else - power:=exp(ln(bas)*expo); - end; - - function power(bas,expo : longint) : longint; - begin - if bas=0 then - begin - if expo<>0 then - power:=0 - else - HandleError(207); - end - else if expo=0 then - power:=1 - else - begin - if bas<0 then - begin - if odd(expo) then - power:=-round(exp(ln(-bas)*expo)) - else - power:=round(exp(ln(-bas)*expo)); - end - else - power:=round(exp(ln(bas)*expo)); - end; - end; - - -{ - $Log$ - Revision 1.3 2002-09-07 16:01:26 peter - * old logs removed and tabs fixed - -}