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