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, _("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, "%s", 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, _("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, _("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, _("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, _("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_GAP, 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", OUTP_PROPORTIONAL);
1485   if (cmd.miss == CRS_REPORT)
1486     c += outp_string_width (d, "M", OUTP_PROPORTIONAL);
1487
1488   /* Set width for header columns. */
1489   if (t->l != 0)
1490     {
1491       size_t i;
1492       int w;
1493
1494       w = d->width - c * (t->nc - t->l);
1495       for (i = 0; i <= t->nc; i++)
1496         w -= t->wrv[i];
1497       w /= t->l;
1498
1499       if (w < d->prop_em_width * 8)
1500         w = d->prop_em_width * 8;
1501
1502       if (w > d->prop_em_width * 15)
1503         w = d->prop_em_width * 15;
1504
1505       for (i = 0; i < t->l; i++)
1506         t->w[i] = w;
1507     }
1508
1509   for (i = t->l; i < t->nc; i++)
1510     t->w[i] = c;
1511
1512   for (i = 0; i < t->nr; i++)
1513     t->h[i] = tab_natural_height (t, d, i);
1514 }
1515
1516 static struct table_entry **find_pivot_extent_general (struct table_entry **tp,
1517                                                 int *cnt, int pivot);
1518 static struct table_entry **find_pivot_extent_integer (struct table_entry **tp,
1519                                                 int *cnt, int pivot);
1520
1521 /* Calls find_pivot_extent_general or find_pivot_extent_integer, as
1522    appropriate. */
1523 static struct table_entry **
1524 find_pivot_extent (struct table_entry **tp, int *cnt, int pivot)
1525 {
1526   return (mode == GENERAL
1527           ? find_pivot_extent_general (tp, cnt, pivot)
1528           : find_pivot_extent_integer (tp, cnt, pivot));
1529 }
1530
1531 /* Find the extent of a region in TP that contains one table.  If
1532    PIVOT != 0 that means a set of table entries with identical table
1533    number; otherwise they also have to have the same values for every
1534    dimension after the row and column dimensions.  The table that is
1535    searched starts at TP and has length CNT.  Returns the first entry
1536    after the last one in the table; sets *CNT to the number of
1537    remaining values.  If there are no entries in TP at all, returns
1538    NULL.  A yucky interface, admittedly, but it works. */
1539 static struct table_entry **
1540 find_pivot_extent_general (struct table_entry **tp, int *cnt, int pivot)
1541 {
1542   struct table_entry *fp = *tp;
1543   struct crosstab *x;
1544
1545   if (*cnt == 0)
1546     return NULL;
1547   x = xtab[(*tp)->table];
1548   for (;;)
1549     {
1550       tp++;
1551       if (--*cnt == 0)
1552         break;
1553       assert (*cnt > 0);
1554
1555       if ((*tp)->table != fp->table)
1556         break;
1557       if (pivot)
1558         continue;
1559
1560       if (memcmp (&(*tp)->values[2], &fp->values[2], sizeof (union value) * (x->nvar - 2)))
1561         break;
1562     }
1563
1564   return tp;
1565 }
1566
1567 /* Integer mode correspondent to find_pivot_extent_general().  This
1568    could be optimized somewhat, but I just don't give a crap about
1569    CROSSTABS performance in integer mode, which is just a
1570    CROSSTABS wart as far as I'm concerned.
1571
1572    That said, feel free to send optimization patches to me. */
1573 static struct table_entry **
1574 find_pivot_extent_integer (struct table_entry **tp, int *cnt, int pivot)
1575 {
1576   struct table_entry *fp = *tp;
1577   struct crosstab *x;
1578
1579   if (*cnt == 0)
1580     return NULL;
1581   x = xtab[(*tp)->table];
1582   for (;;)
1583     {
1584       tp++;
1585       if (--*cnt == 0)
1586         break;
1587       assert (*cnt > 0);
1588
1589       if ((*tp)->table != fp->table)
1590         break;
1591       if (pivot)
1592         continue;
1593
1594       if (memcmp (&(*tp)->values[2], &fp->values[2],
1595                   sizeof (union value) * (x->nvar - 2)))
1596         break;
1597     }
1598
1599   return tp;
1600 }
1601
1602 /* Compares `union value's A_ and B_ and returns a strcmp()-like
1603    result.  WIDTH_ points to an int which is either 0 for a
1604    numeric value or a string width for a string value. */
1605 static int
1606 compare_value (const void *a_, const void *b_, void *width_)
1607 {
1608   const union value *a = a_;
1609   const union value *b = b_;
1610   const int *pwidth = width_;
1611   const int width = *pwidth;
1612
1613   if (width == 0)
1614     return (a->f < b->f) ? -1 : (a->f > b->f);
1615   else
1616     return strncmp (a->s, b->s, width);
1617 }
1618
1619 /* Given an array of ENTRY_CNT table_entry structures starting at
1620    ENTRIES, creates a sorted list of the values that the variable
1621    with index VAR_IDX takes on.  The values are returned as a
1622    malloc()'darray stored in *VALUES, with the number of values
1623    stored in *VALUE_CNT.
1624    */
1625 static void
1626 enum_var_values (struct table_entry **entries, int entry_cnt, int var_idx,
1627                  union value **values, int *value_cnt)
1628 {
1629   struct variable *v = xtab[(*entries)->table]->vars[var_idx];
1630
1631   if (mode == GENERAL)
1632     {
1633       int width = v->width;
1634       int i;
1635
1636       *values = xnmalloc (entry_cnt, sizeof **values);
1637       for (i = 0; i < entry_cnt; i++)
1638         (*values)[i] = entries[i]->values[var_idx];
1639       *value_cnt = sort_unique (*values, entry_cnt, sizeof **values,
1640                                 compare_value, &width);
1641     }
1642   else
1643     {
1644       struct var_range *vr = get_var_range (v);
1645       int i;
1646
1647       assert (mode == INTEGER);
1648       *values = xnmalloc (vr->count, sizeof **values);
1649       for (i = 0; i < vr->count; i++)
1650         (*values)[i].f = i + vr->min;
1651       *value_cnt = vr->count;
1652     }
1653 }
1654
1655 /* Sets cell (C,R) in TABLE, with options OPT, to have a value taken
1656    from V, displayed with print format spec from variable VAR.  When
1657    in REPORT missing-value mode, missing values have an M appended. */
1658 static void
1659 table_value_missing (struct tab_table *table, int c, int r, unsigned char opt,
1660                      const union value *v, const struct variable *var)
1661 {
1662   struct fixed_string s;
1663
1664   const char *label = val_labs_find (var->val_labs, *v);
1665   if (label)
1666     {
1667       tab_text (table, c, r, TAB_LEFT, label);
1668       return;
1669     }
1670
1671   s.string = tab_alloc (table, var->print.w);
1672   format_short (s.string, &var->print, v);
1673   s.length = strlen (s.string);
1674   if (cmd.miss == CRS_REPORT && mv_is_num_user_missing (&var->miss, v->f))
1675     s.string[s.length++] = 'M';
1676   while (s.length && *s.string == ' ')
1677     {
1678       s.length--;
1679       s.string++;
1680     }
1681   tab_raw (table, c, r, opt, &s);
1682 }
1683
1684 /* Draws a line across TABLE at the current row to indicate the most
1685    major dimension variable with index FIRST_DIFFERENCE out of NVAR
1686    that changed, and puts the values that changed into the table.  TB
1687    and X must be the corresponding table_entry and crosstab,
1688    respectively. */
1689 static void
1690 display_dimensions (struct tab_table *table, int first_difference, struct table_entry *tb)
1691 {
1692   tab_hline (table, TAL_1, nvar - first_difference - 1, tab_nc (table) - 1, 0);
1693
1694   for (; first_difference >= 2; first_difference--)
1695     table_value_missing (table, nvar - first_difference - 1, 0,
1696                          TAB_RIGHT, &tb->values[first_difference],
1697                          x->vars[first_difference]);
1698 }
1699
1700 /* Put VALUE into cell (C,R) of TABLE, suffixed with character
1701    SUFFIX if nonzero.  If MARK_MISSING is nonzero the entry is
1702    additionally suffixed with a letter `M'. */
1703 static void
1704 format_cell_entry (struct tab_table *table, int c, int r, double value,
1705                    char suffix, int mark_missing)
1706 {
1707   const struct fmt_spec f = {FMT_F, 10, 1};
1708   union value v;
1709   struct fixed_string s;
1710
1711   s.length = 10;
1712   s.string = tab_alloc (table, 16);
1713   v.f = value;
1714   data_out (s.string, &f, &v);
1715   while (*s.string == ' ')
1716     {
1717       s.length--;
1718       s.string++;
1719     }
1720   if (suffix != 0)
1721     s.string[s.length++] = suffix;
1722   if (mark_missing)
1723     s.string[s.length++] = 'M';
1724
1725   tab_raw (table, c, r, TAB_RIGHT, &s);
1726 }
1727
1728 /* Displays the crosstabulation table. */
1729 static void
1730 display_crosstabulation (void)
1731 {
1732   {
1733     int r;
1734
1735     for (r = 0; r < n_rows; r++)
1736       table_value_missing (table, nvar - 2, r * num_cells,
1737                            TAB_RIGHT, &rows[r], x->vars[ROW_VAR]);
1738   }
1739   tab_text (table, nvar - 2, n_rows * num_cells,
1740             TAB_LEFT, _("Total"));
1741
1742   /* Put in the actual cells. */
1743   {
1744     double *mp = mat;
1745     int r, c, i;
1746
1747     tab_offset (table, nvar - 1, -1);
1748     for (r = 0; r < n_rows; r++)
1749       {
1750         if (num_cells > 1)
1751           tab_hline (table, TAL_1, -1, n_cols, 0);
1752         for (c = 0; c < n_cols; c++)
1753           {
1754             int mark_missing = 0;
1755             double expected_value = row_tot[r] * col_tot[c] / W;
1756             if (cmd.miss == CRS_REPORT
1757                 && (mv_is_num_user_missing (&x->vars[COL_VAR]->miss, cols[c].f)
1758                     || mv_is_num_user_missing (&x->vars[ROW_VAR]->miss,
1759                                                rows[r].f)))
1760               mark_missing = 1;
1761             for (i = 0; i < num_cells; i++)
1762               {
1763                 double v;
1764                 int suffix = 0;
1765
1766                 switch (cells[i])
1767                   {
1768                   case CRS_CL_COUNT:
1769                     v = *mp;
1770                     break;
1771                   case CRS_CL_ROW:
1772                     v = *mp / row_tot[r] * 100.;
1773                     suffix = '%';
1774                     break;
1775                   case CRS_CL_COLUMN:
1776                     v = *mp / col_tot[c] * 100.;
1777                     suffix = '%';
1778                     break;
1779                   case CRS_CL_TOTAL:
1780                     v = *mp / W * 100.;
1781                     suffix = '%';
1782                     break;
1783                   case CRS_CL_EXPECTED:
1784                     v = expected_value;
1785                     break;
1786                   case CRS_CL_RESIDUAL:
1787                     v = *mp - expected_value;
1788                     break;
1789                   case CRS_CL_SRESIDUAL:
1790                     v = (*mp - expected_value) / sqrt (expected_value);
1791                     break;
1792                   case CRS_CL_ASRESIDUAL:
1793                     v = ((*mp - expected_value)
1794                          / sqrt (expected_value
1795                                  * (1. - row_tot[r] / W)
1796                                  * (1. - col_tot[c] / W)));
1797                     break;
1798                   default:
1799                     assert (0);
1800                     abort ();
1801                   }
1802
1803                 format_cell_entry (table, c, i, v, suffix, mark_missing);
1804               }
1805
1806             mp++;
1807           }
1808
1809         tab_offset (table, -1, tab_row (table) + num_cells);
1810       }
1811   }
1812
1813   /* Row totals. */
1814   {
1815     int r, i;
1816
1817     tab_offset (table, -1, tab_row (table) - num_cells * n_rows);
1818     for (r = 0; r < n_rows; r++)
1819       {
1820         char suffix = 0;
1821         int mark_missing = 0;
1822
1823         if (cmd.miss == CRS_REPORT
1824             && mv_is_num_user_missing (&x->vars[ROW_VAR]->miss, rows[r].f))
1825           mark_missing = 1;
1826
1827         for (i = 0; i < num_cells; i++)
1828           {
1829             double v;
1830
1831             switch (cells[i])
1832               {
1833               case CRS_CL_COUNT:
1834                 v = row_tot[r];
1835                 break;
1836               case CRS_CL_ROW:
1837                 v = 100.;
1838                 suffix = '%';
1839                 break;
1840               case CRS_CL_COLUMN:
1841                 v = row_tot[r] / W * 100.;
1842                 suffix = '%';
1843                 break;
1844               case CRS_CL_TOTAL:
1845                 v = row_tot[r] / W * 100.;
1846                 suffix = '%';
1847                 break;
1848               case CRS_CL_EXPECTED:
1849               case CRS_CL_RESIDUAL:
1850               case CRS_CL_SRESIDUAL:
1851               case CRS_CL_ASRESIDUAL:
1852                 v = 0.;
1853                 break;
1854               default:
1855                 assert (0);
1856                 abort ();
1857               }
1858
1859             format_cell_entry (table, n_cols, 0, v, suffix, mark_missing);
1860             tab_next_row (table);
1861           }
1862       }
1863   }
1864
1865   /* Column totals, grand total. */
1866   {
1867     int c;
1868     int last_row = 0;
1869
1870     if (num_cells > 1)
1871       tab_hline (table, TAL_1, -1, n_cols, 0);
1872     for (c = 0; c <= n_cols; c++)
1873       {
1874         double ct = c < n_cols ? col_tot[c] : W;
1875         int mark_missing = 0;
1876         char suffix = 0;
1877         int i;
1878
1879         if (cmd.miss == CRS_REPORT && c < n_cols
1880             && mv_is_num_user_missing (&x->vars[COL_VAR]->miss, cols[c].f))
1881           mark_missing = 1;
1882
1883         for (i = 0; i < num_cells; i++)
1884           {
1885             double v;
1886
1887             switch (cells[i])
1888               {
1889               case CRS_CL_COUNT:
1890                 v = ct;
1891                 suffix = '%';
1892                 break;
1893               case CRS_CL_ROW:
1894                 v = ct / W * 100.;
1895                 suffix = '%';
1896                 break;
1897               case CRS_CL_COLUMN:
1898                 v = 100.;
1899                 suffix = '%';
1900                 break;
1901               case CRS_CL_TOTAL:
1902                 v = ct / W * 100.;
1903                 suffix = '%';
1904                 break;
1905               case CRS_CL_EXPECTED:
1906               case CRS_CL_RESIDUAL:
1907               case CRS_CL_SRESIDUAL:
1908               case CRS_CL_ASRESIDUAL:
1909                 continue;
1910               default:
1911                 assert (0);
1912                 abort ();
1913               }
1914
1915             format_cell_entry (table, c, i, v, suffix, mark_missing);
1916           }
1917         last_row = i;
1918       }
1919
1920     tab_offset (table, -1, tab_row (table) + last_row);
1921   }
1922
1923   tab_offset (table, 0, -1);
1924 }
1925
1926 static void calc_r (double *X, double *Y, double *, double *, double *);
1927 static void calc_chisq (double[N_CHISQ], int[N_CHISQ], double *, double *);
1928
1929 /* Display chi-square statistics. */
1930 static void
1931 display_chisq (void)
1932 {
1933   static const char *chisq_stats[N_CHISQ] =
1934     {
1935       N_("Pearson Chi-Square"),
1936       N_("Likelihood Ratio"),
1937       N_("Fisher's Exact Test"),
1938       N_("Continuity Correction"),
1939       N_("Linear-by-Linear Association"),
1940     };
1941   double chisq_v[N_CHISQ];
1942   double fisher1, fisher2;
1943   int df[N_CHISQ];
1944   int s = 0;
1945
1946   int i;
1947
1948   calc_chisq (chisq_v, df, &fisher1, &fisher2);
1949
1950   tab_offset (chisq, nvar - 2, -1);
1951
1952   for (i = 0; i < N_CHISQ; i++)
1953     {
1954       if ((i != 2 && chisq_v[i] == SYSMIS)
1955           || (i == 2 && fisher1 == SYSMIS))
1956         continue;
1957       s = 1;
1958
1959       tab_text (chisq, 0, 0, TAB_LEFT, gettext (chisq_stats[i]));
1960       if (i != 2)
1961         {
1962           tab_float (chisq, 1, 0, TAB_RIGHT, chisq_v[i], 8, 3);
1963           tab_float (chisq, 2, 0, TAB_RIGHT, df[i], 8, 0);
1964           tab_float (chisq, 3, 0, TAB_RIGHT,
1965                      gsl_cdf_chisq_Q (chisq_v[i], df[i]), 8, 3);
1966         }
1967       else
1968         {
1969           chisq_fisher = 1;
1970           tab_float (chisq, 4, 0, TAB_RIGHT, fisher2, 8, 3);
1971           tab_float (chisq, 5, 0, TAB_RIGHT, fisher1, 8, 3);
1972         }
1973       tab_next_row (chisq);
1974     }
1975
1976   tab_text (chisq, 0, 0, TAB_LEFT, _("N of Valid Cases"));
1977   tab_float (chisq, 1, 0, TAB_RIGHT, W, 8, 0);
1978   tab_next_row (chisq);
1979
1980   tab_offset (chisq, 0, -1);
1981 }
1982
1983 static int calc_symmetric (double[N_SYMMETRIC], double[N_SYMMETRIC],
1984                            double[N_SYMMETRIC]);
1985
1986 /* Display symmetric measures. */
1987 static void
1988 display_symmetric (void)
1989 {
1990   static const char *categories[] =
1991     {
1992       N_("Nominal by Nominal"),
1993       N_("Ordinal by Ordinal"),
1994       N_("Interval by Interval"),
1995       N_("Measure of Agreement"),
1996     };
1997
1998   static const char *stats[N_SYMMETRIC] =
1999     {
2000       N_("Phi"),
2001       N_("Cramer's V"),
2002       N_("Contingency Coefficient"),
2003       N_("Kendall's tau-b"),
2004       N_("Kendall's tau-c"),
2005       N_("Gamma"),
2006       N_("Spearman Correlation"),
2007       N_("Pearson's R"),
2008       N_("Kappa"),
2009     };
2010
2011   static const int stats_categories[N_SYMMETRIC] =
2012     {
2013       0, 0, 0, 1, 1, 1, 1, 2, 3,
2014     };
2015
2016   int last_cat = -1;
2017   double sym_v[N_SYMMETRIC], sym_ase[N_SYMMETRIC], sym_t[N_SYMMETRIC];
2018   int i;
2019
2020   if (!calc_symmetric (sym_v, sym_ase, sym_t))
2021     return;
2022
2023   tab_offset (sym, nvar - 2, -1);
2024
2025   for (i = 0; i < N_SYMMETRIC; i++)
2026     {
2027       if (sym_v[i] == SYSMIS)
2028         continue;
2029
2030       if (stats_categories[i] != last_cat)
2031         {
2032           last_cat = stats_categories[i];
2033           tab_text (sym, 0, 0, TAB_LEFT, gettext (categories[last_cat]));
2034         }
2035
2036       tab_text (sym, 1, 0, TAB_LEFT, gettext (stats[i]));
2037       tab_float (sym, 2, 0, TAB_RIGHT, sym_v[i], 8, 3);
2038       if (sym_ase[i] != SYSMIS)
2039         tab_float (sym, 3, 0, TAB_RIGHT, sym_ase[i], 8, 3);
2040       if (sym_t[i] != SYSMIS)
2041         tab_float (sym, 4, 0, TAB_RIGHT, sym_t[i], 8, 3);
2042       /*tab_float (sym, 5, 0, TAB_RIGHT, normal_sig (sym_v[i]), 8, 3);*/
2043       tab_next_row (sym);
2044     }
2045
2046   tab_text (sym, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2047   tab_float (sym, 2, 0, TAB_RIGHT, W, 8, 0);
2048   tab_next_row (sym);
2049
2050   tab_offset (sym, 0, -1);
2051 }
2052
2053 static int calc_risk (double[], double[], double[], union value *);
2054
2055 /* Display risk estimate. */
2056 static void
2057 display_risk (void)
2058 {
2059   char buf[256];
2060   double risk_v[3], lower[3], upper[3];
2061   union value c[2];
2062   int i;
2063
2064   if (!calc_risk (risk_v, upper, lower, c))
2065     return;
2066
2067   tab_offset (risk, nvar - 2, -1);
2068
2069   for (i = 0; i < 3; i++)
2070     {
2071       if (risk_v[i] == SYSMIS)
2072         continue;
2073
2074       switch (i)
2075         {
2076         case 0:
2077           if (x->vars[COL_VAR]->type == NUMERIC)
2078             sprintf (buf, _("Odds Ratio for %s (%g / %g)"),
2079                      x->vars[COL_VAR]->name, c[0].f, c[1].f);
2080           else
2081             sprintf (buf, _("Odds Ratio for %s (%.*s / %.*s)"),
2082                      x->vars[COL_VAR]->name,
2083                      x->vars[COL_VAR]->width, c[0].s,
2084                      x->vars[COL_VAR]->width, c[1].s);
2085           break;
2086         case 1:
2087         case 2:
2088           if (x->vars[ROW_VAR]->type == NUMERIC)
2089             sprintf (buf, _("For cohort %s = %g"),
2090                      x->vars[ROW_VAR]->name, rows[i - 1].f);
2091           else
2092             sprintf (buf, _("For cohort %s = %.*s"),
2093                      x->vars[ROW_VAR]->name,
2094                      x->vars[ROW_VAR]->width, rows[i - 1].s);
2095           break;
2096         }
2097
2098       tab_text (risk, 0, 0, TAB_LEFT, buf);
2099       tab_float (risk, 1, 0, TAB_RIGHT, risk_v[i], 8, 3);
2100       tab_float (risk, 2, 0, TAB_RIGHT, lower[i], 8, 3);
2101       tab_float (risk, 3, 0, TAB_RIGHT, upper[i], 8, 3);
2102       tab_next_row (risk);
2103     }
2104
2105   tab_text (risk, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2106   tab_float (risk, 1, 0, TAB_RIGHT, W, 8, 0);
2107   tab_next_row (risk);
2108
2109   tab_offset (risk, 0, -1);
2110 }
2111
2112 static int calc_directional (double[N_DIRECTIONAL], double[N_DIRECTIONAL],
2113                              double[N_DIRECTIONAL]);
2114
2115 /* Display directional measures. */
2116 static void
2117 display_directional (void)
2118 {
2119   static const char *categories[] =
2120     {
2121       N_("Nominal by Nominal"),
2122       N_("Ordinal by Ordinal"),
2123       N_("Nominal by Interval"),
2124     };
2125
2126   static const char *stats[] =
2127     {
2128       N_("Lambda"),
2129       N_("Goodman and Kruskal tau"),
2130       N_("Uncertainty Coefficient"),
2131       N_("Somers' d"),
2132       N_("Eta"),
2133     };
2134
2135   static const char *types[] =
2136     {
2137       N_("Symmetric"),
2138       N_("%s Dependent"),
2139       N_("%s Dependent"),
2140     };
2141
2142   static const int stats_categories[N_DIRECTIONAL] =
2143     {
2144       0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2,
2145     };
2146
2147   static const int stats_stats[N_DIRECTIONAL] =
2148     {
2149       0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4,
2150     };
2151
2152   static const int stats_types[N_DIRECTIONAL] =
2153     {
2154       0, 1, 2, 1, 2, 0, 1, 2, 0, 1, 2, 1, 2,
2155     };
2156
2157   static const int *stats_lookup[] =
2158     {
2159       stats_categories,
2160       stats_stats,
2161       stats_types,
2162     };
2163
2164   static const char **stats_names[] =
2165     {
2166       categories,
2167       stats,
2168       types,
2169     };
2170
2171   int last[3] =
2172     {
2173       -1, -1, -1,
2174     };
2175
2176   double direct_v[N_DIRECTIONAL];
2177   double direct_ase[N_DIRECTIONAL];
2178   double direct_t[N_DIRECTIONAL];
2179
2180   int i;
2181
2182   if (!calc_directional (direct_v, direct_ase, direct_t))
2183     return;
2184
2185   tab_offset (direct, nvar - 2, -1);
2186
2187   for (i = 0; i < N_DIRECTIONAL; i++)
2188     {
2189       if (direct_v[i] == SYSMIS)
2190         continue;
2191
2192       {
2193         int j;
2194
2195         for (j = 0; j < 3; j++)
2196           if (last[j] != stats_lookup[j][i])
2197             {
2198               if (j < 2)
2199                 tab_hline (direct, TAL_1, j, 6, 0);
2200
2201               for (; j < 3; j++)
2202                 {
2203                   char *string;
2204                   int k = last[j] = stats_lookup[j][i];
2205
2206                   if (k == 0)
2207                     string = NULL;
2208                   else if (k == 1)
2209                     string = x->vars[0]->name;
2210                   else
2211                     string = x->vars[1]->name;
2212
2213                   tab_text (direct, j, 0, TAB_LEFT | TAT_PRINTF,
2214                             gettext (stats_names[j][k]), string);
2215                 }
2216             }
2217       }
2218
2219       tab_float (direct, 3, 0, TAB_RIGHT, direct_v[i], 8, 3);
2220       if (direct_ase[i] != SYSMIS)
2221         tab_float (direct, 4, 0, TAB_RIGHT, direct_ase[i], 8, 3);
2222       if (direct_t[i] != SYSMIS)
2223         tab_float (direct, 5, 0, TAB_RIGHT, direct_t[i], 8, 3);
2224       /*tab_float (direct, 6, 0, TAB_RIGHT, normal_sig (direct_v[i]), 8, 3);*/
2225       tab_next_row (direct);
2226     }
2227
2228   tab_offset (direct, 0, -1);
2229 }
2230 \f
2231 /* Statistical calculations. */
2232
2233 /* Returns the value of the gamma (factorial) function for an integer
2234    argument X. */
2235 static double
2236 gamma_int (double x)
2237 {
2238   double r = 1;
2239   int i;
2240
2241   for (i = 2; i < x; i++)
2242     r *= i;
2243   return r;
2244 }
2245
2246 /* Calculate P_r as specified in _SPSS Statistical Algorithms_,
2247    Appendix 5. */
2248 static inline double
2249 Pr (int a, int b, int c, int d)
2250 {
2251   return (gamma_int (a + b + 1.) / gamma_int (a + 1.)
2252           * gamma_int (c + d + 1.) / gamma_int (b + 1.)
2253           * gamma_int (a + c + 1.) / gamma_int (c + 1.)
2254           * gamma_int (b + d + 1.) / gamma_int (d + 1.)
2255           / gamma_int (a + b + c + d + 1.));
2256 }
2257
2258 /* Swap the contents of A and B. */
2259 static inline void
2260 swap (int *a, int *b)
2261 {
2262   int t = *a;
2263   *a = *b;
2264   *b = t;
2265 }
2266
2267 /* Calculate significance for Fisher's exact test as specified in
2268    _SPSS Statistical Algorithms_, Appendix 5. */
2269 static void
2270 calc_fisher (int a, int b, int c, int d, double *fisher1, double *fisher2)
2271 {
2272   int x;
2273
2274   if (min (c, d) < min (a, b))
2275     swap (&a, &c), swap (&b, &d);
2276   if (min (b, d) < min (a, c))
2277     swap (&a, &b), swap (&c, &d);
2278   if (b * c < a * d)
2279     {
2280       if (b < c)
2281         swap (&a, &b), swap (&c, &d);
2282       else
2283         swap (&a, &c), swap (&b, &d);
2284     }
2285
2286   *fisher1 = 0.;
2287   for (x = 0; x <= a; x++)
2288     *fisher1 += Pr (a - x, b + x, c + x, d - x);
2289
2290   *fisher2 = *fisher1;
2291   for (x = 1; x <= b; x++)
2292     *fisher2 += Pr (a + x, b - x, c - x, d + x);
2293 }
2294
2295 /* Calculates chi-squares into CHISQ.  MAT is a matrix with N_COLS
2296    columns with values COLS and N_ROWS rows with values ROWS.  Values
2297    in the matrix sum to W. */
2298 static void
2299 calc_chisq (double chisq[N_CHISQ], int df[N_CHISQ],
2300             double *fisher1, double *fisher2)
2301 {
2302   int r, c;
2303
2304   chisq[0] = chisq[1] = 0.;
2305   chisq[2] = chisq[3] = chisq[4] = SYSMIS;
2306   *fisher1 = *fisher2 = SYSMIS;
2307
2308   df[0] = df[1] = (ns_cols - 1) * (ns_rows - 1);
2309
2310   if (ns_rows <= 1 || ns_cols <= 1)
2311     {
2312       chisq[0] = chisq[1] = SYSMIS;
2313       return;
2314     }
2315
2316   for (r = 0; r < n_rows; r++)
2317     for (c = 0; c < n_cols; c++)
2318       {
2319         const double expected = row_tot[r] * col_tot[c] / W;
2320         const double freq = mat[n_cols * r + c];
2321         const double residual = freq - expected;
2322
2323         chisq[0] += residual * residual / expected;
2324         if (freq)
2325           chisq[1] += freq * log (expected / freq);
2326       }
2327
2328   if (chisq[0] == 0.)
2329     chisq[0] = SYSMIS;
2330
2331   if (chisq[1] != 0.)
2332     chisq[1] *= -2.;
2333   else
2334     chisq[1] = SYSMIS;
2335
2336   /* Calculate Yates and Fisher exact test. */
2337   if (ns_cols == 2 && ns_rows == 2)
2338     {
2339       double f11, f12, f21, f22;
2340
2341       {
2342         int nz_cols[2];
2343         int i, j;
2344
2345         for (i = j = 0; i < n_cols; i++)
2346           if (col_tot[i] != 0.)
2347             {
2348               nz_cols[j++] = i;
2349               if (j == 2)
2350                 break;
2351             }
2352
2353         assert (j == 2);
2354
2355         f11 = mat[nz_cols[0]];
2356         f12 = mat[nz_cols[1]];
2357         f21 = mat[nz_cols[0] + n_cols];
2358         f22 = mat[nz_cols[1] + n_cols];
2359       }
2360
2361       /* Yates. */
2362       {
2363         const double x = fabs (f11 * f22 - f12 * f21) - 0.5 * W;
2364
2365         if (x > 0.)
2366           chisq[3] = (W * x * x
2367                       / (f11 + f12) / (f21 + f22)
2368                       / (f11 + f21) / (f12 + f22));
2369         else
2370           chisq[3] = 0.;
2371
2372         df[3] = 1.;
2373       }
2374
2375       /* Fisher. */
2376       if (f11 < 5. || f12 < 5. || f21 < 5. || f22 < 5.)
2377         calc_fisher (f11 + .5, f12 + .5, f21 + .5, f22 + .5, fisher1, fisher2);
2378     }
2379
2380   /* Calculate Mantel-Haenszel. */
2381   if (x->vars[ROW_VAR]->type == NUMERIC && x->vars[COL_VAR]->type == NUMERIC)
2382     {
2383       double r, ase_0, ase_1;
2384       calc_r ((double *) rows, (double *) cols, &r, &ase_0, &ase_1);
2385
2386       chisq[4] = (W - 1.) * r * r;
2387       df[4] = 1;
2388     }
2389 }
2390
2391 /* Calculate the value of Pearson's r.  r is stored into R, ase_1 into
2392    ASE_1, and ase_0 into ASE_0.  The row and column values must be
2393    passed in X and Y. */
2394 static void
2395 calc_r (double *X, double *Y, double *r, double *ase_0, double *ase_1)
2396 {
2397   double SX, SY, S, T;
2398   double Xbar, Ybar;
2399   double sum_XYf, sum_X2Y2f;
2400   double sum_Xr, sum_X2r;
2401   double sum_Yc, sum_Y2c;
2402   int i, j;
2403
2404   for (sum_X2Y2f = sum_XYf = 0., i = 0; i < n_rows; i++)
2405     for (j = 0; j < n_cols; j++)
2406       {
2407         double fij = mat[j + i * n_cols];
2408         double product = X[i] * Y[j];
2409         double temp = fij * product;
2410         sum_XYf += temp;
2411         sum_X2Y2f += temp * product;
2412       }
2413
2414   for (sum_Xr = sum_X2r = 0., i = 0; i < n_rows; i++)
2415     {
2416       sum_Xr += X[i] * row_tot[i];
2417       sum_X2r += X[i] * X[i] * row_tot[i];
2418     }
2419   Xbar = sum_Xr / W;
2420
2421   for (sum_Yc = sum_Y2c = 0., i = 0; i < n_cols; i++)
2422     {
2423       sum_Yc += Y[i] * col_tot[i];
2424       sum_Y2c += Y[i] * Y[i] * col_tot[i];
2425     }
2426   Ybar = sum_Yc / W;
2427
2428   S = sum_XYf - sum_Xr * sum_Yc / W;
2429   SX = sum_X2r - sum_Xr * sum_Xr / W;
2430   SY = sum_Y2c - sum_Yc * sum_Yc / W;
2431   T = sqrt (SX * SY);
2432   *r = S / T;
2433   *ase_0 = sqrt ((sum_X2Y2f - (sum_XYf * sum_XYf) / W) / (sum_X2r * sum_Y2c));
2434
2435   {
2436     double s, c, y, t;
2437
2438     for (s = c = 0., i = 0; i < n_rows; i++)
2439       for (j = 0; j < n_cols; j++)
2440         {
2441           double Xresid, Yresid;
2442           double temp;
2443
2444           Xresid = X[i] - Xbar;
2445           Yresid = Y[j] - Ybar;
2446           temp = (T * Xresid * Yresid
2447                   - ((S / (2. * T))
2448                      * (Xresid * Xresid * SY + Yresid * Yresid * SX)));
2449           y = mat[j + i * n_cols] * temp * temp - c;
2450           t = s + y;
2451           c = (t - s) - y;
2452           s = t;
2453         }
2454     *ase_1 = sqrt (s) / (T * T);
2455   }
2456 }
2457
2458 static double somers_d_v[3];
2459 static double somers_d_ase[3];
2460 static double somers_d_t[3];
2461
2462 /* Calculate symmetric statistics and their asymptotic standard
2463    errors.  Returns 0 if none could be calculated. */
2464 static int
2465 calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
2466                 double t[N_SYMMETRIC])
2467 {
2468   int q = min (ns_rows, ns_cols);
2469
2470   if (q <= 1)
2471     return 0;
2472
2473   {
2474     int i;
2475
2476     if (v)
2477       for (i = 0; i < N_SYMMETRIC; i++)
2478         v[i] = ase[i] = t[i] = SYSMIS;
2479   }
2480
2481   /* Phi, Cramer's V, contingency coefficient. */
2482   if (cmd.a_statistics[CRS_ST_PHI] || cmd.a_statistics[CRS_ST_CC])
2483     {
2484       double Xp = 0.;   /* Pearson chi-square. */
2485
2486       {
2487         int r, c;
2488
2489         for (r = 0; r < n_rows; r++)
2490           for (c = 0; c < n_cols; c++)
2491             {
2492               const double expected = row_tot[r] * col_tot[c] / W;
2493               const double freq = mat[n_cols * r + c];
2494               const double residual = freq - expected;
2495
2496               Xp += residual * residual / expected;
2497             }
2498       }
2499
2500       if (cmd.a_statistics[CRS_ST_PHI])
2501         {
2502           v[0] = sqrt (Xp / W);
2503           v[1] = sqrt (Xp / (W * (q - 1)));
2504         }
2505       if (cmd.a_statistics[CRS_ST_CC])
2506         v[2] = sqrt (Xp / (Xp + W));
2507     }
2508
2509   if (cmd.a_statistics[CRS_ST_BTAU] || cmd.a_statistics[CRS_ST_CTAU]
2510       || cmd.a_statistics[CRS_ST_GAMMA] || cmd.a_statistics[CRS_ST_D])
2511     {
2512       double *cum;
2513       double Dr, Dc;
2514       double P, Q;
2515       double btau_cum, ctau_cum, gamma_cum, d_yx_cum, d_xy_cum;
2516       double btau_var;
2517
2518       {
2519         int r, c;
2520
2521         Dr = Dc = W * W;
2522         for (r = 0; r < n_rows; r++)
2523           Dr -= row_tot[r] * row_tot[r];
2524         for (c = 0; c < n_cols; c++)
2525           Dc -= col_tot[c] * col_tot[c];
2526       }
2527
2528       {
2529         int r, c;
2530
2531         cum = xnmalloc (n_cols * n_rows, sizeof *cum);
2532         for (c = 0; c < n_cols; c++)
2533           {
2534             double ct = 0.;
2535
2536             for (r = 0; r < n_rows; r++)
2537               cum[c + r * n_cols] = ct += mat[c + r * n_cols];
2538           }
2539       }
2540
2541       /* P and Q. */
2542       {
2543         int i, j;
2544         double Cij, Dij;
2545
2546         P = Q = 0.;
2547         for (i = 0; i < n_rows; i++)
2548           {
2549             Cij = Dij = 0.;
2550
2551             for (j = 1; j < n_cols; j++)
2552               Cij += col_tot[j] - cum[j + i * n_cols];
2553
2554             if (i > 0)
2555               for (j = 1; j < n_cols; j++)
2556                 Dij += cum[j + (i - 1) * n_cols];
2557
2558             for (j = 0;;)
2559               {
2560                 double fij = mat[j + i * n_cols];
2561                 P += fij * Cij;
2562                 Q += fij * Dij;
2563
2564                 if (++j == n_cols)
2565                   break;
2566                 assert (j < n_cols);
2567
2568                 Cij -= col_tot[j] - cum[j + i * n_cols];
2569                 Dij += col_tot[j - 1] - cum[j - 1 + i * n_cols];
2570
2571                 if (i > 0)
2572                   {
2573                     Cij += cum[j - 1 + (i - 1) * n_cols];
2574                     Dij -= cum[j + (i - 1) * n_cols];
2575                   }
2576               }
2577           }
2578       }
2579
2580       if (cmd.a_statistics[CRS_ST_BTAU])
2581         v[3] = (P - Q) / sqrt (Dr * Dc);
2582       if (cmd.a_statistics[CRS_ST_CTAU])
2583         v[4] = (q * (P - Q)) / ((W * W) * (q - 1));
2584       if (cmd.a_statistics[CRS_ST_GAMMA])
2585         v[5] = (P - Q) / (P + Q);
2586
2587       /* ASE for tau-b, tau-c, gamma.  Calculations could be
2588          eliminated here, at expense of memory.  */
2589       {
2590         int i, j;
2591         double Cij, Dij;
2592
2593         btau_cum = ctau_cum = gamma_cum = d_yx_cum = d_xy_cum = 0.;
2594         for (i = 0; i < n_rows; i++)
2595           {
2596             Cij = Dij = 0.;
2597
2598             for (j = 1; j < n_cols; j++)
2599               Cij += col_tot[j] - cum[j + i * n_cols];
2600
2601             if (i > 0)
2602               for (j = 1; j < n_cols; j++)
2603                 Dij += cum[j + (i - 1) * n_cols];
2604
2605             for (j = 0;;)
2606               {
2607                 double fij = mat[j + i * n_cols];
2608
2609                 if (cmd.a_statistics[CRS_ST_BTAU])
2610                   {
2611                     const double temp = (2. * sqrt (Dr * Dc) * (Cij - Dij)
2612                                          + v[3] * (row_tot[i] * Dc
2613                                                    + col_tot[j] * Dr));
2614                     btau_cum += fij * temp * temp;
2615                   }
2616
2617                 {
2618                   const double temp = Cij - Dij;
2619                   ctau_cum += fij * temp * temp;
2620                 }
2621
2622                 if (cmd.a_statistics[CRS_ST_GAMMA])
2623                   {
2624                     const double temp = Q * Cij - P * Dij;
2625                     gamma_cum += fij * temp * temp;
2626                   }
2627
2628                 if (cmd.a_statistics[CRS_ST_D])
2629                   {
2630                     d_yx_cum += fij * pow2 (Dr * (Cij - Dij)
2631                                             - (P - Q) * (W - row_tot[i]));
2632                     d_xy_cum += fij * pow2 (Dc * (Dij - Cij)
2633                                             - (Q - P) * (W - col_tot[j]));
2634                   }
2635
2636                 if (++j == n_cols)
2637                   break;
2638                 assert (j < n_cols);
2639
2640                 Cij -= col_tot[j] - cum[j + i * n_cols];
2641                 Dij += col_tot[j - 1] - cum[j - 1 + i * n_cols];
2642
2643                 if (i > 0)
2644                   {
2645                     Cij += cum[j - 1 + (i - 1) * n_cols];
2646                     Dij -= cum[j + (i - 1) * n_cols];
2647                   }
2648               }
2649           }
2650       }
2651
2652       btau_var = ((btau_cum
2653                    - (W * pow2 (W * (P - Q) / sqrt (Dr * Dc) * (Dr + Dc))))
2654                   / pow2 (Dr * Dc));
2655       if (cmd.a_statistics[CRS_ST_BTAU])
2656         {
2657           ase[3] = sqrt (btau_var);
2658           t[3] = v[3] / (2 * sqrt ((ctau_cum - (P - Q) * (P - Q) / W)
2659                                    / (Dr * Dc)));
2660         }
2661       if (cmd.a_statistics[CRS_ST_CTAU])
2662         {
2663           ase[4] = ((2 * q / ((q - 1) * W * W))
2664                     * sqrt (ctau_cum - (P - Q) * (P - Q) / W));
2665           t[4] = v[4] / ase[4];
2666         }
2667       if (cmd.a_statistics[CRS_ST_GAMMA])
2668         {
2669           ase[5] = ((4. / ((P + Q) * (P + Q))) * sqrt (gamma_cum));
2670           t[5] = v[5] / (2. / (P + Q)
2671                          * sqrt (ctau_cum - (P - Q) * (P - Q) / W));
2672         }
2673       if (cmd.a_statistics[CRS_ST_D])
2674         {
2675           somers_d_v[0] = (P - Q) / (.5 * (Dc + Dr));
2676           somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
2677           somers_d_t[0] = (somers_d_v[0]
2678                            / (4 / (Dc + Dr)
2679                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2680           somers_d_v[1] = (P - Q) / Dc;
2681           somers_d_ase[1] = 2. / pow2 (Dc) * sqrt (d_xy_cum);
2682           somers_d_t[1] = (somers_d_v[1]
2683                            / (2. / Dc
2684                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2685           somers_d_v[2] = (P - Q) / Dr;
2686           somers_d_ase[2] = 2. / pow2 (Dr) * sqrt (d_yx_cum);
2687           somers_d_t[2] = (somers_d_v[2]
2688                            / (2. / Dr
2689                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2690         }
2691
2692       free (cum);
2693     }
2694
2695   /* Spearman correlation, Pearson's r. */
2696   if (cmd.a_statistics[CRS_ST_CORR])
2697     {
2698       double *R = local_alloc (sizeof *R * n_rows);
2699       double *C = local_alloc (sizeof *C * n_cols);
2700
2701       {
2702         double y, t, c = 0., s = 0.;
2703         int i = 0;
2704
2705         for (;;)
2706           {
2707             R[i] = s + (row_tot[i] + 1.) / 2.;
2708             y = row_tot[i] - c;
2709             t = s + y;
2710             c = (t - s) - y;
2711             s = t;
2712             if (++i == n_rows)
2713               break;
2714             assert (i < n_rows);
2715           }
2716       }
2717
2718       {
2719         double y, t, c = 0., s = 0.;
2720         int j = 0;
2721
2722         for (;;)
2723           {
2724             C[j] = s + (col_tot[j] + 1.) / 2;
2725             y = col_tot[j] - c;
2726             t = s + y;
2727             c = (t - s) - y;
2728             s = t;
2729             if (++j == n_cols)
2730               break;
2731             assert (j < n_cols);
2732           }
2733       }
2734
2735       calc_r (R, C, &v[6], &t[6], &ase[6]);
2736       t[6] = v[6] / t[6];
2737
2738       local_free (R);
2739       local_free (C);
2740
2741       calc_r ((double *) rows, (double *) cols, &v[7], &t[7], &ase[7]);
2742       t[7] = v[7] / t[7];
2743     }
2744
2745   /* Cohen's kappa. */
2746   if (cmd.a_statistics[CRS_ST_KAPPA] && ns_rows == ns_cols)
2747     {
2748       double sum_fii, sum_rici, sum_fiiri_ci, sum_fijri_ci2, sum_riciri_ci;
2749       int i, j;
2750
2751       for (sum_fii = sum_rici = sum_fiiri_ci = sum_riciri_ci = 0., i = j = 0;
2752            i < ns_rows; i++, j++)
2753         {
2754           double prod, sum;
2755
2756           while (col_tot[j] == 0.)
2757             j++;
2758
2759           prod = row_tot[i] * col_tot[j];
2760           sum = row_tot[i] + col_tot[j];
2761
2762           sum_fii += mat[j + i * n_cols];
2763           sum_rici += prod;
2764           sum_fiiri_ci += mat[j + i * n_cols] * sum;
2765           sum_riciri_ci += prod * sum;
2766         }
2767       for (sum_fijri_ci2 = 0., i = 0; i < ns_rows; i++)
2768         for (j = 0; j < ns_cols; j++)
2769           {
2770             double sum = row_tot[i] + col_tot[j];
2771             sum_fijri_ci2 += mat[j + i * n_cols] * sum * sum;
2772           }
2773
2774       v[8] = (W * sum_fii - sum_rici) / (W * W - sum_rici);
2775
2776       ase[8] = sqrt ((W * W * sum_rici
2777                       + sum_rici * sum_rici
2778                       - W * sum_riciri_ci)
2779                      / (W * (W * W - sum_rici) * (W * W - sum_rici)));
2780 #if 0
2781       t[8] = v[8] / sqrt (W * (((sum_fii * (W - sum_fii))
2782                                 / pow2 (W * W - sum_rici))
2783                                + ((2. * (W - sum_fii)
2784                                    * (2. * sum_fii * sum_rici
2785                                       - W * sum_fiiri_ci))
2786                                   / cube (W * W - sum_rici))
2787                                + (pow2 (W - sum_fii)
2788                                   * (W * sum_fijri_ci2 - 4.
2789                                      * sum_rici * sum_rici)
2790                                   / pow4 (W * W - sum_rici))));
2791 #else
2792       t[8] = v[8] / ase[8];
2793 #endif
2794     }
2795
2796   return 1;
2797 }
2798
2799 /* Calculate risk estimate. */
2800 static int
2801 calc_risk (double *value, double *upper, double *lower, union value *c)
2802 {
2803   double f11, f12, f21, f22;
2804   double v;
2805
2806   {
2807     int i;
2808
2809     for (i = 0; i < 3; i++)
2810       value[i] = upper[i] = lower[i] = SYSMIS;
2811   }
2812
2813   if (ns_rows != 2 || ns_cols != 2)
2814     return 0;
2815
2816   {
2817     int nz_cols[2];
2818     int i, j;
2819
2820     for (i = j = 0; i < n_cols; i++)
2821       if (col_tot[i] != 0.)
2822         {
2823           nz_cols[j++] = i;
2824           if (j == 2)
2825             break;
2826         }
2827
2828     assert (j == 2);
2829
2830     f11 = mat[nz_cols[0]];
2831     f12 = mat[nz_cols[1]];
2832     f21 = mat[nz_cols[0] + n_cols];
2833     f22 = mat[nz_cols[1] + n_cols];
2834
2835     c[0] = cols[nz_cols[0]];
2836     c[1] = cols[nz_cols[1]];
2837   }
2838
2839   value[0] = (f11 * f22) / (f12 * f21);
2840   v = sqrt (1. / f11 + 1. / f12 + 1. / f21 + 1. / f22);
2841   lower[0] = value[0] * exp (-1.960 * v);
2842   upper[0] = value[0] * exp (1.960 * v);
2843
2844   value[1] = (f11 * (f21 + f22)) / (f21 * (f11 + f12));
2845   v = sqrt ((f12 / (f11 * (f11 + f12)))
2846             + (f22 / (f21 * (f21 + f22))));
2847   lower[1] = value[1] * exp (-1.960 * v);
2848   upper[1] = value[1] * exp (1.960 * v);
2849
2850   value[2] = (f12 * (f21 + f22)) / (f22 * (f11 + f12));
2851   v = sqrt ((f11 / (f12 * (f11 + f12)))
2852             + (f21 / (f22 * (f21 + f22))));
2853   lower[2] = value[2] * exp (-1.960 * v);
2854   upper[2] = value[2] * exp (1.960 * v);
2855
2856   return 1;
2857 }
2858
2859 /* Calculate directional measures. */
2860 static int
2861 calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
2862                   double t[N_DIRECTIONAL])
2863 {
2864   {
2865     int i;
2866
2867     for (i = 0; i < N_DIRECTIONAL; i++)
2868       v[i] = ase[i] = t[i] = SYSMIS;
2869   }
2870
2871   /* Lambda. */
2872   if (cmd.a_statistics[CRS_ST_LAMBDA])
2873     {
2874       double *fim = xnmalloc (n_rows, sizeof *fim);
2875       int *fim_index = xnmalloc (n_rows, sizeof *fim_index);
2876       double *fmj = xnmalloc (n_cols, sizeof *fmj);
2877       int *fmj_index = xnmalloc (n_cols, sizeof *fmj_index);
2878       double sum_fim, sum_fmj;
2879       double rm, cm;
2880       int rm_index, cm_index;
2881       int i, j;
2882
2883       /* Find maximum for each row and their sum. */
2884       for (sum_fim = 0., i = 0; i < n_rows; i++)
2885         {
2886           double max = mat[i * n_cols];
2887           int index = 0;
2888
2889           for (j = 1; j < n_cols; j++)
2890             if (mat[j + i * n_cols] > max)
2891               {
2892                 max = mat[j + i * n_cols];
2893                 index = j;
2894               }
2895
2896           sum_fim += fim[i] = max;
2897           fim_index[i] = index;
2898         }
2899
2900       /* Find maximum for each column. */
2901       for (sum_fmj = 0., j = 0; j < n_cols; j++)
2902         {
2903           double max = mat[j];
2904           int index = 0;
2905
2906           for (i = 1; i < n_rows; i++)
2907             if (mat[j + i * n_cols] > max)
2908               {
2909                 max = mat[j + i * n_cols];
2910                 index = i;
2911               }
2912
2913           sum_fmj += fmj[j] = max;
2914           fmj_index[j] = index;
2915         }
2916
2917       /* Find maximum row total. */
2918       rm = row_tot[0];
2919       rm_index = 0;
2920       for (i = 1; i < n_rows; i++)
2921         if (row_tot[i] > rm)
2922           {
2923             rm = row_tot[i];
2924             rm_index = i;
2925           }
2926
2927       /* Find maximum column total. */
2928       cm = col_tot[0];
2929       cm_index = 0;
2930       for (j = 1; j < n_cols; j++)
2931         if (col_tot[j] > cm)
2932           {
2933             cm = col_tot[j];
2934             cm_index = j;
2935           }
2936
2937       v[0] = (sum_fim + sum_fmj - cm - rm) / (2. * W - rm - cm);
2938       v[1] = (sum_fmj - rm) / (W - rm);
2939       v[2] = (sum_fim - cm) / (W - cm);
2940
2941       /* ASE1 for Y given X. */
2942       {
2943         double accum;
2944
2945         for (accum = 0., i = 0; i < n_rows; i++)
2946           for (j = 0; j < n_cols; j++)
2947             {
2948               const int deltaj = j == cm_index;
2949               accum += (mat[j + i * n_cols]
2950                         * pow2 ((j == fim_index[i])
2951                                - deltaj
2952                                + v[0] * deltaj));
2953             }
2954
2955         ase[2] = sqrt (accum - W * v[0]) / (W - cm);
2956       }
2957
2958       /* ASE0 for Y given X. */
2959       {
2960         double accum;
2961
2962         for (accum = 0., i = 0; i < n_rows; i++)
2963           if (cm_index != fim_index[i])
2964             accum += (mat[i * n_cols + fim_index[i]]
2965                       + mat[i * n_cols + cm_index]);
2966         t[2] = v[2] / (sqrt (accum - pow2 (sum_fim - cm) / W) / (W - cm));
2967       }
2968
2969       /* ASE1 for X given Y. */
2970       {
2971         double accum;
2972
2973         for (accum = 0., i = 0; i < n_rows; i++)
2974           for (j = 0; j < n_cols; j++)
2975             {
2976               const int deltaj = i == rm_index;
2977               accum += (mat[j + i * n_cols]
2978                         * pow2 ((i == fmj_index[j])
2979                                - deltaj
2980                                + v[0] * deltaj));
2981             }
2982
2983         ase[1] = sqrt (accum - W * v[0]) / (W - rm);
2984       }
2985
2986       /* ASE0 for X given Y. */
2987       {
2988         double accum;
2989
2990         for (accum = 0., j = 0; j < n_cols; j++)
2991           if (rm_index != fmj_index[j])
2992             accum += (mat[j + n_cols * fmj_index[j]]
2993                       + mat[j + n_cols * rm_index]);
2994         t[1] = v[1] / (sqrt (accum - pow2 (sum_fmj - rm) / W) / (W - rm));
2995       }
2996
2997       /* Symmetric ASE0 and ASE1. */
2998       {
2999         double accum0;
3000         double accum1;
3001
3002         for (accum0 = accum1 = 0., i = 0; i < n_rows; i++)
3003           for (j = 0; j < n_cols; j++)
3004             {
3005               int temp0 = (fmj_index[j] == i) + (fim_index[i] == j);
3006               int temp1 = (i == rm_index) + (j == cm_index);
3007               accum0 += mat[j + i * n_cols] * pow2 (temp0 - temp1);
3008               accum1 += (mat[j + i * n_cols]
3009                          * pow2 (temp0 + (v[0] - 1.) * temp1));
3010             }
3011         ase[0] = sqrt (accum1 - 4. * W * v[0] * v[0]) / (2. * W - rm - cm);
3012         t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / W))
3013                        / (2. * W - rm - cm));
3014       }
3015
3016       free (fim);
3017       free (fim_index);
3018       free (fmj);
3019       free (fmj_index);
3020
3021       {
3022         double sum_fij2_ri, sum_fij2_ci;
3023         double sum_ri2, sum_cj2;
3024
3025         for (sum_fij2_ri = sum_fij2_ci = 0., i = 0; i < n_rows; i++)
3026           for (j = 0; j < n_cols; j++)
3027             {
3028               double temp = pow2 (mat[j + i * n_cols]);
3029               sum_fij2_ri += temp / row_tot[i];
3030               sum_fij2_ci += temp / col_tot[j];
3031             }
3032
3033         for (sum_ri2 = 0., i = 0; i < n_rows; i++)
3034           sum_ri2 += row_tot[i] * row_tot[i];
3035
3036         for (sum_cj2 = 0., j = 0; j < n_cols; j++)
3037           sum_cj2 += col_tot[j] * col_tot[j];
3038
3039         v[3] = (W * sum_fij2_ci - sum_ri2) / (W * W - sum_ri2);
3040         v[4] = (W * sum_fij2_ri - sum_cj2) / (W * W - sum_cj2);
3041       }
3042     }
3043
3044   if (cmd.a_statistics[CRS_ST_UC])
3045     {
3046       double UX, UY, UXY, P;
3047       double ase1_yx, ase1_xy, ase1_sym;
3048       int i, j;
3049
3050       for (UX = 0., i = 0; i < n_rows; i++)
3051         if (row_tot[i] > 0.)
3052           UX -= row_tot[i] / W * log (row_tot[i] / W);
3053
3054       for (UY = 0., j = 0; j < n_cols; j++)
3055         if (col_tot[j] > 0.)
3056           UY -= col_tot[j] / W * log (col_tot[j] / W);
3057
3058       for (UXY = P = 0., i = 0; i < n_rows; i++)
3059         for (j = 0; j < n_cols; j++)
3060           {
3061             double entry = mat[j + i * n_cols];
3062
3063             if (entry <= 0.)
3064               continue;
3065
3066             P += entry * pow2 (log (col_tot[j] * row_tot[i] / (W * entry)));
3067             UXY -= entry / W * log (entry / W);
3068           }
3069
3070       for (ase1_yx = ase1_xy = ase1_sym = 0., i = 0; i < n_rows; i++)
3071         for (j = 0; j < n_cols; j++)
3072           {
3073             double entry = mat[j + i * n_cols];
3074
3075             if (entry <= 0.)
3076               continue;
3077
3078             ase1_yx += entry * pow2 (UY * log (entry / row_tot[i])
3079                                     + (UX - UXY) * log (col_tot[j] / W));
3080             ase1_xy += entry * pow2 (UX * log (entry / col_tot[j])
3081                                     + (UY - UXY) * log (row_tot[i] / W));
3082             ase1_sym += entry * pow2 ((UXY
3083                                       * log (row_tot[i] * col_tot[j] / (W * W)))
3084                                      - (UX + UY) * log (entry / W));
3085           }
3086
3087       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
3088       ase[5] = (2. / (W * pow2 (UX + UY))) * sqrt (ase1_sym);
3089       t[5] = v[5] / ((2. / (W * (UX + UY)))
3090                      * sqrt (P - pow2 (UX + UY - UXY) / W));
3091
3092       v[6] = (UX + UY - UXY) / UX;
3093       ase[6] = sqrt (ase1_xy) / (W * UX * UX);
3094       t[6] = v[6] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UX));
3095
3096       v[7] = (UX + UY - UXY) / UY;
3097       ase[7] = sqrt (ase1_yx) / (W * UY * UY);
3098       t[7] = v[7] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UY));
3099     }
3100
3101   /* Somers' D. */
3102   if (cmd.a_statistics[CRS_ST_D])
3103     {
3104       int i;
3105
3106       if (!sym)
3107         calc_symmetric (NULL, NULL, NULL);
3108       for (i = 0; i < 3; i++)
3109         {
3110           v[8 + i] = somers_d_v[i];
3111           ase[8 + i] = somers_d_ase[i];
3112           t[8 + i] = somers_d_t[i];
3113         }
3114     }
3115
3116   /* Eta. */
3117   if (cmd.a_statistics[CRS_ST_ETA])
3118     {
3119       {
3120         double sum_Xr, sum_X2r;
3121         double SX, SXW;
3122         int i, j;
3123
3124         for (sum_Xr = sum_X2r = 0., i = 0; i < n_rows; i++)
3125           {
3126             sum_Xr += rows[i].f * row_tot[i];
3127             sum_X2r += rows[i].f * rows[i].f * row_tot[i];
3128           }
3129         SX = sum_X2r - sum_Xr * sum_Xr / W;
3130
3131         for (SXW = 0., j = 0; j < n_cols; j++)
3132           {
3133             double cum;
3134
3135             for (cum = 0., i = 0; i < n_rows; i++)
3136               {
3137                 SXW += rows[i].f * rows[i].f * mat[j + i * n_cols];
3138                 cum += rows[i].f * mat[j + i * n_cols];
3139               }
3140
3141             SXW -= cum * cum / col_tot[j];
3142           }
3143         v[11] = sqrt (1. - SXW / SX);
3144       }
3145
3146       {
3147         double sum_Yc, sum_Y2c;
3148         double SY, SYW;
3149         int i, j;
3150
3151         for (sum_Yc = sum_Y2c = 0., i = 0; i < n_cols; i++)
3152           {
3153             sum_Yc += cols[i].f * col_tot[i];
3154             sum_Y2c += cols[i].f * cols[i].f * col_tot[i];
3155           }
3156         SY = sum_Y2c - sum_Yc * sum_Yc / W;
3157
3158         for (SYW = 0., i = 0; i < n_rows; i++)
3159           {
3160             double cum;
3161
3162             for (cum = 0., j = 0; j < n_cols; j++)
3163               {
3164                 SYW += cols[j].f * cols[j].f * mat[j + i * n_cols];
3165                 cum += cols[j].f * mat[j + i * n_cols];
3166               }
3167
3168             SYW -= cum * cum / row_tot[i];
3169           }
3170         v[12] = sqrt (1. - SYW / SY);
3171       }
3172     }
3173
3174   return 1;
3175 }
3176
3177 /* A wrapper around data_out() that limits string output to short
3178    string width and null terminates the result. */
3179 static void
3180 format_short (char *s, const struct fmt_spec *fp, const union value *v)
3181 {
3182   struct fmt_spec fmt_subst;
3183
3184   /* Limit to short string width. */
3185   if (formats[fp->type].cat & FCAT_STRING)
3186     {
3187       fmt_subst = *fp;
3188
3189       assert (fmt_subst.type == FMT_A || fmt_subst.type == FMT_AHEX);
3190       if (fmt_subst.type == FMT_A)
3191         fmt_subst.w = min (8, fmt_subst.w);
3192       else
3193         fmt_subst.w = min (16, fmt_subst.w);
3194
3195       fp = &fmt_subst;
3196     }
3197
3198   /* Format. */
3199   data_out (s, fp, v);
3200
3201   /* Null terminate. */
3202   s[fp->w] = '\0';
3203 }
3204
3205 /*
3206    Local Variables:
3207    mode: c
3208    End:
3209 */