becc2edbd03fdb3299975e36a1ef70e14c174b37
[pspp-builds.git] / src / language / stats / crosstabs.q
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3    Written by Ben Pfaff <blp@gnu.org>.
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18    02110-1301, USA. */
19
20 /* FIXME:
21
22    - Pearson's R (but not Spearman!) is off a little.
23    - T values for Spearman's R and Pearson's R are wrong.
24    - How to calculate significance of symmetric and directional measures?
25    - Asymmetric ASEs and T values for lambda are wrong.
26    - ASE of Goodman and Kruskal's tau is not calculated.
27    - ASE of symmetric somers' d is wrong.
28    - Approx. T of uncertainty coefficient is wrong.
29
30 */
31
32 #include <config.h>
33
34 #include <ctype.h>
35 #include <gsl/gsl_cdf.h>
36 #include <stdlib.h>
37 #include <stdio.h>
38
39 #include <data/case.h>
40 #include <data/dictionary.h>
41 #include <data/procedure.h>
42 #include <data/value-labels.h>
43 #include <data/variable.h>
44 #include <language/command.h>
45 #include <language/dictionary/split-file.h>
46 #include <language/lexer/lexer.h>
47 #include <language/lexer/variable-parser.h>
48 #include <libpspp/alloc.h>
49 #include <libpspp/array.h>
50 #include <libpspp/assertion.h>
51 #include <libpspp/compiler.h>
52 #include <libpspp/hash.h>
53 #include <libpspp/magic.h>
54 #include <libpspp/message.h>
55 #include <libpspp/message.h>
56 #include <libpspp/misc.h>
57 #include <libpspp/pool.h>
58 #include <libpspp/str.h>
59 #include <output/output.h>
60 #include <output/table.h>
61
62 #include "gettext.h"
63 #define _(msgid) gettext (msgid)
64 #define N_(msgid) msgid
65
66 /* (headers) */
67
68 /* (specification)
69    crosstabs (crs_):
70      *^tables=custom;
71      +variables=custom;
72      missing=miss:!table/include/report;
73      +write[wr_]=none,cells,all;
74      +format=fmt:!labels/nolabels/novallabs,
75              val:!avalue/dvalue,
76              indx:!noindex/index,
77              tabl:!tables/notables,
78              box:!box/nobox,
79              pivot:!pivot/nopivot;
80      +cells[cl_]=count,none,expected,row,column,total,residual,sresidual,
81                  asresidual,all;
82      +statistics[st_]=chisq,phi,cc,lambda,uc,none,btau,ctau,risk,gamma,d,
83                       kappa,eta,corr,all.
84 */
85 /* (declarations) */
86 /* (functions) */
87
88 /* Number of chi-square statistics. */
89 #define N_CHISQ 5
90
91 /* Number of symmetric statistics. */
92 #define N_SYMMETRIC 9
93
94 /* Number of directional statistics. */
95 #define N_DIRECTIONAL 13
96
97 /* A single table entry for general mode. */
98 struct table_entry
99   {
100     int table;          /* Flattened table number. */
101     union
102       {
103         double freq;    /* Frequency count. */
104         double *data;   /* Crosstabulation table for integer mode. */
105       }
106     u;
107     union value values[1];      /* Values. */
108   };
109
110 /* A crosstabulation. */
111 struct crosstab
112   {
113     int nvar;                   /* Number of variables. */
114     double missing;             /* Missing cases count. */
115     int ofs;                    /* Integer mode: Offset into sorted_tab[]. */
116     struct variable *vars[2];   /* At least two variables; sorted by
117                                    larger indices first. */
118   };
119
120 /* Integer mode variable info. */
121 struct var_range
122   {
123     int min;                    /* Minimum value. */
124     int max;                    /* Maximum value + 1. */
125     int count;                  /* max - min. */
126   };
127
128 static inline struct var_range *
129 get_var_range (struct variable *v) 
130 {
131   assert (v != NULL);
132   assert (v->aux != NULL);
133   return v->aux;
134 }
135
136 /* Indexes into crosstab.v. */
137 enum
138   {
139     ROW_VAR = 0,
140     COL_VAR = 1
141   };
142
143 /* General mode crosstabulation table. */
144 static struct hsh_table *gen_tab;       /* Hash table. */
145 static int n_sorted_tab;                /* Number of entries in sorted_tab. */
146 static struct table_entry **sorted_tab; /* Sorted table. */
147
148 /* Variables specifies on VARIABLES. */
149 static struct variable **variables;
150 static size_t variables_cnt;
151
152 /* TABLES. */
153 static struct crosstab **xtab;
154 static int nxtab;
155
156 /* Integer or general mode? */
157 enum
158   {
159     INTEGER,
160     GENERAL
161   };
162 static int mode;
163
164 /* CELLS. */
165 static int num_cells;           /* Number of cells requested. */
166 static int cells[8];            /* Cells requested. */
167
168 /* WRITE. */
169 static int write;               /* One of WR_* that specifies the WRITE style. */
170
171 /* Command parsing info. */
172 static struct cmd_crosstabs cmd;
173
174 /* Pools. */
175 static struct pool *pl_tc;      /* For table cells. */
176 static struct pool *pl_col;     /* For column data. */
177
178 static int internal_cmd_crosstabs (struct dataset *ds);
179 static void precalc (const struct ccase *, void *, const struct dataset *);
180 static bool calc_general (const struct ccase *, void *, const struct dataset *);
181 static bool calc_integer (const struct ccase *, void *, const struct dataset *);
182 static bool postcalc (void *, const struct dataset *);
183 static void submit (struct tab_table *);
184
185 static void format_short (char *s, const struct fmt_spec *fp,
186                          const union value *v);
187
188 /* Parse and execute CROSSTABS, then clean up. */
189 int
190 cmd_crosstabs (struct dataset *ds)
191 {
192   int result = internal_cmd_crosstabs (ds);
193
194   free (variables);
195   pool_destroy (pl_tc);
196   pool_destroy (pl_col);
197   
198   return result;
199 }
200
201 /* Parses and executes the CROSSTABS procedure. */
202 static int
203 internal_cmd_crosstabs (struct dataset *ds)
204 {
205   int i;
206   bool ok;
207
208   variables = NULL;
209   variables_cnt = 0;
210   xtab = NULL;
211   nxtab = 0;
212   pl_tc = pool_create ();
213   pl_col = pool_create ();
214
215   if (!parse_crosstabs (ds, &cmd, NULL))
216     return CMD_FAILURE;
217
218   mode = variables ? INTEGER : GENERAL;
219
220   /* CELLS. */
221   if (!cmd.sbc_cells)
222     {
223       cmd.a_cells[CRS_CL_COUNT] = 1;
224     }
225   else 
226     {
227       int count = 0;
228
229       for (i = 0; i < CRS_CL_count; i++)
230         if (cmd.a_cells[i])
231           count++;
232       if (count == 0)
233         {
234           cmd.a_cells[CRS_CL_COUNT] = 1;
235           cmd.a_cells[CRS_CL_ROW] = 1;
236           cmd.a_cells[CRS_CL_COLUMN] = 1;
237           cmd.a_cells[CRS_CL_TOTAL] = 1;
238         }
239       if (cmd.a_cells[CRS_CL_ALL])
240         {
241           for (i = 0; i < CRS_CL_count; i++)
242             cmd.a_cells[i] = 1;
243           cmd.a_cells[CRS_CL_ALL] = 0;
244         }
245       cmd.a_cells[CRS_CL_NONE] = 0;
246     }
247   for (num_cells = i = 0; i < CRS_CL_count; i++)
248     if (cmd.a_cells[i])
249       cells[num_cells++] = i;
250
251   /* STATISTICS. */
252   if (cmd.sbc_statistics)
253     {
254       int i;
255       int count = 0;
256
257       for (i = 0; i < CRS_ST_count; i++)
258         if (cmd.a_statistics[i])
259           count++;
260       if (count == 0)
261         cmd.a_statistics[CRS_ST_CHISQ] = 1;
262       if (cmd.a_statistics[CRS_ST_ALL])
263         for (i = 0; i < CRS_ST_count; i++)
264           cmd.a_statistics[i] = 1;
265     }
266   
267   /* MISSING. */
268   if (cmd.miss == CRS_REPORT && mode == GENERAL)
269     {
270       msg (SE, _("Missing mode REPORT not allowed in general mode.  "
271                  "Assuming MISSING=TABLE."));
272       cmd.miss = CRS_TABLE;
273     }
274
275   /* WRITE. */
276   if (cmd.a_write[CRS_WR_ALL] && cmd.a_write[CRS_WR_CELLS])
277     cmd.a_write[CRS_WR_ALL] = 0;
278   if (cmd.a_write[CRS_WR_ALL] && mode == GENERAL)
279     {
280       msg (SE, _("Write mode ALL not allowed in general mode.  "
281                  "Assuming WRITE=CELLS."));
282       cmd.a_write[CRS_WR_CELLS] = 1;
283     }
284   if (cmd.sbc_write
285       && (cmd.a_write[CRS_WR_NONE]
286           + cmd.a_write[CRS_WR_ALL]
287           + cmd.a_write[CRS_WR_CELLS] == 0))
288     cmd.a_write[CRS_WR_CELLS] = 1;
289   if (cmd.a_write[CRS_WR_CELLS])
290     write = CRS_WR_CELLS;
291   else if (cmd.a_write[CRS_WR_ALL])
292     write = CRS_WR_ALL;
293   else
294     write = CRS_WR_NONE;
295
296   ok = procedure_with_splits (ds, precalc,
297                               mode == GENERAL ? calc_general : calc_integer,
298                               postcalc, NULL);
299
300   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
301 }
302
303 /* Parses the TABLES subcommand. */
304 static int
305 crs_custom_tables (struct dataset *ds, struct cmd_crosstabs *cmd UNUSED, void *aux UNUSED)
306 {
307   struct var_set *var_set;
308   int n_by;
309   struct variable ***by = NULL;
310   size_t *by_nvar = NULL;
311   size_t nx = 1;
312   int success = 0;
313
314   /* Ensure that this is a TABLES subcommand. */
315   if (!lex_match_id ("TABLES")
316       && (token != T_ID || dict_lookup_var (dataset_dict (ds), tokid) == NULL)
317       && token != T_ALL)
318     return 2;
319   lex_match ('=');
320
321   if (variables != NULL)
322     var_set = var_set_create_from_array (variables, variables_cnt);
323   else
324     var_set = var_set_create_from_dict (dataset_dict (ds));
325   assert (var_set != NULL);
326   
327   for (n_by = 0; ;)
328     {
329       by = xnrealloc (by, n_by + 1, sizeof *by);
330       by_nvar = xnrealloc (by_nvar, n_by + 1, sizeof *by_nvar);
331       if (!parse_var_set_vars (var_set, &by[n_by], &by_nvar[n_by],
332                                PV_NO_DUPLICATE | PV_NO_SCRATCH))
333         goto done;
334       if (xalloc_oversized (nx, by_nvar[n_by])) 
335         {
336           msg (SE, _("Too many crosstabulation variables or dimensions."));
337           goto done;
338         }
339       nx *= by_nvar[n_by];
340       n_by++;
341
342       if (!lex_match (T_BY))
343         {
344           if (n_by < 2)
345             {
346               lex_error (_("expecting BY"));
347               goto done;
348             }
349           else 
350             break;
351         }
352     }
353   
354   {
355     int *by_iter = xcalloc (n_by, sizeof *by_iter);
356     int i;
357
358     xtab = xnrealloc (xtab, nxtab + nx, sizeof *xtab);
359     for (i = 0; i < nx; i++)
360       {
361         struct crosstab *x;
362
363         x = xmalloc (sizeof *x + sizeof (struct variable *) * (n_by - 2));
364         x->nvar = n_by;
365         x->missing = 0.;
366
367         {
368           int i;
369
370           for (i = 0; i < n_by; i++)
371             x->vars[i] = by[i][by_iter[i]];
372         }
373         
374         {
375           int i;
376
377           for (i = n_by - 1; i >= 0; i--)
378             {
379               if (++by_iter[i] < by_nvar[i])
380                 break;
381               by_iter[i] = 0;
382             }
383         }
384
385         xtab[nxtab++] = x;
386       }
387     free (by_iter);
388   }
389   success = 1;
390
391  done:
392   /* All return paths lead here. */
393   {
394     int i;
395
396     for (i = 0; i < n_by; i++)
397       free (by[i]);
398     free (by);
399     free (by_nvar);
400   }
401
402   var_set_destroy (var_set);
403
404   return success;
405 }
406
407 /* Parses the VARIABLES subcommand. */
408 static int
409 crs_custom_variables (struct dataset *ds, struct cmd_crosstabs *cmd UNUSED, void *aux UNUSED)
410 {
411   if (nxtab)
412     {
413       msg (SE, _("VARIABLES must be specified before TABLES."));
414       return 0;
415     }
416
417   lex_match ('=');
418   
419   for (;;)
420     {
421       size_t orig_nv = variables_cnt;
422       size_t i;
423
424       long min, max;
425       
426       if (!parse_variables (dataset_dict (ds), &variables, &variables_cnt,
427                             (PV_APPEND | PV_NUMERIC
428                              | PV_NO_DUPLICATE | PV_NO_SCRATCH)))
429         return 0;
430
431       if (token != '(')
432         {
433           lex_error ("expecting `('");
434           goto lossage;
435         }
436       lex_get ();
437
438       if (!lex_force_int ())
439         goto lossage;
440       min = lex_integer ();
441       lex_get ();
442
443       lex_match (',');
444
445       if (!lex_force_int ())
446         goto lossage;
447       max = lex_integer ();
448       if (max < min)
449         {
450           msg (SE, _("Maximum value (%ld) less than minimum value (%ld)."),
451                max, min);
452           goto lossage;
453         }
454       lex_get ();
455
456       if (token != ')')
457         {
458           lex_error ("expecting `)'");
459           goto lossage;
460         }
461       lex_get ();
462       
463       for (i = orig_nv; i < variables_cnt; i++) 
464         {
465           struct var_range *vr = xmalloc (sizeof *vr);
466           vr->min = min;
467           vr->max = max + 1.;
468           vr->count = max - min + 1;
469           var_attach_aux (variables[i], vr, var_dtor_free);
470         }
471       
472       if (token == '/')
473         break;
474     }
475   
476   return 1;
477
478  lossage:
479   free (variables);
480   variables = NULL;
481   return 0;
482 }
483 \f
484 /* Data file processing. */
485
486 static int compare_table_entry (const void *, const void *, const void *);
487 static unsigned hash_table_entry (const void *, const void *);
488
489 /* Set up the crosstabulation tables for processing. */
490 static  void
491 precalc (const struct ccase *first, void *aux UNUSED, const struct dataset *ds)
492 {
493   output_split_file_values (ds, first);
494   if (mode == GENERAL)
495     {
496       gen_tab = hsh_create (512, compare_table_entry, hash_table_entry,
497                             NULL, NULL);
498     }
499   else 
500     {
501       int i;
502
503       sorted_tab = NULL;
504       n_sorted_tab = 0;
505
506       for (i = 0; i < nxtab; i++)
507         {
508           struct crosstab *x = xtab[i];
509           int count = 1;
510           int *v;
511           int j;
512
513           x->ofs = n_sorted_tab;
514
515           for (j = 2; j < x->nvar; j++) 
516             count *= get_var_range (x->vars[j - 2])->count;
517           
518           sorted_tab = xnrealloc (sorted_tab,
519                                   n_sorted_tab + count, sizeof *sorted_tab);
520           v = local_alloc (sizeof *v * x->nvar);
521           for (j = 2; j < x->nvar; j++) 
522             v[j] = get_var_range (x->vars[j])->min; 
523           for (j = 0; j < count; j++)
524             {
525               struct table_entry *te;
526               int k;
527
528               te = sorted_tab[n_sorted_tab++]
529                 = xmalloc (sizeof *te + sizeof (union value) * (x->nvar - 1));
530               te->table = i;
531               
532               {
533                 int row_cnt = get_var_range (x->vars[0])->count;
534                 int col_cnt = get_var_range (x->vars[1])->count;
535                 const int mat_size = row_cnt * col_cnt;
536                 int m;
537                 
538                 te->u.data = xnmalloc (mat_size, sizeof *te->u.data);
539                 for (m = 0; m < mat_size; m++)
540                   te->u.data[m] = 0.;
541               }
542               
543               for (k = 2; k < x->nvar; k++)
544                 te->values[k].f = v[k];
545               for (k = 2; k < x->nvar; k++) 
546                 {
547                   struct var_range *vr = get_var_range (x->vars[k]);
548                   if (++v[k] >= vr->max)
549                     v[k] = vr->min;
550                   else
551                     break; 
552                 }
553             }
554           local_free (v);
555         }
556
557       sorted_tab = xnrealloc (sorted_tab,
558                               n_sorted_tab + 1, sizeof *sorted_tab);
559       sorted_tab[n_sorted_tab] = NULL;
560     }
561
562 }
563
564 /* Form crosstabulations for general mode. */
565 static bool
566 calc_general (const struct ccase *c, void *aux UNUSED, const struct dataset *ds)
567 {
568   bool bad_warn = true;
569
570   /* Case weight. */
571   double weight = dict_get_case_weight (dataset_dict (ds), c, &bad_warn);
572
573   /* Flattened current table index. */
574   int t;
575
576   for (t = 0; t < nxtab; t++)
577     {
578       struct crosstab *x = xtab[t];
579       const size_t entry_size = (sizeof (struct table_entry)
580                                  + sizeof (union value) * (x->nvar - 1));
581       struct table_entry *te = local_alloc (entry_size);
582
583       /* Construct table entry for the current record and table. */
584       te->table = t;
585       {
586         int j;
587
588         assert (x != NULL);
589         for (j = 0; j < x->nvar; j++)
590           {
591             const union value *v = case_data (c, x->vars[j]->fv);
592             const struct missing_values *mv = &x->vars[j]->miss;
593             if ((cmd.miss == CRS_TABLE && mv_is_value_missing (mv, v))
594                 || (cmd.miss == CRS_INCLUDE
595                     && mv_is_value_system_missing (mv, v)))
596               {
597                 x->missing += weight;
598                 goto next_crosstab;
599               }
600               
601             if (x->vars[j]->type == NUMERIC)
602               te->values[j].f = case_num (c, x->vars[j]->fv);
603             else
604               {
605                 memcpy (te->values[j].s, case_str (c, x->vars[j]->fv),
606                         x->vars[j]->width);
607               
608                 /* Necessary in order to simplify comparisons. */
609                 memset (&te->values[j].s[x->vars[j]->width], 0,
610                         sizeof (union value) - x->vars[j]->width);
611               }
612           }
613       }
614
615       /* Add record to hash table. */
616       {
617         struct table_entry **tepp
618           = (struct table_entry **) hsh_probe (gen_tab, te);
619         if (*tepp == NULL)
620           {
621             struct table_entry *tep = pool_alloc (pl_tc, entry_size);
622             
623             te->u.freq = weight;
624             memcpy (tep, te, entry_size);
625             
626             *tepp = tep;
627           }
628         else
629           (*tepp)->u.freq += weight;
630       }
631
632     next_crosstab:
633       local_free (te);
634     }
635   
636   return true;
637 }
638
639 static bool
640 calc_integer (const struct ccase *c, void *aux UNUSED, const struct dataset *ds)
641 {
642   bool bad_warn = true;
643
644   /* Case weight. */
645   double weight = dict_get_case_weight (dataset_dict (ds), c, &bad_warn);
646   
647   /* Flattened current table index. */
648   int t;
649   
650   for (t = 0; t < nxtab; t++)
651     {
652       struct crosstab *x = xtab[t];
653       int i, fact, ofs;
654       
655       fact = i = 1;
656       ofs = x->ofs;
657       for (i = 0; i < x->nvar; i++)
658         {
659           struct variable *const v = x->vars[i];
660           struct var_range *vr = get_var_range (v);
661           double value = case_num (c, v->fv);
662           
663           /* Note that the first test also rules out SYSMIS. */
664           if ((value < vr->min || value >= vr->max)
665               || (cmd.miss == CRS_TABLE
666                   && mv_is_num_user_missing (&v->miss, value)))
667             {
668               x->missing += weight;
669               goto next_crosstab;
670             }
671           
672           if (i > 1)
673             {
674               ofs += fact * ((int) value - vr->min);
675               fact *= vr->count;
676             }
677         }
678       
679       {
680         struct variable *row_var = x->vars[ROW_VAR];
681         const int row = case_num (c, row_var->fv) - get_var_range (row_var)->min;
682
683         struct variable *col_var = x->vars[COL_VAR];
684         const int col = case_num (c, col_var->fv) - get_var_range (col_var)->min;
685
686         const int col_dim = get_var_range (col_var)->count;
687
688         sorted_tab[ofs]->u.data[col + row * col_dim] += weight;
689       }
690       
691     next_crosstab: ;
692     }
693   
694   return true;
695 }
696
697 /* Compare the table_entry's at A and B and return a strcmp()-type
698    result. */
699 static int 
700 compare_table_entry (const void *a_, const void *b_, const void *aux UNUSED)
701 {
702   const struct table_entry *a = a_;
703   const struct table_entry *b = b_;
704
705   if (a->table > b->table)
706     return 1;
707   else if (a->table < b->table)
708     return -1;
709   
710   {
711     const struct crosstab *x = xtab[a->table];
712     int i;
713
714     for (i = x->nvar - 1; i >= 0; i--)
715       if (x->vars[i]->type == NUMERIC)
716         {
717           const double diffnum = a->values[i].f - b->values[i].f;
718           if (diffnum < 0)
719             return -1;
720           else if (diffnum > 0)
721             return 1;
722         }
723       else 
724         {
725           assert (x->vars[i]->type == ALPHA);
726           {
727             const int diffstr = strncmp (a->values[i].s, b->values[i].s,
728                                          x->vars[i]->width);
729             if (diffstr)
730               return diffstr;
731           }
732         }
733   }
734   
735   return 0;
736 }
737
738 /* Calculate a hash value from table_entry A. */
739 static unsigned
740 hash_table_entry (const void *a_, const void *aux UNUSED)
741 {
742   const struct table_entry *a = a_;
743   unsigned long hash;
744   int i;
745
746   hash = a->table;
747   for (i = 0; i < xtab[a->table]->nvar; i++)
748     hash ^= hsh_hash_bytes (&a->values[i], sizeof a->values[i]);
749   
750   return hash;
751 }
752 \f
753 /* Post-data reading calculations. */
754
755 static struct table_entry **find_pivot_extent (struct table_entry **,
756                                                int *cnt, int pivot);
757 static void enum_var_values (struct table_entry **entries, int entry_cnt,
758                              int var_idx,
759                              union value **values, int *value_cnt);
760 static void output_pivot_table (struct table_entry **, struct table_entry **,
761                                 double **, double **, double **,
762                                 int *, int *, int *);
763 static void make_summary_table (void);
764
765 static bool
766 postcalc (void *aux UNUSED, const struct dataset *ds UNUSED)
767 {
768   if (mode == GENERAL)
769     {
770       n_sorted_tab = hsh_count (gen_tab);
771       sorted_tab = (struct table_entry **) hsh_sort (gen_tab);
772     }
773   
774   make_summary_table ();
775   
776   /* Identify all the individual crosstabulation tables, and deal with
777      them. */
778   {
779     struct table_entry **pb = sorted_tab, **pe; /* Pivot begin, pivot end. */
780     int pc = n_sorted_tab;                      /* Pivot count. */
781
782     double *mat = NULL, *row_tot = NULL, *col_tot = NULL;
783     int maxrows = 0, maxcols = 0, maxcells = 0;
784
785     for (;;)
786       {
787         pe = find_pivot_extent (pb, &pc, cmd.pivot == CRS_PIVOT);
788         if (pe == NULL)
789           break;
790         
791         output_pivot_table (pb, pe, &mat, &row_tot, &col_tot,
792                             &maxrows, &maxcols, &maxcells);
793           
794         pb = pe;
795       }
796     free (mat);
797     free (row_tot);
798     free (col_tot);
799   }
800   
801   hsh_destroy (gen_tab);
802
803   return true;
804 }
805
806 static void insert_summary (struct tab_table *, int tab_index, double valid);
807
808 /* Output a table summarizing the cases processed. */
809 static void
810 make_summary_table (void)
811 {
812   struct tab_table *summary;
813   
814   struct table_entry **pb = sorted_tab, **pe;
815   int pc = n_sorted_tab;
816   int cur_tab = 0;
817
818   summary = tab_create (7, 3 + nxtab, 1);
819   tab_title (summary, _("Summary."));
820   tab_headers (summary, 1, 0, 3, 0);
821   tab_joint_text (summary, 1, 0, 6, 0, TAB_CENTER, _("Cases"));
822   tab_joint_text (summary, 1, 1, 2, 1, TAB_CENTER, _("Valid"));
823   tab_joint_text (summary, 3, 1, 4, 1, TAB_CENTER, _("Missing"));
824   tab_joint_text (summary, 5, 1, 6, 1, TAB_CENTER, _("Total"));
825   tab_hline (summary, TAL_1, 1, 6, 1);
826   tab_hline (summary, TAL_1, 1, 6, 2);
827   tab_vline (summary, TAL_1, 3, 1, 1);
828   tab_vline (summary, TAL_1, 5, 1, 1);
829   {
830     int i;
831
832     for (i = 0; i < 3; i++)
833       {
834         tab_text (summary, 1 + i * 2, 2, TAB_RIGHT, _("N"));
835         tab_text (summary, 2 + i * 2, 2, TAB_RIGHT, _("Percent"));
836       }
837   }
838   tab_offset (summary, 0, 3);
839                   
840   for (;;)
841     {
842       double valid;
843       
844       pe = find_pivot_extent (pb, &pc, cmd.pivot == CRS_PIVOT);
845       if (pe == NULL)
846         break;
847
848       while (cur_tab < (*pb)->table)
849         insert_summary (summary, cur_tab++, 0.);
850
851       if (mode == GENERAL)
852         for (valid = 0.; pb < pe; pb++)
853           valid += (*pb)->u.freq;
854       else
855         {
856           const struct crosstab *const x = xtab[(*pb)->table];
857           const int n_cols = get_var_range (x->vars[COL_VAR])->count;
858           const int n_rows = get_var_range (x->vars[ROW_VAR])->count;
859           const int count = n_cols * n_rows;
860             
861           for (valid = 0.; pb < pe; pb++)
862             {
863               const double *data = (*pb)->u.data;
864               int i;
865                 
866               for (i = 0; i < count; i++)
867                 valid += *data++;
868             }
869         }
870       insert_summary (summary, cur_tab++, valid);
871
872       pb = pe;
873     }
874   
875   while (cur_tab < nxtab)
876     insert_summary (summary, cur_tab++, 0.);
877
878   submit (summary);
879 }
880
881 /* Inserts a line into T describing the crosstabulation at index
882    TAB_INDEX, which has VALID valid observations. */
883 static void
884 insert_summary (struct tab_table *t, int tab_index, double valid)
885 {
886   struct crosstab *x = xtab[tab_index];
887
888   tab_hline (t, TAL_1, 0, 6, 0);
889   
890   /* Crosstabulation name. */
891   {
892     char *buf = local_alloc (128 * x->nvar);
893     char *cp = buf;
894     int i;
895
896     for (i = 0; i < x->nvar; i++)
897       {
898         if (i > 0)
899           cp = stpcpy (cp, " * ");
900
901         cp = stpcpy (cp,
902                      x->vars[i]->label ? x->vars[i]->label : x->vars[i]->name);
903       }
904     tab_text (t, 0, 0, TAB_LEFT, buf);
905
906     local_free (buf);
907   }
908     
909   /* Counts and percentages. */
910   {
911     double n[3];
912     int i;
913
914     n[0] = valid;
915     n[1] = x->missing;
916     n[2] = n[0] + n[1];
917
918
919     for (i = 0; i < 3; i++)
920       {
921         tab_float (t, i * 2 + 1, 0, TAB_RIGHT, n[i], 8, 0);
922         tab_text (t, i * 2 + 2, 0, TAB_RIGHT | TAT_PRINTF, "%.1f%%",
923                   n[i] / n[2] * 100.);
924       }
925   }
926   
927   tab_next_row (t);
928 }
929 \f
930 /* Output. */
931
932 /* Tables. */
933 static struct tab_table *table; /* Crosstabulation table. */
934 static struct tab_table *chisq; /* Chi-square table. */
935 static struct tab_table *sym;           /* Symmetric measures table. */
936 static struct tab_table *risk;          /* Risk estimate table. */
937 static struct tab_table *direct;        /* Directional measures table. */
938
939 /* Statistics. */
940 static int chisq_fisher;        /* Did any rows include Fisher's exact test? */
941
942 /* Column values, number of columns. */
943 static union value *cols;
944 static int n_cols;
945
946 /* Row values, number of rows. */
947 static union value *rows;
948 static int n_rows;
949               
950 /* Number of statistically interesting columns/rows (columns/rows with
951    data in them). */
952 static int ns_cols, ns_rows;
953
954 /* Crosstabulation. */
955 static const struct crosstab *x;
956
957 /* Number of variables from the crosstabulation to consider.  This is
958    either x->nvar, if pivoting is on, or 2, if pivoting is off. */
959 static int nvar;
960
961 /* Matrix contents. */
962 static double *mat;             /* Matrix proper. */
963 static double *row_tot;         /* Row totals. */
964 static double *col_tot;         /* Column totals. */
965 static double W;                /* Grand total. */
966
967 static void display_dimensions (struct tab_table *, int first_difference,
968                                 struct table_entry *);
969 static void display_crosstabulation (void);
970 static void display_chisq (void);
971 static void display_symmetric (void);
972 static void display_risk (void);
973 static void display_directional (void);
974 static void crosstabs_dim (struct tab_table *, struct outp_driver *);
975 static void table_value_missing (struct tab_table *table, int c, int r,
976                                  unsigned char opt, const union value *v,
977                                  const struct variable *var);
978 static void delete_missing (void);
979
980 /* Output pivot table beginning at PB and continuing until PE,
981    exclusive.  For efficiency, *MATP is a pointer to a matrix that can
982    hold *MAXROWS entries. */
983 static void
984 output_pivot_table (struct table_entry **pb, struct table_entry **pe,
985                     double **matp, double **row_totp, double **col_totp,
986                     int *maxrows, int *maxcols, int *maxcells)
987 {
988   /* Subtable. */
989   struct table_entry **tb = pb, **te;   /* Table begin, table end. */
990   int tc = pe - pb;             /* Table count. */
991
992   /* Table entry for header comparison. */
993   struct table_entry *cmp = NULL;
994
995   x = xtab[(*pb)->table];
996   enum_var_values (pb, pe - pb, COL_VAR, &cols, &n_cols);
997
998   nvar = cmd.pivot == CRS_PIVOT ? x->nvar : 2;
999
1000   /* Crosstabulation table initialization. */
1001   if (num_cells)
1002     {
1003       table = tab_create (nvar + n_cols,
1004                           (pe - pb) / n_cols * 3 / 2 * num_cells + 10, 1);
1005       tab_headers (table, nvar - 1, 0, 2, 0);
1006
1007       /* First header line. */
1008       tab_joint_text (table, nvar - 1, 0, (nvar - 1) + (n_cols - 1), 0,
1009                       TAB_CENTER | TAT_TITLE, x->vars[COL_VAR]->name);
1010   
1011       tab_hline (table, TAL_1, nvar - 1, nvar + n_cols - 2, 1);
1012              
1013       /* Second header line. */
1014       {
1015         int i;
1016
1017         for (i = 2; i < nvar; i++)
1018           tab_joint_text (table, nvar - i - 1, 0, nvar - i - 1, 1,
1019                           TAB_RIGHT | TAT_TITLE,
1020                           (x->vars[i]->label
1021                            ? x->vars[i]->label : x->vars[i]->name));
1022         tab_text (table, nvar - 2, 1, TAB_RIGHT | TAT_TITLE,
1023                   x->vars[ROW_VAR]->name);
1024         for (i = 0; i < n_cols; i++)
1025           table_value_missing (table, nvar + i - 1, 1, TAB_RIGHT, &cols[i],
1026                                x->vars[COL_VAR]);
1027         tab_text (table, nvar + n_cols - 1, 1, TAB_CENTER, _("Total"));
1028       }
1029
1030       tab_hline (table, TAL_1, 0, nvar + n_cols - 1, 2);
1031       tab_vline (table, TAL_1, nvar + n_cols - 1, 0, 1);
1032
1033       /* Title. */
1034       {
1035         char *title = local_alloc (x->nvar * 64 + 128);
1036         char *cp = title;
1037         int i;
1038     
1039         if (cmd.pivot == CRS_PIVOT)
1040           for (i = 0; i < nvar; i++)
1041             {
1042               if (i)
1043                 cp = stpcpy (cp, " by ");
1044               cp = stpcpy (cp, x->vars[i]->name);
1045             }
1046         else
1047           {
1048             cp = spprintf (cp, "%s by %s for",
1049                            x->vars[0]->name, x->vars[1]->name);
1050             for (i = 2; i < nvar; i++)
1051               {
1052                 char buf[64], *bufp;
1053
1054                 if (i > 2)
1055                   *cp++ = ',';
1056                 *cp++ = ' ';
1057                 cp = stpcpy (cp, x->vars[i]->name);
1058                 *cp++ = '=';
1059                 format_short (buf, &x->vars[i]->print, &(*pb)->values[i]);
1060                 for (bufp = buf; isspace ((unsigned char) *bufp); bufp++)
1061                   ;
1062                 cp = stpcpy (cp, bufp);
1063               }
1064           }
1065
1066         cp = stpcpy (cp, " [");
1067         for (i = 0; i < num_cells; i++)
1068           {
1069             struct tuple
1070               {
1071                 int value;
1072                 const char *name;
1073               };
1074         
1075             static const struct tuple cell_names[] = 
1076               {
1077                 {CRS_CL_COUNT, N_("count")},
1078                 {CRS_CL_ROW, N_("row %")},
1079                 {CRS_CL_COLUMN, N_("column %")},
1080                 {CRS_CL_TOTAL, N_("total %")},
1081                 {CRS_CL_EXPECTED, N_("expected")},
1082                 {CRS_CL_RESIDUAL, N_("residual")},
1083                 {CRS_CL_SRESIDUAL, N_("std. resid.")},
1084                 {CRS_CL_ASRESIDUAL, N_("adj. resid.")},
1085                 {-1, NULL},
1086               };
1087
1088             const struct tuple *t;
1089
1090             for (t = cell_names; t->value != cells[i]; t++)
1091               assert (t->value != -1);
1092             if (i)
1093               cp = stpcpy (cp, ", ");
1094             cp = stpcpy (cp, gettext (t->name));
1095           }
1096         strcpy (cp, "].");
1097
1098         tab_title (table, "%s", title);
1099         local_free (title);
1100       }
1101       
1102       tab_offset (table, 0, 2);
1103     }
1104   else
1105     table = NULL;
1106   
1107   /* Chi-square table initialization. */
1108   if (cmd.a_statistics[CRS_ST_CHISQ])
1109     {
1110       chisq = tab_create (6 + (nvar - 2),
1111                           (pe - pb) / n_cols * 3 / 2 * N_CHISQ + 10, 1);
1112       tab_headers (chisq, 1 + (nvar - 2), 0, 1, 0);
1113
1114       tab_title (chisq, _("Chi-square tests."));
1115       
1116       tab_offset (chisq, nvar - 2, 0);
1117       tab_text (chisq, 0, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1118       tab_text (chisq, 1, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1119       tab_text (chisq, 2, 0, TAB_RIGHT | TAT_TITLE, _("df"));
1120       tab_text (chisq, 3, 0, TAB_RIGHT | TAT_TITLE,
1121                 _("Asymp. Sig. (2-sided)"));
1122       tab_text (chisq, 4, 0, TAB_RIGHT | TAT_TITLE,
1123                 _("Exact. Sig. (2-sided)"));
1124       tab_text (chisq, 5, 0, TAB_RIGHT | TAT_TITLE,
1125                 _("Exact. Sig. (1-sided)"));
1126       chisq_fisher = 0;
1127       tab_offset (chisq, 0, 1);
1128     }
1129   else
1130     chisq = NULL;
1131   
1132   /* Symmetric measures. */
1133   if (cmd.a_statistics[CRS_ST_PHI] || cmd.a_statistics[CRS_ST_CC]
1134       || cmd.a_statistics[CRS_ST_BTAU] || cmd.a_statistics[CRS_ST_CTAU]
1135       || cmd.a_statistics[CRS_ST_GAMMA] || cmd.a_statistics[CRS_ST_CORR]
1136       || cmd.a_statistics[CRS_ST_KAPPA])
1137     {
1138       sym = tab_create (6 + (nvar - 2), (pe - pb) / n_cols * 7 + 10, 1);
1139       tab_headers (sym, 2 + (nvar - 2), 0, 1, 0);
1140       tab_title (sym, _("Symmetric measures."));
1141
1142       tab_offset (sym, nvar - 2, 0);
1143       tab_text (sym, 0, 0, TAB_LEFT | TAT_TITLE, _("Category"));
1144       tab_text (sym, 1, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1145       tab_text (sym, 2, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1146       tab_text (sym, 3, 0, TAB_RIGHT | TAT_TITLE, _("Asymp. Std. Error"));
1147       tab_text (sym, 4, 0, TAB_RIGHT | TAT_TITLE, _("Approx. T"));
1148       tab_text (sym, 5, 0, TAB_RIGHT | TAT_TITLE, _("Approx. Sig."));
1149       tab_offset (sym, 0, 1);
1150     }
1151   else
1152     sym = NULL;
1153
1154   /* Risk estimate. */
1155   if (cmd.a_statistics[CRS_ST_RISK])
1156     {
1157       risk = tab_create (4 + (nvar - 2), (pe - pb) / n_cols * 4 + 10, 1);
1158       tab_headers (risk, 1 + nvar - 2, 0, 2, 0);
1159       tab_title (risk, _("Risk estimate."));
1160
1161       tab_offset (risk, nvar - 2, 0);
1162       tab_joint_text (risk, 2, 0, 3, 0, TAB_CENTER | TAT_TITLE | TAT_PRINTF,
1163                       _("95%% Confidence Interval"));
1164       tab_text (risk, 0, 1, TAB_LEFT | TAT_TITLE, _("Statistic"));
1165       tab_text (risk, 1, 1, TAB_RIGHT | TAT_TITLE, _("Value"));
1166       tab_text (risk, 2, 1, TAB_RIGHT | TAT_TITLE, _("Lower"));
1167       tab_text (risk, 3, 1, TAB_RIGHT | TAT_TITLE, _("Upper"));
1168       tab_hline (risk, TAL_1, 2, 3, 1);
1169       tab_vline (risk, TAL_1, 2, 0, 1);
1170       tab_offset (risk, 0, 2);
1171     }
1172   else
1173     risk = NULL;
1174
1175   /* Directional measures. */
1176   if (cmd.a_statistics[CRS_ST_LAMBDA] || cmd.a_statistics[CRS_ST_UC]
1177       || cmd.a_statistics[CRS_ST_D] || cmd.a_statistics[CRS_ST_ETA])
1178     {
1179       direct = tab_create (7 + (nvar - 2), (pe - pb) / n_cols * 7 + 10, 1);
1180       tab_headers (direct, 3 + (nvar - 2), 0, 1, 0);
1181       tab_title (direct, _("Directional measures."));
1182
1183       tab_offset (direct, nvar - 2, 0);
1184       tab_text (direct, 0, 0, TAB_LEFT | TAT_TITLE, _("Category"));
1185       tab_text (direct, 1, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1186       tab_text (direct, 2, 0, TAB_LEFT | TAT_TITLE, _("Type"));
1187       tab_text (direct, 3, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1188       tab_text (direct, 4, 0, TAB_RIGHT | TAT_TITLE, _("Asymp. Std. Error"));
1189       tab_text (direct, 5, 0, TAB_RIGHT | TAT_TITLE, _("Approx. T"));
1190       tab_text (direct, 6, 0, TAB_RIGHT | TAT_TITLE, _("Approx. Sig."));
1191       tab_offset (direct, 0, 1);
1192     }
1193   else
1194     direct = NULL;
1195
1196   for (;;)
1197     {
1198       /* Find pivot subtable if applicable. */
1199       te = find_pivot_extent (tb, &tc, 0);
1200       if (te == NULL)
1201         break;
1202
1203       /* Find all the row variable values. */
1204       enum_var_values (tb, te - tb, ROW_VAR, &rows, &n_rows);
1205
1206       /* Allocate memory space for the column and row totals. */
1207       if (n_rows > *maxrows)
1208         {
1209           *row_totp = xnrealloc (*row_totp, n_rows, sizeof **row_totp);
1210           row_tot = *row_totp;
1211           *maxrows = n_rows;
1212         }
1213       if (n_cols > *maxcols)
1214         {
1215           *col_totp = xnrealloc (*col_totp, n_cols, sizeof **col_totp);
1216           col_tot = *col_totp;
1217           *maxcols = n_cols;
1218         }
1219       
1220       /* Allocate table space for the matrix. */
1221       if (table && tab_row (table) + (n_rows + 1) * num_cells > tab_nr (table))
1222         tab_realloc (table, -1,
1223                      max (tab_nr (table) + (n_rows + 1) * num_cells,
1224                           tab_nr (table) * (pe - pb) / (te - tb)));
1225
1226       if (mode == GENERAL)
1227         {
1228           /* Allocate memory space for the matrix. */
1229           if (n_cols * n_rows > *maxcells)
1230             {
1231               *matp = xnrealloc (*matp, n_cols * n_rows, sizeof **matp);
1232               *maxcells = n_cols * n_rows;
1233             }
1234           
1235           mat = *matp;
1236
1237           /* Build the matrix and calculate column totals. */
1238           {
1239             union value *cur_col = cols;
1240             union value *cur_row = rows;
1241             double *mp = mat;
1242             double *cp = col_tot;
1243             struct table_entry **p;
1244
1245             *cp = 0.;
1246             for (p = &tb[0]; p < te; p++)
1247               {
1248                 for (; memcmp (cur_col, &(*p)->values[COL_VAR], sizeof *cur_col);
1249                      cur_row = rows)
1250                   {
1251                     *++cp = 0.;
1252                     for (; cur_row < &rows[n_rows]; cur_row++)
1253                       {
1254                         *mp = 0.;
1255                         mp += n_cols;
1256                       }
1257                     cur_col++;
1258                     mp = &mat[cur_col - cols];
1259                   }
1260
1261                 for (; memcmp (cur_row, &(*p)->values[ROW_VAR], sizeof *cur_row);
1262                      cur_row++)
1263                   {
1264                     *mp = 0.;
1265                     mp += n_cols;
1266                   }
1267
1268                 *cp += *mp = (*p)->u.freq;
1269                 mp += n_cols;
1270                 cur_row++;
1271               }
1272
1273             /* Zero out the rest of the matrix. */
1274             for (; cur_row < &rows[n_rows]; cur_row++)
1275               {
1276                 *mp = 0.;
1277                 mp += n_cols;
1278               }
1279             cur_col++;
1280             if (cur_col < &cols[n_cols])
1281               {
1282                 const int rem_cols = n_cols - (cur_col - cols);
1283                 int c, r;
1284
1285                 for (c = 0; c < rem_cols; c++)
1286                   *++cp = 0.;
1287                 mp = &mat[cur_col - cols];
1288                 for (r = 0; r < n_rows; r++)
1289                   {
1290                     for (c = 0; c < rem_cols; c++)
1291                       *mp++ = 0.;
1292                     mp += n_cols - rem_cols;
1293                   }
1294               }
1295           }
1296         }
1297       else
1298         {
1299           int r, c;
1300           double *tp = col_tot;
1301           
1302           assert (mode == INTEGER);
1303           mat = (*tb)->u.data;
1304           ns_cols = n_cols;
1305
1306           /* Calculate column totals. */
1307           for (c = 0; c < n_cols; c++)
1308             {
1309               double cum = 0.;
1310               double *cp = &mat[c];
1311               
1312               for (r = 0; r < n_rows; r++)
1313                 cum += cp[r * n_cols];
1314               *tp++ = cum;
1315             }
1316         }
1317       
1318       {
1319         double *cp;
1320         
1321         for (ns_cols = 0, cp = col_tot; cp < &col_tot[n_cols]; cp++)
1322           ns_cols += *cp != 0.;
1323       }
1324
1325       /* Calculate row totals. */
1326       {
1327         double *mp = mat;
1328         double *rp = row_tot;
1329         int r, c;
1330                 
1331         for (ns_rows = 0, r = n_rows; r--; )
1332           {
1333             double cum = 0.;
1334             for (c = n_cols; c--; )
1335               cum += *mp++;
1336             *rp++ = cum;
1337             if (cum != 0.)
1338               ns_rows++;
1339           }
1340       }
1341
1342       /* Calculate grand total. */
1343       {
1344         double *tp;
1345         double cum = 0.;
1346         int n;
1347
1348         if (n_rows < n_cols)
1349           tp = row_tot, n = n_rows;
1350         else
1351           tp = col_tot, n = n_cols;
1352         while (n--)
1353           cum += *tp++;
1354         W = cum;
1355       }
1356       
1357       /* Find the first variable that differs from the last subtable,
1358          then display the values of the dimensioning variables for
1359          each table that needs it. */
1360       {
1361         int first_difference = nvar - 1;
1362         
1363         if (tb != pb)
1364           for (; ; first_difference--)
1365             {
1366               assert (first_difference >= 2);
1367               if (memcmp (&cmp->values[first_difference],
1368                           &(*tb)->values[first_difference],
1369                           sizeof *cmp->values))
1370                 break;
1371             }
1372         cmp = *tb;
1373             
1374         if (table)
1375           display_dimensions (table, first_difference, *tb);
1376         if (chisq)
1377           display_dimensions (chisq, first_difference, *tb);
1378         if (sym)
1379           display_dimensions (sym, first_difference, *tb);
1380         if (risk)
1381           display_dimensions (risk, first_difference, *tb);
1382         if (direct)
1383           display_dimensions (direct, first_difference, *tb);
1384       }
1385
1386       if (table)
1387         display_crosstabulation ();
1388       if (cmd.miss == CRS_REPORT)
1389         delete_missing ();
1390       if (chisq)
1391         display_chisq ();
1392       if (sym)
1393         display_symmetric ();
1394       if (risk)
1395         display_risk ();
1396       if (direct)
1397         display_directional ();
1398                 
1399       tb = te;
1400       free (rows);
1401     }
1402
1403   submit (table);
1404   
1405   if (chisq)
1406     {
1407       if (!chisq_fisher)
1408         tab_resize (chisq, 4 + (nvar - 2), -1);
1409       submit (chisq);
1410     }
1411
1412   submit (sym);
1413   submit (risk);
1414   submit (direct);
1415
1416   free (cols);
1417 }
1418
1419 /* Delete missing rows and columns for statistical analysis when
1420    /MISSING=REPORT. */
1421 static void
1422 delete_missing (void)
1423 {
1424   {
1425     int r;
1426
1427     for (r = 0; r < n_rows; r++)
1428       if (mv_is_num_user_missing (&x->vars[ROW_VAR]->miss, rows[r].f))
1429         {
1430           int c;
1431
1432           for (c = 0; c < n_cols; c++)
1433             mat[c + r * n_cols] = 0.;
1434           ns_rows--;
1435         }
1436   }
1437   
1438   {
1439     int c;
1440
1441     for (c = 0; c < n_cols; c++)
1442       if (mv_is_num_user_missing (&x->vars[COL_VAR]->miss, cols[c].f))
1443         {
1444           int r;
1445
1446           for (r = 0; r < n_rows; r++)
1447             mat[c + r * n_cols] = 0.;
1448           ns_cols--;
1449         }
1450   }
1451 }
1452
1453 /* Prepare table T for submission, and submit it. */
1454 static void
1455 submit (struct tab_table *t)
1456 {
1457   int i;
1458   
1459   if (t == NULL)
1460     return;
1461   
1462   tab_resize (t, -1, 0);
1463   if (tab_nr (t) == tab_t (t))
1464     {
1465       tab_destroy (t);
1466       return;
1467     }
1468   tab_offset (t, 0, 0);
1469   if (t != table)
1470     for (i = 2; i < nvar; i++)
1471       tab_text (t, nvar - i - 1, 0, TAB_RIGHT | TAT_TITLE,
1472                 x->vars[i]->label ? x->vars[i]->label : x->vars[i]->name);
1473   tab_box (t, TAL_2, TAL_2, -1, -1, 0, 0, tab_nc (t) - 1, tab_nr (t) - 1);
1474   tab_box (t, -1, -1, -1, TAL_1, tab_l (t), tab_t (t) - 1, tab_nc (t) - 1,
1475            tab_nr (t) - 1);
1476   tab_box (t, -1, -1, -1, TAL_GAP, 0, tab_t (t), tab_l (t) - 1,
1477            tab_nr (t) - 1);
1478   tab_vline (t, TAL_2, tab_l (t), 0, tab_nr (t) - 1);
1479   tab_dim (t, crosstabs_dim);
1480   tab_submit (t);
1481 }
1482
1483 /* Sets the widths of all the columns and heights of all the rows in
1484    table T for driver D. */
1485 static void
1486 crosstabs_dim (struct tab_table *t, struct outp_driver *d)
1487 {
1488   int i;
1489   
1490   /* Width of a numerical column. */
1491   int c = outp_string_width (d, "0.000000", OUTP_PROPORTIONAL);
1492   if (cmd.miss == CRS_REPORT)
1493     c += outp_string_width (d, "M", OUTP_PROPORTIONAL);
1494
1495   /* Set width for header columns. */
1496   if (t->l != 0)
1497     {
1498       size_t i;
1499       int w;
1500
1501       w = d->width - c * (t->nc - t->l);
1502       for (i = 0; i <= t->nc; i++)
1503         w -= t->wrv[i];
1504       w /= t->l;
1505       
1506       if (w < d->prop_em_width * 8)
1507         w = d->prop_em_width * 8;
1508
1509       if (w > d->prop_em_width * 15)
1510         w = d->prop_em_width * 15;
1511
1512       for (i = 0; i < t->l; i++)
1513         t->w[i] = w;
1514     }
1515
1516   for (i = t->l; i < t->nc; i++)
1517     t->w[i] = c;
1518
1519   for (i = 0; i < t->nr; i++)
1520     t->h[i] = tab_natural_height (t, d, i);
1521 }
1522
1523 static struct table_entry **find_pivot_extent_general (struct table_entry **tp,
1524                                                 int *cnt, int pivot);
1525 static struct table_entry **find_pivot_extent_integer (struct table_entry **tp,
1526                                                 int *cnt, int pivot);
1527
1528 /* Calls find_pivot_extent_general or find_pivot_extent_integer, as
1529    appropriate. */
1530 static struct table_entry **
1531 find_pivot_extent (struct table_entry **tp, int *cnt, int pivot)
1532 {
1533   return (mode == GENERAL
1534           ? find_pivot_extent_general (tp, cnt, pivot)
1535           : find_pivot_extent_integer (tp, cnt, pivot));
1536 }
1537
1538 /* Find the extent of a region in TP that contains one table.  If
1539    PIVOT != 0 that means a set of table entries with identical table
1540    number; otherwise they also have to have the same values for every
1541    dimension after the row and column dimensions.  The table that is
1542    searched starts at TP and has length CNT.  Returns the first entry
1543    after the last one in the table; sets *CNT to the number of
1544    remaining values.  If there are no entries in TP at all, returns
1545    NULL.  A yucky interface, admittedly, but it works. */
1546 static struct table_entry **
1547 find_pivot_extent_general (struct table_entry **tp, int *cnt, int pivot)
1548 {
1549   struct table_entry *fp = *tp;
1550   struct crosstab *x;
1551
1552   if (*cnt == 0)
1553     return NULL;
1554   x = xtab[(*tp)->table];
1555   for (;;)
1556     {
1557       tp++;
1558       if (--*cnt == 0)
1559         break;
1560       assert (*cnt > 0);
1561
1562       if ((*tp)->table != fp->table)
1563         break;
1564       if (pivot)
1565         continue;
1566
1567       if (memcmp (&(*tp)->values[2], &fp->values[2], sizeof (union value) * (x->nvar - 2)))
1568         break;
1569     }
1570
1571   return tp;
1572 }
1573
1574 /* Integer mode correspondent to find_pivot_extent_general().  This
1575    could be optimized somewhat, but I just don't give a crap about
1576    CROSSTABS performance in integer mode, which is just a
1577    CROSSTABS wart as far as I'm concerned.
1578
1579    That said, feel free to send optimization patches to me. */
1580 static struct table_entry **
1581 find_pivot_extent_integer (struct table_entry **tp, int *cnt, int pivot)
1582 {
1583   struct table_entry *fp = *tp;
1584   struct crosstab *x;
1585
1586   if (*cnt == 0)
1587     return NULL;
1588   x = xtab[(*tp)->table];
1589   for (;;)
1590     {
1591       tp++;
1592       if (--*cnt == 0)
1593         break;
1594       assert (*cnt > 0);
1595
1596       if ((*tp)->table != fp->table)
1597         break;
1598       if (pivot)
1599         continue;
1600       
1601       if (memcmp (&(*tp)->values[2], &fp->values[2],
1602                   sizeof (union value) * (x->nvar - 2)))
1603         break;
1604     }
1605
1606   return tp;
1607 }
1608
1609 /* Compares `union value's A_ and B_ and returns a strcmp()-like
1610    result.  WIDTH_ points to an int which is either 0 for a
1611    numeric value or a string width for a string value. */
1612 static int
1613 compare_value (const void *a_, const void *b_, const void *width_)
1614 {
1615   const union value *a = a_;
1616   const union value *b = b_;
1617   const int *pwidth = width_;
1618   const int width = *pwidth;
1619
1620   if (width == 0)
1621     return (a->f < b->f) ? -1 : (a->f > b->f);
1622   else
1623     return strncmp (a->s, b->s, width);
1624 }
1625
1626 /* Given an array of ENTRY_CNT table_entry structures starting at
1627    ENTRIES, creates a sorted list of the values that the variable
1628    with index VAR_IDX takes on.  The values are returned as a
1629    malloc()'darray stored in *VALUES, with the number of values
1630    stored in *VALUE_CNT.
1631    */
1632 static void 
1633 enum_var_values (struct table_entry **entries, int entry_cnt, int var_idx,
1634                  union value **values, int *value_cnt)
1635 {
1636   struct variable *v = xtab[(*entries)->table]->vars[var_idx];
1637
1638   if (mode == GENERAL)
1639     {
1640       int width = v->width;
1641       int i;
1642
1643       *values = xnmalloc (entry_cnt, sizeof **values);
1644       for (i = 0; i < entry_cnt; i++)
1645         (*values)[i] = entries[i]->values[var_idx];
1646       *value_cnt = sort_unique (*values, entry_cnt, sizeof **values,
1647                                 compare_value, &width);
1648     }
1649   else
1650     {
1651       struct var_range *vr = get_var_range (v);
1652       int i;
1653       
1654       assert (mode == INTEGER);
1655       *values = xnmalloc (vr->count, sizeof **values);
1656       for (i = 0; i < vr->count; i++)
1657         (*values)[i].f = i + vr->min;
1658       *value_cnt = vr->count;
1659     }
1660 }
1661
1662 /* Sets cell (C,R) in TABLE, with options OPT, to have a value taken
1663    from V, displayed with print format spec from variable VAR.  When
1664    in REPORT missing-value mode, missing values have an M appended. */
1665 static void
1666 table_value_missing (struct tab_table *table, int c, int r, unsigned char opt,
1667                      const union value *v, const struct variable *var)
1668 {
1669   struct substring s;
1670
1671   const char *label = val_labs_find (var->val_labs, *v);
1672   if (label) 
1673     {
1674       tab_text (table, c, r, TAB_LEFT, label);
1675       return;
1676     }
1677
1678   s.string = tab_alloc (table, var->print.w);
1679   format_short (s.string, &var->print, v);
1680   s.length = strlen (s.string);
1681   if (cmd.miss == CRS_REPORT && mv_is_num_user_missing (&var->miss, v->f))
1682     s.string[s.length++] = 'M';
1683   while (s.length && *s.string == ' ')
1684     {
1685       s.length--;
1686       s.string++;
1687     }
1688   tab_raw (table, c, r, opt, &s);
1689 }
1690
1691 /* Draws a line across TABLE at the current row to indicate the most
1692    major dimension variable with index FIRST_DIFFERENCE out of NVAR
1693    that changed, and puts the values that changed into the table.  TB
1694    and X must be the corresponding table_entry and crosstab,
1695    respectively. */
1696 static void
1697 display_dimensions (struct tab_table *table, int first_difference, struct table_entry *tb)
1698 {
1699   tab_hline (table, TAL_1, nvar - first_difference - 1, tab_nc (table) - 1, 0);
1700
1701   for (; first_difference >= 2; first_difference--)
1702     table_value_missing (table, nvar - first_difference - 1, 0,
1703                          TAB_RIGHT, &tb->values[first_difference],
1704                          x->vars[first_difference]);
1705 }
1706
1707 /* Put VALUE into cell (C,R) of TABLE, suffixed with character
1708    SUFFIX if nonzero.  If MARK_MISSING is true the entry is
1709    additionally suffixed with a letter `M'. */
1710 static void
1711 format_cell_entry (struct tab_table *table, int c, int r, double value,
1712                    char suffix, bool mark_missing)
1713 {
1714   const struct fmt_spec f = {FMT_F, 10, 1};
1715   union value v;
1716   struct substring s;
1717   
1718   s.length = 10;
1719   s.string = tab_alloc (table, 16);
1720   v.f = value;
1721   data_out (s.string, &f, &v);
1722   while (*s.string == ' ')
1723     {
1724       s.length--;
1725       s.string++;
1726     }
1727   if (suffix != 0)
1728     s.string[s.length++] = suffix;
1729   if (mark_missing)
1730     s.string[s.length++] = 'M';
1731
1732   tab_raw (table, c, r, TAB_RIGHT, &s);
1733 }
1734
1735 /* Displays the crosstabulation table. */
1736 static void
1737 display_crosstabulation (void)
1738 {
1739   {
1740     int r;
1741         
1742     for (r = 0; r < n_rows; r++)
1743       table_value_missing (table, nvar - 2, r * num_cells,
1744                            TAB_RIGHT, &rows[r], x->vars[ROW_VAR]);
1745   }
1746   tab_text (table, nvar - 2, n_rows * num_cells,
1747             TAB_LEFT, _("Total"));
1748       
1749   /* Put in the actual cells. */
1750   {
1751     double *mp = mat;
1752     int r, c, i;
1753
1754     tab_offset (table, nvar - 1, -1);
1755     for (r = 0; r < n_rows; r++)
1756       {
1757         if (num_cells > 1)
1758           tab_hline (table, TAL_1, -1, n_cols, 0);
1759         for (c = 0; c < n_cols; c++)
1760           {
1761             bool mark_missing = false;
1762             double expected_value = row_tot[r] * col_tot[c] / W;
1763             if (cmd.miss == CRS_REPORT
1764                 && (mv_is_num_user_missing (&x->vars[COL_VAR]->miss, cols[c].f)
1765                     || mv_is_num_user_missing (&x->vars[ROW_VAR]->miss,
1766                                                rows[r].f)))
1767               mark_missing = true;
1768             for (i = 0; i < num_cells; i++)
1769               {
1770                 double v;
1771                 int suffix = 0;
1772
1773                 switch (cells[i])
1774                   {
1775                   case CRS_CL_COUNT:
1776                     v = *mp;
1777                     break;
1778                   case CRS_CL_ROW:
1779                     v = *mp / row_tot[r] * 100.;
1780                     suffix = '%';
1781                     break;
1782                   case CRS_CL_COLUMN:
1783                     v = *mp / col_tot[c] * 100.;
1784                     suffix = '%';
1785                     break;
1786                   case CRS_CL_TOTAL:
1787                     v = *mp / W * 100.;
1788                     suffix = '%';
1789                     break;
1790                   case CRS_CL_EXPECTED:
1791                     v = expected_value;
1792                     break;
1793                   case CRS_CL_RESIDUAL:
1794                     v = *mp - expected_value;
1795                     break;
1796                   case CRS_CL_SRESIDUAL:
1797                     v = (*mp - expected_value) / sqrt (expected_value);
1798                     break;
1799                   case CRS_CL_ASRESIDUAL:
1800                     v = ((*mp - expected_value)
1801                          / sqrt (expected_value
1802                                  * (1. - row_tot[r] / W)
1803                                  * (1. - col_tot[c] / W)));
1804                     break;
1805                   default:
1806                     NOT_REACHED ();
1807                   }
1808
1809                 format_cell_entry (table, c, i, v, suffix, mark_missing);
1810               }
1811
1812             mp++;
1813           }
1814
1815         tab_offset (table, -1, tab_row (table) + num_cells);
1816       }
1817   }
1818
1819   /* Row totals. */
1820   {
1821     int r, i;
1822
1823     tab_offset (table, -1, tab_row (table) - num_cells * n_rows);
1824     for (r = 0; r < n_rows; r++) 
1825       {
1826         char suffix = 0;
1827         bool mark_missing = false;
1828
1829         if (cmd.miss == CRS_REPORT
1830             && mv_is_num_user_missing (&x->vars[ROW_VAR]->miss, rows[r].f))
1831           mark_missing = true;
1832
1833         for (i = 0; i < num_cells; i++)
1834           {
1835             double v;
1836
1837             switch (cells[i])
1838               {
1839               case CRS_CL_COUNT:
1840                 v = row_tot[r];
1841                 break;
1842               case CRS_CL_ROW:
1843                 v = 100.;
1844                 suffix = '%';
1845                 break;
1846               case CRS_CL_COLUMN:
1847                 v = row_tot[r] / W * 100.;
1848                 suffix = '%';
1849                 break;
1850               case CRS_CL_TOTAL:
1851                 v = row_tot[r] / W * 100.;
1852                 suffix = '%';
1853                 break;
1854               case CRS_CL_EXPECTED:
1855               case CRS_CL_RESIDUAL:
1856               case CRS_CL_SRESIDUAL:
1857               case CRS_CL_ASRESIDUAL:
1858                 v = 0.;
1859                 break;
1860               default:
1861                 NOT_REACHED ();
1862               }
1863
1864             format_cell_entry (table, n_cols, 0, v, suffix, mark_missing);
1865             tab_next_row (table);
1866           } 
1867       }
1868   }
1869
1870   /* Column totals, grand total. */
1871   {
1872     int c;
1873     int last_row = 0;
1874
1875     if (num_cells > 1)
1876       tab_hline (table, TAL_1, -1, n_cols, 0);
1877     for (c = 0; c <= n_cols; c++)
1878       {
1879         double ct = c < n_cols ? col_tot[c] : W;
1880         bool mark_missing = false;
1881         char suffix = 0;
1882         int i;
1883             
1884         if (cmd.miss == CRS_REPORT && c < n_cols 
1885             && mv_is_num_user_missing (&x->vars[COL_VAR]->miss, cols[c].f))
1886           mark_missing = true;
1887
1888         for (i = 0; i < num_cells; i++)
1889           {
1890             double v;
1891
1892             switch (cells[i])
1893               {
1894               case CRS_CL_COUNT:
1895                 v = ct;
1896                 suffix = '%';
1897                 break;
1898               case CRS_CL_ROW:
1899                 v = ct / W * 100.;
1900                 suffix = '%';
1901                 break;
1902               case CRS_CL_COLUMN:
1903                 v = 100.;
1904                 suffix = '%';
1905                 break;
1906               case CRS_CL_TOTAL:
1907                 v = ct / W * 100.;
1908                 suffix = '%';
1909                 break;
1910               case CRS_CL_EXPECTED:
1911               case CRS_CL_RESIDUAL:
1912               case CRS_CL_SRESIDUAL:
1913               case CRS_CL_ASRESIDUAL:
1914                 continue;
1915               default:
1916                 NOT_REACHED ();
1917               }
1918
1919             format_cell_entry (table, c, i, v, suffix, mark_missing);
1920           }
1921         last_row = i;
1922       }
1923
1924     tab_offset (table, -1, tab_row (table) + last_row);
1925   }
1926   
1927   tab_offset (table, 0, -1);
1928 }
1929
1930 static void calc_r (double *X, double *Y, double *, double *, double *);
1931 static void calc_chisq (double[N_CHISQ], int[N_CHISQ], double *, double *);
1932
1933 /* Display chi-square statistics. */
1934 static void
1935 display_chisq (void)
1936 {
1937   static const char *chisq_stats[N_CHISQ] = 
1938     {
1939       N_("Pearson Chi-Square"),
1940       N_("Likelihood Ratio"),
1941       N_("Fisher's Exact Test"),
1942       N_("Continuity Correction"),
1943       N_("Linear-by-Linear Association"),
1944     };
1945   double chisq_v[N_CHISQ];
1946   double fisher1, fisher2;
1947   int df[N_CHISQ];
1948   int s = 0;
1949
1950   int i;
1951               
1952   calc_chisq (chisq_v, df, &fisher1, &fisher2);
1953
1954   tab_offset (chisq, nvar - 2, -1);
1955   
1956   for (i = 0; i < N_CHISQ; i++)
1957     {
1958       if ((i != 2 && chisq_v[i] == SYSMIS)
1959           || (i == 2 && fisher1 == SYSMIS))
1960         continue;
1961       s = 1;
1962       
1963       tab_text (chisq, 0, 0, TAB_LEFT, gettext (chisq_stats[i]));
1964       if (i != 2)
1965         {
1966           tab_float (chisq, 1, 0, TAB_RIGHT, chisq_v[i], 8, 3);
1967           tab_float (chisq, 2, 0, TAB_RIGHT, df[i], 8, 0);
1968           tab_float (chisq, 3, 0, TAB_RIGHT,
1969                      gsl_cdf_chisq_Q (chisq_v[i], df[i]), 8, 3);
1970         }
1971       else
1972         {
1973           chisq_fisher = 1;
1974           tab_float (chisq, 4, 0, TAB_RIGHT, fisher2, 8, 3);
1975           tab_float (chisq, 5, 0, TAB_RIGHT, fisher1, 8, 3);
1976         }
1977       tab_next_row (chisq);
1978     }
1979
1980   tab_text (chisq, 0, 0, TAB_LEFT, _("N of Valid Cases"));
1981   tab_float (chisq, 1, 0, TAB_RIGHT, W, 8, 0);
1982   tab_next_row (chisq);
1983     
1984   tab_offset (chisq, 0, -1);
1985 }
1986
1987 static int calc_symmetric (double[N_SYMMETRIC], double[N_SYMMETRIC],
1988                            double[N_SYMMETRIC]);
1989
1990 /* Display symmetric measures. */
1991 static void
1992 display_symmetric (void)
1993 {
1994   static const char *categories[] = 
1995     {
1996       N_("Nominal by Nominal"),
1997       N_("Ordinal by Ordinal"),
1998       N_("Interval by Interval"),
1999       N_("Measure of Agreement"),
2000     };
2001
2002   static const char *stats[N_SYMMETRIC] =
2003     {
2004       N_("Phi"),
2005       N_("Cramer's V"),
2006       N_("Contingency Coefficient"),
2007       N_("Kendall's tau-b"),
2008       N_("Kendall's tau-c"),
2009       N_("Gamma"),
2010       N_("Spearman Correlation"),
2011       N_("Pearson's R"),
2012       N_("Kappa"),
2013     };
2014
2015   static const int stats_categories[N_SYMMETRIC] =
2016     {
2017       0, 0, 0, 1, 1, 1, 1, 2, 3,
2018     };
2019
2020   int last_cat = -1;
2021   double sym_v[N_SYMMETRIC], sym_ase[N_SYMMETRIC], sym_t[N_SYMMETRIC];
2022   int i;
2023
2024   if (!calc_symmetric (sym_v, sym_ase, sym_t))
2025     return;
2026
2027   tab_offset (sym, nvar - 2, -1);
2028   
2029   for (i = 0; i < N_SYMMETRIC; i++)
2030     {
2031       if (sym_v[i] == SYSMIS)
2032         continue;
2033
2034       if (stats_categories[i] != last_cat)
2035         {
2036           last_cat = stats_categories[i];
2037           tab_text (sym, 0, 0, TAB_LEFT, gettext (categories[last_cat]));
2038         }
2039       
2040       tab_text (sym, 1, 0, TAB_LEFT, gettext (stats[i]));
2041       tab_float (sym, 2, 0, TAB_RIGHT, sym_v[i], 8, 3);
2042       if (sym_ase[i] != SYSMIS)
2043         tab_float (sym, 3, 0, TAB_RIGHT, sym_ase[i], 8, 3);
2044       if (sym_t[i] != SYSMIS)
2045         tab_float (sym, 4, 0, TAB_RIGHT, sym_t[i], 8, 3);
2046       /*tab_float (sym, 5, 0, TAB_RIGHT, normal_sig (sym_v[i]), 8, 3);*/
2047       tab_next_row (sym);
2048     }
2049
2050   tab_text (sym, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2051   tab_float (sym, 2, 0, TAB_RIGHT, W, 8, 0);
2052   tab_next_row (sym);
2053     
2054   tab_offset (sym, 0, -1);
2055 }
2056
2057 static int calc_risk (double[], double[], double[], union value *);
2058
2059 /* Display risk estimate. */
2060 static void
2061 display_risk (void)
2062 {
2063   char buf[256];
2064   double risk_v[3], lower[3], upper[3];
2065   union value c[2];
2066   int i;
2067   
2068   if (!calc_risk (risk_v, upper, lower, c))
2069     return;
2070   
2071   tab_offset (risk, nvar - 2, -1);
2072   
2073   for (i = 0; i < 3; i++)
2074     {
2075       if (risk_v[i] == SYSMIS)
2076         continue;
2077
2078       switch (i)
2079         {
2080         case 0:
2081           if (x->vars[COL_VAR]->type == NUMERIC)
2082             sprintf (buf, _("Odds Ratio for %s (%g / %g)"),
2083                      x->vars[COL_VAR]->name, c[0].f, c[1].f);
2084           else
2085             sprintf (buf, _("Odds Ratio for %s (%.*s / %.*s)"),
2086                      x->vars[COL_VAR]->name,
2087                      x->vars[COL_VAR]->width, c[0].s,
2088                      x->vars[COL_VAR]->width, c[1].s);
2089           break;
2090         case 1:
2091         case 2:
2092           if (x->vars[ROW_VAR]->type == NUMERIC)
2093             sprintf (buf, _("For cohort %s = %g"),
2094                      x->vars[ROW_VAR]->name, rows[i - 1].f);
2095           else
2096             sprintf (buf, _("For cohort %s = %.*s"),
2097                      x->vars[ROW_VAR]->name,
2098                      x->vars[ROW_VAR]->width, rows[i - 1].s);
2099           break;
2100         }
2101                    
2102       tab_text (risk, 0, 0, TAB_LEFT, buf);
2103       tab_float (risk, 1, 0, TAB_RIGHT, risk_v[i], 8, 3);
2104       tab_float (risk, 2, 0, TAB_RIGHT, lower[i], 8, 3);
2105       tab_float (risk, 3, 0, TAB_RIGHT, upper[i], 8, 3);
2106       tab_next_row (risk);
2107     }
2108
2109   tab_text (risk, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2110   tab_float (risk, 1, 0, TAB_RIGHT, W, 8, 0);
2111   tab_next_row (risk);
2112     
2113   tab_offset (risk, 0, -1);
2114 }
2115
2116 static int calc_directional (double[N_DIRECTIONAL], double[N_DIRECTIONAL],
2117                              double[N_DIRECTIONAL]);
2118
2119 /* Display directional measures. */
2120 static void
2121 display_directional (void)
2122 {
2123   static const char *categories[] = 
2124     {
2125       N_("Nominal by Nominal"),
2126       N_("Ordinal by Ordinal"),
2127       N_("Nominal by Interval"),
2128     };
2129
2130   static const char *stats[] =
2131     {
2132       N_("Lambda"),
2133       N_("Goodman and Kruskal tau"),
2134       N_("Uncertainty Coefficient"),
2135       N_("Somers' d"),
2136       N_("Eta"),
2137     };
2138
2139   static const char *types[] = 
2140     {
2141       N_("Symmetric"),
2142       N_("%s Dependent"),
2143       N_("%s Dependent"),
2144     };
2145
2146   static const int stats_categories[N_DIRECTIONAL] =
2147     {
2148       0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2,
2149     };
2150   
2151   static const int stats_stats[N_DIRECTIONAL] = 
2152     {
2153       0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4,
2154     };
2155
2156   static const int stats_types[N_DIRECTIONAL] = 
2157     {
2158       0, 1, 2, 1, 2, 0, 1, 2, 0, 1, 2, 1, 2,
2159     };
2160
2161   static const int *stats_lookup[] = 
2162     {
2163       stats_categories,
2164       stats_stats,
2165       stats_types,
2166     };
2167
2168   static const char **stats_names[] =
2169     {
2170       categories,
2171       stats,
2172       types,
2173     };
2174
2175   int last[3] =
2176     {
2177       -1, -1, -1,
2178     };
2179     
2180   double direct_v[N_DIRECTIONAL];
2181   double direct_ase[N_DIRECTIONAL];
2182   double direct_t[N_DIRECTIONAL];
2183   
2184   int i;
2185
2186   if (!calc_directional (direct_v, direct_ase, direct_t))
2187     return;
2188
2189   tab_offset (direct, nvar - 2, -1);
2190   
2191   for (i = 0; i < N_DIRECTIONAL; i++)
2192     {
2193       if (direct_v[i] == SYSMIS)
2194         continue;
2195       
2196       {
2197         int j;
2198
2199         for (j = 0; j < 3; j++)
2200           if (last[j] != stats_lookup[j][i])
2201             {
2202               if (j < 2)
2203                 tab_hline (direct, TAL_1, j, 6, 0);
2204               
2205               for (; j < 3; j++)
2206                 {
2207                   char *string;
2208                   int k = last[j] = stats_lookup[j][i];
2209
2210                   if (k == 0)
2211                     string = NULL;
2212                   else if (k == 1)
2213                     string = x->vars[0]->name;
2214                   else
2215                     string = x->vars[1]->name;
2216                   
2217                   tab_text (direct, j, 0, TAB_LEFT | TAT_PRINTF,
2218                             gettext (stats_names[j][k]), string);
2219                 }
2220             }
2221       }
2222       
2223       tab_float (direct, 3, 0, TAB_RIGHT, direct_v[i], 8, 3);
2224       if (direct_ase[i] != SYSMIS)
2225         tab_float (direct, 4, 0, TAB_RIGHT, direct_ase[i], 8, 3);
2226       if (direct_t[i] != SYSMIS)
2227         tab_float (direct, 5, 0, TAB_RIGHT, direct_t[i], 8, 3);
2228       /*tab_float (direct, 6, 0, TAB_RIGHT, normal_sig (direct_v[i]), 8, 3);*/
2229       tab_next_row (direct);
2230     }
2231
2232   tab_offset (direct, 0, -1);
2233 }
2234 \f
2235 /* Statistical calculations. */
2236
2237 /* Returns the value of the gamma (factorial) function for an integer
2238    argument X. */
2239 static double
2240 gamma_int (double x)
2241 {
2242   double r = 1;
2243   int i;
2244   
2245   for (i = 2; i < x; i++)
2246     r *= i;
2247   return r;
2248 }
2249
2250 /* Calculate P_r as specified in _SPSS Statistical Algorithms_,
2251    Appendix 5. */
2252 static inline double
2253 Pr (int a, int b, int c, int d)
2254 {
2255   return (gamma_int (a + b + 1.) / gamma_int (a + 1.)
2256           * gamma_int (c + d + 1.) / gamma_int (b + 1.)
2257           * gamma_int (a + c + 1.) / gamma_int (c + 1.)
2258           * gamma_int (b + d + 1.) / gamma_int (d + 1.)
2259           / gamma_int (a + b + c + d + 1.));
2260 }
2261
2262 /* Swap the contents of A and B. */
2263 static inline void
2264 swap (int *a, int *b)
2265 {
2266   int t = *a;
2267   *a = *b;
2268   *b = t;
2269 }
2270
2271 /* Calculate significance for Fisher's exact test as specified in
2272    _SPSS Statistical Algorithms_, Appendix 5. */
2273 static void
2274 calc_fisher (int a, int b, int c, int d, double *fisher1, double *fisher2)
2275 {
2276   int x;
2277   
2278   if (min (c, d) < min (a, b))
2279     swap (&a, &c), swap (&b, &d);
2280   if (min (b, d) < min (a, c))
2281     swap (&a, &b), swap (&c, &d);
2282   if (b * c < a * d)
2283     {
2284       if (b < c)
2285         swap (&a, &b), swap (&c, &d);
2286       else
2287         swap (&a, &c), swap (&b, &d);
2288     }
2289
2290   *fisher1 = 0.;
2291   for (x = 0; x <= a; x++)
2292     *fisher1 += Pr (a - x, b + x, c + x, d - x);
2293
2294   *fisher2 = *fisher1;
2295   for (x = 1; x <= b; x++)
2296     *fisher2 += Pr (a + x, b - x, c - x, d + x);
2297 }
2298
2299 /* Calculates chi-squares into CHISQ.  MAT is a matrix with N_COLS
2300    columns with values COLS and N_ROWS rows with values ROWS.  Values
2301    in the matrix sum to W. */
2302 static void
2303 calc_chisq (double chisq[N_CHISQ], int df[N_CHISQ],
2304             double *fisher1, double *fisher2)
2305 {
2306   int r, c;
2307
2308   chisq[0] = chisq[1] = 0.;
2309   chisq[2] = chisq[3] = chisq[4] = SYSMIS;
2310   *fisher1 = *fisher2 = SYSMIS;
2311
2312   df[0] = df[1] = (ns_cols - 1) * (ns_rows - 1);
2313
2314   if (ns_rows <= 1 || ns_cols <= 1)
2315     {
2316       chisq[0] = chisq[1] = SYSMIS;
2317       return;
2318     }
2319
2320   for (r = 0; r < n_rows; r++)
2321     for (c = 0; c < n_cols; c++)
2322       {
2323         const double expected = row_tot[r] * col_tot[c] / W;
2324         const double freq = mat[n_cols * r + c];
2325         const double residual = freq - expected;
2326     
2327         chisq[0] += residual * residual / expected;
2328         if (freq)
2329           chisq[1] += freq * log (expected / freq);
2330       }
2331
2332   if (chisq[0] == 0.)
2333     chisq[0] = SYSMIS;
2334
2335   if (chisq[1] != 0.)
2336     chisq[1] *= -2.;
2337   else
2338     chisq[1] = SYSMIS;
2339
2340   /* Calculate Yates and Fisher exact test. */
2341   if (ns_cols == 2 && ns_rows == 2)
2342     {
2343       double f11, f12, f21, f22;
2344       
2345       {
2346         int nz_cols[2];
2347         int i, j;
2348
2349         for (i = j = 0; i < n_cols; i++)
2350           if (col_tot[i] != 0.)
2351             {
2352               nz_cols[j++] = i;
2353               if (j == 2)
2354                 break;
2355             }
2356
2357         assert (j == 2);
2358
2359         f11 = mat[nz_cols[0]];
2360         f12 = mat[nz_cols[1]];
2361         f21 = mat[nz_cols[0] + n_cols];
2362         f22 = mat[nz_cols[1] + n_cols];
2363       }
2364
2365       /* Yates. */
2366       {
2367         const double x = fabs (f11 * f22 - f12 * f21) - 0.5 * W;
2368
2369         if (x > 0.)
2370           chisq[3] = (W * x * x
2371                       / (f11 + f12) / (f21 + f22)
2372                       / (f11 + f21) / (f12 + f22));
2373         else
2374           chisq[3] = 0.;
2375
2376         df[3] = 1.;
2377       }
2378
2379       /* Fisher. */
2380       if (f11 < 5. || f12 < 5. || f21 < 5. || f22 < 5.)
2381         calc_fisher (f11 + .5, f12 + .5, f21 + .5, f22 + .5, fisher1, fisher2);
2382     }
2383
2384   /* Calculate Mantel-Haenszel. */
2385   if (x->vars[ROW_VAR]->type == NUMERIC && x->vars[COL_VAR]->type == NUMERIC)
2386     {
2387       double r, ase_0, ase_1;
2388       calc_r ((double *) rows, (double *) cols, &r, &ase_0, &ase_1);
2389     
2390       chisq[4] = (W - 1.) * r * r;
2391       df[4] = 1;
2392     }
2393 }
2394
2395 /* Calculate the value of Pearson's r.  r is stored into R, ase_1 into
2396    ASE_1, and ase_0 into ASE_0.  The row and column values must be
2397    passed in X and Y. */
2398 static void
2399 calc_r (double *X, double *Y, double *r, double *ase_0, double *ase_1)
2400 {
2401   double SX, SY, S, T;
2402   double Xbar, Ybar;
2403   double sum_XYf, sum_X2Y2f;
2404   double sum_Xr, sum_X2r;
2405   double sum_Yc, sum_Y2c;
2406   int i, j;
2407
2408   for (sum_X2Y2f = sum_XYf = 0., i = 0; i < n_rows; i++)
2409     for (j = 0; j < n_cols; j++)
2410       {
2411         double fij = mat[j + i * n_cols];
2412         double product = X[i] * Y[j];
2413         double temp = fij * product;
2414         sum_XYf += temp;
2415         sum_X2Y2f += temp * product;
2416       }
2417
2418   for (sum_Xr = sum_X2r = 0., i = 0; i < n_rows; i++)
2419     {
2420       sum_Xr += X[i] * row_tot[i];
2421       sum_X2r += X[i] * X[i] * row_tot[i];
2422     }
2423   Xbar = sum_Xr / W;
2424
2425   for (sum_Yc = sum_Y2c = 0., i = 0; i < n_cols; i++)
2426     {
2427       sum_Yc += Y[i] * col_tot[i];
2428       sum_Y2c += Y[i] * Y[i] * col_tot[i];
2429     }
2430   Ybar = sum_Yc / W;
2431
2432   S = sum_XYf - sum_Xr * sum_Yc / W;
2433   SX = sum_X2r - sum_Xr * sum_Xr / W;
2434   SY = sum_Y2c - sum_Yc * sum_Yc / W;
2435   T = sqrt (SX * SY);
2436   *r = S / T;
2437   *ase_0 = sqrt ((sum_X2Y2f - (sum_XYf * sum_XYf) / W) / (sum_X2r * sum_Y2c));
2438   
2439   {
2440     double s, c, y, t;
2441     
2442     for (s = c = 0., i = 0; i < n_rows; i++)
2443       for (j = 0; j < n_cols; j++)
2444         {
2445           double Xresid, Yresid;
2446           double temp;
2447
2448           Xresid = X[i] - Xbar;
2449           Yresid = Y[j] - Ybar;
2450           temp = (T * Xresid * Yresid
2451                   - ((S / (2. * T))
2452                      * (Xresid * Xresid * SY + Yresid * Yresid * SX)));
2453           y = mat[j + i * n_cols] * temp * temp - c;
2454           t = s + y;
2455           c = (t - s) - y;
2456           s = t;
2457         }
2458     *ase_1 = sqrt (s) / (T * T);
2459   }
2460 }
2461
2462 static double somers_d_v[3];
2463 static double somers_d_ase[3];
2464 static double somers_d_t[3];
2465
2466 /* Calculate symmetric statistics and their asymptotic standard
2467    errors.  Returns 0 if none could be calculated. */
2468 static int
2469 calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
2470                 double t[N_SYMMETRIC])
2471 {
2472   int q = min (ns_rows, ns_cols);
2473   
2474   if (q <= 1)
2475     return 0;
2476   
2477   {
2478     int i;
2479
2480     if (v) 
2481       for (i = 0; i < N_SYMMETRIC; i++)
2482         v[i] = ase[i] = t[i] = SYSMIS;
2483   }
2484
2485   /* Phi, Cramer's V, contingency coefficient. */
2486   if (cmd.a_statistics[CRS_ST_PHI] || cmd.a_statistics[CRS_ST_CC])
2487     {
2488       double Xp = 0.;   /* Pearson chi-square. */
2489
2490       {
2491         int r, c;
2492     
2493         for (r = 0; r < n_rows; r++)
2494           for (c = 0; c < n_cols; c++)
2495             {
2496               const double expected = row_tot[r] * col_tot[c] / W;
2497               const double freq = mat[n_cols * r + c];
2498               const double residual = freq - expected;
2499     
2500               Xp += residual * residual / expected;
2501             }
2502       }
2503
2504       if (cmd.a_statistics[CRS_ST_PHI])
2505         {
2506           v[0] = sqrt (Xp / W);
2507           v[1] = sqrt (Xp / (W * (q - 1)));
2508         }
2509       if (cmd.a_statistics[CRS_ST_CC])
2510         v[2] = sqrt (Xp / (Xp + W));
2511     }
2512   
2513   if (cmd.a_statistics[CRS_ST_BTAU] || cmd.a_statistics[CRS_ST_CTAU]
2514       || cmd.a_statistics[CRS_ST_GAMMA] || cmd.a_statistics[CRS_ST_D])
2515     {
2516       double *cum;
2517       double Dr, Dc;
2518       double P, Q;
2519       double btau_cum, ctau_cum, gamma_cum, d_yx_cum, d_xy_cum;
2520       double btau_var;
2521       
2522       {
2523         int r, c;
2524         
2525         Dr = Dc = W * W;
2526         for (r = 0; r < n_rows; r++)
2527           Dr -= row_tot[r] * row_tot[r];
2528         for (c = 0; c < n_cols; c++)
2529           Dc -= col_tot[c] * col_tot[c];
2530       }
2531       
2532       {
2533         int r, c;
2534
2535         cum = xnmalloc (n_cols * n_rows, sizeof *cum);
2536         for (c = 0; c < n_cols; c++)
2537           {
2538             double ct = 0.;
2539             
2540             for (r = 0; r < n_rows; r++)
2541               cum[c + r * n_cols] = ct += mat[c + r * n_cols];
2542           }
2543       }
2544       
2545       /* P and Q. */
2546       {
2547         int i, j;
2548         double Cij, Dij;
2549
2550         P = Q = 0.;
2551         for (i = 0; i < n_rows; i++)
2552           {
2553             Cij = Dij = 0.;
2554
2555             for (j = 1; j < n_cols; j++)
2556               Cij += col_tot[j] - cum[j + i * n_cols];
2557
2558             if (i > 0)
2559               for (j = 1; j < n_cols; j++)
2560                 Dij += cum[j + (i - 1) * n_cols];
2561
2562             for (j = 0;;)
2563               {
2564                 double fij = mat[j + i * n_cols];
2565                 P += fij * Cij;
2566                 Q += fij * Dij;
2567                 
2568                 if (++j == n_cols)
2569                   break;
2570                 assert (j < n_cols);
2571
2572                 Cij -= col_tot[j] - cum[j + i * n_cols];
2573                 Dij += col_tot[j - 1] - cum[j - 1 + i * n_cols];
2574                 
2575                 if (i > 0)
2576                   {
2577                     Cij += cum[j - 1 + (i - 1) * n_cols];
2578                     Dij -= cum[j + (i - 1) * n_cols];
2579                   }
2580               }
2581           }
2582       }
2583
2584       if (cmd.a_statistics[CRS_ST_BTAU])
2585         v[3] = (P - Q) / sqrt (Dr * Dc);
2586       if (cmd.a_statistics[CRS_ST_CTAU])
2587         v[4] = (q * (P - Q)) / ((W * W) * (q - 1));
2588       if (cmd.a_statistics[CRS_ST_GAMMA])
2589         v[5] = (P - Q) / (P + Q);
2590
2591       /* ASE for tau-b, tau-c, gamma.  Calculations could be
2592          eliminated here, at expense of memory.  */
2593       {
2594         int i, j;
2595         double Cij, Dij;
2596
2597         btau_cum = ctau_cum = gamma_cum = d_yx_cum = d_xy_cum = 0.;
2598         for (i = 0; i < n_rows; i++)
2599           {
2600             Cij = Dij = 0.;
2601
2602             for (j = 1; j < n_cols; j++)
2603               Cij += col_tot[j] - cum[j + i * n_cols];
2604
2605             if (i > 0)
2606               for (j = 1; j < n_cols; j++)
2607                 Dij += cum[j + (i - 1) * n_cols];
2608
2609             for (j = 0;;)
2610               {
2611                 double fij = mat[j + i * n_cols];
2612
2613                 if (cmd.a_statistics[CRS_ST_BTAU])
2614                   {
2615                     const double temp = (2. * sqrt (Dr * Dc) * (Cij - Dij)
2616                                          + v[3] * (row_tot[i] * Dc
2617                                                    + col_tot[j] * Dr));
2618                     btau_cum += fij * temp * temp;
2619                   }
2620                 
2621                 {
2622                   const double temp = Cij - Dij;
2623                   ctau_cum += fij * temp * temp;
2624                 }
2625
2626                 if (cmd.a_statistics[CRS_ST_GAMMA])
2627                   {
2628                     const double temp = Q * Cij - P * Dij;
2629                     gamma_cum += fij * temp * temp;
2630                   }
2631
2632                 if (cmd.a_statistics[CRS_ST_D])
2633                   {
2634                     d_yx_cum += fij * pow2 (Dr * (Cij - Dij)
2635                                             - (P - Q) * (W - row_tot[i]));
2636                     d_xy_cum += fij * pow2 (Dc * (Dij - Cij)
2637                                             - (Q - P) * (W - col_tot[j]));
2638                   }
2639                 
2640                 if (++j == n_cols)
2641                   break;
2642                 assert (j < n_cols);
2643
2644                 Cij -= col_tot[j] - cum[j + i * n_cols];
2645                 Dij += col_tot[j - 1] - cum[j - 1 + i * n_cols];
2646                 
2647                 if (i > 0)
2648                   {
2649                     Cij += cum[j - 1 + (i - 1) * n_cols];
2650                     Dij -= cum[j + (i - 1) * n_cols];
2651                   }
2652               }
2653           }
2654       }
2655
2656       btau_var = ((btau_cum
2657                    - (W * pow2 (W * (P - Q) / sqrt (Dr * Dc) * (Dr + Dc))))
2658                   / pow2 (Dr * Dc));
2659       if (cmd.a_statistics[CRS_ST_BTAU])
2660         {
2661           ase[3] = sqrt (btau_var);
2662           t[3] = v[3] / (2 * sqrt ((ctau_cum - (P - Q) * (P - Q) / W)
2663                                    / (Dr * Dc)));
2664         }
2665       if (cmd.a_statistics[CRS_ST_CTAU])
2666         {
2667           ase[4] = ((2 * q / ((q - 1) * W * W))
2668                     * sqrt (ctau_cum - (P - Q) * (P - Q) / W));
2669           t[4] = v[4] / ase[4];
2670         }
2671       if (cmd.a_statistics[CRS_ST_GAMMA])
2672         {
2673           ase[5] = ((4. / ((P + Q) * (P + Q))) * sqrt (gamma_cum));
2674           t[5] = v[5] / (2. / (P + Q)
2675                          * sqrt (ctau_cum - (P - Q) * (P - Q) / W));
2676         }
2677       if (cmd.a_statistics[CRS_ST_D])
2678         {
2679           somers_d_v[0] = (P - Q) / (.5 * (Dc + Dr));
2680           somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
2681           somers_d_t[0] = (somers_d_v[0]
2682                            / (4 / (Dc + Dr)
2683                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2684           somers_d_v[1] = (P - Q) / Dc;
2685           somers_d_ase[1] = 2. / pow2 (Dc) * sqrt (d_xy_cum);
2686           somers_d_t[1] = (somers_d_v[1]
2687                            / (2. / Dc
2688                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2689           somers_d_v[2] = (P - Q) / Dr;
2690           somers_d_ase[2] = 2. / pow2 (Dr) * sqrt (d_yx_cum);
2691           somers_d_t[2] = (somers_d_v[2]
2692                            / (2. / Dr
2693                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2694         }
2695
2696       free (cum);
2697     }
2698
2699   /* Spearman correlation, Pearson's r. */
2700   if (cmd.a_statistics[CRS_ST_CORR])
2701     {
2702       double *R = local_alloc (sizeof *R * n_rows);
2703       double *C = local_alloc (sizeof *C * n_cols);
2704       
2705       {
2706         double y, t, c = 0., s = 0.;
2707         int i = 0;
2708         
2709         for (;;)
2710           {
2711             R[i] = s + (row_tot[i] + 1.) / 2.;
2712             y = row_tot[i] - c;
2713             t = s + y;
2714             c = (t - s) - y;
2715             s = t;
2716             if (++i == n_rows)
2717               break;
2718             assert (i < n_rows);
2719           }
2720       }
2721       
2722       {
2723         double y, t, c = 0., s = 0.;
2724         int j = 0;
2725         
2726         for (;;)
2727           {
2728             C[j] = s + (col_tot[j] + 1.) / 2;
2729             y = col_tot[j] - c;
2730             t = s + y;
2731             c = (t - s) - y;
2732             s = t;
2733             if (++j == n_cols)
2734               break;
2735             assert (j < n_cols);
2736           }
2737       }
2738       
2739       calc_r (R, C, &v[6], &t[6], &ase[6]);
2740       t[6] = v[6] / t[6];
2741
2742       local_free (R);
2743       local_free (C);
2744
2745       calc_r ((double *) rows, (double *) cols, &v[7], &t[7], &ase[7]);
2746       t[7] = v[7] / t[7];
2747     }
2748
2749   /* Cohen's kappa. */
2750   if (cmd.a_statistics[CRS_ST_KAPPA] && ns_rows == ns_cols)
2751     {
2752       double sum_fii, sum_rici, sum_fiiri_ci, sum_fijri_ci2, sum_riciri_ci;
2753       int i, j;
2754       
2755       for (sum_fii = sum_rici = sum_fiiri_ci = sum_riciri_ci = 0., i = j = 0;
2756            i < ns_rows; i++, j++)
2757         {
2758           double prod, sum;
2759           
2760           while (col_tot[j] == 0.)
2761             j++;
2762           
2763           prod = row_tot[i] * col_tot[j];
2764           sum = row_tot[i] + col_tot[j];
2765           
2766           sum_fii += mat[j + i * n_cols];
2767           sum_rici += prod;
2768           sum_fiiri_ci += mat[j + i * n_cols] * sum;
2769           sum_riciri_ci += prod * sum;
2770         }
2771       for (sum_fijri_ci2 = 0., i = 0; i < ns_rows; i++)
2772         for (j = 0; j < ns_cols; j++)
2773           {
2774             double sum = row_tot[i] + col_tot[j];
2775             sum_fijri_ci2 += mat[j + i * n_cols] * sum * sum;
2776           }
2777       
2778       v[8] = (W * sum_fii - sum_rici) / (W * W - sum_rici);
2779
2780       ase[8] = sqrt ((W * W * sum_rici
2781                       + sum_rici * sum_rici
2782                       - W * sum_riciri_ci)
2783                      / (W * (W * W - sum_rici) * (W * W - sum_rici)));
2784 #if 0
2785       t[8] = v[8] / sqrt (W * (((sum_fii * (W - sum_fii))
2786                                 / pow2 (W * W - sum_rici))
2787                                + ((2. * (W - sum_fii)
2788                                    * (2. * sum_fii * sum_rici
2789                                       - W * sum_fiiri_ci))
2790                                   / cube (W * W - sum_rici))
2791                                + (pow2 (W - sum_fii)
2792                                   * (W * sum_fijri_ci2 - 4.
2793                                      * sum_rici * sum_rici)
2794                                   / pow4 (W * W - sum_rici))));
2795 #else
2796       t[8] = v[8] / ase[8];
2797 #endif
2798     }
2799
2800   return 1;
2801 }
2802
2803 /* Calculate risk estimate. */
2804 static int
2805 calc_risk (double *value, double *upper, double *lower, union value *c)
2806 {
2807   double f11, f12, f21, f22;
2808   double v;
2809
2810   {
2811     int i;
2812       
2813     for (i = 0; i < 3; i++)
2814       value[i] = upper[i] = lower[i] = SYSMIS;
2815   }
2816     
2817   if (ns_rows != 2 || ns_cols != 2)
2818     return 0;
2819   
2820   {
2821     int nz_cols[2];
2822     int i, j;
2823
2824     for (i = j = 0; i < n_cols; i++)
2825       if (col_tot[i] != 0.)
2826         {
2827           nz_cols[j++] = i;
2828           if (j == 2)
2829             break;
2830         }
2831
2832     assert (j == 2);
2833
2834     f11 = mat[nz_cols[0]];
2835     f12 = mat[nz_cols[1]];
2836     f21 = mat[nz_cols[0] + n_cols];
2837     f22 = mat[nz_cols[1] + n_cols];
2838
2839     c[0] = cols[nz_cols[0]];
2840     c[1] = cols[nz_cols[1]];
2841   }
2842
2843   value[0] = (f11 * f22) / (f12 * f21);
2844   v = sqrt (1. / f11 + 1. / f12 + 1. / f21 + 1. / f22);
2845   lower[0] = value[0] * exp (-1.960 * v);
2846   upper[0] = value[0] * exp (1.960 * v);
2847
2848   value[1] = (f11 * (f21 + f22)) / (f21 * (f11 + f12));
2849   v = sqrt ((f12 / (f11 * (f11 + f12)))
2850             + (f22 / (f21 * (f21 + f22))));
2851   lower[1] = value[1] * exp (-1.960 * v);
2852   upper[1] = value[1] * exp (1.960 * v);
2853     
2854   value[2] = (f12 * (f21 + f22)) / (f22 * (f11 + f12));
2855   v = sqrt ((f11 / (f12 * (f11 + f12)))
2856             + (f21 / (f22 * (f21 + f22))));
2857   lower[2] = value[2] * exp (-1.960 * v);
2858   upper[2] = value[2] * exp (1.960 * v);
2859
2860   return 1;
2861 }
2862
2863 /* Calculate directional measures. */
2864 static int
2865 calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
2866                   double t[N_DIRECTIONAL])
2867 {
2868   {
2869     int i;
2870
2871     for (i = 0; i < N_DIRECTIONAL; i++)
2872       v[i] = ase[i] = t[i] = SYSMIS;
2873   }
2874
2875   /* Lambda. */
2876   if (cmd.a_statistics[CRS_ST_LAMBDA])
2877     {
2878       double *fim = xnmalloc (n_rows, sizeof *fim);
2879       int *fim_index = xnmalloc (n_rows, sizeof *fim_index);
2880       double *fmj = xnmalloc (n_cols, sizeof *fmj);
2881       int *fmj_index = xnmalloc (n_cols, sizeof *fmj_index);
2882       double sum_fim, sum_fmj;
2883       double rm, cm;
2884       int rm_index, cm_index;
2885       int i, j;
2886
2887       /* Find maximum for each row and their sum. */
2888       for (sum_fim = 0., i = 0; i < n_rows; i++)
2889         {
2890           double max = mat[i * n_cols];
2891           int index = 0;
2892
2893           for (j = 1; j < n_cols; j++)
2894             if (mat[j + i * n_cols] > max)
2895               {
2896                 max = mat[j + i * n_cols];
2897                 index = j;
2898               }
2899         
2900           sum_fim += fim[i] = max;
2901           fim_index[i] = index;
2902         }
2903
2904       /* Find maximum for each column. */
2905       for (sum_fmj = 0., j = 0; j < n_cols; j++)
2906         {
2907           double max = mat[j];
2908           int index = 0;
2909
2910           for (i = 1; i < n_rows; i++)
2911             if (mat[j + i * n_cols] > max)
2912               {
2913                 max = mat[j + i * n_cols];
2914                 index = i;
2915               }
2916         
2917           sum_fmj += fmj[j] = max;
2918           fmj_index[j] = index;
2919         }
2920
2921       /* Find maximum row total. */
2922       rm = row_tot[0];
2923       rm_index = 0;
2924       for (i = 1; i < n_rows; i++)
2925         if (row_tot[i] > rm)
2926           {
2927             rm = row_tot[i];
2928             rm_index = i;
2929           }
2930
2931       /* Find maximum column total. */
2932       cm = col_tot[0];
2933       cm_index = 0;
2934       for (j = 1; j < n_cols; j++)
2935         if (col_tot[j] > cm)
2936           {
2937             cm = col_tot[j];
2938             cm_index = j;
2939           }
2940
2941       v[0] = (sum_fim + sum_fmj - cm - rm) / (2. * W - rm - cm);
2942       v[1] = (sum_fmj - rm) / (W - rm);
2943       v[2] = (sum_fim - cm) / (W - cm);
2944
2945       /* ASE1 for Y given X. */
2946       {
2947         double accum;
2948
2949         for (accum = 0., i = 0; i < n_rows; i++)
2950           for (j = 0; j < n_cols; j++)
2951             {
2952               const int deltaj = j == cm_index;
2953               accum += (mat[j + i * n_cols]
2954                         * pow2 ((j == fim_index[i])
2955                                - deltaj
2956                                + v[0] * deltaj));
2957             }
2958       
2959         ase[2] = sqrt (accum - W * v[0]) / (W - cm);
2960       }
2961
2962       /* ASE0 for Y given X. */
2963       {
2964         double accum;
2965       
2966         for (accum = 0., i = 0; i < n_rows; i++)
2967           if (cm_index != fim_index[i])
2968             accum += (mat[i * n_cols + fim_index[i]]
2969                       + mat[i * n_cols + cm_index]);
2970         t[2] = v[2] / (sqrt (accum - pow2 (sum_fim - cm) / W) / (W - cm));
2971       }
2972
2973       /* ASE1 for X given Y. */
2974       {
2975         double accum;
2976
2977         for (accum = 0., i = 0; i < n_rows; i++)
2978           for (j = 0; j < n_cols; j++)
2979             {
2980               const int deltaj = i == rm_index;
2981               accum += (mat[j + i * n_cols]
2982                         * pow2 ((i == fmj_index[j])
2983                                - deltaj
2984                                + v[0] * deltaj));
2985             }
2986       
2987         ase[1] = sqrt (accum - W * v[0]) / (W - rm);
2988       }
2989
2990       /* ASE0 for X given Y. */
2991       {
2992         double accum;
2993       
2994         for (accum = 0., j = 0; j < n_cols; j++)
2995           if (rm_index != fmj_index[j])
2996             accum += (mat[j + n_cols * fmj_index[j]]
2997                       + mat[j + n_cols * rm_index]);
2998         t[1] = v[1] / (sqrt (accum - pow2 (sum_fmj - rm) / W) / (W - rm));
2999       }
3000
3001       /* Symmetric ASE0 and ASE1. */
3002       {
3003         double accum0;
3004         double accum1;
3005
3006         for (accum0 = accum1 = 0., i = 0; i < n_rows; i++)
3007           for (j = 0; j < n_cols; j++)
3008             {
3009               int temp0 = (fmj_index[j] == i) + (fim_index[i] == j);
3010               int temp1 = (i == rm_index) + (j == cm_index);
3011               accum0 += mat[j + i * n_cols] * pow2 (temp0 - temp1);
3012               accum1 += (mat[j + i * n_cols]
3013                          * pow2 (temp0 + (v[0] - 1.) * temp1));
3014             }
3015         ase[0] = sqrt (accum1 - 4. * W * v[0] * v[0]) / (2. * W - rm - cm);
3016         t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / W))
3017                        / (2. * W - rm - cm));
3018       }
3019
3020       free (fim);
3021       free (fim_index);
3022       free (fmj);
3023       free (fmj_index);
3024       
3025       {
3026         double sum_fij2_ri, sum_fij2_ci;
3027         double sum_ri2, sum_cj2;
3028
3029         for (sum_fij2_ri = sum_fij2_ci = 0., i = 0; i < n_rows; i++)
3030           for (j = 0; j < n_cols; j++)
3031             {
3032               double temp = pow2 (mat[j + i * n_cols]);
3033               sum_fij2_ri += temp / row_tot[i];
3034               sum_fij2_ci += temp / col_tot[j];
3035             }
3036
3037         for (sum_ri2 = 0., i = 0; i < n_rows; i++)
3038           sum_ri2 += row_tot[i] * row_tot[i];
3039
3040         for (sum_cj2 = 0., j = 0; j < n_cols; j++)
3041           sum_cj2 += col_tot[j] * col_tot[j];
3042
3043         v[3] = (W * sum_fij2_ci - sum_ri2) / (W * W - sum_ri2);
3044         v[4] = (W * sum_fij2_ri - sum_cj2) / (W * W - sum_cj2);
3045       }
3046     }
3047
3048   if (cmd.a_statistics[CRS_ST_UC])
3049     {
3050       double UX, UY, UXY, P;
3051       double ase1_yx, ase1_xy, ase1_sym;
3052       int i, j;
3053
3054       for (UX = 0., i = 0; i < n_rows; i++)
3055         if (row_tot[i] > 0.)
3056           UX -= row_tot[i] / W * log (row_tot[i] / W);
3057       
3058       for (UY = 0., j = 0; j < n_cols; j++)
3059         if (col_tot[j] > 0.)
3060           UY -= col_tot[j] / W * log (col_tot[j] / W);
3061
3062       for (UXY = P = 0., i = 0; i < n_rows; i++)
3063         for (j = 0; j < n_cols; j++)
3064           {
3065             double entry = mat[j + i * n_cols];
3066
3067             if (entry <= 0.)
3068               continue;
3069             
3070             P += entry * pow2 (log (col_tot[j] * row_tot[i] / (W * entry)));
3071             UXY -= entry / W * log (entry / W);
3072           }
3073
3074       for (ase1_yx = ase1_xy = ase1_sym = 0., i = 0; i < n_rows; i++)
3075         for (j = 0; j < n_cols; j++)
3076           {
3077             double entry = mat[j + i * n_cols];
3078
3079             if (entry <= 0.)
3080               continue;
3081             
3082             ase1_yx += entry * pow2 (UY * log (entry / row_tot[i])
3083                                     + (UX - UXY) * log (col_tot[j] / W));
3084             ase1_xy += entry * pow2 (UX * log (entry / col_tot[j])
3085                                     + (UY - UXY) * log (row_tot[i] / W));
3086             ase1_sym += entry * pow2 ((UXY
3087                                       * log (row_tot[i] * col_tot[j] / (W * W)))
3088                                      - (UX + UY) * log (entry / W));
3089           }
3090       
3091       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
3092       ase[5] = (2. / (W * pow2 (UX + UY))) * sqrt (ase1_sym);
3093       t[5] = v[5] / ((2. / (W * (UX + UY)))
3094                      * sqrt (P - pow2 (UX + UY - UXY) / W));
3095                     
3096       v[6] = (UX + UY - UXY) / UX;
3097       ase[6] = sqrt (ase1_xy) / (W * UX * UX);
3098       t[6] = v[6] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UX));
3099       
3100       v[7] = (UX + UY - UXY) / UY;
3101       ase[7] = sqrt (ase1_yx) / (W * UY * UY);
3102       t[7] = v[7] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UY));
3103     }
3104
3105   /* Somers' D. */
3106   if (cmd.a_statistics[CRS_ST_D])
3107     {
3108       int i;
3109       
3110       if (!sym)
3111         calc_symmetric (NULL, NULL, NULL);
3112       for (i = 0; i < 3; i++)
3113         {
3114           v[8 + i] = somers_d_v[i];
3115           ase[8 + i] = somers_d_ase[i];
3116           t[8 + i] = somers_d_t[i];
3117         }
3118     }
3119
3120   /* Eta. */
3121   if (cmd.a_statistics[CRS_ST_ETA])
3122     {
3123       {
3124         double sum_Xr, sum_X2r;
3125         double SX, SXW;
3126         int i, j;
3127       
3128         for (sum_Xr = sum_X2r = 0., i = 0; i < n_rows; i++)
3129           {
3130             sum_Xr += rows[i].f * row_tot[i];
3131             sum_X2r += rows[i].f * rows[i].f * row_tot[i];
3132           }
3133         SX = sum_X2r - sum_Xr * sum_Xr / W;
3134       
3135         for (SXW = 0., j = 0; j < n_cols; j++)
3136           {
3137             double cum;
3138
3139             for (cum = 0., i = 0; i < n_rows; i++)
3140               {
3141                 SXW += rows[i].f * rows[i].f * mat[j + i * n_cols];
3142                 cum += rows[i].f * mat[j + i * n_cols];
3143               }
3144
3145             SXW -= cum * cum / col_tot[j];
3146           }
3147         v[11] = sqrt (1. - SXW / SX);
3148       }
3149
3150       {
3151         double sum_Yc, sum_Y2c;
3152         double SY, SYW;
3153         int i, j;
3154
3155         for (sum_Yc = sum_Y2c = 0., i = 0; i < n_cols; i++)
3156           {
3157             sum_Yc += cols[i].f * col_tot[i];
3158             sum_Y2c += cols[i].f * cols[i].f * col_tot[i];
3159           }
3160         SY = sum_Y2c - sum_Yc * sum_Yc / W;
3161
3162         for (SYW = 0., i = 0; i < n_rows; i++)
3163           {
3164             double cum;
3165
3166             for (cum = 0., j = 0; j < n_cols; j++)
3167               {
3168                 SYW += cols[j].f * cols[j].f * mat[j + i * n_cols];
3169                 cum += cols[j].f * mat[j + i * n_cols];
3170               }
3171           
3172             SYW -= cum * cum / row_tot[i];
3173           }
3174         v[12] = sqrt (1. - SYW / SY);
3175       }
3176     }
3177
3178   return 1;
3179 }
3180
3181 /* A wrapper around data_out() that limits string output to short
3182    string width and null terminates the result. */
3183 static void
3184 format_short (char *s, const struct fmt_spec *fp, const union value *v)
3185 {
3186   struct fmt_spec fmt_subst;
3187
3188   /* Limit to short string width. */
3189   if (formats[fp->type].cat & FCAT_STRING) 
3190     {
3191       fmt_subst = *fp;
3192
3193       assert (fmt_subst.type == FMT_A || fmt_subst.type == FMT_AHEX);
3194       if (fmt_subst.type == FMT_A)
3195         fmt_subst.w = min (8, fmt_subst.w);
3196       else
3197         fmt_subst.w = min (16, fmt_subst.w);
3198
3199       fp = &fmt_subst;
3200     }
3201
3202   /* Format. */
3203   data_out (s, fp, v);
3204   
3205   /* Null terminate. */
3206   s[fp->w] = '\0';
3207 }
3208
3209 /* 
3210    Local Variables:
3211    mode: c
3212    End:
3213 */