paweld/fpc
1
0
Fork 0
mirror of https://gitlab.com/freepascal.org/fpc/source.git synced 2025-04-13 01:49:35 +02:00
Code Issues Packages Projects Releases Wiki Activity
15d85c7e3b
fpc/rtl/template/math.inc
root d1a1f501c8 Initial revision
1998-03-25 11:18:12 +00:00

876 lines
31 KiB
PHP
Raw Blame History

{
$Id$
This file is part of the Free Pascal run time library.
Copyright (c) 1993,97 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
power:=exp(ln(bas)*expo);
end;
function power(bas,expo : longint) : longint;
begin
power:=round(exp(ln(bas)*expo));
end;
{
$Log$
Revision 1.1 1998-03-25 11:18:46 root
Initial revision
Revision 1.4 1998/01/26 12:02:04 michael
+ Added log at the end
Working file: rtl/template/math.inc
description:
----------------------------
revision 1.3
date: 1998/01/05 00:41:29; author: carl; state: Exp; lines: +835 -836
* Licence problems fixed
----------------------------
revision 1.2
date: 1997/12/01 12:45:48; author: michael; state: Exp; lines: +14 -0
+ added copyright reference in header.
----------------------------
revision 1.1
date: 1997/11/28 16:48:02; author: carl; state: Exp;
math.inc template.
=============================================================================
}
Reference in New Issue View Git Blame Copy Permalink