/[pspp]/pspp/src/pfm-write.c
ViewVC logotype

Diff of /pspp/src/pfm-write.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.12 by blp, Mon May 2 06:21:20 2005 UTC revision 1.13 by blp, Mon Jul 25 03:32:18 2005 UTC
# Line 32  Line 32 
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"
# Line 67  static int write_version_data (struct pf Line 67  static int write_version_data (struct pf
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. */
# Line 162  buf_write (struct pfm_writer *w, const v Line 165  buf_write (struct pfm_writer *w, const v
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. */
# Line 497  pfm_close_writer (struct pfm_writer *w) Line 411  pfm_close_writer (struct pfm_writer *w)
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    }

Legend:
Removed from v.1.12  
changed lines
  Added in v.1.13

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26