mirror of
https://gitlab.com/freepascal.org/fpc/source.git
synced 2025-04-21 00:29:24 +02:00
524 lines
15 KiB
PHP
524 lines
15 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 mathematical 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.
|
|
|
|
**********************************************************************}
|
|
|
|
|
|
const
|
|
longint_to_real_helper: int64 = $4330000080000000;
|
|
cardinal_to_real_helper: int64 = $4330000000000000;
|
|
int_to_real_factor: double = double(high(cardinal))+1.0;
|
|
|
|
|
|
{****************************************************************************
|
|
EXTENDED data type routines
|
|
****************************************************************************}
|
|
|
|
{$define FPC_SYSTEM_HAS_PI}
|
|
function pi : double;[internproc:in_pi];
|
|
|
|
{$define FPC_SYSTEM_HAS_ABS}
|
|
function abs(d : extended) : extended;[internproc:in_abs_extended];
|
|
|
|
{$define FPC_SYSTEM_HAS_SQR}
|
|
function sqr(d : extended) : extended;[internproc:in_sqr_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;[internconst:in_const_exp];
|
|
begin
|
|
runerror(207);
|
|
end;
|
|
|
|
}
|
|
|
|
const
|
|
factor: double = double(int64(1) shl 32);
|
|
factor2: double = double(int64(1) shl 31);
|
|
|
|
{$define FPC_SYSTEM_HAS_TRUNC}
|
|
function trunc(d : extended) : int64;assembler;[internconst:in_const_trunc];
|
|
{ input: d in fr1 }
|
|
{ output: result in r3 }
|
|
assembler;
|
|
var
|
|
temp: packed record
|
|
case byte of
|
|
0: (l1,l2: longint);
|
|
1: (d: double);
|
|
end;
|
|
asm
|
|
// store d in temp
|
|
stfd f1,temp
|
|
// extract sign bit (record in cr0)
|
|
lwz r3,temp
|
|
rlwinm. r3,r3,1,31,31
|
|
// make d positive
|
|
fabs f1,f1
|
|
// load 2^32 in f2
|
|
{$ifndef macos}
|
|
lis r4,factor@ha
|
|
lfd f2,factor@l(r4)
|
|
{$else}
|
|
lwz r4,factor[TC](r2)
|
|
lfd f2,0(r4)
|
|
{$endif}
|
|
// check if value is < 0
|
|
// f3 := d / 2^32;
|
|
fdiv f3,f1,f2
|
|
// round
|
|
fctiwz f4,f3
|
|
// store
|
|
stfd f4,temp
|
|
// and load into r4
|
|
lwz r3,temp+4
|
|
// convert back to float
|
|
lis r0,0x4330
|
|
stw r0,temp
|
|
xoris r0,r3,0x8000
|
|
stw r0,temp+4
|
|
{$ifndef macos}
|
|
lis r4,longint_to_real_helper@ha
|
|
lfd f0,longint_to_real_helper@l(r4)
|
|
{$else}
|
|
lwz r4,longint_to_real_helper[TC](r2)
|
|
lfd f0,0(r4)
|
|
{$endif}
|
|
lfd f3,temp
|
|
fsub f3,f3,f0
|
|
|
|
|
|
// f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
|
|
fnmsub f4,f3,f2,f1
|
|
|
|
// now, convert to unsigned 32 bit
|
|
|
|
// load 2^31 in f2
|
|
{$ifndef macos}
|
|
lis r4,factor2@ha
|
|
lfd f2,factor2@l(r4)
|
|
{$else}
|
|
lwz r4,factor2[TC](r2)
|
|
lfd f2,0(r4)
|
|
{$endif}
|
|
|
|
// subtract 2^31
|
|
fsub f3,f4,f2
|
|
// was the value > 2^31?
|
|
fcmpu cr1,f4,f2
|
|
// use diff if >= 2^31
|
|
fsel f4,f3,f3,f4
|
|
|
|
// next part same as conversion to signed integer word
|
|
fctiwz f4,f4
|
|
stfd f4,temp
|
|
lwz r4,temp+4
|
|
// add 2^31 if value was >=2^31
|
|
blt cr1, .LTruncNoAdd
|
|
xoris r4,r4,0x8000
|
|
.LTruncNoAdd:
|
|
// negate value if it was negative to start with
|
|
beq cr0,.LTruncPositive
|
|
subfic r4,r4,0
|
|
subfze r3,r3
|
|
.LTruncPositive:
|
|
end;
|
|
|
|
|
|
{$define FPC_SYSTEM_HAS_ROUND}
|
|
{$ifdef hascompilerproc}
|
|
function round(d : extended) : int64;[internconst:in_const_round, external name 'FPC_ROUND'];
|
|
|
|
function fpc_round(d : extended) : int64;assembler;[public, alias:'FPC_ROUND'];{$ifdef hascompilerproc}compilerproc;{$endif hascompilerproc}
|
|
{$else}
|
|
function round(d : extended) : int64;assembler;[internconst:in_const_round];
|
|
{$endif hascompilerproc}
|
|
{ exactly the same as trunc, except that one fctiwz has become fctiw }
|
|
{ input: d in fr1 }
|
|
{ output: result in r3 }
|
|
assembler;
|
|
var
|
|
temp: packed record
|
|
case byte of
|
|
0: (l1,l2: longint);
|
|
1: (d: double);
|
|
end;
|
|
asm
|
|
// store d in temp
|
|
stfd f1, temp
|
|
// extract sign bit (record in cr0)
|
|
lwz r4,temp
|
|
rlwinm. r4,r4,1,31,31
|
|
// make d positive
|
|
fabs f1,f1
|
|
// load 2^32 in f2
|
|
{$ifndef macos}
|
|
lis r4,factor@ha
|
|
lfd f2,factor@l(r4)
|
|
{$else}
|
|
lwz r4,factor[TC](r2)
|
|
lfd f2,0(r4)
|
|
{$endif}
|
|
// check if value is < 0
|
|
// f3 := d / 2^32;
|
|
fdiv f3,f1,f2
|
|
// round
|
|
fctiwz f4,f3
|
|
// store
|
|
stfd f4,temp
|
|
// and load into r4
|
|
lwz r3,temp+4
|
|
// convert back to float
|
|
lis r0,0x4330
|
|
stw r0,temp
|
|
xoris r0,r3,0x8000
|
|
stw r0,temp+4
|
|
{$ifndef macos}
|
|
lis r4,longint_to_real_helper@ha
|
|
lfd f0,longint_to_real_helper@l(r4)
|
|
{$else}
|
|
lwz r4,longint_to_real_helper[TC](r2)
|
|
lfd f0,0(r4)
|
|
{$endif}
|
|
lfd f3,temp
|
|
fsub f3,f3,f0
|
|
|
|
|
|
// f4 := d "mod" 2^32 ( = d - ((d / 2^32) * 2^32))
|
|
fnmsub f4,f3,f2,f1
|
|
|
|
// now, convert to unsigned 32 bit
|
|
|
|
// load 2^31 in f2
|
|
{$ifndef macos}
|
|
lis r4,factor2@ha
|
|
lfd f2,factor2@l(r4)
|
|
{$else}
|
|
lwz r4,factor2[TC](r2)
|
|
lfd f2,0(r4)
|
|
{$endif}
|
|
|
|
// subtract 2^31
|
|
fsub f3,f4,f2
|
|
// was the value > 2^31?
|
|
fcmpu cr1,f4,f2
|
|
// use diff if >= 2^31
|
|
fsel f4,f3,f3,f4
|
|
|
|
// next part same as conversion to signed integer word
|
|
fctiw f4,f4
|
|
stfd f4,temp
|
|
lwz r4,temp+4
|
|
// add 2^31 if value was >=2^31
|
|
blt cr1, .LRoundNoAdd
|
|
xoris r4,r4,0x8000
|
|
.LRoundNoAdd:
|
|
// negate value if it was negative to start with
|
|
beq cr0,.LRoundPositive
|
|
subfic r4,r4,0
|
|
subfze r3,r3
|
|
.LRoundPositive:
|
|
end;
|
|
|
|
|
|
{$define FPC_SYSTEM_HAS_POWER}
|
|
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
|
|
****************************************************************************}
|
|
|
|
{$define FPC_SYSTEM_HAS_POWER_INT64}
|
|
function power(bas,expo : int64) : int64;
|
|
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! }
|
|
{$define FPC_SYSTEM_HAS_REAL2DOUBLE}
|
|
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
|
|
****************************************************************************}
|
|
|
|
{$define FPC_SYSTEM_HAS_INT64_TO_DOUBLE}
|
|
function fpc_int64_to_double(i: int64): double; compilerproc;
|
|
assembler;
|
|
{ input: high(i) in r4, low(i) in r3 }
|
|
{ output: double(i) in f0 }
|
|
var
|
|
temp: packed record
|
|
case byte of
|
|
0: (l1,l2: cardinal);
|
|
1: (d: double);
|
|
end;
|
|
asm
|
|
lis r0,0x4330
|
|
stw r0,temp
|
|
xoris r3,r3,0x8000
|
|
stw r3,temp+4
|
|
{$ifndef macos}
|
|
lis r3,longint_to_real_helper@ha
|
|
lfd f1,longint_to_real_helper@l(r3)
|
|
{$else}
|
|
lwz r3,longint_to_real_helper[TC](r2)
|
|
lfd f1,0(r3)
|
|
{$endif}
|
|
lfd f0,temp
|
|
stw r4,temp+4
|
|
fsub f0,f0,f1
|
|
{$ifndef macos}
|
|
lis r4,cardinal_to_real_helper@ha
|
|
lfd f1,cardinal_to_real_helper@l(r4)
|
|
lis r4,int_to_real_factor@ha
|
|
lfd f3,temp
|
|
lfd f2,int_to_real_factor@l(r4)
|
|
{$else}
|
|
lwz r4,cardinal_to_real_helper[TC](r2)
|
|
lwz r3,int_to_real_factor[TC](r2)
|
|
lfd f3,temp
|
|
lfd f1,0(r4)
|
|
lfd f2,0(r3)
|
|
{$endif}
|
|
fsub f3,f3,f1
|
|
fmadd f1,f0,f2,f3
|
|
end;
|
|
|
|
|
|
{$define FPC_SYSTEM_HAS_QWORD_TO_DOUBLE}
|
|
function fpc_qword_to_double(q: qword): double; compilerproc;
|
|
assembler;
|
|
{ input: high(q) in r4, low(q) in r3 }
|
|
{ output: double(q) in f0 }
|
|
var
|
|
temp: packed record
|
|
case byte of
|
|
0: (l1,l2: cardinal);
|
|
1: (d: double);
|
|
end;
|
|
asm
|
|
lis r0,0x4330
|
|
stw r0,temp
|
|
stw r3,temp+4
|
|
lfd f0,temp
|
|
{$ifndef macos}
|
|
lis r3,cardinal_to_real_helper@ha
|
|
lfd f1,cardinal_to_real_helper@l(r3)
|
|
{$else}
|
|
lwz r3,longint_to_real_helper[TC](r2)
|
|
lfd f1,0(r3)
|
|
{$endif}
|
|
stw r4,temp+4
|
|
fsub f0,f0,f1
|
|
lfd f3,temp
|
|
{$ifndef macos}
|
|
lis r4,int_to_real_factor@ha
|
|
lfd f2,int_to_real_factor@l(r4)
|
|
{$else}
|
|
lwz r4,int_to_real_factor[TC](r2)
|
|
lfd f2,0(r4)
|
|
{$endif}
|
|
fsub f3,f3,f1
|
|
fmadd f1,f0,f2,f3
|
|
end;
|
|
|
|
|
|
{
|
|
$Log$
|
|
Revision 1.32 2003-12-07 19:55:37 jonas
|
|
- reverted previous patch, solved with the new assembler reader
|
|
(which didn't understand the new syntax)
|
|
|
|
Revision 1.30 2003/11/15 19:01:27 florian
|
|
* fixed rtl to work with the integrated fpc ppc assembler reader
|
|
|
|
Revision 1.29 2003/09/04 16:07:31 florian
|
|
* fixed qword_to_double conversion on powerpc
|
|
|
|
Revision 1.28 2003/09/03 14:09:37 florian
|
|
* arm fixes to the common rtl code
|
|
* some generic math code fixed
|
|
* ...
|
|
|
|
Revision 1.27 2003/08/08 22:02:05 olle
|
|
* small bugfix macos
|
|
|
|
Revision 1.26 2003/06/14 12:41:08 jonas
|
|
* fixed compilation problems (removed unnecessary modified registers
|
|
lists from procedures)
|
|
|
|
Revision 1.25 2003/05/31 20:22:06 jonas
|
|
* fixed 64 bit results of trunc and round
|
|
|
|
Revision 1.24 2003/05/30 23:56:41 florian
|
|
* fixed parameter passing for int64
|
|
|
|
Revision 1.23 2003/05/24 13:39:32 jonas
|
|
* fsqrt is an optional instruction in the ppc architecture and isn't
|
|
implemented by any current ppc afaik, so use the generic sqrt routine
|
|
instead (adapted so it works with compilerproc)
|
|
|
|
Revision 1.22 2003/05/16 16:04:33 jonas
|
|
* fixed round() (almost the same as trunc)
|
|
|
|
Revision 1.21 2003/05/11 18:09:45 jonas
|
|
* fixed qword and int64 to double conversion
|
|
|
|
Revision 1.20 2003/05/02 15:12:19 jonas
|
|
- removed empty ppc-specific frac()
|
|
+ added correct generic frac() implementation for doubles (translated
|
|
from glibc code)
|
|
|
|
Revision 1.19 2003/04/26 20:36:24 jonas
|
|
* trunc now also supports int64 (no NaN's etc though)
|
|
|
|
Revision 1.18 2003/04/26 17:20:16 florian
|
|
* fixed trunc, now it's working at least for longint range
|
|
|
|
Revision 1.17 2003/04/23 21:28:21 peter
|
|
* fpc_round added, needed for int64 currency
|
|
|
|
Revision 1.16 2003/01/16 11:29:11 olle
|
|
* changed access of globals to be indirect via TOC
|
|
|
|
Revision 1.15 2003/01/15 01:09:04 florian
|
|
* changed power(...) prototype to int64
|
|
|
|
Revision 1.14 2002/11/28 11:04:16 olle
|
|
* macos: refs to globals in asm adapted to macos
|
|
|
|
Revision 1.13 2002/10/21 18:08:28 jonas
|
|
* round has int64 instead of longint result
|
|
|
|
Revision 1.12 2002/09/08 13:00:21 jonas
|
|
* made pi an internproc instead of internconst
|
|
|
|
Revision 1.11 2002/09/07 16:01:26 peter
|
|
* old logs removed and tabs fixed
|
|
|
|
Revision 1.10 2002/08/18 22:11:10 florian
|
|
* fixed remaining assembler errors
|
|
|
|
Revision 1.9 2002/08/18 21:37:48 florian
|
|
* several errors in inline assembler fixed
|
|
|
|
Revision 1.8 2002/08/10 17:14:36 jonas
|
|
* various fixes, mostly changing the names of the modifies registers to
|
|
upper case since that seems to be required by the compiler
|
|
|
|
Revision 1.7 2002/07/31 16:58:12 jonas
|
|
* fixed conversion from int64/qword to double errors
|
|
|
|
Revision 1.6 2002/07/29 21:28:17 florian
|
|
* several fixes to get further with linux/ppc system unit compilation
|
|
|
|
Revision 1.5 2002/07/28 21:39:29 florian
|
|
* made abs a compiler proc if it is generic
|
|
|
|
Revision 1.4 2002/07/28 20:43:49 florian
|
|
* several fixes for linux/powerpc
|
|
* several fixes to MT
|
|
|
|
}
|