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