1{
2    Implementation of mathematical routines for x86_64
3
4    This file is part of the Free Pascal run time library.
5    Copyright (c) 1999-2005 by the Free Pascal development team
6
7    See the file COPYING.FPC, included in this distribution,
8    for details about the copyright.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13
14 **********************************************************************}
15{-------------------------------------------------------------------------
16 Using functions from AMath/DAMath libraries, which are covered by the
17 following license:
18
19 (C) Copyright 2009-2013 Wolfgang Ehrhardt
20
21 This software is provided 'as-is', without any express or implied warranty.
22 In no event will the authors be held liable for any damages arising from
23 the use of this software.
24
25 Permission is granted to anyone to use this software for any purpose,
26 including commercial applications, and to alter it and redistribute it
27 freely, subject to the following restrictions:
28
29 1. The origin of this software must not be misrepresented; you must not
30    claim that you wrote the original software. If you use this software in
31    a product, an acknowledgment in the product documentation would be
32    appreciated but is not required.
33
34 2. Altered source versions must be plainly marked as such, and must not be
35    misrepresented as being the original software.
36
37 3. This notice may not be removed or altered from any source distribution.
38----------------------------------------------------------------------------}
39
40{$push}
41{$codealign constmin=16}
42const
43  FPC_ABSMASK_SINGLE: array[0..1] of qword=($7fffffff7fffffff,$7fffffff7fffffff); cvar; public;
44  FPC_ABSMASK_DOUBLE: array[0..1] of qword=($7fffffffffffffff,$7fffffffffffffff); cvar; public;
45{$pop}
46
47{****************************************************************************
48                            FPU Control word
49 ****************************************************************************}
50
51    procedure Set8087CW(cw:word);
52      begin
53         default8087cw:=cw;
54         asm
55           fnclex
56           fldcw cw
57         end;
58      end;
59
60
61    function Get8087CW:word;assembler;
62      var
63        tmp: word;
64      asm
65        fnstcw tmp
66        movw   tmp,%ax
67        andl   $0xffff,%eax  { clears bits 32-63 }
68      end;
69
70
71    procedure SetMXCSR(w : dword);
72      begin
73        defaultmxcsr:=w;
74        asm
75          ldmxcsr w
76        end;
77      end;
78
79
80    function GetMXCSR : dword;assembler;
81      var
82        _w : dword;
83      asm
84        stmxcsr _w
85        movl    _w,%eax
86      end;
87
88
89    procedure SetSSECSR(w : dword);
90      begin
91        SetMXCSR(w);
92      end;
93
94
95    function GetSSECSR: dword;
96      begin
97        result:=GetMXCSR;
98      end;
99
100{****************************************************************************
101                       EXTENDED data type routines
102 ****************************************************************************}
103
104    {$ifndef FPC_SYSTEM_HAS_ABS}
105    {$define FPC_SYSTEM_HAS_ABS}
106    function fpc_abs_real(d : ValReal) : ValReal;compilerproc;
107    begin
108      { Function is handled internal in the compiler }
109      runerror(207);
110      result:=0;
111    end;
112    {$endif FPC_SYSTEM_HAS_ABS}
113    {$ifndef FPC_SYSTEM_HAS_SQR}
114    {$define FPC_SYSTEM_HAS_SQR}
115    function fpc_sqr_real(d : ValReal) : ValReal;compilerproc;
116    begin
117      { Function is handled internal in the compiler }
118      runerror(207);
119      result:=0;
120    end;
121    {$endif FPC_SYSTEM_HAS_SQR}
122    {$ifndef FPC_SYSTEM_HAS_SQRT}
123    {$define FPC_SYSTEM_HAS_SQRT}
124    function fpc_sqrt_real(d : ValReal) : ValReal;compilerproc;
125{$ifndef cpullvm}
126    begin
127      { Function is handled internal in the compiler }
128      runerror(207);
129      result:=0;
130{$else not cpullvm}
131    assembler;
132    asm
133      fldt d
134      fsqrt
135{$endif not cpullvm}
136    end;
137    {$endif FPC_SYSTEM_HAS_SQRT}
138
139{$ifdef FPC_HAS_TYPE_EXTENDED}
140
141    {$ifndef FPC_SYSTEM_HAS_ARCTAN}
142    {$define FPC_SYSTEM_HAS_ARCTAN}
143    function fpc_arctan_real(d : ValReal) : ValReal;compilerproc;
144{$ifndef cpullvm}
145    begin
146      { Function is handled internal in the compiler }
147      runerror(207);
148      result:=0;
149{$else not cpullvm}
150    assembler;
151    asm
152      fldt d
153      fld1
154      fpatan
155{$endif not cpullvm}
156    end;
157    {$endif FPC_SYSTEM_HAS_ARCTAN}
158    {$ifndef FPC_SYSTEM_HAS_LN}
159    {$define FPC_SYSTEM_HAS_LN}
160    function fpc_ln_real(d : ValReal) : ValReal;compilerproc;
161{$ifndef cpullvm}
162    begin
163      { Function is handled internal in the compiler }
164      runerror(207);
165      result:=0;
166{$else not cpullvm}
167    assembler;
168    asm
169      fldln2
170      fldt d
171      fyl2x
172{$endif not cpullvm}
173    end;
174    {$endif FPC_SYSTEM_HAS_LN}
175    {$ifndef FPC_SYSTEM_HAS_SIN}
176    {$define FPC_SYSTEM_HAS_SIN}
177    function fpc_sin_real(d : ValReal) : ValReal;compilerproc;
178{$ifndef cpullvm}
179    begin
180      { Function is handled internal in the compiler }
181      runerror(207);
182      result:=0;
183{$else not cpullvm}
184    assembler;
185    asm
186      fldt d
187      fsin
188{$endif not cpullvm}
189    end;
190    {$endif FPC_SYSTEM_HAS_SIN}
191    {$ifndef FPC_SYSTEM_HAS_COS}
192    {$define FPC_SYSTEM_HAS_COS}
193    function fpc_cos_real(d : ValReal) : ValReal;compilerproc;
194{$ifndef cpullvm}
195    begin
196      { Function is handled internal in the compiler }
197      runerror(207);
198      result:=0;
199{$else not cpullvm}
200    assembler;
201    asm
202      fldt d
203      fcos
204{$endif not cpullvm}
205    end;
206    {$endif FPC_SYSTEM_HAS_COS}
207
208    {$ifndef FPC_SYSTEM_HAS_EXP}
209    {$define FPC_SYSTEM_HAS_EXP}
210    { exp function adapted from AMath library (C) Copyright 2009-2013 Wolfgang Ehrhardt
211      * translated into AT&T syntax
212      + PIC support
213      * return +Inf/0 for +Inf/-Inf input, instead of NaN }
214    function fpc_exp_real(d : ValReal) : ValReal;assembler;compilerproc;
215      const
216        ln2hi: double=6.9314718036912382E-001;
217        ln2lo: double=1.9082149292705877E-010;
218        large: single=24576.0;
219        two:   single=2.0;
220        half:  single=0.5;
221      asm
222            fldt        d
223            fldl2e
224            fmul        %st(1),%st        { z = d * log2(e) }
225            frndint
226          { Calculate frac(z) using modular arithmetic to avoid precision loss }
227            fldl        ln2hi(%rip)
228            fmul        %st(1),%st
229            fsubrp      %st,%st(2)
230            fldl        ln2lo(%rip)
231            fmul        %st(1),%st
232            fsubrp      %st,%st(2)
233            fxch        %st(1)            { (d-int(z)*ln2_hi)-int(z)*ln2_lo }
234            fldl2e
235            fmulp       %st,%st(1)        { frac(z) }
236
237          { Above calculation can yield |frac(z)|>1, particularly when rounding mode
238            is not "round to nearest". f2xm1 is undefined in that case, so it's
239            necessary to check }
240            fld         %st
241            fabs
242            fld1
243            fcomip      %st(1),%st(0)
244            fstp        %st
245            jp          .L3               { NaN }
246            jae         .L1               { |frac(z)| <= 1, good }
247
248            fld         %st(1)
249            fabs
250            flds        large(%rip)
251            fcomip      %st(1),%st(0)
252            fstp        %st
253            jb          .L3               { int(z) >= 24576 }
254   .L0:
255            { Calculate 2**frac(z)-1 as N*(N+2), where N=2**(frac(z)/2)-1 }
256            fmuls       half(%rip)
257            f2xm1
258            fld         %st
259            fadds       two(%rip)
260            fmulp       %st,%st(1)
261            jmp         .L2
262   .L3:
263            fstp        %st               { pop frac(z) and load 0 }
264            fldz
265   .L1:
266            f2xm1
267   .L2:
268            fld1
269            faddp       %st,%st(1)
270            fscale
271            fstp        %st(1)
272      end;
273    {$endif FPC_SYSTEM_HAS_EXP}
274
275
276    {$ifndef FPC_SYSTEM_HAS_FRAC}
277    {$define FPC_SYSTEM_HAS_FRAC}
278    function fpc_frac_real(d : ValReal) : ValReal;assembler;compilerproc;
279      var
280        oldcw,newcw: word;
281      asm
282            fnstcw oldcw
283            fwait
284            movw oldcw,%cx
285            orw $0x0c3f,%cx
286            movw %cx,newcw
287            fldcw newcw
288            fwait
289            fldt d
290            frndint
291            fldt d
292            fsub %st(1),%st
293            fstp %st(1)
294            fnclex
295            fldcw oldcw
296      end;
297    {$endif FPC_SYSTEM_HAS_FRAC}
298
299
300    {$ifndef FPC_SYSTEM_HAS_INT}
301    {$define FPC_SYSTEM_HAS_INT}
302    function fpc_int_real(d : ValReal) : ValReal;assembler;compilerproc;
303      var
304        oldcw,newcw: word;
305      asm
306            fnstcw oldcw
307            fwait
308            movw oldcw,%cx
309            orw $0x0c3f,%cx
310            movw %cx,newcw
311            fldcw newcw
312            fwait
313            fldt d
314            frndint
315            fwait
316            fldcw oldcw
317      end;
318    {$endif FPC_SYSTEM_HAS_INT}
319
320
321    {$ifndef FPC_SYSTEM_HAS_TRUNC}
322    {$define FPC_SYSTEM_HAS_TRUNC}
323    function fpc_trunc_real(d : ValReal) : int64;assembler;compilerproc;
324      var
325        oldcw,
326        newcw : word;
327        res   : int64;
328      asm
329        fnstcw oldcw
330        fwait
331        movw oldcw,%cx
332        orw $0x0c3f,%cx
333        movw %cx,newcw
334        fldcw newcw
335        fldt d
336        fistpq res
337        fwait
338        movq res,%rax
339        fldcw oldcw
340      end;
341    {$endif FPC_SYSTEM_HAS_TRUNC}
342
343
344    {$ifndef FPC_SYSTEM_HAS_ROUND}
345    {$define FPC_SYSTEM_HAS_ROUND}
346    function fpc_round_real(d : ValReal) : int64;assembler;compilerproc;
347      var
348        res   : int64;
349      asm
350        fldt d
351        fistpq res
352        fwait
353        movq res,%rax
354      end;
355    {$endif FPC_SYSTEM_HAS_ROUND}
356
357{$else FPC_HAS_TYPE_EXTENDED}
358
359    {$ifndef FPC_SYSTEM_HAS_INT}
360    {$define FPC_SYSTEM_HAS_INT}
361    function fpc_int_real(d : ValReal) : ValReal;compilerproc; assembler; nostackframe;
362      asm
363        movq      %xmm0,  %rax
364        shr       $48,    %rax
365        and       $0x7ff0,%ax
366        cmp       $0x4330,%ax
367        jge       .L0
368        cvttsd2si %xmm0,  %rax
369        cvtsi2sd  %rax,   %xmm0
370    .L0:
371      end;
372    {$endif FPC_SYSTEM_HAS_INT}
373
374    {$ifndef FPC_SYSTEM_HAS_TRUNC}
375    {$define FPC_SYSTEM_HAS_TRUNC}
376    function fpc_trunc_real(d : ValReal) : int64;compilerproc; assembler; nostackframe;
377      asm
378        cvttsd2si   %xmm0,%rax;
379      end;
380    {$endif FPC_SYSTEM_HAS_TRUNC}
381
382    {$ifndef FPC_SYSTEM_HAS_ROUND}
383    {$define FPC_SYSTEM_HAS_ROUND}
384    function fpc_round_real(d : ValReal) : int64;compilerproc; assembler; nostackframe;
385      asm
386        cvtsd2si   %xmm0,%rax;
387      end;
388    {$endif FPC_SYSTEM_HAS_ROUND}
389
390    {$ifndef FPC_SYSTEM_HAS_FRAC}
391    {$define FPC_SYSTEM_HAS_FRAC}
392    function fpc_frac_real(d: ValReal) : ValReal;compilerproc; assembler; nostackframe;
393      asm
394        { Windows defines %xmm4 and %xmm5 as first non-parameter volatile registers;
395          on SYSV systems all are considered as such, so use %xmm4 }
396        movq      %xmm0,  %rax
397        movapd    %xmm0,  %xmm4
398        shr       $48,    %rax
399        and       $0x7ff0,%ax
400        cmp       $0x4330,%ax
401        jge       .L0
402        cvttsd2si %xmm0,  %rax
403        cvtsi2sd  %rax,   %xmm4
404  .L0:
405        subsd     %xmm4,  %xmm0
406      end;
407    {$endif FPC_SYSTEM_HAS_FRAC}
408
409{$endif FPC_HAS_TYPE_EXTENDED}
410