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
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
|
; $Id: powcore.asm $
;; @file
; IPRT - No-CRT common pow code - AMD64 & X86.
;
;
; Copyright (C) 2006-2023 Oracle and/or its affiliates.
;
; This file is part of VirtualBox base platform packages, as
; available from https://www.virtualbox.org.
;
; This program is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public License
; as published by the Free Software Foundation, in version 3 of the
; License.
;
; 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. See the GNU
; General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with this program; if not, see <https://www.gnu.org/licenses>.
;
; The contents of this file may alternatively be used under the terms
; of the Common Development and Distribution License Version 1.0
; (CDDL), a copy of it is provided in the "COPYING.CDDL" file included
; in the VirtualBox distribution, in which case the provisions of the
; CDDL are applicable instead of those of the GPL.
;
; You may elect to license modified versions of this file under the
; terms and conditions of either the GPL or the CDDL or both.
;
; SPDX-License-Identifier: GPL-3.0-only OR CDDL-1.0
;
%define RT_ASM_WITH_SEH64
%include "iprt/asmdefs.mac"
%include "iprt/x86.mac"
BEGINCODE
extern NAME(RT_NOCRT(feraiseexcept))
;;
; Call feraiseexcept(%1)
%macro CALL_feraiseexcept_WITH 1
%ifdef RT_ARCH_X86
mov dword [xSP], X86_FSW_IE
%elifdef ASM_CALL64_GCC
mov edi, X86_FSW_IE
%elifdef ASM_CALL64_MSC
mov ecx, X86_FSW_IE
%else
%error calling conv.
%endif
call NAME(RT_NOCRT(feraiseexcept))
%endmacro
;;
; Compute the st1 to the power of st0.
;
; @returns st(0) = result
; eax = what's being returned:
; 0 - Just a value.
; 1 - The rBase value. Caller may take steps to ensure it's exactly the same.
; 2 - The rExp value. Caller may take steps to ensure it's exactly the same.
; @param rBase/st1 The base.
; @param rExp/st0 The exponent
; @param fFxamBase/dx The status flags after fxam(rBase).
; @param enmType/ebx The original parameter and return types:
; 0 - 32-bit / float
; 1 - 64-bit / double
; 2 - 80-bit / long double
;
BEGINPROC rtNoCrtMathPowCore
push xBP
SEH64_PUSH_xBP
mov xBP, xSP
SEH64_SET_FRAME_xBP 0
sub xSP, 30h
SEH64_ALLOCATE_STACK 30h
SEH64_END_PROLOGUE
;
; Weed out special values, starting with the exponent.
;
fxam
fnstsw ax
mov cx, ax ; cx=fxam(exp)
and ax, X86_FSW_C3 | X86_FSW_C2 | X86_FSW_C0
cmp ax, X86_FSW_C2 ; Normal finite number (excluding zero)
je .exp_finite
cmp ax, X86_FSW_C3 ; Zero
je .exp_zero
cmp ax, X86_FSW_C3 | X86_FSW_C2 ; Denormals
je .exp_finite
cmp ax, X86_FSW_C0 | X86_FSW_C2 ; Infinity.
je .exp_inf
jmp .exp_nan
.exp_finite:
;
; Detect special base values.
;
mov ax, dx ; ax=fxam(base)
and ax, X86_FSW_C3 | X86_FSW_C2 | X86_FSW_C0
cmp ax, X86_FSW_C2 ; Normal finite number (excluding zero)
je .base_finite
cmp ax, X86_FSW_C3 ; Zero
je .base_zero
cmp ax, X86_FSW_C3 | X86_FSW_C2 ; Denormals
je .base_finite
cmp ax, X86_FSW_C0 | X86_FSW_C2 ; Infinity.
je .base_inf
jmp .base_nan
.base_finite:
;
; 1 in the base is also special.
; Rule 6 (see below): base == +1 and exponent = whatever: Return +1.0
;
fld1
fcomip st0, st2
je .return_base_value
;
; Check if the exponent is an integer value we can handle in a 64-bit
; GRP as that is simpler to handle accurately.
;
; In 64-bit integer range?
fld tword [.s_r80MaxInt xWrtRIP]
fcomip st0, st1
jb .not_integer_exp
fld tword [.s_r80MinInt xWrtRIP]
fcomip st0, st1
ja .not_integer_exp
; Convert it to integer.
fld st0 ; -> st0=exp; st1=exp; st2=base
fistp qword [xBP - 8] ; Save and pop 64-bit int (no non-popping version of this instruction).
fild qword [xBP - 8] ; Load it again for comparison.
fucomip st0, st1 ; Compare integer exp and floating point exp to see if they are the same. Pop.
jne .not_integer_exp
;
;
; Ok, we've got an integer exponent value in that fits into a 64-bit.
; We'll multiply the base exponention bit by exponention bit, applying
; it as a factor for bits that are set.
;
;
.integer_exp:
; Load the integer value into edx:exx / rdx and ditch the floating point exponent.
mov xDX, [xBP - 8]
%ifdef RT_ARCH_X86
mov eax, [xBP - 8 + 4]
%endif
ffreep st0 ; -> st0=base;
; Load a 1 onto the stack, we'll need it below as well as for converting
; a negative exponent to a positive one.
fld1 ; -> st0=1.0; st1=base;
; If the exponent is negative, negate it and change base to 1/base.
or xDX, xDX
jns .integer_exp_positive
neg xDX
%ifdef RT_ARCH_X86
neg eax
sbb edx, 0
%endif
fdivr st1, st0 ; -> st0=1.0; st1=1/base
.integer_exp_positive:
;
; We'll process edx:eax / rdx bit by bit till it's zero, using st0 for
; the multiplication factor corresponding to the current exponent bit
; and st1 as the result.
;
fxch ; -> st0=base; st1=1.0;
.integer_exp_loop:
%ifdef RT_ARCH_X86
shrd eax, edx, 1
%else
shr rdx, 1
%endif
jnc .integer_exp_loop_advance
fmul st1, st0
.integer_exp_loop_advance:
; Check if we're done.
%ifdef RT_ARCH_AMD64
jz .integer_exp_return ; (we will have the flags for the shr rdx above)
%else
shr edx, 1 ; complete the above shift operation
mov ecx, edx ; check if edx:eax is zero.
or ecx, eax
jz .integer_exp_return
%endif
; Calculate the factor for the next bit.
fmul st0, st0
jmp .integer_exp_loop
.integer_exp_return:
ffreep st0 ; drop the factor -> st0=result; no st1.
jmp .return_val
;
;
; Non-integer or value was out of range for an int64_t.
;
; The approach here is the same as in exp.asm, only we have to do the
; log2(base) calculation first as it's a parameter and not a constant.
;
;
.not_integer_exp:
; First reject negative numbers. We still have the fxam(base) status in dx.
test dx, X86_FSW_C1
jnz .base_negative_non_integer_exp
; Swap the items on the stack, so we can process the base first.
fxch st0, st1 ; -> st0=base; st1=exponent;
;
; From log2.asm:
;
; The fyl2xp1 instruction (ST1=ST1*log2(ST0+1.0), popping ST0) has a
; valid ST0 range of 1(1-sqrt(0.5)) (approx 0.29289321881) on both
; sides of zero. We try use it if we can.
;
.above_one:
; For both fyl2xp1 and fyl2xp1 we need st1=1.0.
fld1
fxch st0, st1 ; -> st0=base; st1=1.0; st2=exponent
; Check if the input is within the fyl2xp1 range.
fld qword [.s_r64AbsFyL2xP1InputMax xWrtRIP]
fcomip st0, st1
jbe .cannot_use_fyl2xp1
fld qword [.s_r64AbsFyL2xP1InputMin xWrtRIP]
fcomip st0, st1
jae .cannot_use_fyl2xp1
; Do the calculation.
.use_fyl2xp1:
fsub st0, st1 ; -> st0=base-1; st1=1.0; st2=exponent
fyl2xp1 ; -> st0=1.0*log2(base-1.0+1.0); st1=exponent
jmp .done_log2
.cannot_use_fyl2xp1:
fyl2x ; -> st0=1.0*log2(base); st1=exponent
.done_log2:
;
; From exp.asm:
;
; Convert to power of 2 and it'll be the same as exp2.
;
fmulp ; st0=log2(base); st1=exponent -> st0=pow2exp
;
; Split the job in two on the fraction and integer l2base parts.
;
fld st0 ; Push a copy of the pow2exp on the stack.
frndint ; st0 = (int)pow2exp
fsub st1, st0 ; st1 = pow2exp - (int)pow2exp; i.e. st1 = fraction, st0 = integer.
fxch ; st0 = fraction, st1 = integer.
; 1. Calculate on the fraction.
f2xm1 ; st0 = 2**fraction - 1.0
fld1
faddp ; st0 = 2**fraction
; 2. Apply the integer power of two.
fscale ; st0 = result; st1 = integer part of pow2exp.
fstp st1 ; st0 = result; no st1.
;
; Return st0.
;
.return_val:
xor eax, eax
.return:
leave
ret
;
;
; pow() has a lot of defined behavior for special values, which is why
; this is the largest and most difficult part of the code. :-)
;
; On https://pubs.opengroup.org/onlinepubs/9699919799/functions/pow.html
; there are 21 error conditions listed in the return value section.
; The code below refers to this by number.
;
; When we get here:
; dx=fxam(base)
; cx=fxam(exponent)
; st1=base
; st0=exponent
;
;
; 1. Finit base < 0 and finit non-interger exponent: -> domain error (#IE) + NaN.
;
; The non-integer exponent claim might be wrong, as we only check if it
; fits into a int64_t register. But, I don't see how we can calculate
; it right now.
;
.base_negative_non_integer_exp:
CALL_feraiseexcept_WITH X86_FSW_IE
jmp .return_nan
;
; 7. Exponent = +/-0.0, any base value including NaN: return +1.0
; Note! According to https://en.cppreference.com/w/c/numeric/math/pow a
; domain error (#IE) occur if base=+/-0. Not implemented.
.exp_zero:
.return_plus_one:
fld1
jmp .return_pop_pop_val
;
; 6. Exponent = whatever and base = 1: Return 1.0
; 10. Exponent = +/-Inf and base = -1: Return 1.0
;6+10 => Exponent = +/-Inf and |base| = 1: Return 1.0
; 11. Exponent = -Inf and |base| < 1: Return +Inf
; 12. Exponent = -Inf and |base| > 1: Return +0
; 13. Exponent = +Inf and |base| < 1: Return +0
; 14. Exponent = +Inf and |base| > 1: Return +Inf
;
; Note! Rule 4 would trigger for the same conditions as 11 when base == 0,
; but it's optional to raise div/0 and it's apparently marked as
; obsolete in C23, so not implemented.
;
.exp_inf:
; Check if base is NaN or unsupported.
and dx, X86_FSW_C3 | X86_FSW_C2 | X86_FSW_C0 ; fxam(base)
cmp dx, X86_FSW_C0
jbe .return_base_nan
; Calc fabs(base) and replace the exponent with 1.0 as we're very likely to need this here.
ffreep st0
fabs
fld1 ; st0=1.0; st1=|rdBase|
fcomi st0, st1
je .return_plus_one ; Matches rule 6 + 10 (base is +/-1).
ja .exp_inf_base_smaller_than_one
.exp_inf_base_larger_than_one:
test cx, X86_FSW_C1 ; cx=faxm(exponent); C1=sign
jz .return_plus_inf ; Matches rule 14 (exponent is +Inf).
jmp .return_plus_zero ; Matches rule 12 (exponent is -Inf).
.exp_inf_base_smaller_than_one:
test cx, X86_FSW_C1 ; cx=faxm(exponent); C1=sign
jnz .return_plus_inf ; Matches rule 11 (exponent is -Inf).
jmp .return_plus_zero ; Matches rule 13 (exponent is +Inf).
;
; 6. Exponent = whatever and base = 1: Return 1.0
; 5. Unless specified elsewhere, return NaN if any of the parameters are NaN.
;
.exp_nan:
; Check if base is a number and possible 1.
test dx, X86_FSW_C2 ; dx=fxam(base); C2 is set for finite number, infinity and denormals.
jz .return_exp_nan
fld1
fcomip st0, st2
jne .return_exp_nan
jmp .return_plus_one
;
; 4a. base == +/-0.0 and exp < 0 and exp is odd integer: Return +/-Inf, raise div/0.
; 4b. base == +/-0.0 and exp < 0 and exp is not odd int: Return +Inf, raise div/0.
; 8. base == +/-0.0 and exp > 0 and exp is odd integer: Return +/-0.0
; 9. base == +/-0.0 and exp > 0 and exp is not odd int: Return +0
;
; Note! Exponent must be finite and non-zero if we get here.
;
.base_zero:
fldz
fcomip st0, st1
jbe .base_zero_plus_exp
.base_zero_minus_exp:
mov cx, dx ; stashing fxam(base) in CX because EDX is trashed by .is_exp_odd_integer
call .is_exp_odd_integer ; trashes EDX but no ECX.
or eax, eax
jz .base_zero_minus_exp_not_odd_int
; Matching 4a.
.base_zero_minus_exp_odd_int:
test cx, X86_FSW_C1 ; base sign
jz .raise_de_and_return_plus_inf
.raise_de_and_return_minus_inf:
CALL_feraiseexcept_WITH X86_FSW_DE
jmp .return_minus_inf
.raise_de_and_return_plus_inf:
CALL_feraiseexcept_WITH X86_FSW_DE
jmp .return_plus_inf
; Matching 4b.
.base_zero_minus_exp_not_odd_int:
CALL_feraiseexcept_WITH X86_FSW_DE
jmp .return_plus_inf
.base_zero_plus_exp:
call .is_exp_odd_integer
or eax, eax
jnz .return_base_value ; Matching 8
.return_plus_zero: ; Matching 9
fldz
jmp .return_pop_pop_val
;
; 15. base == -Inf and exp < 0 and exp is odd integer: Return -0
; 16. base == -Inf and exp < 0 and exp is not odd int: Return +0
; 17. base == -Inf and exp > 0 and exp is odd integer: Return -Inf
; 18. base == -Inf and exp > 0 and exp is not odd int: Return +Inf
; 19. base == +Inf and exp < 0: Return +0
; 20. base == +Inf and exp > 0: Return +Inf
;
; Note! Exponent must be finite and non-zero if we get here.
;
.base_inf:
fldz
fcomip st0, st1
jbe .base_inf_plus_exp
.base_inf_minus_exp:
test dx, X86_FSW_C1
jz .return_plus_zero ; Matches 19 (base == +Inf).
.base_minus_inf_minus_exp:
call .is_exp_odd_integer
or eax, eax
jz .return_plus_zero ; Matches 16 (exp not odd and < 0, base == -Inf)
.return_minus_zero: ; Matches 15 (exp is odd and < 0, base == -Inf)
fldz
fchs
jmp .return_pop_pop_val
.base_inf_plus_exp:
test dx, X86_FSW_C1
jz .return_plus_inf ; Matches 20 (base == +Inf).
.base_minus_inf_plus_exp:
call .is_exp_odd_integer
or eax, eax
jnz .return_minus_inf ; Matches 17 (exp is odd and > 0, base == +Inf)
jmp .return_plus_inf ; Matches 18 (exp not odd and > 0, base == +Inf)
;
; Return the exponent NaN (or whatever) value.
;
.return_exp_nan:
fld st0
mov eax, 2 ; return param 2
jmp .return_pop_pop_val_with_eax
;
; Return the base NaN (or whatever) value.
;
.return_base_nan:
.return_base_value:
.base_nan: ; 5. Unless specified elsewhere, return NaN if any of the parameters are NaN.
fld st1
mov eax, 1 ; return param 1
jmp .return_pop_pop_val_with_eax
;
; Pops the two values off the FPU stack and returns NaN.
;
.return_nan:
fld qword [.s_r64QNan xWrtRIP]
jmp .return_pop_pop_val
;
; Pops the two values off the FPU stack and returns +Inf.
;
.return_plus_inf:
fld qword [.s_r64PlusInf xWrtRIP]
jmp .return_pop_pop_val
;
; Pops the two values off the FPU stack and returns -Inf.
;
.return_minus_inf:
fld qword [.s_r64MinusInf xWrtRIP]
jmp .return_pop_pop_val
;
; Return st0, remove st1 and st2.
;
.return_pop_pop_val:
xor eax, eax
.return_pop_pop_val_with_eax:
fstp st2
ffreep st0
jmp .return
ALIGNCODE(8)
.s_r80MaxInt:
dt +9223372036854775807.0
ALIGNCODE(8)
.s_r80MinInt:
dt -9223372036854775807.0
ALIGNCODE(8)
;; The fyl2xp1 instruction only works between +/-1(1-sqrt(0.5)).
; These two variables is that range + 1.0, so we can compare directly
; with the input w/o any extra fsub and fabs work.
.s_r64AbsFyL2xP1InputMin:
dq 0.708 ; -0.292 + 1.0
.s_r64AbsFyL2xP1InputMax:
dq 1.292
.s_r64QNan:
dq RTFLOAT64U_QNAN_MINUS
.s_r64PlusInf:
dq RTFLOAT64U_INF_PLUS
.s_r64MinusInf:
dq RTFLOAT64U_INF_MINUS
;;
; Sub-function that checks if the exponent (st0) is an odd integer or not.
;
; @returns eax = 1 if odd, 0 if even or not integer.
; @uses eax, edx, eflags.
;
.is_exp_odd_integer:
;
; Save the FPU enviornment and mask all exceptions.
;
fnstenv [xBP - 30h]
mov ax, [xBP - 30h + X86FSTENV32P.FCW]
or word [xBP - 30h + X86FSTENV32P.FCW], X86_FCW_MASK_ALL
fldcw [xBP - 30h + X86FSTENV32P.FCW]
mov [xBP - 30h + X86FSTENV32P.FCW], ax
;
; Convert to 64-bit integer (probably not 100% correct).
;
fld st0 ; -> st0=exponent st1=exponent; st2=base;
fistp qword [xBP - 10h]
fild qword [xBP - 10h] ; -> st0=int(exponent) st1=exponent; st2=base;
fcomip st0, st1 ; -> st0=exponent; st1=base;
jne .is_exp_odd_integer__return_false ; jump if not integer.
mov xAX, [xBP - 10h]
%ifdef
mov edx, [xBP - 10h + 4]
%endif
;
; Check the lowest bit if it might be odd.
; This works both for positive and negative numbers.
;
test al, 1
jz .is_exp_odd_integer__return_false ; jump if even.
;
; If the result is negative, convert to positive.
;
%ifdef RT_ARCH_AMD64
bt rax, 63
%else
bt edx, 31
%endif
jnc .is_exp_odd_integer__positive
%ifdef RT_ARCH_AMD64
neg xAX
%else
neg edx
neg eax
sbb edx, 0
%endif
.is_exp_odd_integer__positive:
;
; Now find the most significant bit in the value so we can verify that
; the odd bit was part of the mantissa/fraction of the input.
;
cmp bl, 3 ; Skip if 80-bit input, as it has a 64-bit mantissa which
je .is_exp_odd_integer__return_true ; makes it a 1 bit more precision than out integer reg(s).
%ifdef RT_ARCH_AMD64
bsr rax, rax
%else
bsr edx, edx
jnz .is_exp_odd_integer__high_dword_is_zero
lea eax, [edx + 20h]
jmp .is_exp_odd_integer__first_bit_in_eax
.is_exp_odd_integer__high_dword_is_zero:
bsr eax, eax
.is_exp_odd_integer__first_bit_in_eax:
%endif
;
; The limit is 53 for double precision (one implicit bit + 52 bits fraction),
; and 24 for single precision types.
;
mov ah, 53 ; RTFLOAT64U_FRACTION_BITS + 1
cmp bl, 0
jne .is_exp_odd_integer__is_double_limit
mov ah, 24 ; RTFLOAT32U_FRACTION_BITS + 1
.is_exp_odd_integer__is_double_limit:
cmp al, ah
jae .is_exp_odd_integer__return_false
mov eax, 1
; Return.
.is_exp_odd_integer__return_true:
jmp .is_exp_odd_integer__return
.is_exp_odd_integer__return_false:
xor eax, eax
.is_exp_odd_integer__return:
ffreep st0
fldenv [xBP - 30h]
ret
ENDPROC rtNoCrtMathPowCore
|