-
Notifications
You must be signed in to change notification settings - Fork 17
/
FastMath.pas
312 lines (254 loc) · 6.98 KB
/
FastMath.pas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
unit FastMath;
//note that this unit only contains optimized versions for X64.
//The other versions are already optimized.
interface
{ MinValue: Returns the smallest signed value in the data array (MIN) }
function MinValue(const Data: array of Single): Single; overload;
function MinValue(const Data: array of Double): Double; overload;
function MinValue(const Data: array of Extended): Extended; overload;
function MinIntValue(const Data: array of Integer): Integer;
function Min(const A, B: Integer): Integer; overload; {$ifdef PurePascal}inline;{$endif}
function Min(const A, B: Int64): Int64; overload; {$if defined(PurePascal)}inline;{$endif}
function Min(const A, B: UInt64): UInt64; overload; inline;
function Min(const A, B: Single): Single; overload; inline;
function Min(const A, B: Double): Double; overload; inline;
function Min(const A, B: Extended): Extended; overload; inline;
{ MaxValue: Returns the largest signed value in the data array (MAX) }
function MaxValue(const Data: array of Single): Single; overload;
function MaxValue(const Data: array of Double): Double; overload;
function MaxValue(const Data: array of Extended): Extended; overload;
function MaxIntValue(const Data: array of Integer): Integer;
function Max(const A, B: Integer): Integer; overload; inline;
function Max(const A, B: Int64): Int64; overload; inline;
function Max(const A, B: UInt64): UInt64; overload; inline;
function Max(const A, B: Single): Single; overload; inline;
function Max(const A, B: Double): Double; overload; inline;
function Max(const A, B: Extended): Extended; overload; inline;
implementation
uses Math;
{ MinValue: Returns the smallest signed value in the data array (MIN) }
function MinValue(const Data: array of Single): Single; overload;
{$ifdef cpux64}
asm
or rcx,rcx
jz @error
mov eax,[rcx-4] //get length
neg eax
sub rcx,rax
add eax,32
jns @tail
movups xmm0,[rcx+rax-32]
@loop:
minps xmm0,[rcx+rax-32+16]
add eax,16
js @loop
movaps xmm1,xmm0
shufps xmm1,xmm0,3+(2 shl 2)+(1 shl 4)+(0 shl 6)
minps xmm0,xmm1
shufps xmm1,xmm1,2+(3 shl 2)+(0 shl 4)+(1 shl 6)
minps xmm0,xmm1
movd eax,xmm0
@tail:
@error:
{$else}
begin
Result:= Math.MinValue(Data)
{$endif}
end;
function MinValue(const Data: array of Double): Double; overload;
begin
end;
function MinValue(const Data: array of Extended): Extended; overload;
begin
end;
function MinIntValue(const Data: array of Integer): Integer;
begin
end;
function Min(const A, B: Integer): Integer; overload;
{$if defined(PurePascal)}
inline;
begin
Result:= ((a-b)*(a>b))+b;
end;
{$elseif defined(CPUX86)}
asm
sub a,b
sbb ecx,ecx
and a,ecx
add a,b
end;
{$else}
asm
sub a,b
sbb eax,eax //RAX = 0 if a>b else RAX=-1
and eax,a //RAX = 0 if a>b else RAX=(a-b)
add eax,b //rax = a-b+b = a if b>a else rax =0+b = b
end;
{$endif}
function Min(const A, B: Int64): Int64; overload;
{$if defined(PurePascal)}
inline;
begin
Result:= ((a-b)*(a>b))+b;
end;
{$else}
{$if defined(CPUX86)}
asm
sub eax,edx
sbb ecx,ecx
and eax,ecx
add eax,edx
end;
{$else}
asm
sub a,b
sbb rax,rax //RAX = 0 if a>b else RAX=-1
and rax,a //RAX = 0 if a>b else RAX=(a-b)
add rax,b //rax = a-b+b = a if b>a else rax =0+b = b
end;
{$endif}
{$endif}
function Min(const A, B: UInt64): UInt64; overload; inline;
begin
end;
function Min(const A, B: Single): Single; overload; inline;
begin
end;
function Min(const A, B: Double): Double; overload; inline;
begin
end;
function Min(const A, B: Extended): Extended; overload; inline;
begin
end;
{ MaxValue: Returns the largest signed value in the data array (MAX) }
function MaxValue(const Data: array of Single): Single; overload;
begin
end;
function MaxValue(const Data: array of Double): Double; overload;
begin
end;
function MaxValue(const Data: array of Extended): Extended; overload;
begin
end;
function MaxIntValue(const Data: array of Integer): Integer;
begin
end;
function Max(const A, B: Integer): Integer; overload; inline;
begin
end;
function Max(const A, B: Int64): Int64; overload; inline;
begin
end;
function Max(const A, B: UInt64): UInt64; overload; inline;
begin
end;
function Max(const A, B: Single): Single; overload; inline;
begin
end;
function Max(const A, B: Double): Double; overload; inline;
begin
end;
function Max(const A, B: Extended): Extended; overload; inline;
begin
end;
//// func Abs(x float64) float64
//TEXT ·Abs(SB),NOSPLIT,$0
// MOVQ $(1<<63), BX
// MOVQ BX, X0 // movsd $(-0.0), x0
// MOVSD x+0(FP), X1
// ANDNPD X1, X0
// MOVSD X0, ret+8(FP)
// RET
//
//// func Frexp(f float64) (frac float64, exp int)
//TEXT ·Frexp(SB),NOSPLIT,$0
// FMOVD f+0(FP), F0 // F0=f
// FXAM
// FSTSW AX
// SAHF
// JNP nan_zero_inf
// JCS nan_zero_inf
// FXTRACT // F0=f (0<=f<1), F1=e
// FMULD $(0.5), F0 // F0=f (0.5<=f<1), F1=e
// FMOVDP F0, frac+8(FP) // F0=e
// FLD1 // F0=1, F1=e
// FADDDP F0, F1 // F0=e+1
// FMOVLP F0, exp+16(FP) // (int=int32)
// RET
//nan_zero_inf:
// FMOVDP F0, frac+8(FP) // F0=e
// MOVL $0, exp+16(FP) // exp=0
// RET
//
//
//#define Big 0x4330000000000000 // 2**52
//
//// func Floor(x float64) float64
//TEXT ·Floor(SB),NOSPLIT,$0
// MOVQ x+0(FP), AX
// MOVQ $~(1<<63), DX // sign bit mask
// ANDQ AX,DX // DX = |x|
// SUBQ $1,DX
// MOVQ $(Big - 1), CX // if |x| >= 2**52-1 or IsNaN(x) or |x| == 0, return x
// CMPQ DX,CX
// JAE isBig_floor
// MOVQ AX, X0 // X0 = x
// CVTTSD2SQ X0, AX
// CVTSQ2SD AX, X1 // X1 = float(int(x))
// CMPSD X1, X0, 1 // compare LT; X0 = 0xffffffffffffffff or 0
// MOVSD $(-1.0), X2
// ANDPD X2, X0 // if x < float(int(x)) {X0 = -1} else {X0 = 0}
// ADDSD X1, X0
// MOVSD X0, ret+8(FP)
// RET
//isBig_floor:
// MOVQ AX, ret+8(FP) // return x
// RET
//
//// func Ceil(x float64) float64
//TEXT ·Ceil(SB),NOSPLIT,$0
// MOVQ x+0(FP), AX
// MOVQ $~(1<<63), DX // sign bit mask
// MOVQ AX, BX // BX = copy of x
// ANDQ DX, BX // BX = |x|
// MOVQ $Big, CX // if |x| >= 2**52 or IsNaN(x), return x
// CMPQ BX, CX
// JAE isBig_ceil
// MOVQ AX, X0 // X0 = x
// MOVQ DX, X2 // X2 = sign bit mask
// CVTTSD2SQ X0, AX
// ANDNPD X0, X2 // X2 = sign
// CVTSQ2SD AX, X1 // X1 = float(int(x))
// CMPSD X1, X0, 2 // compare LE; X0 = 0xffffffffffffffff or 0
// ORPD X2, X1 // if X1 = 0.0, incorporate sign
// MOVSD $1.0, X3
// ANDNPD X3, X0
// ORPD X2, X0 // if float(int(x)) <= x {X0 = 1} else {X0 = -0}
// ADDSD X1, X0
// MOVSD X0, ret+8(FP)
// RET
//isBig_ceil:
// MOVQ AX, ret+8(FP)
// RET
//
//// func Trunc(x float64) float64
//TEXT ·Trunc(SB),NOSPLIT,$0
// MOVQ x+0(FP), AX
// MOVQ $~(1<<63), DX // sign bit mask
// MOVQ AX, BX // BX = copy of x
// ANDQ DX, BX // BX = |x|
// MOVQ $Big, CX // if |x| >= 2**52 or IsNaN(x), return x
// CMPQ BX, CX
// JAE isBig_trunc
// MOVQ AX, X0
// MOVQ DX, X2 // X2 = sign bit mask
// CVTTSD2SQ X0, AX
// ANDNPD X0, X2 // X2 = sign
// CVTSQ2SD AX, X0 // X0 = float(int(x))
// ORPD X2, X0 // if X0 = 0.0, incorporate sign
// MOVSD X0, ret+8(FP)
// RET
//isBig_trunc:
// MOVQ AX, ret+8(FP) // return x
// RET
end.