32 |
#include "dictionary.h" |
#include "dictionary.h" |
33 |
#include "error.h" |
#include "error.h" |
34 |
#include "file-handle.h" |
#include "file-handle.h" |
|
#include "gmp.h" |
|
35 |
#include "hash.h" |
#include "hash.h" |
36 |
#include "magic.h" |
#include "magic.h" |
37 |
|
#include "misc.h" |
38 |
#include "str.h" |
#include "str.h" |
39 |
#include "value-labels.h" |
#include "value-labels.h" |
40 |
#include "var.h" |
#include "var.h" |
67 |
static int write_variables (struct pfm_writer *, struct dictionary *); |
static int write_variables (struct pfm_writer *, struct dictionary *); |
68 |
static int write_value_labels (struct pfm_writer *, const struct dictionary *); |
static int write_value_labels (struct pfm_writer *, const struct dictionary *); |
69 |
|
|
70 |
|
static void format_trig_double (long double, int base_10_precision, char[]); |
71 |
|
static char *format_trig_int (int, bool force_sign, char[]); |
72 |
|
|
73 |
/* Writes the dictionary DICT to portable file HANDLE. Returns |
/* Writes the dictionary DICT to portable file HANDLE. Returns |
74 |
nonzero only if successful. DICT will not be modified, except |
nonzero only if successful. DICT will not be modified, except |
75 |
to assign short names. */ |
to assign short names. */ |
165 |
static int |
static int |
166 |
write_float (struct pfm_writer *w, double d) |
write_float (struct pfm_writer *w, double d) |
167 |
{ |
{ |
168 |
int neg = 0; |
char buffer[64]; |
169 |
char *mantissa; |
format_trig_double (d, DBL_DIG, buffer); |
170 |
int mantissa_len; |
return buf_write (w, buffer, strlen (buffer)) && buf_write (w, "/", 1); |
|
mp_exp_t exponent; |
|
|
char *buf, *cp; |
|
|
int success; |
|
|
|
|
|
if (d < 0.) |
|
|
{ |
|
|
d = -d; |
|
|
neg = 1; |
|
|
} |
|
|
|
|
|
if (d == fabs (SYSMIS) || d == HUGE_VAL) |
|
|
return buf_write (w, "*.", 2); |
|
|
|
|
|
/* Use GNU libgmp2 to convert D into base-30. */ |
|
|
{ |
|
|
mpf_t f; |
|
|
|
|
|
mpf_init_set_d (f, d); |
|
|
mantissa = mpf_get_str (NULL, &exponent, 30, 0, f); |
|
|
mpf_clear (f); |
|
|
|
|
|
for (cp = mantissa; *cp; cp++) |
|
|
*cp = toupper (*cp); |
|
|
} |
|
|
|
|
|
/* Choose standard or scientific notation. */ |
|
|
mantissa_len = (int) strlen (mantissa); |
|
|
cp = buf = local_alloc (mantissa_len + 32); |
|
|
if (neg) |
|
|
*cp++ = '-'; |
|
|
if (mantissa_len == 0) |
|
|
*cp++ = '0'; |
|
|
else if (exponent < -4 || exponent > (mp_exp_t) mantissa_len) |
|
|
{ |
|
|
/* Scientific notation. */ |
|
|
*cp++ = mantissa[0]; |
|
|
*cp++ = '.'; |
|
|
cp = stpcpy (cp, &mantissa[1]); |
|
|
cp = spprintf (cp, "%+ld", (long) (exponent - 1)); |
|
|
} |
|
|
else if (exponent <= 0) |
|
|
{ |
|
|
/* Standard notation, D <= 1. */ |
|
|
*cp++ = '.'; |
|
|
memset (cp, '0', -exponent); |
|
|
cp += -exponent; |
|
|
cp = stpcpy (cp, mantissa); |
|
|
} |
|
|
else |
|
|
{ |
|
|
/* Standard notation, D > 1. */ |
|
|
memcpy (cp, mantissa, exponent); |
|
|
cp += exponent; |
|
|
*cp++ = '.'; |
|
|
cp = stpcpy (cp, &mantissa[exponent]); |
|
|
} |
|
|
*cp++ = '/'; |
|
|
|
|
|
success = buf_write (w, buf, cp - buf); |
|
|
local_free (buf); |
|
|
free (mantissa); |
|
|
return success; |
|
171 |
} |
} |
172 |
|
|
173 |
/* Write N to the portable file as an integer field, and return success. */ |
/* Write N to the portable file as an integer field, and return success. */ |
174 |
static int |
static int |
175 |
write_int (struct pfm_writer *w, int n) |
write_int (struct pfm_writer *w, int n) |
176 |
{ |
{ |
177 |
char buf[64]; |
char buffer[64]; |
178 |
char *bp = &buf[64]; |
format_trig_int (n, false, buffer); |
179 |
int neg = 0; |
return buf_write (w, buffer, strlen (buffer)) && buf_write (w, "/", 1); |
|
|
|
|
*--bp = '/'; |
|
|
|
|
|
if (n < 0) |
|
|
{ |
|
|
n = -n; |
|
|
neg = 1; |
|
|
} |
|
|
|
|
|
do |
|
|
{ |
|
|
int r = n % 30; |
|
|
|
|
|
/* PORTME: character codes. */ |
|
|
if (r < 10) |
|
|
*--bp = r + '0'; |
|
|
else |
|
|
*--bp = r - 10 + 'A'; |
|
|
|
|
|
n /= 30; |
|
|
} |
|
|
while (n > 0); |
|
|
|
|
|
if (neg) |
|
|
*--bp = '-'; |
|
|
|
|
|
return buf_write (w, bp, &buf[64] - bp); |
|
180 |
} |
} |
181 |
|
|
182 |
/* Write S to the portable file as a string field. */ |
/* Write S to the portable file as a string field. */ |
411 |
free (w->vars); |
free (w->vars); |
412 |
free (w); |
free (w); |
413 |
} |
} |
414 |
|
|
415 |
|
/* Base-30 conversion. */ |
416 |
|
|
417 |
|
/* Conversion base. */ |
418 |
|
#define BASE 30 /* As an integer. */ |
419 |
|
#define LDBASE ((long double) BASE) /* As a long double. */ |
420 |
|
|
421 |
|
/* This is floor(log30(2**31)), the minimum number of trigesimal |
422 |
|
digits that a `long int' can hold. */ |
423 |
|
#define CHUNK_SIZE 6 |
424 |
|
|
425 |
|
/* Yields the square of X. */ |
426 |
|
#define Q(X) ((X) * (X)) |
427 |
|
|
428 |
|
/* Returns 30**EXPONENT, for 0 <= EXPONENT <= log30(DBL_MAX). */ |
429 |
|
static long double |
430 |
|
pow30_nonnegative (int exponent) |
431 |
|
{ |
432 |
|
/* pow_tab[i] = pow (30, pow (2, i)) */ |
433 |
|
static const long double pow_tab[] = |
434 |
|
{ |
435 |
|
LDBASE, |
436 |
|
Q (LDBASE), |
437 |
|
Q (Q (LDBASE)), |
438 |
|
Q (Q (Q (LDBASE))), |
439 |
|
Q (Q (Q (Q (LDBASE)))), |
440 |
|
Q (Q (Q (Q (Q (LDBASE))))), |
441 |
|
Q (Q (Q (Q (Q (Q (LDBASE)))))), |
442 |
|
Q (Q (Q (Q (Q (Q (Q (LDBASE))))))), |
443 |
|
Q (Q (Q (Q (Q (Q (Q (Q (LDBASE)))))))), |
444 |
|
Q (Q (Q (Q (Q (Q (Q (Q (Q (LDBASE))))))))), |
445 |
|
Q (Q (Q (Q (Q (Q (Q (Q (Q (Q (LDBASE)))))))))), |
446 |
|
Q (Q (Q (Q (Q (Q (Q (Q (Q (Q (Q (LDBASE))))))))))), |
447 |
|
}; |
448 |
|
|
449 |
|
long double power; |
450 |
|
int i; |
451 |
|
|
452 |
|
assert (exponent >= 0); |
453 |
|
assert (exponent < 1L << (sizeof pow_tab / sizeof *pow_tab)); |
454 |
|
|
455 |
|
power = 1.L; |
456 |
|
for (i = 0; exponent > 0; exponent >>= 1, i++) |
457 |
|
if (exponent & 1) |
458 |
|
power *= pow_tab[i]; |
459 |
|
|
460 |
|
return power; |
461 |
|
} |
462 |
|
|
463 |
|
/* Returns 30**EXPONENT, for log30(DBL_MIN) <= EXPONENT <= |
464 |
|
log30(DBL_MAX). */ |
465 |
|
static long double |
466 |
|
pow30 (int exponent) |
467 |
|
{ |
468 |
|
if (exponent >= 0) |
469 |
|
return pow30_nonnegative (exponent); |
470 |
|
else |
471 |
|
return 1.L / pow30_nonnegative (-exponent); |
472 |
|
} |
473 |
|
|
474 |
|
/* Returns the character corresponding to TRIG. */ |
475 |
|
static int |
476 |
|
trig_to_char (int trig) |
477 |
|
{ |
478 |
|
assert (trig >= 0 && trig < 30); |
479 |
|
return "0123456789ABCDEFGHIJKLMNOPQRST"[trig]; |
480 |
|
} |
481 |
|
|
482 |
|
/* Formats the TRIG_CNT trigs in TRIGS[], writing them as |
483 |
|
null-terminated STRING. The trigesimal point is inserted |
484 |
|
after TRIG_PLACES characters have been printed, if necessary |
485 |
|
adding extra zeros at either end for correctness. Returns the |
486 |
|
character after the formatted number. */ |
487 |
|
static char * |
488 |
|
format_trig_digits (char *string, |
489 |
|
const char trigs[], int trig_cnt, int trig_places) |
490 |
|
{ |
491 |
|
if (trig_places < 0) |
492 |
|
{ |
493 |
|
*string++ = '.'; |
494 |
|
while (trig_places++ < 0) |
495 |
|
*string++ = '0'; |
496 |
|
trig_places = -1; |
497 |
|
} |
498 |
|
while (trig_cnt-- > 0) |
499 |
|
{ |
500 |
|
if (trig_places-- == 0) |
501 |
|
*string++ = '.'; |
502 |
|
*string++ = trig_to_char (*trigs++); |
503 |
|
} |
504 |
|
while (trig_places-- > 0) |
505 |
|
*string++ = '0'; |
506 |
|
*string = '\0'; |
507 |
|
return string; |
508 |
|
} |
509 |
|
|
510 |
|
/* Helper function for format_trig_int() that formats VALUE as a |
511 |
|
trigesimal integer at CP. VALUE must be nonnegative. |
512 |
|
Returns the character following the formatted integer. */ |
513 |
|
static char * |
514 |
|
recurse_format_trig_int (char *cp, int value) |
515 |
|
{ |
516 |
|
int trig = value % BASE; |
517 |
|
value /= BASE; |
518 |
|
if (value > 0) |
519 |
|
cp = recurse_format_trig_int (cp, value); |
520 |
|
*cp++ = trig_to_char (trig); |
521 |
|
return cp; |
522 |
|
} |
523 |
|
|
524 |
|
/* Formats VALUE as a trigesimal integer in null-terminated |
525 |
|
STRING[]. VALUE must be in the range -DBL_MAX...DBL_MAX. If |
526 |
|
FORCE_SIGN is true, a sign is always inserted; otherwise, a |
527 |
|
sign is only inserted if VALUE is negative. */ |
528 |
|
static char * |
529 |
|
format_trig_int (int value, bool force_sign, char string[]) |
530 |
|
{ |
531 |
|
/* Insert sign. */ |
532 |
|
if (value < 0) |
533 |
|
{ |
534 |
|
*string++ = '-'; |
535 |
|
value = -value; |
536 |
|
} |
537 |
|
else if (force_sign) |
538 |
|
*string++ = '+'; |
539 |
|
|
540 |
|
/* Format integer. */ |
541 |
|
string = recurse_format_trig_int (string, value); |
542 |
|
*string = '\0'; |
543 |
|
return string; |
544 |
|
} |
545 |
|
|
546 |
|
/* Determines whether the TRIG_CNT trigesimals in TRIGS[] warrant |
547 |
|
rounding up or down. Returns true if TRIGS[] represents a |
548 |
|
value greater than half, false if less than half. If TRIGS[] |
549 |
|
is exactly half, examines TRIGS[-1] and returns true if odd, |
550 |
|
false if even ("round to even"). */ |
551 |
|
static bool |
552 |
|
should_round_up (const char trigs[], int trig_cnt) |
553 |
|
{ |
554 |
|
assert (trig_cnt > 0); |
555 |
|
|
556 |
|
if (*trigs < BASE / 2) |
557 |
|
{ |
558 |
|
/* Less than half: round down. */ |
559 |
|
return false; |
560 |
|
} |
561 |
|
else if (*trigs > BASE / 2) |
562 |
|
{ |
563 |
|
/* Greater than half: round up. */ |
564 |
|
return true; |
565 |
|
} |
566 |
|
else |
567 |
|
{ |
568 |
|
/* Approximately half: look more closely. */ |
569 |
|
int i; |
570 |
|
for (i = 1; i < trig_cnt; i++) |
571 |
|
if (trigs[i] > 0) |
572 |
|
{ |
573 |
|
/* Slightly greater than half: round up. */ |
574 |
|
return true; |
575 |
|
} |
576 |
|
|
577 |
|
/* Exactly half: round to even. */ |
578 |
|
return trigs[-1] % 2; |
579 |
|
} |
580 |
|
} |
581 |
|
|
582 |
|
/* Rounds up the rightmost trig in the TRIG_CNT trigs in TRIGS[], |
583 |
|
carrying to the left as necessary. Returns true if |
584 |
|
successful, false on failure (due to a carry out of the |
585 |
|
leftmost position). */ |
586 |
|
static bool |
587 |
|
try_round_up (char *trigs, int trig_cnt) |
588 |
|
{ |
589 |
|
while (trig_cnt > 0) |
590 |
|
{ |
591 |
|
char *round_trig = trigs + --trig_cnt; |
592 |
|
if (*round_trig != BASE - 1) |
593 |
|
{ |
594 |
|
/* Round this trig up to the next value. */ |
595 |
|
++*round_trig; |
596 |
|
return true; |
597 |
|
} |
598 |
|
|
599 |
|
/* Carry over to the next trig to the left. */ |
600 |
|
*round_trig = 0; |
601 |
|
} |
602 |
|
|
603 |
|
/* Ran out of trigs to carry. */ |
604 |
|
return false; |
605 |
|
} |
606 |
|
|
607 |
|
/* Converts VALUE to trigesimal format in string OUTPUT[] with the |
608 |
|
equivalent of at least BASE_10_PRECISION decimal digits of |
609 |
|
precision. The output format may use conventional or |
610 |
|
scientific notation. Missing, infinite, and extreme values |
611 |
|
are represented with "*.". */ |
612 |
|
static void |
613 |
|
format_trig_double (long double value, int base_10_precision, char output[]) |
614 |
|
{ |
615 |
|
/* Original VALUE was negative? */ |
616 |
|
bool negative; |
617 |
|
|
618 |
|
/* Number of significant trigesimals. */ |
619 |
|
int base_30_precision; |
620 |
|
|
621 |
|
/* Base-2 significand and exponent for original VALUE. */ |
622 |
|
double base_2_sig; |
623 |
|
int base_2_exp; |
624 |
|
|
625 |
|
/* VALUE as a set of trigesimals. */ |
626 |
|
char buffer[DBL_DIG + 16]; |
627 |
|
char *trigs; |
628 |
|
int trig_cnt; |
629 |
|
|
630 |
|
/* Number of trigesimal places for trigs. |
631 |
|
trigs[0] has coefficient 30**(trig_places - 1), |
632 |
|
trigs[1] has coefficient 30**(trig_places - 2), |
633 |
|
and so on. |
634 |
|
In other words, the trigesimal point is just before trigs[0]. |
635 |
|
*/ |
636 |
|
int trig_places; |
637 |
|
|
638 |
|
/* Number of trigesimal places left to write into BUFFER. */ |
639 |
|
int trigs_to_output; |
640 |
|
|
641 |
|
/* Handle special cases. */ |
642 |
|
if (value == SYSMIS) |
643 |
|
goto missing_value; |
644 |
|
if (value == 0.) |
645 |
|
goto zero; |
646 |
|
|
647 |
|
/* Make VALUE positive. */ |
648 |
|
if (value < 0) |
649 |
|
{ |
650 |
|
value = -value; |
651 |
|
negative = true; |
652 |
|
} |
653 |
|
else |
654 |
|
negative = false; |
655 |
|
|
656 |
|
/* Adjust VALUE to roughly 30**3, by shifting the trigesimal |
657 |
|
point left or right as necessary. We approximate the |
658 |
|
base-30 exponent by obtaining the base-2 exponent, then |
659 |
|
multiplying by log30(2). This approximation is sufficient |
660 |
|
to ensure that the adjusted VALUE is always in the range |
661 |
|
0...30**6, an invariant of the loop below. */ |
662 |
|
errno = 0; |
663 |
|
base_2_sig = frexp (value, &base_2_exp); |
664 |
|
if (errno != 0 || !finite (base_2_sig)) |
665 |
|
goto missing_value; |
666 |
|
if (base_2_exp == 0 && base_2_sig == 0.) |
667 |
|
goto zero; |
668 |
|
if (base_2_exp <= INT_MIN / 20379L || base_2_exp >= INT_MAX / 20379L) |
669 |
|
goto missing_value; |
670 |
|
trig_places = (base_2_exp * 20379L / 100000L) + CHUNK_SIZE / 2; |
671 |
|
value *= pow30 (CHUNK_SIZE - trig_places); |
672 |
|
|
673 |
|
/* Dump all the trigs to buffer[], CHUNK_SIZE at a time. */ |
674 |
|
trigs = buffer; |
675 |
|
trig_cnt = 0; |
676 |
|
for (trigs_to_output = DIV_RND_UP (DBL_DIG * 2, 3) + 1 + (CHUNK_SIZE / 2); |
677 |
|
trigs_to_output > 0; |
678 |
|
trigs_to_output -= CHUNK_SIZE) |
679 |
|
{ |
680 |
|
long chunk; |
681 |
|
int trigs_left; |
682 |
|
|
683 |
|
/* The current chunk is just the integer part of VALUE, |
684 |
|
truncated to the nearest integer. The chunk fits in a |
685 |
|
long. */ |
686 |
|
chunk = value; |
687 |
|
assert (pow30 (CHUNK_SIZE) <= LONG_MAX); |
688 |
|
assert (chunk >= 0 && chunk < pow30 (CHUNK_SIZE)); |
689 |
|
|
690 |
|
value -= chunk; |
691 |
|
|
692 |
|
/* Append the chunk, in base 30, to trigs[]. */ |
693 |
|
for (trigs_left = CHUNK_SIZE; chunk > 0 && trigs_left > 0; ) |
694 |
|
{ |
695 |
|
trigs[trig_cnt + --trigs_left] = chunk % 30; |
696 |
|
chunk /= 30; |
697 |
|
} |
698 |
|
while (trigs_left > 0) |
699 |
|
trigs[trig_cnt + --trigs_left] = 0; |
700 |
|
trig_cnt += CHUNK_SIZE; |
701 |
|
|
702 |
|
/* Proceed to the next chunk. */ |
703 |
|
if (value == 0.) |
704 |
|
break; |
705 |
|
value *= pow (LDBASE, CHUNK_SIZE); |
706 |
|
} |
707 |
|
|
708 |
|
/* Strip leading zeros. */ |
709 |
|
while (trig_cnt > 1 && *trigs == 0) |
710 |
|
{ |
711 |
|
trigs++; |
712 |
|
trig_cnt--; |
713 |
|
trig_places--; |
714 |
|
} |
715 |
|
|
716 |
|
/* Round to requested precision, conservatively estimating the |
717 |
|
required base-30 precision as 2/3 of the base-10 precision |
718 |
|
(log30(10) = .68). */ |
719 |
|
assert (base_10_precision > 0); |
720 |
|
base_30_precision = DIV_RND_UP (base_10_precision * 2, 3); |
721 |
|
if (trig_cnt > base_30_precision) |
722 |
|
{ |
723 |
|
if (should_round_up (trigs + base_30_precision, |
724 |
|
trig_cnt - base_30_precision)) |
725 |
|
{ |
726 |
|
/* Try to round up. */ |
727 |
|
if (try_round_up (trigs, base_30_precision)) |
728 |
|
{ |
729 |
|
/* Rounding up worked. */ |
730 |
|
trig_cnt = base_30_precision; |
731 |
|
} |
732 |
|
else |
733 |
|
{ |
734 |
|
/* Couldn't round up because we ran out of trigs to |
735 |
|
carry into. Do the carry here instead. */ |
736 |
|
*trigs = 1; |
737 |
|
trig_cnt = 1; |
738 |
|
trig_places++; |
739 |
|
} |
740 |
|
} |
741 |
|
else |
742 |
|
{ |
743 |
|
/* Round down. */ |
744 |
|
trig_cnt = base_30_precision; |
745 |
|
} |
746 |
|
} |
747 |
|
else |
748 |
|
{ |
749 |
|
/* No rounding required: fewer digits available than |
750 |
|
requested. */ |
751 |
|
} |
752 |
|
|
753 |
|
/* Strip trailing zeros. */ |
754 |
|
while (trig_cnt > 1 && trigs[trig_cnt - 1] == 0) |
755 |
|
trig_cnt--; |
756 |
|
|
757 |
|
/* Write output. */ |
758 |
|
if (negative) |
759 |
|
*output++ = '-'; |
760 |
|
if (trig_places >= -1 && trig_places < trig_cnt + 3) |
761 |
|
{ |
762 |
|
/* Use conventional notation. */ |
763 |
|
format_trig_digits (output, trigs, trig_cnt, trig_places); |
764 |
|
} |
765 |
|
else |
766 |
|
{ |
767 |
|
/* Use scientific notation. */ |
768 |
|
char *op; |
769 |
|
op = format_trig_digits (output, trigs, trig_cnt, trig_cnt); |
770 |
|
op = format_trig_int (trig_places - trig_cnt, true, op); |
771 |
|
} |
772 |
|
return; |
773 |
|
|
774 |
|
zero: |
775 |
|
strcpy (output, "0"); |
776 |
|
return; |
777 |
|
|
778 |
|
missing_value: |
779 |
|
strcpy (output, "*."); |
780 |
|
return; |
781 |
|
} |