forked from douglascrockford/DEC64
-
Notifications
You must be signed in to change notification settings - Fork 0
/
dec64.asm
1581 lines (1219 loc) · 55.7 KB
/
dec64.asm
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
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
title dec64.asm for x64.
; dec64.com
; 2015-12-25
; Public Domain
; No warranty expressed or implied. Use at your own risk. You have been warned.
; This file implements the elementary arithmetic operations for DEC64, a decimal
; floating point type. DEC64 uses 64 bits to represent a number. The low order
; 8 bits contain an exponent that ranges from -127 to 127. The exponent is not
; biased. The exponent -128 is reserved for nan (not a number). The remaining
; 56 bits, including the sign bit, are the coefficient in the range
; -36028797018963968 thru 36028797018963967. The exponent and the coefficient
; are both two's complement signed numbers.
;
; The value of any non-nan DEC64 number is coefficient * (10 ** exponent).
;
; Rounding is to the nearest value. Ties are rounded away from zero. Integer
; division is floored. The result of modulo has the sign of the divisor.
; There is no negative zero.
;
; These operations will produce a result of nan:
;
; dec64_abs(nan)
; dec64_ceiling(nan)
; dec64_dec(nan)
; dec64_floor(nan)
; dec64_half(nan)
; dec64_inc(nan)
; dec64_int(nan)
; dec64_neg(nan)
; dec64_normal(nan)
; dec64_not(nan)
; dec64_signum(nan)
;
; These operations will produce a result of zero for all values of n,
; even if n is nan:
;
; dec64_divide(0, n)
; dec64_integer_divide(0, n)
; dec64_modulo(0, n)
; dec64_multiply(0, n)
; dec64_multiply(n, 0)
;
; These operations will produce a result of nan for all values of n
; except zero:
;
; dec64_divide(n, 0)
; dec64_divide(n, nan)
; dec64_integer_divide(n, 0)
; dec64_integer_divide(n, nan)
; dec64_modulo(n, 0)
; dec64_modulo(n, nan)
; dec64_multiply(n, nan)
; dec64_multiply(nan, n)
;
; These operations will produce a result of nan for all values of n:
;
; dec64_add(n, nan)
; dec64_add(nan, n)
; dec64_divide(nan, n)
; dec64_integer_divide(nan, n)
; dec64_less(n, nan)
; dec64_less(nan, n)
; dec64_modulo(nan, n)
; dec64_round(n, nan)
; dec64_round(nan, n)
; dec64_subtract(n, nan)
; dec64_subtract(nan, n)
;
; This file can be processed with Microsoft's ML64.exe. There might be other
; assemblers that can process this file, but that has not been tested.
;
; You know what goes great with those defective pentium chips?
; Defective pentium salsa! (David Letterman)
; All public symbols have a dec64_ prefix. All other symbols are private.
; When these functions return nan, they will always return 0x80, the normal
; nan.
; There are 72057594037927936 possible nan values. Two are reserved to
; represent true and false. The comparison functions will return a boolean
; encoded as a DEC64 nan. Be care when using these comparisons with C
; because C thinks that the DEC64 false is truthy.
true equ 380h
false equ 280h
public dec64_abs;(number: dec64)
; returns absolution: dec64
public dec64_add;(augend: dec64, addend: dec64)
; returns sum: dec64
public dec64_ceiling;(number: dec64)
; returns integer: dec64
public dec64_coefficient;(number: dec64)
; returns coefficient: int64
public dec64_dec;(augend: dec64)
; returns decrementation: dec64
public dec64_divide;(dividend: dec64, divisor: dec64)
; returns quotient: dec64
public dec64_equal;(comparahend: dec64, comparator: dec64)
; returns comparison: dec64
public dec64_exponent;(number: dec64)
; returns exponent: int64
public dec64_floor;(number: dec64)
; returns integer: dec64
public dec64_half;(dividend: dec64)
; returns quotient: dec64
public dec64_inc;(augend: dec64)
; returns incrementation: dec64
public dec64_int;(number: dec64)
; returns integer: dec64
public dec64_integer_divide;(dividend: dec64, divisor: dec64)
; returns quotient: dec64
public dec64_is_any_nan;(number: dec64)
; returns comparison: dec64
public dec64_is_integer;(number: dec64)
; returns comparison: dec64
public dec64_is_zero;(number: dec64)
; returns comparison: dec64
public dec64_less;(comparahend: dec64, comparator: dec64)
; returns comparison: dec64
public dec64_modulo;(dividend: dec64, divisor: dec64)
; returns modulus: dec64
public dec64_multiply;(multiplicand: dec64, multiplier: dec64)
; returns product: dec64
public dec64_neg;(number: dec64)
; returns negation: dec64
public dec64_new;(coefficient: int64, exponent: int64)
; returns number: dec64
public dec64_normal;(number: dec64)
; returns normalization: dec64
public dec64_not;(boolean: dec64)
; returns notation: dec64
public dec64_round;(number: dec64, exponent: dec64)
; returns quantization: dec64
public dec64_signum;(number: dec64)
; returns signature: dec64
public dec64_subtract;(minuend: dec64, subtrahend: dec64)
; returns difference: dec64
; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
; Repair the register names. Over the long and twisted evolution of x86, the
; register names have picked up some weird, inconsistent conventions. We can
; simplify, naming them r0 thru r15, with a _b suffix to indicate the low
; order byte, an _h suffix to indicate the so-called high byte and a _w suffix
; to indicate the lower 16 bits. We would use _d to indicate the lower 32 bits,
; but that was not needed here.
r0 equ rax
r1 equ rcx
r2 equ rdx
r6 equ rsi
r7 equ rdi
r0_b equ al
r1_b equ cl
r2_b equ dl
r8_b equ r8b
r9_b equ r9b
r10_b equ r10b
r11_b equ r11b
r0_h equ ah
r1_h equ ch
r2_h equ dh
r0_w equ ax
r1_w equ cx
; All of the public functions in this file accept up to two arguments, which
; are passed in registers (either r1, r2 or r7, r6), returning a result in r0.
; Registers r1, r2, r8, r9, r10, and r11 are clobbered. Register r0 is the
; return value. The other registers are not disturbed.
; There is painfully inadequate standardization around x64 calling conventions.
; On Win64, the first two arguments are passed in r1 and r2. On Unix, the first
; two arguments are passed in r7 and r6. We try to hide this behind macros. The
; two systems also have different conventions about which registers may be
; clobbered and which must be preserved. This code lives in the intersection.
; This has not yet been tested on Unix.
UNIX equ 0 ; calling convention: 0 for Windows, 1 for Unix
function_with_one_parameter macro
if UNIX
mov r1,r7 ;; UNIX
endif
endm
function_with_two_parameters macro
if UNIX
mov r1,r7 ;; UNIX
mov r2,r6 ;; UNIX
endif
endm
call_with_one_parameter macro function
if UNIX
mov r7,r1 ;; UNIX
endif
call function
endm
call_with_two_parameters macro function
if UNIX
mov r7,r1 ;; UNIX
mov r6,r2 ;; UNIX
endif
call function
endm
tail_with_one_parameter macro function
if UNIX
mov r7,r1 ;; UNIX
endif
jmp function
endm
tail_with_two_parameters macro function
if UNIX
mov r7,r1 ;; UNIX
mov r6,r2 ;; UNIX
endif
jmp function
endm
; There may be a performance benefit in padding programs so that most jump
; destinations are aligned on 16 byte boundaries.
pad macro
align 16
endm
; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_data segment para read
power: ; the powers of 10
qword 1 ; 0
qword 10 ; 1
qword 100 ; 2
qword 1000 ; 3
qword 10000 ; 4
qword 100000 ; 5
qword 1000000 ; 6
qword 10000000 ; 7
qword 100000000 ; 8
qword 1000000000 ; 9
qword 10000000000 ; 10
qword 100000000000 ; 11
qword 1000000000000 ; 12
qword 10000000000000 ; 13
qword 100000000000000 ; 14
qword 1000000000000000 ; 15
qword 10000000000000000 ; 16
qword 100000000000000000 ; 17
qword 1000000000000000000 ; 18
qword 10000000000000000000 ; 19
dec64_data ends
; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_code segment para execute
dec64_coefficient: function_with_one_parameter
;(number: dec64) returns coefficient: int64
; Return the coefficient part from a dec64 number.
mov r0,r1
sar r0,8
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_exponent: function_with_one_parameter
;(number: dec64) returns exponent: int64
; Return the exponent part, sign extended to 64 bits, from a dec64 number.
; dec64_exponent(nan) returns -128.
movsx r0,r1_b ; r0 is the exponent
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_new: function_with_two_parameters
;(coefficient: int64, exponent: int64) returns number: dec64
; Construct a new dec64 number with a coefficient and an exponent.
mov r0,r1 ; r0 is the coefficient
test r0,r0 ; is the coefficient zero?
jz return_zero
mov r8,r2 ; r8 is the exponent
; Fall into pack.
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
pack:
; The pack function will combine the coefficient and exponent into a dec64.
; Numbers that are too huge to be contained in this format become nan or -nan.
; Numbers that are too tiny to be contained in this format become zero.
; The coefficient is in r0.
; The exponent is in r8.
; If the exponent is greater than 127, then the number is too big.
; But it might still be possible to salvage a value.
cmp r8,127 ; compare exponent with 127
jg pack_decrease ; it might be possible to decrease it
; If the exponent is too small, or if the coefficient is too large, then some
; division is necessary. The absolute value of the coefficient is off by one
; for the negative because
; negative_extreme_coefficent = -(extreme_coefficent + 1)
mov r10,r0 ; r10 is the coefficient
mov r1,3602879701896396800 ; the ultimate coefficient * 100
not r10 ; r10 is -coefficient
xor r11,r11 ; r11 is zero
test r10,r10 ; look at the sign bit
cmovs r10,r0 ; r10 is the absolute value of the coefficient
cmp r10,r1 ; compare with the actual coefficient
jae pack_large ; deal with the very large coefficient
mov r1,36028797018963967 ; the ultimate coefficient - 1
mov r9,-127 ; r9 is the ultimate exponent
cmp r1,r10 ; compare with the actual coefficient
adc r11,0 ; add 1 to r11 if 1 digit too big
mov r1,360287970189639679 ; the ultimate coefficient * 10 - 1
sub r9,r8 ; r9 is the difference from the actual exponent
cmp r1,r10 ; compare with the actual coefficient
adc r11,0 ; add 1 to r11 if 2 digits too big
cmp r9,r11 ; which excess is larger?
cmovl r9,r11 ; take the max
test r9,r9 ; if neither was zero
jnz pack_increase ; then increase the exponent by the excess
shl r0,8 ; shift the exponent into position
cmovz r8,r0 ; if the coefficient is zero, also zero the exp
movzx r8,r8_b ; zero out all but the bottom 8 bits of the exp
or r0,r8 ; mix in the exponent
ret ; whew
pad
pack_large:
mov r10,10 ; r10 is 10
cqo ; sign extend r0 into r2
idiv r10 ; divide r0 by ten
add r8,1 ; add 1 to the exponent
jmp pack ; start over
pad
pack_increase:
mov r10,power
mov r10,[r10][r9*8] ; r10 is 10^r9
cmp r9,20 ; is the difference more than 20?
jae return_zero ; if so, the result is zero (rare)
mov r11,r10 ; r11 is the power of ten
neg r11 ; r11 is the negation of the power of ten
test r0,r0 ; examine the sign of the coefficient
cmovns r11,r10 ; r11 has the sign of the coefficient
sar r11,1 ; r11 is half the signed power of ten
add r0,r11 ; r0 is the coefficient plus the rounding fudge
cqo ; sign extend r0 into r2
idiv r10 ; divide by the power of ten
add r8,r9 ; increase the exponent
jmp pack ; start over
pad
pack_decrease:
; The exponent is too big. We can attempt to reduce it by scaling back.
; This can salvage values in a small set of cases.
imul r0,10 ; try multiplying the coefficient by 10
jo return_nan ; if it overflows, we failed to salvage
sub r8,1 ; decrement the exponent
cmp r8,127 ; test the exponent
jg pack_decrease ; until the exponent is less than 127
mov r9,r0 ; r9 is the coefficient
xor r1,r1 ; the rounding flag is zero
sar r9,56 ; r9 is top 8 bits of the coefficient
adc r9,0 ; add the ninth bit
jnz return_nan ; the number is still too large
shl r0,8 ; shift the exponent into position
cmovz r8,r0 ; if the coefficient is zero, also zero the exponent
movzx r8,r8_b ; zero out all but the bottom 8 bits of the exponent
or r0,r8 ; mix in the exponent
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_round: function_with_two_parameters
;(number: dec64, place: dec64) returns quantization: dec64
; The place argument indicates at what decimal place to round.
; -2 nearest cent
; 0 nearest integer
; 3 nearest thousand
; 6 nearest million
; 9 nearest billion
; The place should be between -16 and 16.
mov r11,r1 ; preserve the number
mov r1,r2 ; pass the place
call_with_one_parameter dec64_int
mov r9,r0 ; r9 is the normalized place
cmp r11_b,128 ; is the first operand not nan?
setne r1_b
sar r9,8 ; r9 is the place as int64
cmp r0_b,128 ; is the second operand nan?
setne r1_h
movsx r8,r11_b ; r8 is the current exponent
mov r0,r11 ; r0 is the number
test r1_b,r1_h ; is either nan?
jz return_nan ; if so, the result is nan
sar r0,8 ; r0 is the coefficient
jz return ; if the coefficient is zero, the result is zero
cmp r8,r9 ; compare the exponents
jge pack ; no rounding required
mov r10,10 ; r10 is the divisor
pad
round_loop:
; Increment the exponent and divide the coefficient by 10 until the target
; exponent is reached.
cqo ; sign extend r0 into r2
idiv r10 ; divide r2:r0 by 10
add r8,1 ; increment the exponent
cmp r8,r9 ; compare the exponents
jne round_loop ; loop if the exponent has not reached the target
; Round if necessary and return the result.
xor r1,r1 ; r1 is zero
cmp r2_b,-5 ; compare the remainder to -5
setle r1_b ; r1 is 1 if remainder <= -5 else 0
cmp r2_b,5 ; compare the remainder to 5
setge r2_h ; if the remainder is greater than or equal to 5
neg r1 ; r1 is -1 or 0
or r1_b,r2_h ; r1 is the rounding flag: -1, 0, or 1
add r0,r1 ; r0 is the rounded coefficient
jmp pack ; pack it up
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_add: function_with_two_parameters
;(augend: dec64, addend: dec64) returns sum: dec64
; Add two dec64 numbers together.
; Registers:
; r0 result, imul/idiv
; r1 first argument
; r2 second argument, imul/idiv
; r8 exponent of first argument
; r9 exponent of second argument
; r10 coefficient of first argument
; r11 coefficient of second argument
; If the two exponents are both zero (which is usually the case for integers)
; we can take the fast path. Since the exponents are both zero, we can simply
; add the numbers together and check for overflow.
pad
add_begin:
mov r0,r1 ; r0 is the first number
or r0_b,r2_b ; r0_b is the two exponents or'd together
jnz add_slow ; if either exponent is not zero, take the slow path
add r0,r2 ; add the coefficients together
jo add_overflow ; if there was no overflow, we are done
ret ; no need to pack
pad
add_overflow:
; If there was an overflow (extremely unlikely) then we must make it fit.
; pack knows how to do that.
rcr r0,1 ; divide the sum by 2 and repair its sign
movsx r8,r1_b ; r8 is the exponent
sar r0,7 ; r0 is the coefficient of the sum
jmp pack ; pack it up
pad
add_slow:
; The slow path is taken if the two operands do not both have zero exponents.
cmp r1_b,128 ; Is the first operand nan?
je return_nan
; Are the two exponents the same? This will happen often, especially with
; money values.
cmp r1_b,r2_b ; compare the two exponents
jne add_slower ; if not equal, take the slower path
; The exponents match so we may add now. Zero out the exponents so there
; will be no carry into the coefficients when the coefficients are added.
; If the result is zero, then return the normal zero.
xor r0_b,r0_b ; remove the first exponent
xor r2_b,r2_b ; remove the second exponent
add r0,r2 ; add the coefficients
jo add_overflow ; if it overflows, it must be repaired
cmovz r1,r0 ; if the coefficient is zero, the exponent is zero
mov r0_b,r1_b ; mix in the exponent
ret ; no need to pack
pad
add_slower:
; The slower path is taken when neither operand is nan, and their
; exponents are different. Before addition can take place, the exponents
; must be made to match. Swap the numbers if the second exponent is greater
; than the first.
cmp r2_b,128 ; Is the second operand nan?
je return_nan
cmp r1_b,r2_b ; compare the exponents
mov r0,r1 ; r0 is the first number
cmovl r1,r2 ; r1 is the number with the larger exponent
cmovl r2,r0 ; r2 is the number with the smaller exponent
; Shift the coefficients of r1 and r2 into r10 and r11 and unpack the exponents.
mov r10,r1 ; r10 is the first number
mov r11,r2 ; r11 is the second number
movsx r8,r1_b ; r8 is the first exponent
movsx r9,r2_b ; r9 is the second exponent
sar r10,8 ; r10 is the first coefficient
sar r11,8 ; r11 is the second coefficient
mov r0,r10 ; r0 is the first coefficient
pad
add_slower_decrease:
; The coefficients are not the same. Before we can add, they must be the same.
; We will try to decrease the first exponent. When we decrease the exponent
; by 1, we must also multiply the coefficient by 10. We can do this as long as
; there is no overflow. We have 8 extra bits to work with, so we can do this
; at least twice, possibly more.
imul r0,10 ; before decrementing the exponent, multiply
jo add_slower_increase
sub r8,1 ; decrease the first exponent
mov r10,r0 ; r10 is the enlarged first coefficient
cmp r8,r9 ; are the exponents equal yet?
jg add_slower_decrease
mov r0,r11 ; r0 is the second coefficient
; The exponents are now equal, so the coefficients may be added.
add r0,r10 ; add the two coefficients
jmp pack ; pack it up
pad
add_slower_increase:
; We cannot decrease the first exponent any more, so we must instead try to
; increase the second exponent, which will result in a loss of significance.
; That is the heartbreak of floating point.
; Determine how many places need to be shifted. If it is more than 17, there is
; nothing more to add.
mov r2,r8 ; r2 is the first exponent
sub r2,r9 ; r2 is the remaining exponent difference
mov r0,r11 ; r0 is the second coefficient
cmp r2,17 ; 17 is the max digits in a packed coefficient
ja return_r1 ; too small to matter
mov r9,power
mov r9,[r9][r2*8] ; r9 is the power of ten
cqo ; sign extend r0 into r2
idiv r9 ; divide the second coefficient by the power of 10
test r0,r0 ; examine the scaled coefficient
jz return_r1 ; too insignificant to add?
; The exponents are now equal, so the coefficients may be added.
add r0,r10 ; add the two coefficients
jmp pack
pad
return_r1:
mov r0,r1 ; r0 is the original number
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_inc: function_with_one_parameter
;(augend: dec64) returns incrementation: dec64
; Increment a number. In most cases, this will be a faster way to add one than
; dec64_add.
test r1_b,r1_b ; what is the exponent?
jnz inc_not_integer
; The number is an integer. This might be easy.
mov r0,100h ; r0 is one
add r0,r1 ; r0 is the incrementation
jo inc_hardway ; overflow (very rare)
ret
pad
inc_not_integer:
js inc_negative_exponent
test r1,-256 ; is the coefficient zero?
jz return_one ; if so, the result is one
cmp r1_b,17 ; is the number too enormous to increment?
jge return_r1 ; then return the number
inc_hardway:
mov r2,100h ; r2 is one
tail_with_two_parameters dec64_add
pad
inc_negative_exponent:
cmp r1_b,128 ; is the number nan?
je return_nan
cmp r1_b,-17 ; is the number too small to increment?
jle return_one ; then return one
movsx r8,r1_b ; r8 is the negative exponent
neg r8 ; flip the sign
mov r9,power
mov r0,[r9][r8*8] ; r0 is 10^(-exponent)
shl r0,8 ; convert to dec64
add r0,r1 ; now we add
jo inc_hardway ; if it overflows, do it the hard way
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_dec: function_with_one_parameter
;(augend: dec64) returns decrementation: dec64
; Increment a number. In most cases, this will be a faster way to subtract one
; than dec64_subtract.
test r1_b,r1_b ; what is the exponent?
jnz dec_not_integer
; The number is an integer. This might be easy.
mov r0,-256 ; r0 is negative one
add r0,r1 ; r0 is the decrementation
jo dec_hardway ; overflow (very rare)
ret
pad
dec_not_integer:
js dec_negative_exponent
test r1,-256 ; is the coefficient zero?
jz dec_neg_one ; if so, the result is negative one
cmp r1_b,17 ; is the number too enormous to decrement?
jge return_r1 ; then return the number
dec_hardway:
mov r2,100h ; r2 is one
tail_with_two_parameters dec64_subtract
pad
dec_negative_exponent:
cmp r1_b,128 ; is the number nan?
je return_nan
cmp r1_b,-17 ; is the number too small to decrement?
jle dec_neg_one ; then return negative one
movsx r8,r1_b ; r8 is the negative exponent
neg r8 ; flip the sign
mov r9,power
mov r0,[r9][r8*8] ; r0 is 10^(-exponent)
neg r0 ; go negative
shl r0,8 ; convert to dec64
add r0,r1 ; now we subtract
jo dec_hardway ; if it overflows, do it the hard way
ret
pad
dec_neg_one:
mov r0,-256 ; r0 is -1
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_ceiling: function_with_one_parameter
;(number: dec64) returns integer: dec64
; Produce the smallest integer that is greater than or equal to the number. In
; the result, the exponent will be greater than or equal to zero unless it is
; nan. Numbers with positive exponents will not be modified, even if the
; numbers are outside of the safe integer range.
; Preserved: r11.
mov r9,1 ; r9 is the round up flag
jmp floor_begin
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_floor: function_with_one_parameter
;(number: dec64) returns integer: dec64
; Produce the largest integer that is less than or equal to the number. This
; is sometimes called the entier function. In the result, the exponent will be
; greater than or equal to zero unless it is nan. Numbers with positive
; exponents will not be modified, even if the numbers are outside of the safe
; integer range.
; Preserved: r11.
mov r9,-1 ; r9 is the round down flag
pad
floor_begin:
cmp r1_b,128 ; compare the exponent to nan
je return_nan ; if the exponent is nan, the result is nan
mov r0,r1 ; r0 is the number
movsx r8,r1_b ; r8 is the exponent
sar r0,8 ; r0 is the coefficient
cmovz r1,r0 ; if the coefficient is zero, the number is zero
neg r8 ; r8 is the negated exponent
test r1_b,r1_b ; examine the exponent
jns return_r1 ; nothing to do unless the exponent was negative
cmp r8,17 ; is the exponent is too extreme?
jae floor_micro ; deal with a micro number
mov r10,power
mov r10,[r10][r8*8] ; r10 is the power of ten
cqo ; sign extend r0 into r2
idiv r10 ; divide r2:r0 by 10
test r2,r2 ; examine the remainder
jnz floor_remains ; deal with the remainder
shl r0,8 ; pack the coefficient
ret
pad
floor_micro:
mov r2,r0 ; r2 is the coefficient
xor r0,r0 ; r0 is zero
pad
floor_remains:
; If the remainder is negative and the rounding flag is negative, then we need
; to decrement r0. But if the remainder and the rounding flag are both
; positive, then we need to increment r0.
xor r10,r10 ; r10 is zero
xor r2,r9 ; xor the remainder and the rounding
cmovs r9,r10 ; if they had different signs, clear the rounding
add r0,r9 ; add the rounding to the coefficient
shl r0,8 ; pack the coefficient
ret
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_subtract: function_with_two_parameters
;(minuend: dec64, subtrahend: dec64) returns difference: dec64
; Subtract the dec64 number in r2 from the dec64 number in r1.
; The result is in r0.
; This is the same as dec64_add, except that the operand in r2 has its
; coefficient complemented first.
xor r2,-256 ; not the coefficient only
add r2,256 ; the two's complement increment of the coefficient
jno add_begin ; if there is no overflow, begin the beguine
; The subtrahend coefficient is -36028797018963968. This value cannot easily be
; complemented, so take the slower path. This should be extremely rare.
cmp r1_b,128 ; is the first operand nan
sete r0_b
cmp r2_b,128 ; is the second operand nan?
sete r0_h
or r0_b,r0_h ; is either nan?
jnz return_nan
mov r10,r1 ; r10 is the first coefficient
movsx r8,r1_b ; r8 is the first exponent
sar r10,8
movsx r9,r2_b ; r9 is the second exponent
mov r11,80000000000000H ; r11 is 36028797018963968
mov r0,r10 ; r0 is the first coefficient
cmp r8,r9 ; if the second exponent is larger, swap
jge subtract_slower_decrease_compare
mov r0,r11 ; r0 is the second coefficient
xchg r8,r9 ; swap the exponents
xchg r10,r11 ; swap the coefficients
jmp subtract_slower_decrease_compare
subtract_slower_decrease:
; The coefficients are not the same. Before we can add, they must be the same.
; We will try to decrease the first exponent. When we decrease the exponent
; by 1, we must also multiply the coefficient by 10. We can do this as long as
; there is no overflow. We have 8 extra bits to work with, so we can do this
; at least twice, possibly more.
imul r0,10 ; before decrementing the exponent, multiply
jo subtract_slower_increase
sub r8,1 ; decrease the first exponent
mov r10,r0 ; r10 is the enlarged first coefficient
subtract_slower_decrease_compare:
cmp r8,r9 ; are the exponents equal yet?
jg subtract_slower_decrease
mov r0,r11 ; r0 is the second coefficient
; The exponents are now equal, so the coefficients may be added.
add r0,r10 ; add the two coefficients
jmp pack ; pack it up
subtract_slower_increase:
; We cannot decrease the first exponent any more, so we must instead try to
; increase the second exponent, which will result in a loss of significance.
; That is the heartbreak of floating point.
; Determine how many places need to be shifted. If it is more than 17, there is
; nothing more to add.
mov r2,r8 ; r2 is the first exponent
sub r2,r9 ; r2 is the remaining exponent difference
mov r0,r11 ; r0 is the second coefficient
cmp r2,17 ; 17 is the max digits in a packed coefficient
ja subtract_underflow ; too small to matter
mov r9,power
mov r9,[r9][r2*8] ; r9 is the power of ten
cqo ; sign extend r0 into r2
idiv r9 ; divide the second coefficient by the power of 10
; The exponents are now equal, so the coefficients may be added.
add r0,r10 ; add the two coefficients
jmp pack
subtract_underflow:
mov r0,r10 ; r0 is the first coefficient
jmp pack
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_multiply: function_with_two_parameters
;(multiplicand: dec64, multiplier: dec64) returns product: dec64
; Multiply two dec64 numbers together.
; Unpack the exponents into r8 and r9.
movsx r8,r1_b ; r8 is the first exponent
movsx r9,r2_b ; r9 is the second exponent
; Set flags in r0 if either operand is nan.
cmp r1_b,128 ; is the first operand nan?
sete r0_b ; r0_b is 1 if the first operand is nan
cmp r2_b,128 ; is the second operand nan?
sete r0_h ; r0_h is 1 if the second operand is nan
; Unpack the coefficients. Set flags in r1 if either is not zero.
sar r1,8 ; r1 is the first coefficient
mov r10,r1 ; r10 is the first coefficient
setnz r1_b ; r1_b is 1 if the first coefficient is not zero
sar r2,8 ; r2 is the second coefficient
setnz r1_h ; r1_h is 1 if the second coefficient is not zero
; The result is nan if one or both of the operands is nan and neither of the
; operands is zero.
or r1_w,r0_w ; is either coefficient zero and not nan?
xchg r1_b,r1_h
test r0_w,r1_w
jnz return_nan
mov r0,r10 ; r0 is the first coefficient
add r8,r9 ; r8 is the product exponent
imul r2 ; r2:r0 is r1 * r2
jno pack ; if the product fits in 64 bits, start packing
; There was overflow.
; Make the 110 bit coefficient in r2:r0Er8 all fit. Estimate the number of
; digits of excess, and increase the exponent by that many digits.
; We use 77/256 to convert log2 to log10.
mov r9,r2 ; r9 is the excess significance
xor r1,r1 ; r1 is zero anticipating bsr
neg r9
cmovs r9,r2 ; r9 is absolute value of the excess significance
bsr r1,r9 ; find the position of the most significant bit
imul r1,77 ; multiply the bit number by 77/256 to
shr r1,8 ; convert a bit number to a digit number
add r1,2 ; add two extra digits to the scale
add r8,r1 ; increase the exponent
mov r9,power
idiv qword ptr [r9][r1*8] ; divide by the power of ten
jmp pack
pad; -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
dec64_divide: function_with_two_parameters
;(dividend: dec64, divisor: dec64) returns quotient: dec64
; Divide a dec64 number by another.
; Begin unpacking the components.
movsx r8,r1_b ; r8 is the first exponent
movsx r9,r2_b ; r9 is the second exponent
mov r10,r1 ; r10 is the first number
mov r11,r2 ; r11 is the second number
; Set nan flags in r0.
cmp r1_b,128 ; is the first operand nan?
sete r0_b ; r0_b is 1 if the first operand is nan
cmp r2_b,128 ; is the second operand nan?
sete r0_h ; r0_h is 1 if the second operand is nan
sar r10,8 ; r10 is the dividend coefficient
setnz r1_b ; r1_b is 1 if the dividend coefficient is zero
sar r11,8 ; r11 is the divisor coefficient
setz r1_h ; r1_h is 1 if dividing by zero
or r0_h,r0_b ; r0_h is 1 if either is nan
or r1_b,r0_b ; r1_b is zero if the dividend is zero and not nan
jz return_zero ; if the dividend is zero, the quotient is zero
sub r8,r9 ; r8 is the quotient exponent
or r0_b,r1_h ; r0_b is 1 if either is nan or the divisor is zero
jnz return_nan
pad
divide_measure: