Get rid of dependency on libgmp by writing our own routine for
[pspp-builds.git] / src / pfm-write.c
index e229072ea5578678e0bfe51c3ac2042e7872c2e5..a4b80bd079321d4b669809ec6cbd7e66d8c273f7 100644 (file)
 
    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
-   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
-   02111-1307, USA. */
+   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+   02110-1301, USA. */
 
 #include <config.h>
-#include "pfm.h"
-#include <assert.h>
+#include "pfm-write.h"
+#include "error.h"
 #include <ctype.h>
 #include <errno.h>
 #include <float.h>
 #include <stdlib.h>
 #include <time.h>
 #include "alloc.h"
+#include "case.h"
+#include "dictionary.h"
 #include "error.h"
 #include "file-handle.h"
-#include "gmp.h"
 #include "hash.h"
 #include "magic.h"
+#include "misc.h"
 #include "str.h"
 #include "value-labels.h"
 #include "var.h"
 
 #include "debug-print.h"
 
-/* pfm writer file_handle extension. */
-struct pfm_fhuser_ext
+/* Portable file writer. */
+struct pfm_writer
   {
-    FILE *file;                        /* Actual file. */
+    struct file_handle *fh;     /* File handle. */
+    FILE *file;                        /* File stream. */
 
     int lc;                    /* Number of characters on this line so far. */
 
-    int nvars;                 /* Number of variables. */
-    int *vars;                 /* Variable widths. */
+    size_t var_cnt;             /* Number of variables. */
+    struct pfm_var *vars;       /* Variables. */
   };
 
-static struct fh_ext_class pfm_w_class;
+/* A variable to write to the portable file. */
+struct pfm_var 
+  {
+    int width;                  /* 0=numeric, otherwise string var width. */
+    int fv;                     /* Starting case index. */
+  };
 
-static int bufwrite (struct file_handle *h, const void *buf, size_t nbytes);
-static int write_header (struct file_handle *h);
-static int write_version_data (struct file_handle *h);
-static int write_variables (struct file_handle *h, struct dictionary *d);
-static int write_value_labels (struct file_handle *h, struct dictionary *d);
+static int buf_write (struct pfm_writer *, const void *, size_t);
+static int write_header (struct pfm_writer *);
+static int write_version_data (struct pfm_writer *);
+static int write_variables (struct pfm_writer *, struct dictionary *);
+static int write_value_labels (struct pfm_writer *, const struct dictionary *);
+
+static void format_trig_double (long double, int base_10_precision, char[]);
+static char *format_trig_int (int, bool force_sign, char[]);
 
 /* Writes the dictionary DICT to portable file HANDLE.  Returns
-   nonzero only if successful. */
-int
-pfm_write_dictionary (struct file_handle *handle, struct dictionary *dict)
+   nonzero only if successful.  DICT will not be modified, except
+   to assign short names. */
+struct pfm_writer *
+pfm_open_writer (struct file_handle *fh, struct dictionary *dict)
 {
-  struct pfm_fhuser_ext *ext;
-  
-  if (handle->class != NULL)
-    {
-      msg (ME, _("Cannot write file %s as portable file: already opened "
-                "for %s."),
-          handle_get_name (handle), handle->class->name);
-      return 0;
-    }
+  struct pfm_writer *w = NULL;
+  size_t i;
 
-  msg (VM (1), _("%s: Opening portable-file handle %s for writing."),
-       handle_get_filename (handle), handle_get_name (handle));
+  if (!fh_open (fh, "portable file", "we"))
+    goto error;
   
   /* Open the physical disk file. */
-  handle->class = &pfm_w_class;
-  handle->ext = ext = xmalloc (sizeof (struct pfm_fhuser_ext));
-  ext->file = fopen (handle_get_filename (handle), "wb");
-  ext->lc = 0;
-  if (ext->file == NULL)
+  w = xmalloc (sizeof *w);
+  w->fh = fh;
+  w->file = fopen (handle_get_filename (fh), "wb");
+  w->lc = 0;
+  w->var_cnt = 0;
+  w->vars = NULL;
+  
+  /* Check that file create succeeded. */
+  if (w->file == NULL)
     {
       msg (ME, _("An error occurred while opening \"%s\" for writing "
           "as a portable file: %s."),
-           handle_get_filename (handle), strerror (errno));
+           handle_get_filename (fh), strerror (errno));
       err_cond_fail ();
-      free (ext);
-      return 0;
+      goto error;
     }
   
-  {
-    int i;
-
-    ext->nvars = dict_get_var_cnt (dict);
-    ext->vars = xmalloc (sizeof *ext->vars * ext->nvars);
-    for (i = 0; i < ext->nvars; i++)
-      ext->vars[i] = dict_get_var (dict, i)->width;
-  }
-
-  /* Write the file header. */
-  if (!write_header (handle))
-    goto lossage;
-
-  /* Write version data. */
-  if (!write_version_data (handle))
-    goto lossage;
-
-  /* Write variables. */
-  if (!write_variables (handle, dict))
-    goto lossage;
-
-  /* Write value labels. */
-  if (!write_value_labels (handle, dict))
-    goto lossage;
+  w->var_cnt = dict_get_var_cnt (dict);
+  w->vars = xmalloc (sizeof *w->vars * w->var_cnt);
+  for (i = 0; i < w->var_cnt; i++) 
+    {
+      const struct variable *dv = dict_get_var (dict, i);
+      struct pfm_var *pv = &w->vars[i];
+      pv->width = dv->width;
+      pv->fv = dv->fv;
+    }
 
-  /* Write beginning of data marker. */
-  if (!bufwrite (handle, "F", 1))
-    goto lossage;
+  /* Write file header. */
+  if (!write_header (w)
+      || !write_version_data (w)
+      || !write_variables (w, dict)
+      || !write_value_labels (w, dict)
+      || !buf_write (w, "F", 1))
+    goto error;
 
-  msg (VM (2), _("Wrote portable-file header successfully."));
+  return w;
 
-  return 1;
-
-lossage:
-  msg (VM (1), _("Error writing portable-file header."));
-  fclose (ext->file);
-  free (ext->vars);
-  handle->class = NULL;
-  handle->ext = NULL;
-  return 0;
+error:
+  pfm_close_writer (w);
+  return NULL;
 }
 \f  
 /* Write NBYTES starting at BUF to the portable file represented by
    H.  Break lines properly every 80 characters.  */
 static int
-bufwrite (struct file_handle *h, const void *buf_, size_t nbytes)
+buf_write (struct pfm_writer *w, const void *buf_, size_t nbytes)
 {
   const char *buf = buf_;
-  struct pfm_fhuser_ext *ext = h->ext;
 
   assert (buf != NULL);
-  while (nbytes + ext->lc >= 80)
+  while (nbytes + w->lc >= 80)
     {
-      size_t n = 80 - ext->lc;
+      size_t n = 80 - w->lc;
       
-      if (n && 1 != fwrite (buf, n, 1, ext->file))
-       goto lossage;
+      if (n && fwrite (buf, n, 1, w->file) != 1)
+       goto error;
       
-      /* PORTME: line ends. */
-      if (1 != fwrite ("\r\n", 2, 1, ext->file))
-       goto lossage;
+      if (fwrite ("\r\n", 2, 1, w->file) != 1)
+       goto error;
 
       nbytes -= n;
       buf += n;
-      ext->lc = 0;
+      w->lc = 0;
     }
 
-  if (nbytes && 1 != fwrite (buf, nbytes, 1, ext->file))
-    goto lossage;
-  ext->lc += nbytes;
+  if (nbytes && 1 != fwrite (buf, nbytes, 1, w->file))
+    goto error;
+  w->lc += nbytes;
   
   return 1;
 
- lossage:
-  abort ();
+ error:
   msg (ME, _("%s: Writing portable file: %s."),
-       handle_get_filename (h), strerror (errno));
+       handle_get_filename (w->fh), strerror (errno));
   return 0;
 }
 
 /* Write D to the portable file as a floating-point field, and return
    success. */
 static int
-write_float (struct file_handle *h, double d)
+write_float (struct pfm_writer *w, double d)
 {
-  int neg = 0;
-  char *mantissa;
-  int mantissa_len;
-  mp_exp_t exponent;
-  char *buf, *cp;
-  int success;
-
-  if (d < 0.)
-    {
-      d = -d;
-      neg = 1;
-    }
-  
-  if (d == fabs (SYSMIS) || d == HUGE_VAL)
-    return bufwrite (h, "*.", 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 = bufwrite (h, buf, cp - buf);
-  local_free (buf);
-  free (mantissa);
-  return success;
+  char buffer[64];
+  format_trig_double (d, DBL_DIG, buffer);
+  return buf_write (w, buffer, strlen (buffer)) && buf_write (w, "/", 1);
 }
 
 /* Write N to the portable file as an integer field, and return success. */
 static int
-write_int (struct file_handle *h, int n)
+write_int (struct pfm_writer *w, int n)
 {
-  char buf[64];
-  char *bp = &buf[64];
-  int neg = 0;
-
-  *--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 bufwrite (h, bp, &buf[64] - bp);
+  char buffer[64];
+  format_trig_int (n, false, buffer);
+  return buf_write (w, buffer, strlen (buffer)) && buf_write (w, "/", 1);
 }
 
 /* Write S to the portable file as a string field. */
 static int
-write_string (struct file_handle *h, const char *s)
+write_string (struct pfm_writer *w, const char *s)
 {
   size_t n = strlen (s);
-  return write_int (h, (int) n) && bufwrite (h, s, n);
+  return write_int (w, (int) n) && buf_write (w, s, n);
 }
 \f
 /* Write file header. */
 static int
-write_header (struct file_handle *h)
+write_header (struct pfm_writer *w)
 {
   /* PORTME. */
   {
     int i;
 
     for (i = 0; i < 5; i++)
-      if (!bufwrite (h, "ASCII SPSS PORT FILE                    ", 40))
+      if (!buf_write (w, "ASCII SPSS PORT FILE                    ", 40))
        return 0;
   }
   
@@ -312,11 +211,11 @@ write_header (struct file_handle *h)
        "0000000000000000000000000000000000000000000000000000000000000000"
       };
 
-    if (!bufwrite (h, spss2ascii, 256))
+    if (!buf_write (w, spss2ascii, 256))
       return 0;
   }
 
-  if (!bufwrite (h, "SPSSPORT", 8))
+  if (!buf_write (w, "SPSSPORT", 8))
     return 0;
 
   return 1;
@@ -324,9 +223,9 @@ write_header (struct file_handle *h)
 
 /* Writes version, date, and identification records. */
 static int
-write_version_data (struct file_handle *h)
+write_version_data (struct pfm_writer *w)
 {
-  if (!bufwrite (h, "A", 1))
+  if (!buf_write (w, "A", 1))
     return 0;
   
   {
@@ -348,16 +247,16 @@ write_version_data (struct file_handle *h)
     sprintf (date_str, "%04d%02d%02d",
             tmp->tm_year + 1900, tmp->tm_mon + 1, tmp->tm_mday);
     sprintf (time_str, "%02d%02d%02d", tmp->tm_hour, tmp->tm_min, tmp->tm_sec);
-    if (!write_string (h, date_str) || !write_string (h, time_str))
+    if (!write_string (w, date_str) || !write_string (w, time_str))
       return 0;
   }
 
   /* Product identification. */
-  if (!bufwrite (h, "1", 1) || !write_string (h, version))
+  if (!buf_write (w, "1", 1) || !write_string (w, version))
     return 0;
 
   /* Subproduct identification. */
-  if (!bufwrite (h, "3", 1) || !write_string (h, host_system))
+  if (!buf_write (w, "3", 1) || !write_string (w, host_system))
     return 0;
 
   return 1;
@@ -365,31 +264,33 @@ write_version_data (struct file_handle *h)
 
 /* Write format F to file H, and return success. */
 static int
-write_format (struct file_handle *h, struct fmt_spec *f)
+write_format (struct pfm_writer *w, struct fmt_spec *f)
 {
-  return (write_int (h, formats[f->type].spss)
-         && write_int (h, f->w)
-         && write_int (h, f->d));
+  return (write_int (w, formats[f->type].spss)
+         && write_int (w, f->w)
+         && write_int (w, f->d));
 }
 
 /* Write value V for variable VV to file H, and return success. */
 static int
-write_value (struct file_handle *h, union value *v, struct variable *vv)
+write_value (struct pfm_writer *w, union value *v, struct variable *vv)
 {
   if (vv->type == NUMERIC)
-    return write_float (h, v->f);
+    return write_float (w, v->f);
   else
-    return write_int (h, vv->width) && bufwrite (h, v->s, vv->width);
+    return write_int (w, vv->width) && buf_write (w, v->s, vv->width);
 }
 
 /* Write variable records, and return success. */
 static int
-write_variables (struct file_handle *h, struct dictionary *dict)
+write_variables (struct pfm_writer *w, struct dictionary *dict)
 {
   int i;
+
+  dict_assign_short_names (dict);
   
-  if (!bufwrite (h, "4", 1) || !write_int (h, dict_get_var_cnt (dict))
-      || !write_int (h, 161))
+  if (!buf_write (w, "4", 1) || !write_int (w, dict_get_var_cnt (dict))
+      || !write_int (w, 161))
     return 0;
 
   for (i = 0; i < dict_get_var_cnt (dict); i++)
@@ -404,17 +305,17 @@ write_variables (struct file_handle *h, struct dictionary *dict)
 
       struct variable *v = dict_get_var (dict, i);
       
-      if (!bufwrite (h, "7", 1) || !write_int (h, v->width)
-         || !write_string (h, v->name)
-         || !write_format (h, &v->print) || !write_format (h, &v->write))
+      if (!buf_write (w, "7", 1) || !write_int (w, v->width)
+         || !write_string (w, v->short_name)
+         || !write_format (w, &v->print) || !write_format (w, &v->write))
        return 0;
 
       for (m = miss_types[v->miss_type], j = 0; j < (int) strlen (m); j++)
-       if ((m[j] != ' ' && !bufwrite (h, &m[j], 1))
-           || !write_value (h, &v->missing[j], v))
+       if ((m[j] != ' ' && !buf_write (w, &m[j], 1))
+           || !write_value (w, &v->missing[j], v))
          return 0;
 
-      if (v->label && (!bufwrite (h, "C", 1) || !write_string (h, v->label)))
+      if (v->label && (!buf_write (w, "C", 1) || !write_string (w, v->label)))
        return 0;
     }
 
@@ -423,7 +324,7 @@ write_variables (struct file_handle *h, struct dictionary *dict)
 
 /* Write value labels to disk.  FIXME: Inefficient. */
 static int
-write_value_labels (struct file_handle *h, struct dictionary *dict)
+write_value_labels (struct pfm_writer *w, const struct dictionary *dict)
 {
   int i;
 
@@ -436,16 +337,16 @@ write_value_labels (struct file_handle *h, struct dictionary *dict)
       if (!val_labs_count (v->val_labs))
        continue;
 
-      if (!bufwrite (h, "D", 1)
-         || !write_int (h, 1)
-         || !write_string (h, v->name)
-         || !write_int (h, val_labs_count (v->val_labs)))
+      if (!buf_write (w, "D", 1)
+         || !write_int (w, 1)
+         || !write_string (w, v->short_name)
+         || !write_int (w, val_labs_count (v->val_labs)))
        return 0;
 
       for (vl = val_labs_first_sorted (v->val_labs, &j); vl != NULL;
            vl = val_labs_next (v->val_labs, &j)) 
-       if (!write_value (h, &vl->value, v)
-           || !write_string (h, vl->label)) 
+       if (!write_value (w, &vl->value, v)
+           || !write_string (w, vl->label)) 
           {
             val_labs_done (&j);
             return 0; 
@@ -458,24 +359,23 @@ write_value_labels (struct file_handle *h, struct dictionary *dict)
 /* Writes case ELEM to the portable file represented by H.  Returns
    success. */
 int 
-pfm_write_case (struct file_handle *h, const union value *elem)
+pfm_write_case (struct pfm_writer *w, struct ccase *c)
 {
-  struct pfm_fhuser_ext *ext = h->ext;
-  
   int i;
   
-  for (i = 0; i < ext->nvars; i++)
+  for (i = 0; i < w->var_cnt; i++)
     {
-      const int width = ext->vars[i];
+      struct pfm_var *v = &w->vars[i];
       
-      if (width == 0)
+      if (v->width == 0)
        {
-         if (!write_float (h, elem[i].f))
+         if (!write_float (w, case_num (c, v->fv)))
            return 0;
        }
       else
        {
-         if (!write_int (h, width) || !bufwrite (h, elem[i].c, width))
+         if (!write_int (w, v->width)
+              || !buf_write (w, case_str (c, v->fv), v->width))
            return 0;
        }
     }
@@ -484,33 +384,398 @@ pfm_write_case (struct file_handle *h, const union value *elem)
 }
 
 /* Closes a portable file after we're done with it. */
-static void
-pfm_close (struct file_handle *h)
+void
+pfm_close_writer (struct pfm_writer *w)
 {
-  struct pfm_fhuser_ext *ext = h->ext;
+  if (w == NULL)
+    return;
+
+  fh_close (w->fh, "portable file", "we");
   
-  {
-    char buf[80];
+  if (w->file != NULL)
+    {
+      char buf[80];
     
-    int n = 80 - ext->lc;
-    if (n == 0)
-      n = 80;
+      int n = 80 - w->lc;
+      if (n == 0)
+        n = 80;
 
-    memset (buf, 'Z', n);
-    bufwrite (h, buf, n);
-  }
+      memset (buf, 'Z', n);
+      buf_write (w, buf, n);
+
+      if (fclose (w->file) == EOF)
+        msg (ME, _("%s: Closing portable file: %s."),
+             handle_get_filename (w->fh), strerror (errno));
+    }
+
+  free (w->vars);
+  free (w);
+}
+\f
+/* Base-30 conversion. */
+
+/* Conversion base. */
+#define BASE 30                         /* As an integer. */
+#define LDBASE ((long double) BASE)     /* As a long double. */
+
+/* This is floor(log30(2**31)), the minimum number of trigesimal
+   digits that a `long int' can hold. */
+#define CHUNK_SIZE 6                    
+
+/* Yields the square of X. */
+#define Q(X) ((X) * (X))
+
+/* Returns 30**EXPONENT, for 0 <= EXPONENT <= log30(DBL_MAX). */
+static long double
+pow30_nonnegative (int exponent)
+{
+  /* pow_tab[i] = pow (30, pow (2, i)) */
+  static const long double pow_tab[] =
+    {
+      LDBASE,
+      Q (LDBASE),
+      Q (Q (LDBASE)),
+      Q (Q (Q (LDBASE))),
+      Q (Q (Q (Q (LDBASE)))),
+      Q (Q (Q (Q (Q (LDBASE))))),
+      Q (Q (Q (Q (Q (Q (LDBASE)))))),
+      Q (Q (Q (Q (Q (Q (Q (LDBASE))))))),
+      Q (Q (Q (Q (Q (Q (Q (Q (LDBASE)))))))),
+      Q (Q (Q (Q (Q (Q (Q (Q (Q (LDBASE))))))))),
+      Q (Q (Q (Q (Q (Q (Q (Q (Q (Q (LDBASE)))))))))),
+      Q (Q (Q (Q (Q (Q (Q (Q (Q (Q (Q (LDBASE))))))))))),
+    };
+
+  long double power;
+  int i;
+
+  assert (exponent >= 0);
+  assert (exponent < 1L << (sizeof pow_tab / sizeof *pow_tab));
 
-  if (EOF == fclose (ext->file))
-    msg (ME, _("%s: Closing portable file: %s."),
-         handle_get_filename (h), strerror (errno));
+  power = 1.L;
+  for (i = 0; exponent > 0; exponent >>= 1, i++)
+    if (exponent & 1)
+      power *= pow_tab[i];
 
-  free (ext->vars);
-  free (ext);
+  return power;
 }
 
-static struct fh_ext_class pfm_w_class =
+/* Returns 30**EXPONENT, for log30(DBL_MIN) <= EXPONENT <=
+   log30(DBL_MAX). */
+static long double
+pow30 (int exponent)
 {
-  6,
-  N_("writing as a portable file"),
-  pfm_close,
-};
+  if (exponent >= 0)
+    return pow30_nonnegative (exponent);
+  else
+    return 1.L / pow30_nonnegative (-exponent);
+}
+
+/* Returns the character corresponding to TRIG. */
+static int
+trig_to_char (int trig)
+{
+  assert (trig >= 0 && trig < 30);
+  return "0123456789ABCDEFGHIJKLMNOPQRST"[trig];
+}
+
+/* Formats the TRIG_CNT trigs in TRIGS[], writing them as
+   null-terminated STRING.  The trigesimal point is inserted
+   after TRIG_PLACES characters have been printed, if necessary
+   adding extra zeros at either end for correctness.  Returns the
+   character after the formatted number. */
+static char *
+format_trig_digits (char *string,
+                    const char trigs[], int trig_cnt, int trig_places)
+{
+  if (trig_places < 0)
+    {
+      *string++ = '.';
+      while (trig_places++ < 0)
+        *string++ = '0';
+      trig_places = -1;
+    }
+  while (trig_cnt-- > 0)
+    {
+      if (trig_places-- == 0)
+        *string++ = '.';
+      *string++ = trig_to_char (*trigs++);
+    }
+  while (trig_places-- > 0)
+    *string++ = '0';
+  *string = '\0';
+  return string;
+}
+
+/* Helper function for format_trig_int() that formats VALUE as a
+   trigesimal integer at CP.  VALUE must be nonnegative.
+   Returns the character following the formatted integer. */
+static char *
+recurse_format_trig_int (char *cp, int value)
+{
+  int trig = value % BASE;
+  value /= BASE;
+  if (value > 0)
+    cp = recurse_format_trig_int (cp, value);
+  *cp++ = trig_to_char (trig);
+  return cp;
+}
+
+/* Formats VALUE as a trigesimal integer in null-terminated
+   STRING[].  VALUE must be in the range -DBL_MAX...DBL_MAX.  If
+   FORCE_SIGN is true, a sign is always inserted; otherwise, a
+   sign is only inserted if VALUE is negative. */
+static char *
+format_trig_int (int value, bool force_sign, char string[])
+{
+  /* Insert sign. */
+  if (value < 0)
+    {
+      *string++ = '-';
+      value = -value;
+    }
+  else if (force_sign)
+    *string++ = '+';
+
+  /* Format integer. */
+  string = recurse_format_trig_int (string, value);
+  *string = '\0';
+  return string;
+}
+
+/* Determines whether the TRIG_CNT trigesimals in TRIGS[] warrant
+   rounding up or down.  Returns true if TRIGS[] represents a
+   value greater than half, false if less than half.  If TRIGS[]
+   is exactly half, examines TRIGS[-1] and returns true if odd,
+   false if even ("round to even"). */
+static bool
+should_round_up (const char trigs[], int trig_cnt)
+{
+  assert (trig_cnt > 0);
+
+  if (*trigs < BASE / 2)
+    {
+      /* Less than half: round down. */
+      return false;
+    }
+  else if (*trigs > BASE / 2)
+    {
+      /* Greater than half: round up. */
+      return true;
+    }
+  else
+    {
+      /* Approximately half: look more closely. */
+      int i;
+      for (i = 1; i < trig_cnt; i++)
+        if (trigs[i] > 0)
+          {
+            /* Slightly greater than half: round up. */
+            return true;
+          }
+
+      /* Exactly half: round to even. */
+      return trigs[-1] % 2;
+    }
+}
+
+/* Rounds up the rightmost trig in the TRIG_CNT trigs in TRIGS[],
+   carrying to the left as necessary.  Returns true if
+   successful, false on failure (due to a carry out of the
+   leftmost position). */
+static bool
+try_round_up (char *trigs, int trig_cnt)
+{
+  while (trig_cnt > 0)
+    {
+      char *round_trig = trigs + --trig_cnt;
+      if (*round_trig != BASE - 1)
+        {
+          /* Round this trig up to the next value. */
+          ++*round_trig;
+          return true;
+        }
+
+      /* Carry over to the next trig to the left. */
+      *round_trig = 0;
+    }
+
+  /* Ran out of trigs to carry. */
+  return false;
+}
+
+/* Converts VALUE to trigesimal format in string OUTPUT[] with the
+   equivalent of at least BASE_10_PRECISION decimal digits of
+   precision.  The output format may use conventional or
+   scientific notation.  Missing, infinite, and extreme values
+   are represented with "*.". */
+static void
+format_trig_double (long double value, int base_10_precision, char output[])
+{
+  /* Original VALUE was negative? */
+  bool negative;
+
+  /* Number of significant trigesimals. */
+  int base_30_precision;
+
+  /* Base-2 significand and exponent for original VALUE. */
+  double base_2_sig;
+  int base_2_exp;
+
+  /* VALUE as a set of trigesimals. */
+  char buffer[DBL_DIG + 16];
+  char *trigs;
+  int trig_cnt;
+
+  /* Number of trigesimal places for trigs.
+     trigs[0] has coefficient 30**(trig_places - 1),
+     trigs[1] has coefficient 30**(trig_places - 2),
+     and so on.
+     In other words, the trigesimal point is just before trigs[0].
+   */
+  int trig_places;
+
+  /* Number of trigesimal places left to write into BUFFER. */
+  int trigs_to_output;
+
+  /* Handle special cases. */
+  if (value == SYSMIS)
+    goto missing_value;
+  if (value == 0.)
+    goto zero;
+
+  /* Make VALUE positive. */
+  if (value < 0)
+    {
+      value = -value;
+      negative = true;
+    }
+  else
+    negative = false;
+
+  /* Adjust VALUE to roughly 30**3, by shifting the trigesimal
+     point left or right as necessary.  We approximate the
+     base-30 exponent by obtaining the base-2 exponent, then
+     multiplying by log30(2).  This approximation is sufficient
+     to ensure that the adjusted VALUE is always in the range
+     0...30**6, an invariant of the loop below. */
+  errno = 0;
+  base_2_sig = frexp (value, &base_2_exp);
+  if (errno != 0 || !finite (base_2_sig))
+    goto missing_value;
+  if (base_2_exp == 0 && base_2_sig == 0.)
+    goto zero;
+  if (base_2_exp <= INT_MIN / 20379L || base_2_exp >= INT_MAX / 20379L)
+    goto missing_value;
+  trig_places = (base_2_exp * 20379L / 100000L) + CHUNK_SIZE / 2;
+  value *= pow30 (CHUNK_SIZE - trig_places);
+
+  /* Dump all the trigs to buffer[], CHUNK_SIZE at a time. */
+  trigs = buffer;
+  trig_cnt = 0;
+  for (trigs_to_output = DIV_RND_UP (DBL_DIG * 2, 3) + 1 + (CHUNK_SIZE / 2);
+       trigs_to_output > 0;
+       trigs_to_output -= CHUNK_SIZE)
+    {
+      long chunk;
+      int trigs_left;
+
+      /* The current chunk is just the integer part of VALUE,
+         truncated to the nearest integer.  The chunk fits in a
+         long. */
+      chunk = value;
+      assert (pow30 (CHUNK_SIZE) <= LONG_MAX);
+      assert (chunk >= 0 && chunk < pow30 (CHUNK_SIZE));
+
+      value -= chunk;
+
+      /* Append the chunk, in base 30, to trigs[]. */
+      for (trigs_left = CHUNK_SIZE; chunk > 0 && trigs_left > 0; )
+        {
+          trigs[trig_cnt + --trigs_left] = chunk % 30;
+          chunk /= 30;
+        }
+      while (trigs_left > 0)
+        trigs[trig_cnt + --trigs_left] = 0;
+      trig_cnt += CHUNK_SIZE;
+
+      /* Proceed to the next chunk. */
+      if (value == 0.)
+        break;
+      value *= pow (LDBASE, CHUNK_SIZE);
+    }
+
+  /* Strip leading zeros. */
+  while (trig_cnt > 1 && *trigs == 0)
+    {
+      trigs++;
+      trig_cnt--;
+      trig_places--;
+    }
+
+  /* Round to requested precision, conservatively estimating the
+     required base-30 precision as 2/3 of the base-10 precision
+     (log30(10) = .68). */
+  assert (base_10_precision > 0);
+  base_30_precision = DIV_RND_UP (base_10_precision * 2, 3);
+  if (trig_cnt > base_30_precision)
+    {
+      if (should_round_up (trigs + base_30_precision,
+                           trig_cnt - base_30_precision))
+        {
+          /* Try to round up. */
+          if (try_round_up (trigs, base_30_precision))
+            {
+              /* Rounding up worked. */
+              trig_cnt = base_30_precision;
+            }
+          else
+            {
+              /* Couldn't round up because we ran out of trigs to
+                 carry into.  Do the carry here instead. */
+              *trigs = 1;
+              trig_cnt = 1;
+              trig_places++;
+            }
+        }
+      else
+        {
+          /* Round down. */
+          trig_cnt = base_30_precision;
+        }
+    }
+  else
+    {
+      /* No rounding required: fewer digits available than
+         requested. */
+    }
+
+  /* Strip trailing zeros. */
+  while (trig_cnt > 1 && trigs[trig_cnt - 1] == 0)
+    trig_cnt--;
+
+  /* Write output. */
+  if (negative)
+    *output++ = '-';
+  if (trig_places >= -1 && trig_places < trig_cnt + 3)
+    {
+      /* Use conventional notation. */
+      format_trig_digits (output, trigs, trig_cnt, trig_places);
+    }
+  else
+    {
+      /* Use scientific notation. */
+      char *op;
+      op = format_trig_digits (output, trigs, trig_cnt, trig_cnt);
+      op = format_trig_int (trig_places - trig_cnt, true, op);
+    }
+  return;
+
+ zero:
+  strcpy (output, "0");
+  return;
+
+ missing_value:
+  strcpy (output, "*.");
+  return;
+}