PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 1/7
TITLE : TRANSCENDENTAL FUNCTIONS

The following example routines are public domain programs that
have been uploaded to our Forum on CompuServe. As a courtesy to
Technical Support distributes these routines free of charge.

However, because these routines are public domain programs, not
developed by Borland International, we are unable to provide any
technical support or assistance using these routines. If you need
assistance using these routines, or are experiencing
difficulties, we recommend that you log onto CompuServe and
request assistance from the Forum members that developed these
routines.

Written by Randall A. Gacek

This is a first approximation of a set of routines to do the
transcendental functions LOG, LN, SQRT, ARCTAN, SIN, COS and EXP
in the BCD version of Turbo Pascal.

WARNING: The following code is specific to the implementation of
Turbo Pascal with BCD support. These functions should only be
used with this implementation.

program checkfuncs;

function sqrt(x:real):real;

var
n,i,m :integer;
f,y :real;
v :record case boolean of
true:(y:real);
false:(z:array[1..10] of byte)
end;
begin
if x = 0.0 then
sqrt:=0.0
else if x < 0.0 then
halt
else begin
v.y:=x;
n:=v.z[1]-63;
v.z[1]:=63;

PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 2/7
TITLE : TRANSCENDENTAL FUNCTIONS

f:=v.y;
y:=0.580661+f/2.0-0.086462/(f+0.175241);
for i:=1 to 2 do
y:=0.5*(y+f/y);
y:=y+0.5*(f/y-y);
if odd(n) then
begin
y:=y*0.316227766016837933;
n:=n+1;
end;

checkfuncs Con't.

m:=n div 2;
v.y:=y;
v.z[1]:=v.z[1]+m;
sqrt:=v.y;
end;
end; { sqrt }

function log(x:real):real;
const
c0= 0.316227766016837933;
a0=-0.260447002405557636E+2;
a1= 0.554085912041205931E+2;
a2=-0.392737410203156250E+2;
a3= 0.103338571514793865E+2;
a4=-0.741010784161919239E+0;
b0=-0.899552077881033117E+2;
b1= 0.245347618868489348E+3;
b2=-0.244303035341829542E+3;
b3= 0.107109789115668009E+3;
b4=-0.193732345832854786E+2;
c= 0.868588963806503655;

var
n:integer;
xn,f,s,w,aw,bw,rs2,rs:real;
v:record case boolean of
true:(y:real);
false:(z:array[1..10] of byte)
end;

PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 3/7
TITLE : TRANSCENDENTAL FUNCTIONS

begin
if x <= 0.0 then
halt;
v.y:=x;
n:=v.z[1]-63;
v.z[1]:=63;
f:=v.y;
if f <= c0 then
begin
n:=n-1;
f:=f*10.0;
end;

checkfuncs Con't.

s:=((f-0.5)-0.5)/(f+1.0);
w:=sqr(s);
aw:= (((a4*w+a3)*w+a2)*w+a1)*w+a0;
bw:=((((w+b4)*w+b3)*w+b2)*w+b1)*w+b0;
rs2:=w*aw/bw;
rs:=s*(c+rs2);
xn:=n;
log:=xn+rs;
end; { log }

function ln(x:real):real;
const
c3=2.30258509299404568;

begin
ln:=log(x)*c3;
end;

function exp(x:real):real;
const
bigx=147.365445951618923;
smallx=-145.062860858624878;
eps=5.0e-19;
onelnsqrt10=0.868588963806503655;
c1=1.151;
c2=2.92546497022842009e-4;
p0=0.333267029226801611e+6;

PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 4/7
TITLE : TRANSCENDENTAL FUNCTIONS

p1=0.100974148724273918E+5;
p2=0.420414268137450315E+2;
q0=0.666534058453603223E+6;
q1=0.757393346159883444E+5;
q2=0.841243584514154545E+3;
sqrt10=3.16227766016837933;
var
n:integer;
xn,g,z,gpz,qz,rg:real;
v:record case boolean of
true:(y:real);
false:(z:array[1..10] of byte)
end;

checkfuncs Con't.

begin
if x > bigx then
halt;
if x < smallx then
halt;
if abs(x) < eps then
exp:=1.0
else begin
n:=round(x*onelnsqrt10);
xn:=n;
g:=(x-xn*c1)-xn*c2;
z:=sqr(g);
gpz:=((p2*z+p1)*z+p0)*g;
qz:= ((z+q2)*z+q1)*z+q0;
rg:=(0.5+gpz/(qz-gpz))*2.0;
if odd(n) then
if n >= 0 then
rg:=sqrt10*rg
else
rg:=rg/sqrt10;
n:=n div 2;
v.y:=rg;
v.z[1]:=v.z[1]+n;
exp:=v.y;
end;
end; { exp }

PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 5/7
TITLE : TRANSCENDENTAL FUNCTIONS

function sincos(x,y,sgn:real):real;
const
ymax=3141592654.0;
onepi=0.318309886183790672;
c1= 3.141; { pi to 22 digits }
c2= 0.000592653589793238463;
eps=1.0e-9;
r1=-0.166666666666666651e+0;
r2= 0.833333333333316503E-2;
r3=-0.198412698412018405E-3;
r4= 0.275573192101527561E-5;
r5=-0.250521067982745845E-7;
r6= 0.160589364903715891E-9;
r7=-0.764291780689104677E-12;
r8= 0.272047909578888462E-14;

checkfuncs Con't.

var
xn,f,t,g,rg:real;
begin
if y >= ymax then
halt;
xn:=y*onepi;
xn:=int(xn+0.5);
if frac(xn / 2.0) <> 0.0 then
sgn:=-sgn;
if abs(x) <> y then { cos wanted }
xn:=xn-0.5;
f:=(abs(x)-xn*c1)-xn*c2;
if abs(f) < eps then
t:=f
else begin
g:=sqr(f);
rg:=(((((((r8*g+r7)*g+r6)*g+r5)*g+r4)*g+r3)*g+r2)*g+r1)*g;
t:=f+f*rg;
end;
sincos:=sgn*t;
end; { sincos }

function sin(x:real):real;
begin
if x < 0.0 then

PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 6/7
TITLE : TRANSCENDENTAL FUNCTIONS

sin:=sincos(x,-x,-1.0)
else
sin:=sincos(x,x,1.0);
end; {sin}

function cos(x:real):real;
begin
cos:=sincos(x,abs(x)+1.57079632679489662,1.0);
end; {cos}

checkfuncs Con't.

function arctan(x:real):real;
const
twomsqrt3=0.267949192431122706;
sqrt3=1.73205080756887729;
a=0.732050807568877294;
eps=1e-9;
p0=-0.136887688941919269e+2;
p1=-0.205058551958616520e+2;
p2=-0.849462403513206835e+1;
p3=-0.837582993681500593e+0;
q0= 0.410663066825757813e+2;
q1= 0.861573495971302425e+2;
q2= 0.595784361425973445e+2;
q3= 0.150240011600285761e+1;

var
n:integer;
f,result,g,gpg,qg,r:real;
begin
f:=abs(x);
if f > 1.0 then
begin
f:=1.0/f;
n:=2;
end
else
n:=0;
if f > twomsqrt3 then
begin
f:=(((a*f-0.5)-0.5)+f)/(sqrt3+f);
n:=n+1;

PRODUCT : TURBO PASCAL WITH BCD SUPPORT NUMBER : 227
VERSION : 3.01
OS : MS-DOS, PC-DOS, CP/M-86
DATE : August 4, 1986 PAGE : 7/7
TITLE : TRANSCENDENTAL FUNCTIONS

end;
if abs(f) < eps then
result:=f
else begin
g:=sqr(f);
gpg:=(((p3*g+p2)*g+p1)*g+p0)*g;
qg:=(((g+q3)*g+q2)*g+q1)*g+q0;
r:=gpg/qg;
result:=f+f*r;
end;
if n > 1 then
result:=-result;

checkfuncs Con't.

case n of
0:;
1:result:=0.523598775598298873+result;
2:result:=1.57079632679489662+result;
3:result:=1.04719755119659775+result;
end;
if x < 0.0 then
result:=-result;
arctan:=result;
end; { arctan }

begin
writeln('sqrt= ',sqrt(25));
writeln('ln = ',ln(25));
writeln('exp = ',exp(25));
writeln('cos = ',cos(25));
writeln('sin = ',sin(25));
writeln('log = ',log(25));
writeln('arctan = ',arctan(25));
end.