-
Notifications
You must be signed in to change notification settings - Fork 137
/
Copy pathfloatconv-submodule-code.c
1956 lines (1785 loc) · 66.2 KB
/
floatconv-submodule-code.c
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
// After editing this file, run "go generate" in the ../data directory.
// Copyright 2020 The Wuffs Authors.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// https://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// ---------------- IEEE 754 Floating Point
WUFFS_BASE__MAYBE_STATIC wuffs_base__lossy_value_u16 //
wuffs_base__ieee_754_bit_representation__from_f64_to_u16_truncate(double f) {
uint64_t u = 0;
if (sizeof(uint64_t) == sizeof(double)) {
memcpy(&u, &f, sizeof(uint64_t));
}
uint16_t neg = ((uint16_t)((u >> 63) << 15));
u &= 0x7FFFFFFFFFFFFFFF;
uint64_t exp = u >> 52;
uint64_t man = u & 0x000FFFFFFFFFFFFF;
if (exp == 0x7FF) {
if (man == 0) { // Infinity.
wuffs_base__lossy_value_u16 ret;
ret.value = neg | 0x7C00;
ret.lossy = false;
return ret;
}
// NaN. Shift the 52 mantissa bits to 10 mantissa bits, keeping the most
// significant mantissa bit (quiet vs signaling NaNs). Also set the low 9
// bits of ret.value so that the 10-bit mantissa is non-zero.
wuffs_base__lossy_value_u16 ret;
ret.value = neg | 0x7DFF | ((uint16_t)(man >> 42));
ret.lossy = false;
return ret;
} else if (exp > 0x40E) { // Truncate to the largest finite f16.
wuffs_base__lossy_value_u16 ret;
ret.value = neg | 0x7BFF;
ret.lossy = true;
return ret;
} else if (exp <= 0x3E6) { // Truncate to zero.
wuffs_base__lossy_value_u16 ret;
ret.value = neg;
ret.lossy = (u != 0);
return ret;
} else if (exp <= 0x3F0) { // Normal f64, subnormal f16.
// Convert from a 53-bit mantissa (after realizing the implicit bit) to a
// 10-bit mantissa and then adjust for the exponent.
man |= 0x0010000000000000;
uint32_t shift = ((uint32_t)(1051 - exp)); // 1051 = 0x3F0 + 53 - 10.
uint64_t shifted_man = man >> shift;
wuffs_base__lossy_value_u16 ret;
ret.value = neg | ((uint16_t)shifted_man);
ret.lossy = (shifted_man << shift) != man;
return ret;
}
// Normal f64, normal f16.
// Re-bias from 1023 to 15 and shift above f16's 10 mantissa bits.
exp = (exp - 1008) << 10; // 1008 = 1023 - 15 = 0x3FF - 0xF.
// Convert from a 52-bit mantissa (excluding the implicit bit) to a 10-bit
// mantissa (again excluding the implicit bit). We lose some information if
// any of the bottom 42 bits are non-zero.
wuffs_base__lossy_value_u16 ret;
ret.value = neg | ((uint16_t)exp) | ((uint16_t)(man >> 42));
ret.lossy = (man << 22) != 0;
return ret;
}
WUFFS_BASE__MAYBE_STATIC wuffs_base__lossy_value_u32 //
wuffs_base__ieee_754_bit_representation__from_f64_to_u32_truncate(double f) {
uint64_t u = 0;
if (sizeof(uint64_t) == sizeof(double)) {
memcpy(&u, &f, sizeof(uint64_t));
}
uint32_t neg = ((uint32_t)(u >> 63)) << 31;
u &= 0x7FFFFFFFFFFFFFFF;
uint64_t exp = u >> 52;
uint64_t man = u & 0x000FFFFFFFFFFFFF;
if (exp == 0x7FF) {
if (man == 0) { // Infinity.
wuffs_base__lossy_value_u32 ret;
ret.value = neg | 0x7F800000;
ret.lossy = false;
return ret;
}
// NaN. Shift the 52 mantissa bits to 23 mantissa bits, keeping the most
// significant mantissa bit (quiet vs signaling NaNs). Also set the low 22
// bits of ret.value so that the 23-bit mantissa is non-zero.
wuffs_base__lossy_value_u32 ret;
ret.value = neg | 0x7FBFFFFF | ((uint32_t)(man >> 29));
ret.lossy = false;
return ret;
} else if (exp > 0x47E) { // Truncate to the largest finite f32.
wuffs_base__lossy_value_u32 ret;
ret.value = neg | 0x7F7FFFFF;
ret.lossy = true;
return ret;
} else if (exp <= 0x369) { // Truncate to zero.
wuffs_base__lossy_value_u32 ret;
ret.value = neg;
ret.lossy = (u != 0);
return ret;
} else if (exp <= 0x380) { // Normal f64, subnormal f32.
// Convert from a 53-bit mantissa (after realizing the implicit bit) to a
// 23-bit mantissa and then adjust for the exponent.
man |= 0x0010000000000000;
uint32_t shift = ((uint32_t)(926 - exp)); // 926 = 0x380 + 53 - 23.
uint64_t shifted_man = man >> shift;
wuffs_base__lossy_value_u32 ret;
ret.value = neg | ((uint32_t)shifted_man);
ret.lossy = (shifted_man << shift) != man;
return ret;
}
// Normal f64, normal f32.
// Re-bias from 1023 to 127 and shift above f32's 23 mantissa bits.
exp = (exp - 896) << 23; // 896 = 1023 - 127 = 0x3FF - 0x7F.
// Convert from a 52-bit mantissa (excluding the implicit bit) to a 23-bit
// mantissa (again excluding the implicit bit). We lose some information if
// any of the bottom 29 bits are non-zero.
wuffs_base__lossy_value_u32 ret;
ret.value = neg | ((uint32_t)exp) | ((uint32_t)(man >> 29));
ret.lossy = (man << 35) != 0;
return ret;
}
// --------
#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE 2047
#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION 800
// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL is the largest N
// such that ((10 << N) < (1 << 64)).
#define WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL 60
// wuffs_base__private_implementation__high_prec_dec (abbreviated as HPD) is a
// fixed precision floating point decimal number, augmented with ±infinity
// values, but it cannot represent NaN (Not a Number).
//
// "High precision" means that the mantissa holds 800 decimal digits. 800 is
// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION.
//
// An HPD isn't for general purpose arithmetic, only for conversions to and
// from IEEE 754 double-precision floating point, where the largest and
// smallest positive, finite values are approximately 1.8e+308 and 4.9e-324.
// HPD exponents above +2047 mean infinity, below -2047 mean zero. The ±2047
// bounds are further away from zero than ±(324 + 800), where 800 and 2047 is
// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION and
// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
//
// digits[.. num_digits] are the number's digits in big-endian order. The
// uint8_t values are in the range [0 ..= 9], not ['0' ..= '9'], where e.g. '7'
// is the ASCII value 0x37.
//
// decimal_point is the index (within digits) of the decimal point. It may be
// negative or be larger than num_digits, in which case the explicit digits are
// padded with implicit zeroes.
//
// For example, if num_digits is 3 and digits is "\x07\x08\x09":
// - A decimal_point of -2 means ".00789"
// - A decimal_point of -1 means ".0789"
// - A decimal_point of +0 means ".789"
// - A decimal_point of +1 means "7.89"
// - A decimal_point of +2 means "78.9"
// - A decimal_point of +3 means "789."
// - A decimal_point of +4 means "7890."
// - A decimal_point of +5 means "78900."
//
// As above, a decimal_point higher than +2047 means that the overall value is
// infinity, lower than -2047 means zero.
//
// negative is a sign bit. An HPD can distinguish positive and negative zero.
//
// truncated is whether there are more than
// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION digits, and at
// least one of those extra digits are non-zero. The existence of long-tail
// digits can affect rounding.
//
// The "all fields are zero" value is valid, and represents the number +0.
typedef struct wuffs_base__private_implementation__high_prec_dec__struct {
uint32_t num_digits;
int32_t decimal_point;
bool negative;
bool truncated;
uint8_t digits[WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION];
} wuffs_base__private_implementation__high_prec_dec;
// wuffs_base__private_implementation__high_prec_dec__trim trims trailing
// zeroes from the h->digits[.. h->num_digits] slice. They have no benefit,
// since we explicitly track h->decimal_point.
//
// Preconditions:
// - h is non-NULL.
static inline void //
wuffs_base__private_implementation__high_prec_dec__trim(
wuffs_base__private_implementation__high_prec_dec* h) {
while ((h->num_digits > 0) && (h->digits[h->num_digits - 1] == 0)) {
h->num_digits--;
}
}
// wuffs_base__private_implementation__high_prec_dec__assign sets h to
// represent the number x.
//
// Preconditions:
// - h is non-NULL.
static void //
wuffs_base__private_implementation__high_prec_dec__assign(
wuffs_base__private_implementation__high_prec_dec* h,
uint64_t x,
bool negative) {
uint32_t n = 0;
// Set h->digits.
if (x > 0) {
// Calculate the digits, working right-to-left. After we determine n (how
// many digits there are), copy from buf to h->digits.
//
// UINT64_MAX, 18446744073709551615, is 20 digits long. It can be faster to
// copy a constant number of bytes than a variable number (20 instead of
// n). Make buf large enough (and start writing to it from the middle) so
// that can we always copy 20 bytes: the slice buf[(20-n) .. (40-n)].
uint8_t buf[40] = {0};
uint8_t* ptr = &buf[20];
do {
uint64_t remaining = x / 10;
x -= remaining * 10;
ptr--;
*ptr = (uint8_t)x;
n++;
x = remaining;
} while (x > 0);
memcpy(h->digits, ptr, 20);
}
// Set h's other fields.
h->num_digits = n;
h->decimal_point = (int32_t)n;
h->negative = negative;
h->truncated = false;
wuffs_base__private_implementation__high_prec_dec__trim(h);
}
static wuffs_base__status //
wuffs_base__private_implementation__high_prec_dec__parse(
wuffs_base__private_implementation__high_prec_dec* h,
wuffs_base__slice_u8 s,
uint32_t options) {
if (!h) {
return wuffs_base__make_status(wuffs_base__error__bad_receiver);
}
h->num_digits = 0;
h->decimal_point = 0;
h->negative = false;
h->truncated = false;
uint8_t* p = s.ptr;
uint8_t* q = s.ptr + s.len;
if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) {
for (;; p++) {
if (p >= q) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
} else if (*p != '_') {
break;
}
}
}
// Parse sign.
do {
if (*p == '+') {
p++;
} else if (*p == '-') {
h->negative = true;
p++;
} else {
break;
}
if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) {
for (;; p++) {
if (p >= q) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
} else if (*p != '_') {
break;
}
}
}
} while (0);
// Parse digits, up to (and including) a '.', 'E' or 'e'. Examples for each
// limb in this if-else chain:
// - "0.789"
// - "1002.789"
// - ".789"
// - Other (invalid input).
uint32_t nd = 0;
int32_t dp = 0;
bool no_digits_before_separator = false;
if (('0' == *p) &&
!(options &
WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_MULTIPLE_LEADING_ZEROES)) {
p++;
for (;; p++) {
if (p >= q) {
goto after_all;
} else if (*p ==
((options &
WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
? ','
: '.')) {
p++;
goto after_sep;
} else if ((*p == 'E') || (*p == 'e')) {
p++;
goto after_exp;
} else if ((*p != '_') ||
!(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
}
} else if (('0' <= *p) && (*p <= '9')) {
if (*p == '0') {
for (; (p < q) && (*p == '0'); p++) {
}
} else {
h->digits[nd++] = (uint8_t)(*p - '0');
dp = (int32_t)nd;
p++;
}
for (;; p++) {
if (p >= q) {
goto after_all;
} else if (('0' <= *p) && (*p <= '9')) {
if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->digits[nd++] = (uint8_t)(*p - '0');
dp = (int32_t)nd;
} else if ('0' != *p) {
// Long-tail non-zeroes set the truncated bit.
h->truncated = true;
}
} else if (*p ==
((options &
WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
? ','
: '.')) {
p++;
goto after_sep;
} else if ((*p == 'E') || (*p == 'e')) {
p++;
goto after_exp;
} else if ((*p != '_') ||
!(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
}
} else if (*p == ((options &
WUFFS_BASE__PARSE_NUMBER_FXX__DECIMAL_SEPARATOR_IS_A_COMMA)
? ','
: '.')) {
p++;
no_digits_before_separator = true;
} else {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
after_sep:
for (;; p++) {
if (p >= q) {
goto after_all;
} else if ('0' == *p) {
if (nd == 0) {
// Track leading zeroes implicitly.
dp--;
} else if (nd <
WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->digits[nd++] = (uint8_t)(*p - '0');
}
} else if (('0' < *p) && (*p <= '9')) {
if (nd < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->digits[nd++] = (uint8_t)(*p - '0');
} else {
// Long-tail non-zeroes set the truncated bit.
h->truncated = true;
}
} else if ((*p == 'E') || (*p == 'e')) {
p++;
goto after_exp;
} else if ((*p != '_') ||
!(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
}
after_exp:
do {
if (options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES) {
for (;; p++) {
if (p >= q) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
} else if (*p != '_') {
break;
}
}
}
int32_t exp_sign = +1;
if (*p == '+') {
p++;
} else if (*p == '-') {
exp_sign = -1;
p++;
}
int32_t exp = 0;
const int32_t exp_large =
WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE +
WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION;
bool saw_exp_digits = false;
for (; p < q; p++) {
if ((*p == '_') &&
(options & WUFFS_BASE__PARSE_NUMBER_XXX__ALLOW_UNDERSCORES)) {
// No-op.
} else if (('0' <= *p) && (*p <= '9')) {
saw_exp_digits = true;
if (exp < exp_large) {
exp = (10 * exp) + ((int32_t)(*p - '0'));
}
} else {
break;
}
}
if (!saw_exp_digits) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
dp += exp_sign * exp;
} while (0);
after_all:
if (p != q) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
h->num_digits = nd;
if (nd == 0) {
if (no_digits_before_separator) {
return wuffs_base__make_status(wuffs_base__error__bad_argument);
}
h->decimal_point = 0;
} else if (dp <
-WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
h->decimal_point =
-WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE - 1;
} else if (dp >
+WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
h->decimal_point =
+WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE + 1;
} else {
h->decimal_point = dp;
}
wuffs_base__private_implementation__high_prec_dec__trim(h);
return wuffs_base__make_status(NULL);
}
// --------
// wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits
// returns the number of additional decimal digits when left-shifting by shift.
//
// See below for preconditions.
static uint32_t //
wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits(
wuffs_base__private_implementation__high_prec_dec* h,
uint32_t shift) {
// Masking with 0x3F should be unnecessary (assuming the preconditions) but
// it's cheap and ensures that we don't overflow the
// wuffs_base__private_implementation__hpd_left_shift array.
shift &= 63;
uint32_t x_a = wuffs_base__private_implementation__hpd_left_shift[shift];
uint32_t x_b = wuffs_base__private_implementation__hpd_left_shift[shift + 1];
uint32_t num_new_digits = x_a >> 11;
uint32_t pow5_a = 0x7FF & x_a;
uint32_t pow5_b = 0x7FF & x_b;
const uint8_t* pow5 =
&wuffs_base__private_implementation__powers_of_5[pow5_a];
uint32_t i = 0;
uint32_t n = pow5_b - pow5_a;
for (; i < n; i++) {
if (i >= h->num_digits) {
return num_new_digits - 1;
} else if (h->digits[i] == pow5[i]) {
continue;
} else if (h->digits[i] < pow5[i]) {
return num_new_digits - 1;
} else {
return num_new_digits;
}
}
return num_new_digits;
}
// --------
// wuffs_base__private_implementation__high_prec_dec__rounded_integer returns
// the integral (non-fractional) part of h, provided that it is 18 or fewer
// decimal digits. For 19 or more digits, it returns UINT64_MAX. Note that:
// - (1 << 53) is 9007199254740992, which has 16 decimal digits.
// - (1 << 56) is 72057594037927936, which has 17 decimal digits.
// - (1 << 59) is 576460752303423488, which has 18 decimal digits.
// - (1 << 63) is 9223372036854775808, which has 19 decimal digits.
// and that IEEE 754 double precision has 52 mantissa bits.
//
// That integral part is rounded-to-even: rounding 7.5 or 8.5 both give 8.
//
// h's negative bit is ignored: rounding -8.6 returns 9.
//
// See below for preconditions.
static uint64_t //
wuffs_base__private_implementation__high_prec_dec__rounded_integer(
wuffs_base__private_implementation__high_prec_dec* h) {
if ((h->num_digits == 0) || (h->decimal_point < 0)) {
return 0;
} else if (h->decimal_point > 18) {
return UINT64_MAX;
}
uint32_t dp = (uint32_t)(h->decimal_point);
uint64_t n = 0;
uint32_t i = 0;
for (; i < dp; i++) {
n = (10 * n) + ((i < h->num_digits) ? h->digits[i] : 0);
}
bool round_up = false;
if (dp < h->num_digits) {
round_up = h->digits[dp] >= 5;
if ((h->digits[dp] == 5) && (dp + 1 == h->num_digits)) {
// We are exactly halfway. If we're truncated, round up, otherwise round
// to even.
round_up = h->truncated || //
((dp > 0) && (1 & h->digits[dp - 1]));
}
}
if (round_up) {
n++;
}
return n;
}
// wuffs_base__private_implementation__high_prec_dec__small_xshift shifts h's
// number (where 'x' is 'l' or 'r' for left or right) by a small shift value.
//
// Preconditions:
// - h is non-NULL.
// - h->decimal_point is "not extreme".
// - shift is non-zero.
// - shift is "a small shift".
//
// "Not extreme" means within
// ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
//
// "A small shift" means not more than
// WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL.
//
// wuffs_base__private_implementation__high_prec_dec__rounded_integer and
// wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits
// have the same preconditions.
//
// wuffs_base__private_implementation__high_prec_dec__lshift keeps the first
// two preconditions but not the last two. Its shift argument is signed and
// does not need to be "small": zero is a no-op, positive means left shift and
// negative means right shift.
static void //
wuffs_base__private_implementation__high_prec_dec__small_lshift(
wuffs_base__private_implementation__high_prec_dec* h,
uint32_t shift) {
if (h->num_digits == 0) {
return;
}
uint32_t num_new_digits =
wuffs_base__private_implementation__high_prec_dec__lshift_num_new_digits(
h, shift);
uint32_t rx = h->num_digits - 1; // Read index.
uint32_t wx = h->num_digits - 1 + num_new_digits; // Write index.
uint64_t n = 0;
// Repeat: pick up a digit, put down a digit, right to left.
while (((int32_t)rx) >= 0) {
n += ((uint64_t)(h->digits[rx])) << shift;
uint64_t quo = n / 10;
uint64_t rem = n - (10 * quo);
if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->digits[wx] = (uint8_t)rem;
} else if (rem > 0) {
h->truncated = true;
}
n = quo;
wx--;
rx--;
}
// Put down leading digits, right to left.
while (n > 0) {
uint64_t quo = n / 10;
uint64_t rem = n - (10 * quo);
if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->digits[wx] = (uint8_t)rem;
} else if (rem > 0) {
h->truncated = true;
}
n = quo;
wx--;
}
// Finish.
h->num_digits += num_new_digits;
if (h->num_digits >
WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->num_digits = WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION;
}
h->decimal_point += (int32_t)num_new_digits;
wuffs_base__private_implementation__high_prec_dec__trim(h);
}
static void //
wuffs_base__private_implementation__high_prec_dec__small_rshift(
wuffs_base__private_implementation__high_prec_dec* h,
uint32_t shift) {
uint32_t rx = 0; // Read index.
uint32_t wx = 0; // Write index.
uint64_t n = 0;
// Pick up enough leading digits to cover the first shift.
while ((n >> shift) == 0) {
if (rx < h->num_digits) {
// Read a digit.
n = (10 * n) + h->digits[rx++];
} else if (n == 0) {
// h's number used to be zero and remains zero.
return;
} else {
// Read sufficient implicit trailing zeroes.
while ((n >> shift) == 0) {
n = 10 * n;
rx++;
}
break;
}
}
h->decimal_point -= ((int32_t)(rx - 1));
if (h->decimal_point <
-WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE) {
// After the shift, h's number is effectively zero.
h->num_digits = 0;
h->decimal_point = 0;
h->negative = false;
h->truncated = false;
return;
}
// Repeat: pick up a digit, put down a digit, left to right.
uint64_t mask = (((uint64_t)(1)) << shift) - 1;
while (rx < h->num_digits) {
uint8_t new_digit = ((uint8_t)(n >> shift));
n = (10 * (n & mask)) + h->digits[rx++];
h->digits[wx++] = new_digit;
}
// Put down trailing digits, left to right.
while (n > 0) {
uint8_t new_digit = ((uint8_t)(n >> shift));
n = 10 * (n & mask);
if (wx < WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DIGITS_PRECISION) {
h->digits[wx++] = new_digit;
} else if (new_digit > 0) {
h->truncated = true;
}
}
// Finish.
h->num_digits = wx;
wuffs_base__private_implementation__high_prec_dec__trim(h);
}
static void //
wuffs_base__private_implementation__high_prec_dec__lshift(
wuffs_base__private_implementation__high_prec_dec* h,
int32_t shift) {
if (shift > 0) {
while (shift > +WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
wuffs_base__private_implementation__high_prec_dec__small_lshift(
h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL);
shift -= WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
}
wuffs_base__private_implementation__high_prec_dec__small_lshift(
h, ((uint32_t)(+shift)));
} else if (shift < 0) {
while (shift < -WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL) {
wuffs_base__private_implementation__high_prec_dec__small_rshift(
h, WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL);
shift += WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__SHIFT__MAX_INCL;
}
wuffs_base__private_implementation__high_prec_dec__small_rshift(
h, ((uint32_t)(-shift)));
}
}
// --------
// wuffs_base__private_implementation__high_prec_dec__round_etc rounds h's
// number. For those functions that take an n argument, rounding produces at
// most n digits (which is not necessarily at most n decimal places). Negative
// n values are ignored, as well as any n greater than or equal to h's number
// of digits. The etc__round_just_enough function implicitly chooses an n to
// implement WUFFS_BASE__RENDER_NUMBER_FXX__JUST_ENOUGH_PRECISION.
//
// Preconditions:
// - h is non-NULL.
// - h->decimal_point is "not extreme".
//
// "Not extreme" means within
// ±WUFFS_BASE__PRIVATE_IMPLEMENTATION__HPD__DECIMAL_POINT__RANGE.
static void //
wuffs_base__private_implementation__high_prec_dec__round_down(
wuffs_base__private_implementation__high_prec_dec* h,
int32_t n) {
if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
return;
}
h->num_digits = (uint32_t)(n);
wuffs_base__private_implementation__high_prec_dec__trim(h);
}
static void //
wuffs_base__private_implementation__high_prec_dec__round_up(
wuffs_base__private_implementation__high_prec_dec* h,
int32_t n) {
if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
return;
}
for (n--; n >= 0; n--) {
if (h->digits[n] < 9) {
h->digits[n]++;
h->num_digits = (uint32_t)(n + 1);
return;
}
}
// The number is all 9s. Change to a single 1 and adjust the decimal point.
h->digits[0] = 1;
h->num_digits = 1;
h->decimal_point++;
}
static void //
wuffs_base__private_implementation__high_prec_dec__round_nearest(
wuffs_base__private_implementation__high_prec_dec* h,
int32_t n) {
if ((n < 0) || (h->num_digits <= (uint32_t)n)) {
return;
}
bool up = h->digits[n] >= 5;
if ((h->digits[n] == 5) && ((n + 1) == ((int32_t)(h->num_digits)))) {
up = h->truncated || //
((n > 0) && ((h->digits[n - 1] & 1) != 0));
}
if (up) {
wuffs_base__private_implementation__high_prec_dec__round_up(h, n);
} else {
wuffs_base__private_implementation__high_prec_dec__round_down(h, n);
}
}
static void //
wuffs_base__private_implementation__high_prec_dec__round_just_enough(
wuffs_base__private_implementation__high_prec_dec* h,
int32_t exp2,
uint64_t mantissa) {
// The magic numbers 52 and 53 in this function are because IEEE 754 double
// precision has 52 mantissa bits.
//
// Let f be the floating point number represented by exp2 and mantissa (and
// also the number in h): the number (mantissa * (2 ** (exp2 - 52))).
//
// If f is zero or a small integer, we can return early.
if ((mantissa == 0) ||
((exp2 < 53) && (h->decimal_point >= ((int32_t)(h->num_digits))))) {
return;
}
// The smallest normal f has an exp2 of -1022 and a mantissa of (1 << 52).
// Subnormal numbers have the same exp2 but a smaller mantissa.
static const int32_t min_incl_normal_exp2 = -1022;
static const uint64_t min_incl_normal_mantissa = 0x0010000000000000ul;
// Compute lower and upper bounds such that any number between them (possibly
// inclusive) will round to f. First, the lower bound. Our number f is:
// ((mantissa + 0) * (2 ** ( exp2 - 52)))
//
// The next lowest floating point number is:
// ((mantissa - 1) * (2 ** ( exp2 - 52)))
// unless (mantissa - 1) drops the (1 << 52) bit and exp2 is not the
// min_incl_normal_exp2. Either way, call it:
// ((l_mantissa) * (2 ** (l_exp2 - 52)))
//
// The lower bound is halfway between them (noting that 52 became 53):
// (((2 * l_mantissa) + 1) * (2 ** (l_exp2 - 53)))
int32_t l_exp2 = exp2;
uint64_t l_mantissa = mantissa - 1;
if ((exp2 > min_incl_normal_exp2) && (mantissa <= min_incl_normal_mantissa)) {
l_exp2 = exp2 - 1;
l_mantissa = (2 * mantissa) - 1;
}
wuffs_base__private_implementation__high_prec_dec lower;
wuffs_base__private_implementation__high_prec_dec__assign(
&lower, (2 * l_mantissa) + 1, false);
wuffs_base__private_implementation__high_prec_dec__lshift(&lower,
l_exp2 - 53);
// Next, the upper bound. Our number f is:
// ((mantissa + 0) * (2 ** (exp2 - 52)))
//
// The next highest floating point number is:
// ((mantissa + 1) * (2 ** (exp2 - 52)))
//
// The upper bound is halfway between them (noting that 52 became 53):
// (((2 * mantissa) + 1) * (2 ** (exp2 - 53)))
wuffs_base__private_implementation__high_prec_dec upper;
wuffs_base__private_implementation__high_prec_dec__assign(
&upper, (2 * mantissa) + 1, false);
wuffs_base__private_implementation__high_prec_dec__lshift(&upper, exp2 - 53);
// The lower and upper bounds are possible outputs only if the original
// mantissa is even, so that IEEE round-to-even would round to the original
// mantissa and not its neighbors.
bool inclusive = (mantissa & 1) == 0;
// As we walk the digits, we want to know whether rounding up would fall
// within the upper bound. This is tracked by upper_delta:
// - When -1, the digits of h and upper are the same so far.
// - When +0, we saw a difference of 1 between h and upper on a previous
// digit and subsequently only 9s for h and 0s for upper. Thus, rounding
// up may fall outside of the bound if !inclusive.
// - When +1, the difference is greater than 1 and we know that rounding up
// falls within the bound.
//
// This is a state machine with three states. The numerical value for each
// state (-1, +0 or +1) isn't important, other than their order.
int upper_delta = -1;
// We can now figure out the shortest number of digits required. Walk the
// digits until h has distinguished itself from lower or upper.
//
// The zi and zd variables are indexes and digits, for z in l (lower), h (the
// number) and u (upper).
//
// The lower, h and upper numbers may have their decimal points at different
// places. In this case, upper is the longest, so we iterate ui starting from
// 0 and iterate li and hi starting from either 0 or -1.
int32_t ui = 0;
for (;; ui++) {
// Calculate hd, the middle number's digit.
int32_t hi = ui - upper.decimal_point + h->decimal_point;
if (hi >= ((int32_t)(h->num_digits))) {
break;
}
uint8_t hd = (((uint32_t)hi) < h->num_digits) ? h->digits[hi] : 0;
// Calculate ld, the lower bound's digit.
int32_t li = ui - upper.decimal_point + lower.decimal_point;
uint8_t ld = (((uint32_t)li) < lower.num_digits) ? lower.digits[li] : 0;
// We can round down (truncate) if lower has a different digit than h or if
// lower is inclusive and is exactly the result of rounding down (i.e. we
// have reached the final digit of lower).
bool can_round_down =
(ld != hd) || //
(inclusive && ((li + 1) == ((int32_t)(lower.num_digits))));
// Calculate ud, the upper bound's digit, and update upper_delta.
uint8_t ud = (((uint32_t)ui) < upper.num_digits) ? upper.digits[ui] : 0;
if (upper_delta < 0) {
if ((hd + 1) < ud) {
// For example:
// h = 12345???
// upper = 12347???
upper_delta = +1;
} else if (hd != ud) {
// For example:
// h = 12345???
// upper = 12346???
upper_delta = +0;
}
} else if (upper_delta == 0) {
if ((hd != 9) || (ud != 0)) {
// For example:
// h = 1234598?
// upper = 1234600?
upper_delta = +1;
}
}
// We can round up if upper has a different digit than h and either upper
// is inclusive or upper is bigger than the result of rounding up.
bool can_round_up =
(upper_delta > 0) || //
((upper_delta == 0) && //
(inclusive || ((ui + 1) < ((int32_t)(upper.num_digits)))));
// If we can round either way, round to nearest. If we can round only one
// way, do it. If we can't round, continue the loop.
if (can_round_down) {
if (can_round_up) {
wuffs_base__private_implementation__high_prec_dec__round_nearest(
h, hi + 1);
return;
} else {
wuffs_base__private_implementation__high_prec_dec__round_down(h,
hi + 1);
return;
}
} else {
if (can_round_up) {
wuffs_base__private_implementation__high_prec_dec__round_up(h, hi + 1);
return;
}
}
}
}
// --------
// wuffs_base__private_implementation__parse_number_f64_eisel_lemire produces
// the IEEE 754 double-precision value for an exact mantissa and base-10
// exponent. For example:
// - when parsing "12345.678e+02", man is 12345678 and exp10 is -1.
// - when parsing "-12", man is 12 and exp10 is 0. Processing the leading
// minus sign is the responsibility of the caller, not this function.
//
// On success, it returns a non-negative int64_t such that the low 63 bits hold
// the 11-bit exponent and 52-bit mantissa.
//
// On failure, it returns a negative value.
//
// The algorithm is based on an original idea by Michael Eisel that was refined
// by Daniel Lemire. See
// https://lemire.me/blog/2020/03/10/fast-float-parsing-in-practice/
// and
// https://nigeltao.github.io/blog/2020/eisel-lemire.html
//
// Preconditions:
// - man is non-zero.
// - exp10 is in the range [-307 ..= 288], the same range of the
// wuffs_base__private_implementation__powers_of_10 array.
//
// The exp10 range (and the fact that man is in the range [1 ..= UINT64_MAX],
// approximately [1 ..= 1.85e+19]) means that (man * (10 ** exp10)) is in the
// range [1e-307 ..= 1.85e+307]. This is entirely within the range of normal
// (neither subnormal nor non-finite) f64 values: DBL_MIN and DBL_MAX are
// approximately 2.23e–308 and 1.80e+308.
static int64_t //
wuffs_base__private_implementation__parse_number_f64_eisel_lemire(
uint64_t man,
int32_t exp10) {
// Look up the (possibly truncated) base-2 representation of (10 ** exp10).
// The look-up table was constructed so that it is already normalized: the
// table entry's mantissa's MSB (most significant bit) is on.
const uint64_t* po10 =
&wuffs_base__private_implementation__powers_of_10[exp10 + 307][0];