1 /* PSPP - computes sample statistics.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3 Written by Ben Pfaff <blp@gnu.org>.
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.
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.
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
26 /* Kahan summation formula, Thm. 8, _What Every Computer Scientist
27 Should Know About Floating-Point Arithmetic_, David Goldberg,
28 orig. March 1991 issue of Computing Surveys, also at
29 <URL:http://www.wam.umd.edu/whats_new/workshop3.0/common-tools/numerical_comp_guide/goldberg1.doc.html>.
30 Hopefully your compiler won't try to optimize the code below too
31 much, because that will ruin the precision. */
32 #define KAHAN_SUMMATION_FORMULA(S) \
38 S = SUMMATION_ELEMENT (0); \
40 for (S_j = 1; S_j < SUMMATION_COUNT; S_j++) \
42 double S_y = SUMMATION_ELEMENT (S_j) - S_c; \
43 double S_t = S + S_y; \
44 S_c = (S_t - S) - S_y; \
53 /* Allocate a new vector of length N. */
57 struct vector *vec = xmalloc (sizeof *vec);
58 vec->data = xmalloc (sizeof *vec->data * n);
63 /* Change the length of VEC to N. The amount of space allocated will
64 not be lowered, but may be enlarged. */
66 vec_realloc (struct vector *vec, int n)
71 vec->data = xrealloc (vec->data, sizeof *vec->data * n);
76 /* Free vector VEC. */
78 vec_free (struct vector *vec)
84 /* Set the values in vector VEC to constant VALUE. */
87 vec_init (struct vector *vec, double value)
93 for (i = 0; i < vec->n; i++)
98 /* Print out vector VEC to stdout for debugging purposes. */
101 #include "settings.h"
104 vec_print (const struct vector *vec)
108 for (i = 0; i < vec->n; i++)
110 if (i % ((get_viewwidth() - 4) / 8) == 0)
117 printf ("%8g", vec_elem (vec, i));
122 /* Return the sum of the values in VEC. */
124 vec_total (const struct vector *vec)
128 #define SUMMATION_COUNT (vec->n)
129 #define SUMMATION_ELEMENT(INDEX) (vec_elem (vec, (INDEX)))
130 KAHAN_SUMMATION_FORMULA (sum);
131 #undef SUMMATION_COUNT
132 #undef SUMMATION_ELEMENT
139 /* Allocate a new matrix with NR rows and NC columns. */
141 mat_alloc (int nr, int nc)
143 struct matrix *mat = xmalloc (sizeof *mat);
147 mat->data = xmalloc (sizeof *mat->data * nr * nc);
151 /* Set the size of matrix MAT to NR rows and NC columns. The matrix
152 data array will be enlarged if necessary but will not be shrunk. */
154 mat_realloc (struct matrix *mat, int nr, int nc)
156 if (nc * nr > mat->m)
159 mat->data = xrealloc (mat->data, sizeof *mat->data * mat->m);
165 /* Free matrix MAT. */
167 mat_free (struct matrix *mat)
173 /* Set all matrix MAT entries to VALUE. */
175 mat_init (struct matrix *mat, double value)
181 for (i = 0; i < mat->nr * mat->nc; i++)
185 /* Set all MAT entries in row R to VALUE. */
187 mat_init_row (struct matrix *mat, int r, double value)
192 p = &mat_elem (mat, r, 0);
193 for (i = 0; i < mat->nc; i++)
197 /* Set all MAT entries in column C to VALUE. */
199 mat_init_col (struct matrix *mat, int c, double value)
204 p = &mat_elem (mat, 0, c);
205 for (i = 0; i < mat->nr; i++)
212 /* Print out MAT entries to stdout, optionally with row and column
213 labels ROW_LABELS and COL_LABELS. */
216 mat_print (const struct matrix *mat,
217 const struct vector *row_labels,
218 const struct vector *col_labels)
222 assert (!row_labels || row_labels->n == mat->nr);
227 assert (col_labels->n == mat->nc);
230 for (c = 0; c < mat->nc; c++)
231 printf ("%8g", vec_elem (col_labels, c));
234 for (r = 0; r < mat->nr; r++)
237 printf ("%8g:", vec_elem (row_labels, r));
238 for (c = 0; c < mat->nc; c++)
239 printf ("%8g", mat_elem (mat, r, c));
243 #endif /* GLOBAL_DEBUGGING */
245 /* Calculate row totals for matrix MAT into vector ROW_TOTS. */
247 mat_row_totals (const struct matrix *mat, struct vector *row_tots)
251 vec_realloc (row_tots, mat->nr);
252 for (r = 0; r < mat->nr; r++)
256 #define SUMMATION_COUNT (mat->nc)
257 #define SUMMATION_ELEMENT(INDEX) (mat_elem (mat, r, INDEX))
258 KAHAN_SUMMATION_FORMULA (sum);
259 #undef SUMMATION_COUNT
260 #undef SUMMATION_ELEMENT
262 vec_elem (row_tots, r) = sum;
266 /* Calculate column totals for matrix MAT into vector COL_TOTS. */
268 mat_col_totals (const struct matrix *mat, struct vector *col_tots)
272 vec_realloc (col_tots, mat->nc);
273 for (c = 0; c < mat->nc; c++)
277 #define SUMMATION_COUNT (mat->nr)
278 #define SUMMATION_ELEMENT(INDEX) (mat_elem (mat, INDEX, c))
279 KAHAN_SUMMATION_FORMULA (sum);
280 #undef SUMMATION_COUNT
281 #undef SUMMATION_ELEMENT
283 vec_elem (col_tots, c) = sum;
287 /* Return the grand total for matrix MAT. Of course, if you're also
288 calculating column or row totals, it would be faster to use
289 vec_total on one of those sets of totals. */
291 mat_grand_total (const struct matrix *mat)
295 #define SUMMATION_COUNT (mat->nr * mat->nc)
296 #define SUMMATION_ELEMENT(INDEX) (mat->data[INDEX])
297 KAHAN_SUMMATION_FORMULA (sum);
298 #undef SUMMATION_COUNT
299 #undef SUMMATION_ELEMENT