fpc/rtl/powerpc/math.inc
2001-12-02 16:19:45 +00:00

268 lines
7.3 KiB
PHP

{
$Id$
This file is part of the Free Pascal run time library.
Copyright (c) 2000 by Jonas Maebe and other members of the
Free Pascal development team
Implementation of mathamatical Routines (only for real)
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.
**********************************************************************}
{****************************************************************************
EXTENDED data type routines
****************************************************************************}
function pi : double;[internconst:in_pi];
begin
pi := 3.14159265358979320;
end;
function abs(d : extended) : extended;[internproc:in_abs_extended];
function sqr(d : extended) : extended;[internproc:in_sqr_extended];
function sqrt(d : extended) : extended;[internproc:in_sqrt_extended];
function arctan(d : extended) : extended;[internconst:in_arctan_extended];
begin
runerror(207);
end;
function ln(d : extended) : extended;[internconst:in_ln_extended];
begin
runerror(207);
end;
function sin(d : extended) : extended;[internconst: in_sin_extended];
begin
runerror(207);
end;
function cos(d : extended) : extended;[internconst:in_cos_extended];
begin
runerror(207);
end;
function exp(d : extended) : extended;assembler;[internconst:in_const_exp];
begin
runerror(207);
end;
function frac(d : extended) : extended;assembler;[internconst:in_const_frac];
begin
runerror(207);
end;
function int(d : extended) : extended;assembler;[internconst:in_const_int];
begin
runerror(207);
end;
function trunc(d : extended) : longint;assembler;[internconst:in_const_trunc];
{ input: d in fr1 }
{ output: result in r3 }
assembler;
var
temp:
record
case byte of
0: (l1,l2: longint);
1: (d: double);
end;
asm
fctiwz fr1,fr1
stfd fr1,temp.d
lwz r3,temp.l2
end ['r3','fr1'];
function round(d : extended) : longint;assembler;[internconst:in_const_round];
{ input: d in fr1 }
{ output: result in r3 }
assembler;
var
temp:
record
case byte of
0: (l1,l2: longint);
1: (d: double);
end;
asm
fctiw fr1,fr1
stfd fr1,temp.d
lwz r3,temp.l2
end ['r3','fr1'];
function power(bas,expo : extended) : extended;
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;
{****************************************************************************
Longint data type routines
****************************************************************************}
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;
{****************************************************************************
Helper routines to support old TP styled reals
****************************************************************************}
{ warning: the following converts a little-endian TP-style real }
{ to a big-endian double. So don't byte-swap the TP real! }
function real2double(r : real48) : double;
var
res : array[0..7] of byte;
exponent : word;
begin
{ copy mantissa }
res[6]:=0;
res[5]:=r[1] shl 5;
res[4]:=(r[1] shr 3) or (r[2] shl 5);
res[3]:=(r[2] shr 3) or (r[3] shl 5);
res[2]:=(r[3] shr 3) or (r[4] shl 5);
res[1]:=(r[4] shr 3) or (r[5] and $7f) shl 5;
res[0]:=(r[5] and $7f) shr 3;
{ copy exponent }
{ correct exponent: }
exponent:=(word(r[0])+(1023-129));
res[1]:=res[1] or ((exponent and $f) shl 4);
res[0]:=exponent shr 4;
{ set sign }
res[0]:=res[0] or (r[5] and $80);
real2double:=double(res);
end;
{****************************************************************************
Int to real helpers
****************************************************************************}
function fpc_int64_to_real(i: int64): double; compilerproc;
assembler;
{ input: high(i) in r3, low(i) in r4 }
{ output: double(i) in f0 }
var
temp:
record
case byte of
0: (l1,l2: cardinal);
1: (d: double);
end;
asm
lis r0,$4330
stw r0,temp.l1
xoris r3,r3,$8000
stw r3,temp.l2
lis r3,longint_to_real_helper@ha
lfd fr1,longint_to_real_helper@l(r3)
lfd fr0,temp.d
stw r4,temp.l2
fsub fr0,fr0,fr1
lis r4,cardinal_to_real_helper@ha
lfd fr1,cardinal_to_real_helper@l(r4)
lis r3,int_to_real_factor@ha
lfd fr3,temp
lfd fr2,int_to_real_factor@l(r3)
fsub fr3,fr3,fr1
fmadd fr1,fr0,fr2,fr3
end ['r0','r3','r4','fr0','fr1','fr2','fr3'];
function fpc_qword_to_real(q: qword): double; compilerproc;
assembler;
{ input: high(q) in r3, low(q) in r4 }
{ output: double(q) in f0 }
var
temp:
record
case byte of
0: (l1,l2: cardinal);
1: (d: double);
end;
asm
lis r0,$4330
stw r0,temp.l1
stw r3,temp.l2
lfd fr0,temp.d
lis r3,cardinal_to_real_helper@ha
lfd fr1,cardinal_to_real_helper@l(r3)
stw r4,temp.l2
fsub fr0,fr0,fr1
lfd fr3,temp
lis r3,int_to_real_factor@ha
lfd fr2,int_to_real_factor@l(r3)
fsub fr3,fr3,fr1
fmadd fr1,fr0,fr2,fr3
end ['r0','r3','fr0','fr1','fr2','fr3'];
{
$Log$
Revision 1.3 2001-12-02 16:19:45 jonas
* fpu results are returned in fr1, not fr0
Revision 1.2 2001/10/30 17:18:14 jonas
* fixed fpc_int64_to_double and fpc_int64_to_double (fpc_int64_to_double
is now mostly tested and should work fine, fpc_qword_to_double should
work too since it's almost the same)
Revision 1.1 2001/10/28 14:09:13 jonas
+ initial implementation, lots of things still missing
}