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