sstest: add summary to random numbers test section
[gnumeric.git] / src / sstest.c
blob38f9162bf047f10c09da9e015db8d13a820fa7f7
1 /*
2 * sstest.c: Test code for Gnumeric
4 * Copyright (C) 2009 Morten Welinder (terra@gnome.org)
5 */
6 #include <gnumeric-config.h>
7 #include "gnumeric.h"
8 #include "libgnumeric.h"
9 #include <goffice/goffice.h>
10 #include "command-context-stderr.h"
11 #include "workbook-view.h"
12 #include "workbook.h"
13 #include "gutils.h"
14 #include "gnm-plugin.h"
15 #include "parse-util.h"
16 #include "expr-name.h"
17 #include "expr.h"
18 #include "search.h"
19 #include "sheet.h"
20 #include "cell.h"
21 #include "value.h"
22 #include "func.h"
23 #include "parse-util.h"
24 #include "sheet-object-cell-comment.h"
25 #include "mathfunc.h"
26 #include "gnm-random.h"
27 #include "sf-dpq.h"
28 #include "sf-gamma.h"
29 #include "rangefunc.h"
31 #include <gsf/gsf-input-stdio.h>
32 #include <gsf/gsf-input-textline.h>
33 #include <glib/gi18n.h>
34 #include <string.h>
35 #include <errno.h>
37 static gboolean sstest_show_version = FALSE;
38 static gboolean sstest_fast = FALSE;
40 static GOptionEntry const sstest_options [] = {
42 "fast", 'f',
43 0, G_OPTION_ARG_NONE, &sstest_fast,
44 N_("Run fewer iterations"),
45 NULL
49 "version", 'V',
50 0, G_OPTION_ARG_NONE, &sstest_show_version,
51 N_("Display program version"),
52 NULL
55 { NULL }
58 static void
59 mark_test_start (const char *name)
61 g_printerr ("-----------------------------------------------------------------------------\nStart: %s\n-----------------------------------------------------------------------------\n\n", name);
64 static void
65 mark_test_end (const char *name)
67 g_printerr ("End: %s\n\n", name);
70 static void
71 cb_collect_names (G_GNUC_UNUSED const char *name, GnmNamedExpr *nexpr, GSList **names)
73 *names = g_slist_prepend (*names, nexpr);
76 static GnmCell *
77 fetch_cell (Sheet *sheet, const char *where)
79 GnmCellPos cp;
80 gboolean ok = cellpos_parse (where,
81 gnm_sheet_get_size (sheet),
82 &cp, TRUE) != NULL;
83 g_return_val_if_fail (ok, NULL);
84 return sheet_cell_fetch (sheet, cp.col, cp.row);
87 static void
88 set_cell (Sheet *sheet, const char *where, const char *what)
90 GnmCell *cell = fetch_cell (sheet, where);
91 if (cell)
92 gnm_cell_set_text (cell, what);
95 static void
96 dump_sheet (Sheet *sheet, const char *header)
98 GPtrArray *cells = sheet_cells (sheet, NULL);
99 unsigned ui;
101 if (header)
102 g_printerr ("# %s\n", header);
103 for (ui = 0; ui < cells->len; ui++) {
104 GnmCell *cell = g_ptr_array_index (cells, ui);
105 char *txt = gnm_cell_get_entered_text (cell);
106 g_printerr ("%s: %s\n",
107 cellpos_as_string (&cell->pos), txt);
108 g_free (txt);
110 g_ptr_array_free (cells, TRUE);
114 static void
115 dump_names (Workbook *wb)
117 GSList *l, *names = NULL;
119 workbook_foreach_name (wb, FALSE, (GHFunc)cb_collect_names, &names);
120 names = g_slist_sort (names, (GCompareFunc)expr_name_cmp_by_name);
122 g_printerr ("Dumping names...\n");
123 for (l = names; l; l = l->next) {
124 GnmNamedExpr *nexpr = l->data;
125 GnmConventionsOut out;
127 out.accum = g_string_new (NULL);
128 out.pp = &nexpr->pos;
129 out.convs = gnm_conventions_default;
131 g_string_append (out.accum, "Scope=");
132 if (out.pp->sheet)
133 g_string_append (out.accum, out.pp->sheet->name_quoted);
134 else
135 g_string_append (out.accum, "Global");
137 g_string_append (out.accum, " Name=");
138 go_strescape (out.accum, expr_name_name (nexpr));
140 g_string_append (out.accum, " Expr=");
141 gnm_expr_top_as_gstring (nexpr->texpr, &out);
143 g_printerr ("%s\n", out.accum->str);
144 g_string_free (out.accum, TRUE);
146 g_printerr ("Dumping names... Done\n");
148 g_slist_free (names);
151 static void
152 define_name (const char *name, const char *expr_txt, gpointer scope)
154 GnmParsePos pos;
155 GnmExprTop const *texpr;
156 GnmNamedExpr const *nexpr;
157 GnmConventions const *convs;
159 if (IS_SHEET (scope)) {
160 parse_pos_init_sheet (&pos, scope);
161 convs = sheet_get_conventions (pos.sheet);
162 } else {
163 parse_pos_init (&pos, WORKBOOK (scope), NULL, 0, 0);
164 convs = gnm_conventions_default;
167 texpr = gnm_expr_parse_str (expr_txt, &pos,
168 GNM_EXPR_PARSE_DEFAULT,
169 convs, NULL);
170 if (!texpr) {
171 g_printerr ("Failed to parse %s for name %s\n",
172 expr_txt, name);
173 return;
176 nexpr = expr_name_add (&pos, name, texpr, NULL, TRUE, NULL);
177 if (!nexpr)
178 g_printerr ("Failed to add name %s\n", name);
181 static void
182 test_insdel_rowcol_names (void)
184 Workbook *wb;
185 Sheet *sheet1,*sheet2;
186 const char *test_name = "test_insdel_rowcol_names";
187 GOUndo *undo;
188 int i;
190 mark_test_start (test_name);
192 wb = workbook_new ();
193 sheet1 = workbook_sheet_add (wb, -1,
194 GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
195 sheet2 = workbook_sheet_add (wb, -1,
196 GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
198 define_name ("Print_Area", "Sheet1!$A$1:$IV$65536", sheet1);
199 define_name ("Print_Area", "Sheet2!$A$1:$IV$65536", sheet2);
201 define_name ("NAMEGA1", "A1", wb);
202 define_name ("NAMEG2", "$A$14+Sheet1!$A$14+Sheet2!$A$14", wb);
204 define_name ("NAMEA1", "A1", sheet1);
205 define_name ("NAMEA2", "A2", sheet1);
206 define_name ("NAMEA1ABS", "$A$1", sheet1);
207 define_name ("NAMEA2ABS", "$A$2", sheet1);
209 dump_names (wb);
211 for (i = 3; i >= 0; i--) {
212 g_printerr ("About to insert before column %s on %s\n",
213 col_name (i), sheet1->name_unquoted);
214 sheet_insert_cols (sheet1, i, 12, &undo, NULL);
215 dump_names (wb);
216 g_printerr ("Undoing.\n");
217 go_undo_undo_with_data (undo, NULL);
218 g_object_unref (undo);
219 g_printerr ("Done.\n");
222 for (i = 3; i >= 0; i--) {
223 g_printerr ("About to insert before column %s on %s\n",
224 col_name (i), sheet2->name_unquoted);
225 sheet_insert_cols (sheet2, i, 12, &undo, NULL);
226 dump_names (wb);
227 g_printerr ("Undoing.\n");
228 go_undo_undo_with_data (undo, NULL);
229 g_object_unref (undo);
230 g_printerr ("Done.\n");
233 for (i = 3; i >= 0; i--) {
234 g_printerr ("About to delete column %s on %s\n",
235 col_name (i), sheet1->name_unquoted);
236 sheet_delete_cols (sheet1, i, 1, &undo, NULL);
237 dump_names (wb);
238 g_printerr ("Undoing.\n");
239 go_undo_undo_with_data (undo, NULL);
240 g_object_unref (undo);
241 g_printerr ("Done.\n");
244 g_object_unref (wb);
246 mark_test_end (test_name);
249 /*-------------------------------------------------------------------------- */
251 static void
252 test_insert_delete (void)
254 const char *test_name = "test_insert_delete";
255 Workbook *wb;
256 Sheet *sheet1;
257 int i;
258 GOUndo *u = NULL, *u1;
260 mark_test_start (test_name);
262 wb = workbook_new ();
263 sheet1 = workbook_sheet_add (wb, -1,
264 GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
265 set_cell (sheet1, "B2", "=D4+1");
266 set_cell (sheet1, "D2", "=if(TRUE,B2,2)");
268 dump_sheet (sheet1, "Init");
270 for (i = 5; i >= 0; i--) {
271 g_printerr ("# About to insert column before %s\n",
272 col_name (i));
273 sheet_insert_cols (sheet1, i, 1, &u1, NULL);
274 u = go_undo_combine (u, u1);
275 dump_sheet (sheet1, NULL);
278 for (i = 5; i >= 0; i--) {
279 g_printerr ("# About to insert row before %s\n",
280 row_name (i));
281 sheet_insert_rows (sheet1, i, 1, &u1, NULL);
282 u = go_undo_combine (u, u1);
283 dump_sheet (sheet1, NULL);
286 go_undo_undo (u);
287 g_object_unref (u);
288 u = NULL;
289 dump_sheet (sheet1, "Undo the lot");
291 for (i = 5; i >= 0; i--) {
292 g_printerr ("# About to delete column %s\n",
293 col_name (i));
294 sheet_delete_cols (sheet1, i, 1, &u1, NULL);
295 u = go_undo_combine (u, u1);
296 dump_sheet (sheet1, NULL);
299 for (i = 5; i >= 0; i--) {
300 g_printerr ("# About to delete row %s\n",
301 row_name (i));
302 sheet_delete_rows (sheet1, i, 1, &u1, NULL);
303 u = go_undo_combine (u, u1);
304 dump_sheet (sheet1, NULL);
307 go_undo_undo (u);
308 g_object_unref (u);
309 u = NULL;
310 dump_sheet (sheet1, "Undo the lot");
312 g_object_unref (wb);
314 mark_test_end (test_name);
317 /*-------------------------------------------------------------------------- */
319 static void
320 test_func_help (void)
322 const char *test_name = "test_func_help";
323 int res;
325 mark_test_start (test_name);
327 res = gnm_func_sanity_check ();
328 g_printerr ("Result = %d\n", res);
330 mark_test_end (test_name);
333 /*-------------------------------------------------------------------------- */
335 static int
336 test_strtol_ok (const char *s, long l, size_t expected_len)
338 long l2;
339 char *end;
340 int save_errno;
342 l2 = gnm_utf8_strtol (s, &end);
343 save_errno = errno;
345 if (end != s + expected_len) {
346 g_printerr ("Unexpect conversion end of [%s]\n", s);
347 return 1;
349 if (l != l2) {
350 g_printerr ("Unexpect conversion result of [%s]\n", s);
351 return 1;
353 if (save_errno != 0) {
354 g_printerr ("Unexpect conversion errno of [%s]\n", s);
355 return 1;
358 return 0;
361 static int
362 test_strtol_noconv (const char *s)
364 long l;
365 char *end;
366 int save_errno;
368 l = gnm_utf8_strtol (s, &end);
369 save_errno = errno;
371 if (end != s) {
372 g_printerr ("Unexpect conversion end of [%s]\n", s);
373 return 1;
375 if (l != 0) {
376 g_printerr ("Unexpect conversion result of [%s]\n", s);
377 return 1;
379 if (save_errno != 0) {
380 g_printerr ("Unexpect conversion errno of [%s]\n", s);
381 return 1;
384 return 0;
387 static int
388 test_strtol_overflow (const char *s, gboolean pos)
390 long l;
391 char *end;
392 int save_errno;
393 size_t expected_len = strlen (s);
395 l = gnm_utf8_strtol (s, &end);
396 save_errno = errno;
398 if (end != s + expected_len) {
399 g_printerr ("Unexpect conversion end of [%s]\n", s);
400 return 1;
402 if (l != (pos ? LONG_MAX : LONG_MIN)) {
403 g_printerr ("Unexpect conversion result of [%s]\n", s);
404 return 1;
406 if (save_errno != ERANGE) {
407 g_printerr ("Unexpect conversion errno of [%s]\n", s);
408 return 1;
411 return 0;
414 static int
415 test_strtol_reverse (long l)
417 char buffer[4*sizeof(l) + 4];
418 int res = 0;
420 sprintf(buffer, "%ld", l);
421 res |= test_strtol_ok (buffer, l, strlen (buffer));
423 sprintf(buffer, " %ld", l);
424 res |= test_strtol_ok (buffer, l, strlen (buffer));
426 sprintf(buffer, "\xc2\xa0\n\t%ld", l);
427 res |= test_strtol_ok (buffer, l, strlen (buffer));
429 sprintf(buffer, " \t%ldx", l);
430 res |= test_strtol_ok (buffer, l, strlen (buffer) - 1);
432 return res;
435 static int
436 test_strtod_ok (const char *s, double d, size_t expected_len)
438 gnm_float d2;
439 char *end;
440 int save_errno;
442 d2 = gnm_utf8_strto (s, &end);
443 save_errno = errno;
445 if (end != s + expected_len) {
446 g_printerr ("Unexpect conversion end of [%s]\n", s);
447 return 1;
449 if (d != d2) {
450 g_printerr ("Unexpect conversion result of [%s]\n", s);
451 return 1;
453 if (save_errno != 0) {
454 g_printerr ("Unexpect conversion errno of [%s]\n", s);
455 return 1;
458 return 0;
461 static void
462 test_nonascii_numbers (void)
464 const char *test_name = "test_nonascii_numbers";
465 int res = 0;
467 mark_test_start (test_name);
469 res |= test_strtol_reverse (0);
470 res |= test_strtol_reverse (1);
471 res |= test_strtol_reverse (-1);
472 res |= test_strtol_reverse (LONG_MIN);
473 res |= test_strtol_reverse (LONG_MIN + 1);
474 res |= test_strtol_reverse (LONG_MAX - 1);
476 res |= test_strtol_ok ("\xef\xbc\x8d\xef\xbc\x91", -1, 6);
477 res |= test_strtol_ok ("\xc2\xa0+1", 1, 4);
479 res |= test_strtol_ok ("000000000000000000000000000000", 0, 30);
481 res |= test_strtol_noconv ("");
482 res |= test_strtol_noconv (" ");
483 res |= test_strtol_noconv (" +");
484 res |= test_strtol_noconv (" -");
485 res |= test_strtol_noconv (" .00");
486 res |= test_strtol_noconv (" e0");
487 res |= test_strtol_noconv ("--0");
488 res |= test_strtol_noconv ("+-0");
489 res |= test_strtol_noconv ("+ 0");
490 res |= test_strtol_noconv ("- 0");
493 char buffer[4 * sizeof (long) + 2];
495 sprintf (buffer, "-%lu", 1 + (unsigned long)LONG_MIN);
496 res |= test_strtol_overflow (buffer, FALSE);
497 sprintf (buffer, "-%lu", 10 + (unsigned long)LONG_MIN);
498 res |= test_strtol_overflow (buffer, FALSE);
500 sprintf (buffer, "%lu", 1 + (unsigned long)LONG_MAX);
501 res |= test_strtol_overflow (buffer, TRUE);
502 sprintf (buffer, "%lu", 10 + (unsigned long)LONG_MAX);
503 res |= test_strtol_overflow (buffer, TRUE);
506 /* -------------------- */
508 res |= test_strtod_ok ("0", 0, 1);
509 res |= test_strtod_ok ("1", 1, 1);
510 res |= test_strtod_ok ("-1", -1, 2);
511 res |= test_strtod_ok ("+1", 1, 2);
512 res |= test_strtod_ok (" +1", 1, 3);
513 res |= test_strtod_ok ("\xc2\xa0+1", 1, 4);
514 res |= test_strtod_ok ("\xc2\xa0+1x", 1, 4);
515 res |= test_strtod_ok ("\xc2\xa0+1e", 1, 4);
516 res |= test_strtod_ok ("\xc2\xa0+1e+", 1, 4);
517 res |= test_strtod_ok ("\xc2\xa0+1e+0", 1, 7);
518 res |= test_strtod_ok ("-1e1", -10, 4);
519 res |= test_strtod_ok ("100e-2", 1, 6);
520 res |= test_strtod_ok ("100e+2", 10000, 6);
521 res |= test_strtod_ok ("1x0p0", 1, 1);
522 res |= test_strtod_ok ("+inf", gnm_pinf, 4);
523 res |= test_strtod_ok ("-inf", gnm_ninf, 4);
524 res |= test_strtod_ok ("1.25", 1.25, 4);
525 res |= test_strtod_ok ("1.25e1", 12.5, 6);
526 res |= test_strtod_ok ("12.5e-1", 1.25, 7);
528 g_printerr ("Result = %d\n", res);
530 mark_test_end (test_name);
533 /*-------------------------------------------------------------------------- */
535 static char *random_summary = NULL;
537 static void
538 add_random_fail (const char *s)
540 if (random_summary) {
541 char *t = g_strconcat (random_summary, ", ", s, NULL);
542 g_free (random_summary);
543 random_summary = t;
544 } else
545 random_summary = g_strdup (s);
548 static void
549 define_cell (Sheet *sheet, int c, int r, const char *expr)
551 GnmCell *cell = sheet_cell_fetch (sheet, c, r);
552 sheet_cell_set_text (cell, expr, NULL);
555 #define GET_PROB(i_) ((i_) <= 0 ? 0 : ((i_) >= nf ? 1 : probs[(i_)]))
557 static gboolean
558 rand_fractile_test (gnm_float const *vals, int N, int nf,
559 gnm_float const *fractiles, gnm_float const *probs)
561 gnm_float f = 1.0 / nf;
562 int *fractilecount = g_new (int, nf + 1);
563 int *expected = g_new (int, nf + 1);
564 int i;
565 gboolean ok = TRUE;
566 gboolean debug = TRUE;
568 if (debug) {
569 g_printerr ("Bin upper limit:");
570 for (i = 1; i <= nf; i++) {
571 gnm_float U = (i == nf) ? gnm_pinf : fractiles[i];
572 g_printerr ("%s%" GNM_FORMAT_g,
573 (i == 1) ? " " : ", ",
576 g_printerr (".\n");
579 if (debug && probs) {
580 g_printerr ("Cumulative probabilities:");
581 for (i = 1; i <= nf; i++)
582 g_printerr ("%s%.1" GNM_FORMAT_f "%%",
583 (i == 1) ? " " : ", ", 100 * GET_PROB (i));
584 g_printerr (".\n");
587 for (i = 1; i < nf - 1; i++) {
588 if (!(fractiles[i] <= fractiles[i + 1])) {
589 g_printerr ("Severe fractile ordering problem.\n");
590 return FALSE;
593 if (probs && !(probs[i] <= probs[i + 1])) {
594 g_printerr ("Severe cumulative probabilities ordering problem.\n");
595 return FALSE;
598 if (probs && (probs[1] < 0 || probs[nf - 1] > 1)) {
599 g_printerr ("Severe cumulative probabilities range problem.\n");
600 return FALSE;
603 for (i = 0; i <= nf; i++)
604 fractilecount[i] = 0;
606 for (i = 0; i < N; i++) {
607 gnm_float r = vals[i];
608 int j;
609 for (j = 1; j < nf; j++)
610 if (r <= fractiles[j])
611 break;
612 fractilecount[j]++;
614 g_printerr ("Fractile counts:");
615 for (i = 1; i <= nf; i++)
616 g_printerr ("%s%d", (i == 1) ? " " : ", ", fractilecount[i]);
617 g_printerr (".\n");
619 if (probs) {
620 g_printerr ("Expected counts:");
621 for (i = 1; i <= nf; i++) {
622 gnm_float p = GET_PROB (i) - GET_PROB (i-1);
623 expected[i] = gnm_floor (p * N + 0.5);
624 g_printerr ("%s%d", (i == 1) ? " " : ", ", expected[i]);
626 g_printerr (".\n");
627 } else {
628 gnm_float T = f * N;
629 g_printerr ("Expected count in each fractile: %.10" GNM_FORMAT_g "\n", T);
630 for (i = 0; i <= nf; i++)
631 expected[i] = T;
634 for (i = 1; i <= nf; i++) {
635 gnm_float T = expected[i];
636 if (!(gnm_abs (fractilecount[i] - T) <= 3 * gnm_sqrt (T))) {
637 g_printerr ("Fractile test failure for bin %d.\n", i);
638 ok = FALSE;
642 g_free (fractilecount);
643 g_free (expected);
645 return ok;
648 #undef GET_PROB
650 static gnm_float *
651 test_random_1 (int N, const char *expr,
652 gnm_float *mean, gnm_float *var,
653 gnm_float *skew, gnm_float *kurt)
655 Workbook *wb = workbook_new ();
656 Sheet *sheet;
657 gnm_float *res = g_new (gnm_float, N);
658 int i;
659 char *s;
660 int cols = 2, rows = N;
662 g_printerr ("Testing %s\n", expr);
664 gnm_sheet_suggest_size (&cols, &rows);
665 sheet = workbook_sheet_add (wb, -1, cols, rows);
667 for (i = 0; i < N; i++)
668 define_cell (sheet, 0, i, expr);
670 s = g_strdup_printf ("=average(a1:a%d)", N);
671 define_cell (sheet, 1, 0, s);
672 g_free (s);
674 s = g_strdup_printf ("=var(a1:a%d)", N);
675 define_cell (sheet, 1, 1, s);
676 g_free (s);
678 s = g_strdup_printf ("=skew(a1:a%d)", N);
679 define_cell (sheet, 1, 2, s);
680 g_free (s);
682 s = g_strdup_printf ("=kurt(a1:a%d)", N);
683 define_cell (sheet, 1, 3, s);
684 g_free (s);
686 /* Force recalc of all dirty cells even in manual mode. */
687 workbook_recalc (sheet->workbook);
689 for (i = 0; i < N; i++)
690 res[i] = value_get_as_float (sheet_cell_get (sheet, 0, i)->value);
691 *mean = value_get_as_float (sheet_cell_get (sheet, 1, 0)->value);
692 g_printerr ("Mean: %.10" GNM_FORMAT_g "\n", *mean);
694 *var = value_get_as_float (sheet_cell_get (sheet, 1, 1)->value);
695 g_printerr ("Var: %.10" GNM_FORMAT_g "\n", *var);
697 *skew = value_get_as_float (sheet_cell_get (sheet, 1, 2)->value);
698 g_printerr ("Skew: %.10" GNM_FORMAT_g "\n", *skew);
700 *kurt = value_get_as_float (sheet_cell_get (sheet, 1, 3)->value);
701 g_printerr ("Kurt: %.10" GNM_FORMAT_g "\n", *kurt);
703 g_object_unref (wb);
704 return res;
707 static gnm_float *
708 test_random_normality (int N, const char *expr,
709 gnm_float *mean, gnm_float *var,
710 gnm_float *adtest, gnm_float *cvmtest,
711 gnm_float *lkstest, gnm_float *sftest)
713 Workbook *wb = workbook_new ();
714 Sheet *sheet;
715 gnm_float *res = g_new (gnm_float, N);
716 int i;
717 char *s;
718 int cols = 2, rows = N;
720 g_printerr ("Testing %s\n", expr);
722 gnm_sheet_suggest_size (&cols, &rows);
723 sheet = workbook_sheet_add (wb, -1, cols, rows);
725 for (i = 0; i < N; i++)
726 define_cell (sheet, 0, i, expr);
728 s = g_strdup_printf ("=average(a1:a%d)", N);
729 define_cell (sheet, 1, 0, s);
730 g_free (s);
732 s = g_strdup_printf ("=var(a1:a%d)", N);
733 define_cell (sheet, 1, 1, s);
734 g_free (s);
736 s = g_strdup_printf ("=adtest(a1:a%d)", N);
737 define_cell (sheet, 1, 2, s);
738 g_free (s);
740 s = g_strdup_printf ("=cvmtest(a1:a%d)", N);
741 define_cell (sheet, 1, 3, s);
742 g_free (s);
744 s = g_strdup_printf ("=lkstest(a1:a%d)", N);
745 define_cell (sheet, 1, 4, s);
746 g_free (s);
748 s = g_strdup_printf ("=sftest(a1:a%d)", N > 5000 ? 5000 : N);
749 define_cell (sheet, 1, 5, s);
750 g_free (s);
752 /* Force recalc of all dirty cells even in manual mode. */
753 workbook_recalc (sheet->workbook);
755 for (i = 0; i < N; i++)
756 res[i] = value_get_as_float (sheet_cell_get (sheet, 0, i)->value);
757 *mean = value_get_as_float (sheet_cell_get (sheet, 1, 0)->value);
758 g_printerr ("Mean: %.10" GNM_FORMAT_g "\n", *mean);
760 *var = value_get_as_float (sheet_cell_get (sheet, 1, 1)->value);
761 g_printerr ("Var: %.10" GNM_FORMAT_g "\n", *var);
763 *adtest = value_get_as_float (sheet_cell_get (sheet, 1, 2)->value);
764 g_printerr ("ADTest: %.10" GNM_FORMAT_g "\n", *adtest);
766 *cvmtest = value_get_as_float (sheet_cell_get (sheet, 1, 3)->value);
767 g_printerr ("CVMTest: %.10" GNM_FORMAT_g "\n", *cvmtest);
769 *lkstest = value_get_as_float (sheet_cell_get (sheet, 1, 4)->value);
770 g_printerr ("LKSTest: %.10" GNM_FORMAT_g "\n", *lkstest);
772 *sftest = value_get_as_float (sheet_cell_get (sheet, 1, 5)->value);
773 g_printerr ("SFTest: %.10" GNM_FORMAT_g "\n", *sftest);
775 g_object_unref (wb);
776 return res;
779 static void
780 test_random_rand (int N)
782 gnm_float mean, var, skew, kurt;
783 gnm_float *vals;
784 int i;
785 gboolean ok;
786 gnm_float T;
787 gnm_float fractiles[10];
788 const int nf = G_N_ELEMENTS (fractiles);
790 vals = test_random_1 (N, "=RAND()", &mean, &var, &skew, &kurt);
791 ok = TRUE;
792 for (i = 0; i < N; i++) {
793 gnm_float r = vals[i];
794 if (!(r >= 0 && r < 1)) {
795 g_printerr ("Range failure.\n");
796 ok = FALSE;
797 break;
801 T = 0.5;
802 if (gnm_abs (mean - T) > 0.01) {
803 g_printerr ("Mean failure [%.10" GNM_FORMAT_g "]\n", T);
804 ok = FALSE;
806 T = 1.0 / 12;
807 if (gnm_abs (var - T) > 0.01) {
808 g_printerr ("Var failure [%.10" GNM_FORMAT_g "]\n", T);
809 ok = FALSE;
811 T = 0;
812 if (gnm_abs (skew - T) > 0.05) {
813 g_printerr ("Skew failure [%.10" GNM_FORMAT_g "]\n", T);
814 ok = FALSE;
816 T = -6.0 / 5;
817 if (gnm_abs (kurt - T) > 0.05) {
818 g_printerr ("Kurt failure [%.10" GNM_FORMAT_g "]\n", T);
819 ok = FALSE;
822 /* Fractile test */
823 for (i = 1; i < nf; i++)
824 fractiles[i] = i / (double)nf;
825 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
826 ok = FALSE;
828 if (ok)
829 g_printerr ("OK\n");
830 else
831 add_random_fail ("RAND");
832 g_printerr ("\n");
834 g_free (vals);
837 static void
838 test_random_randuniform (int N)
840 gnm_float mean, var, skew, kurt;
841 gnm_float *vals;
842 gboolean ok;
843 gnm_float lsign = (random_01 () > 0.75 ? 1 : -1);
844 gnm_float param_l = lsign * gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
845 gnm_float param_h = param_l + gnm_floor (1 / (0.0001 + gnm_pow (random_01 () / 2, 4)));
846 gnm_float n = param_h - param_l;
847 gnm_float mean_target = (param_l + param_h) / 2;
848 gnm_float var_target = (n * n) / 12;
849 gnm_float skew_target = 0;
850 gnm_float kurt_target = -6 / 5.0;
851 char *expr;
852 gnm_float T;
853 int i;
854 gnm_float fractiles[10];
855 const int nf = G_N_ELEMENTS (fractiles);
857 expr = g_strdup_printf ("=RANDUNIFORM(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")", param_l, param_h);
858 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
859 g_free (expr);
861 ok = TRUE;
862 for (i = 0; i < N; i++) {
863 gnm_float r = vals[i];
864 if (!(r >= param_l && r < param_h)) {
865 g_printerr ("Range failure.\n");
866 ok = FALSE;
867 break;
871 T = mean_target;
872 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
873 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
874 g_printerr ("Mean failure.\n");
875 ok = FALSE;
878 T = var_target;
879 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
880 if (!(var >= 0 && gnm_finite (var))) {
881 /* That is a very simplistic test! */
882 g_printerr ("Var failure.\n");
883 ok = FALSE;
886 T = skew_target;
887 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
888 if (!gnm_finite (skew)) {
889 /* That is a very simplistic test! */
890 g_printerr ("Skew failure.\n");
891 ok = FALSE;
894 T = kurt_target;
895 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
896 if (!(kurt >= -3 && gnm_finite (kurt))) {
897 /* That is a very simplistic test! */
898 g_printerr ("Kurt failure.\n");
899 ok = FALSE;
902 /* Fractile test */
903 for (i = 1; i < nf; i++)
904 fractiles[i] = param_l + n * i / (double)nf;
905 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
906 ok = FALSE;
908 if (ok)
909 g_printerr ("OK\n");
910 else
911 add_random_fail ("RANDUNIFORM");
912 g_printerr ("\n");
914 g_free (vals);
917 static void
918 test_random_randbernoulli (int N)
920 gnm_float mean, var, skew, kurt;
921 gnm_float *vals;
922 int i;
923 gboolean ok;
924 gnm_float p = 0.3;
925 gnm_float q = 1 - p;
926 char *expr;
927 gnm_float T;
929 expr = g_strdup_printf ("=RANDBERNOULLI(%.10" GNM_FORMAT_g ")", p);
930 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
931 g_free (expr);
933 ok = TRUE;
934 for (i = 0; i < N; i++) {
935 gnm_float r = vals[i];
936 if (!(r == 0 || r == 1)) {
937 g_printerr ("Range failure.\n");
938 ok = FALSE;
939 break;
943 T = p;
944 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
945 if (gnm_abs (mean - p) > 0.01) {
946 g_printerr ("Mean failure [%.10" GNM_FORMAT_g "]\n", T);
947 ok = FALSE;
950 T = p * (1 - p);
951 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
952 if (gnm_abs (var - T) > 0.01) {
953 g_printerr ("Var failure [%.10" GNM_FORMAT_g "]\n", T);
954 ok = FALSE;
957 T = (q - p) / gnm_sqrt (p * q);
958 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
959 if (gnm_abs (skew - T) > 0.05) {
960 g_printerr ("Skew failure [%.10" GNM_FORMAT_g "]\n", T);
961 ok = FALSE;
964 T = (1 - 6 * p * q) / (p * q);
965 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
966 if (gnm_abs (kurt - T) > 0.10) {
967 g_printerr ("Kurt failure [%.10" GNM_FORMAT_g "]\n", T);
968 ok = FALSE;
970 if (ok)
971 g_printerr ("OK\n");
972 else
973 add_random_fail ("RANDBERNOULLI");
974 g_printerr ("\n");
976 g_free (vals);
979 static void
980 test_random_randdiscrete (int N)
982 gnm_float mean, var, skew, kurt;
983 gnm_float *vals;
984 int i;
985 gboolean ok;
986 gnm_float mean_target = 13;
987 gnm_float var_target = 156;
988 gnm_float skew_target = 0.6748;
989 gnm_float kurt_target = -0.9057;
990 char *expr;
991 gnm_float T;
993 expr = g_strdup_printf ("=RANDDISCRETE({0;1;4;9;16;25;36})");
994 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
995 g_free (expr);
997 ok = TRUE;
998 for (i = 0; i < N; i++) {
999 gnm_float r = vals[i];
1000 if (!(r >= 0 && r <= 36 && gnm_sqrt (r) == gnm_floor (gnm_sqrt (r)))) {
1001 g_printerr ("Range failure.\n");
1002 ok = FALSE;
1003 break;
1007 T = mean_target;
1008 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1009 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1010 g_printerr ("Mean failure.\n");
1011 ok = FALSE;
1014 T = var_target;
1015 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1016 if (!(var >= 0 && gnm_finite (var))) {
1017 /* That is a very simplistic test! */
1018 g_printerr ("Var failure.\n");
1019 ok = FALSE;
1022 T = skew_target;
1023 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1024 if (!gnm_finite (skew)) {
1025 /* That is a very simplistic test! */
1026 g_printerr ("Skew failure.\n");
1027 ok = FALSE;
1030 T = kurt_target;
1031 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1032 if (!(kurt >= -3 && gnm_finite (kurt))) {
1033 /* That is a very simplistic test! */
1034 g_printerr ("Kurt failure.\n");
1035 ok = FALSE;
1038 if (ok)
1039 g_printerr ("OK\n");
1040 else
1041 add_random_fail ("RANDDISCRETE");
1042 g_printerr ("\n");
1044 g_free (vals);
1047 static void
1048 test_random_randnorm (int N)
1050 gnm_float mean, var, adtest, cvmtest, lkstest, sftest;
1051 gnm_float mean_target = 0, var_target = 1;
1052 gnm_float *vals;
1053 gboolean ok;
1054 char *expr;
1055 gnm_float T;
1056 int i;
1057 gnm_float fractiles[10];
1058 const int nf = G_N_ELEMENTS (fractiles);
1060 expr = g_strdup_printf ("=RANDNORM(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")",
1061 mean_target, var_target);
1062 vals = test_random_normality (N, expr, &mean, &var, &adtest, &cvmtest, &lkstest, &sftest);
1063 g_free (expr);
1065 ok = TRUE;
1066 for (i = 0; i < N; i++) {
1067 gnm_float r = vals[i];
1068 if (!gnm_finite (r)) {
1069 g_printerr ("Range failure.\n");
1070 ok = FALSE;
1071 break;
1075 T = mean_target;
1076 if (gnm_abs (mean - T) > 0.02) {
1077 g_printerr ("Mean failure [%.10" GNM_FORMAT_g "]\n", T);
1078 ok = FALSE;
1080 T = var_target;
1081 if (gnm_abs (var - T) > 0.02) {
1082 g_printerr ("Var failure [%.10" GNM_FORMAT_g "]\n", T);
1083 ok = FALSE;
1086 /* Fractile test */
1087 for (i = 1; i < nf; i++)
1088 fractiles[i] = qnorm (i / (double)nf, mean_target, gnm_sqrt (var_target), TRUE, FALSE);
1089 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1090 ok = FALSE;
1092 if (adtest < 0.05) {
1093 g_printerr ("Anderson Darling Test rejected [%.10" GNM_FORMAT_g "]\n", adtest);
1094 ok = FALSE;
1096 if (cvmtest < 0.05) {
1097 g_printerr ("Cramér-von Mises Test rejected [%.10" GNM_FORMAT_g "]\n", cvmtest);
1098 ok = FALSE;
1100 if (lkstest < 0.01) {
1101 g_printerr ("Lilliefors (Kolmogorov-Smirnov) Test rejected [%.10" GNM_FORMAT_g "]\n",
1102 lkstest);
1103 ok = FALSE;
1105 if (sftest < 0.05) {
1106 g_printerr ("Shapiro-Francia Test rejected [%.10" GNM_FORMAT_g "]\n", sftest);
1107 ok = FALSE;
1110 if (ok)
1111 g_printerr ("OK\n");
1112 else
1113 add_random_fail ("RANDNORM");
1114 g_printerr ("\n");
1116 g_free (vals);
1119 static void
1120 test_random_randsnorm (int N)
1122 gnm_float mean, var, skew, kurt;
1123 gnm_float *vals;
1124 gboolean ok;
1125 gnm_float alpha = 5;
1126 gnm_float delta = alpha/gnm_sqrt(1+alpha*alpha);
1127 gnm_float mean_target = delta * gnm_sqrt (2/M_PIgnum);
1128 gnm_float var_target = 1-mean_target*mean_target;
1129 char *expr;
1130 gnm_float T;
1131 int i;
1133 expr = g_strdup_printf ("=RANDSNORM(%.10" GNM_FORMAT_g ")", alpha);
1134 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1135 g_free (expr);
1137 ok = TRUE;
1138 for (i = 0; i < N; i++) {
1139 gnm_float r = vals[i];
1140 if (!gnm_finite (r)) {
1141 g_printerr ("Range failure.\n");
1142 ok = FALSE;
1143 break;
1147 T = mean_target;
1148 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1149 if (gnm_abs (mean - T) > 0.01) {
1150 g_printerr ("Mean failure [%.10" GNM_FORMAT_g "]\n", T);
1151 ok = FALSE;
1154 T = var_target;
1155 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1156 if (gnm_abs (var - T) > 0.01) {
1157 g_printerr ("Var failure [%.10" GNM_FORMAT_g "]\n", T);
1158 ok = FALSE;
1161 T = mean_target/gnm_sqrt(var_target);
1162 T = T*T*T*(4-M_PIgnum)/2;
1163 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1164 if (gnm_abs (skew - T) > 0.05) {
1165 g_printerr ("Skew failure [%.10" GNM_FORMAT_g "]\n", T);
1166 ok = FALSE;
1169 T = 2*(M_PIgnum - 3)*mean_target*mean_target*mean_target*mean_target/(var_target*var_target);
1170 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1171 if (gnm_abs (kurt - T) > 0.15) {
1172 g_printerr ("Kurt failure [%.10" GNM_FORMAT_g "]\n", T);
1173 ok = FALSE;
1176 if (ok)
1177 g_printerr ("OK\n");
1178 else
1179 add_random_fail ("RANDSNORM");
1180 g_printerr ("\n");
1182 g_free (vals);
1185 static void
1186 test_random_randexp (int N)
1188 gnm_float mean, var, skew, kurt;
1189 gnm_float *vals;
1190 gboolean ok;
1191 gnm_float param_l = 1 / (0.0001 + gnm_pow (random_01 () / 2, 4));
1192 gnm_float mean_target = param_l;
1193 gnm_float var_target = mean_target * mean_target;
1194 gnm_float skew_target = 2;
1195 gnm_float kurt_target = 6;
1196 char *expr;
1197 gnm_float T;
1198 int i;
1199 gnm_float fractiles[10];
1200 const int nf = G_N_ELEMENTS (fractiles);
1202 expr = g_strdup_printf ("=RANDEXP(%.10" GNM_FORMAT_g ")", param_l);
1203 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1204 g_free (expr);
1206 ok = TRUE;
1207 for (i = 0; i < N; i++) {
1208 gnm_float r = vals[i];
1209 if (!(r >= 0 && gnm_finite (r))) {
1210 g_printerr ("Range failure.\n");
1211 ok = FALSE;
1212 break;
1216 T = mean_target;
1217 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1218 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1219 g_printerr ("Mean failure.\n");
1220 ok = FALSE;
1223 T = var_target;
1224 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1225 if (!(var >= 0 && gnm_finite (var))) {
1226 /* That is a very simplistic test! */
1227 g_printerr ("Var failure.\n");
1228 ok = FALSE;
1231 T = skew_target;
1232 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1233 if (!gnm_finite (skew)) {
1234 /* That is a very simplistic test! */
1235 g_printerr ("Skew failure.\n");
1236 ok = FALSE;
1239 T = kurt_target;
1240 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1241 if (!(kurt >= -3 && gnm_finite (kurt))) {
1242 /* That is a very simplistic test! */
1243 g_printerr ("Kurt failure.\n");
1244 ok = FALSE;
1247 /* Fractile test */
1248 for (i = 1; i < nf; i++)
1249 fractiles[i] = qexp (i / (double)nf, param_l, TRUE, FALSE);
1250 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1251 ok = FALSE;
1253 if (ok)
1254 g_printerr ("OK\n");
1255 else
1256 add_random_fail ("RANDEXP");
1257 g_printerr ("\n");
1259 g_free (vals);
1262 static void
1263 test_random_randgamma (int N)
1265 gnm_float mean, var, skew, kurt;
1266 gnm_float *vals;
1267 gboolean ok;
1268 gnm_float param_shape = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 6)));
1269 gnm_float param_scale = 0.001 + gnm_pow (random_01 (), 4) * 1000;
1270 gnm_float mean_target = param_shape * param_scale;
1271 gnm_float var_target = mean_target * param_scale;
1272 gnm_float skew_target = 2 / gnm_sqrt (param_shape);
1273 gnm_float kurt_target = 6 / param_shape;
1274 char *expr;
1275 gnm_float T;
1276 int i;
1277 gnm_float fractiles[10];
1278 const int nf = G_N_ELEMENTS (fractiles);
1280 expr = g_strdup_printf ("=RANDGAMMA(%.0" GNM_FORMAT_f ",%.10" GNM_FORMAT_g ")", param_shape, param_scale);
1281 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1282 g_free (expr);
1284 ok = TRUE;
1285 for (i = 0; i < N; i++) {
1286 gnm_float r = vals[i];
1287 if (!(r > 0 && gnm_finite (r))) {
1288 g_printerr ("Range failure.\n");
1289 ok = FALSE;
1290 break;
1294 T = mean_target;
1295 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1296 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1297 g_printerr ("Mean failure.\n");
1298 ok = FALSE;
1301 T = var_target;
1302 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1303 if (!(var >= 0 && gnm_finite (var))) {
1304 /* That is a very simplistic test! */
1305 g_printerr ("Var failure.\n");
1306 ok = FALSE;
1309 T = skew_target;
1310 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1311 if (!gnm_finite (skew)) {
1312 /* That is a very simplistic test! */
1313 g_printerr ("Skew failure.\n");
1314 ok = FALSE;
1317 T = kurt_target;
1318 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1319 if (!(kurt >= -3 && gnm_finite (kurt))) {
1320 /* That is a very simplistic test! */
1321 g_printerr ("Kurt failure.\n");
1322 ok = FALSE;
1325 /* Fractile test */
1326 for (i = 1; i < nf; i++)
1327 fractiles[i] = qgamma (i / (double)nf, param_shape, param_scale, TRUE, FALSE);
1328 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1329 ok = FALSE;
1331 if (ok)
1332 g_printerr ("OK\n");
1333 else
1334 add_random_fail ("RANDGAMMA");
1335 g_printerr ("\n");
1337 g_free (vals);
1340 static void
1341 test_random_randbeta (int N)
1343 gnm_float mean, var, skew, kurt;
1344 gnm_float *vals;
1345 gboolean ok;
1346 gnm_float param_a = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 6))); gnm_float param_b = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 6)));
1347 gnm_float s = param_a + param_b;
1348 gnm_float mean_target = param_a / s;
1349 gnm_float var_target = mean_target * param_b / (s * (s + 1));
1350 gnm_float skew_target =
1351 (2 * (param_b - param_a) * gnm_sqrt (s + 1))/
1352 ((s + 2) * gnm_sqrt (param_a * param_b));
1353 gnm_float kurt_target = gnm_nan; /* Complicated */
1354 char *expr;
1355 gnm_float T;
1356 int i;
1357 gnm_float fractiles[10];
1358 const int nf = G_N_ELEMENTS (fractiles);
1360 expr = g_strdup_printf ("=RANDBETA(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")", param_a, param_b);
1361 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1362 g_free (expr);
1364 ok = TRUE;
1365 for (i = 0; i < N; i++) {
1366 gnm_float r = vals[i];
1367 if (!(r >= 0 && r <= 1)) {
1368 g_printerr ("Range failure.\n");
1369 ok = FALSE;
1370 break;
1374 T = mean_target;
1375 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1376 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1377 g_printerr ("Mean failure.\n");
1378 ok = FALSE;
1381 T = var_target;
1382 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1383 if (!(var >= 0 && gnm_finite (var))) {
1384 /* That is a very simplistic test! */
1385 g_printerr ("Var failure.\n");
1386 ok = FALSE;
1389 T = skew_target;
1390 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1391 if (!gnm_finite (skew)) {
1392 /* That is a very simplistic test! */
1393 g_printerr ("Skew failure.\n");
1394 ok = FALSE;
1397 T = kurt_target;
1398 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1399 if (!(kurt >= -3 && gnm_finite (kurt))) {
1400 /* That is a very simplistic test! */
1401 g_printerr ("Kurt failure.\n");
1402 ok = FALSE;
1405 /* Fractile test */
1406 for (i = 1; i < nf; i++)
1407 fractiles[i] = qbeta (i / (double)nf, param_a, param_b, TRUE, FALSE);
1408 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1409 ok = FALSE;
1411 if (ok)
1412 g_printerr ("OK\n");
1413 else
1414 add_random_fail ("RANDBETA");
1415 g_printerr ("\n");
1417 g_free (vals);
1420 static void
1421 test_random_randtdist (int N)
1423 gnm_float mean, var, skew, kurt;
1424 gnm_float *vals;
1425 gboolean ok;
1426 gnm_float param_df = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1427 gnm_float mean_target = 0;
1428 gnm_float var_target = param_df > 2 ? param_df / (param_df - 2) : gnm_nan;
1429 gnm_float skew_target = param_df > 3 ? 0 : gnm_nan;
1430 gnm_float kurt_target = param_df > 4 ? 6 / (param_df - 4) : gnm_nan;
1431 char *expr;
1432 gnm_float T;
1433 int i;
1434 gnm_float fractiles[10];
1435 const int nf = G_N_ELEMENTS (fractiles);
1437 expr = g_strdup_printf ("=RANDTDIST(%.0" GNM_FORMAT_f ")", param_df);
1438 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1439 g_free (expr);
1441 ok = TRUE;
1442 for (i = 0; i < N; i++) {
1443 gnm_float r = vals[i];
1444 if (!(gnm_finite (r))) {
1445 g_printerr ("Range failure.\n");
1446 ok = FALSE;
1447 break;
1451 T = mean_target;
1452 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1453 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1454 g_printerr ("Mean failure.\n");
1455 ok = FALSE;
1458 T = var_target;
1459 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1460 if (!(var >= 0 && gnm_finite (var))) {
1461 /* That is a very simplistic test! */
1462 g_printerr ("Var failure.\n");
1463 ok = FALSE;
1466 T = skew_target;
1467 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1468 if (!gnm_finite (skew)) {
1469 /* That is a very simplistic test! */
1470 g_printerr ("Skew failure.\n");
1471 ok = FALSE;
1474 T = kurt_target;
1475 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1476 if (!(kurt >= -3 && gnm_finite (kurt))) {
1477 /* That is a very simplistic test! */
1478 g_printerr ("Kurt failure.\n");
1479 ok = FALSE;
1482 /* Fractile test */
1483 for (i = 1; i < nf; i++)
1484 fractiles[i] = qt (i / (double)nf, param_df, TRUE, FALSE);
1485 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1486 ok = FALSE;
1488 if (ok)
1489 g_printerr ("OK\n");
1490 else
1491 add_random_fail ("RANDTDIST");
1492 g_printerr ("\n");
1494 g_free (vals);
1497 static void
1498 test_random_randfdist (int N)
1500 gnm_float mean, var, skew, kurt;
1501 gnm_float *vals;
1502 gboolean ok;
1503 gnm_float param_df1 = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1504 gnm_float param_df2 = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1505 gnm_float mean_target = param_df2 > 2 ? param_df2 / (param_df2 - 2) : gnm_nan;
1506 gnm_float var_target = param_df2 > 4
1507 ? (2 * param_df2 * param_df2 * (param_df1 + param_df2 - 2) /
1508 (param_df1 * (param_df2 - 2) * (param_df2 - 2) * (param_df2 - 4)))
1509 : gnm_nan;
1510 gnm_float skew_target = gnm_nan; /* Complicated */
1511 gnm_float kurt_target = gnm_nan; /* Complicated */
1512 char *expr;
1513 gnm_float T;
1514 int i;
1515 gnm_float fractiles[10];
1516 const int nf = G_N_ELEMENTS (fractiles);
1518 expr = g_strdup_printf ("=RANDFDIST(%.0" GNM_FORMAT_f ",%.0" GNM_FORMAT_f ")", param_df1, param_df2);
1519 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1520 g_free (expr);
1522 ok = TRUE;
1523 for (i = 0; i < N; i++) {
1524 gnm_float r = vals[i];
1525 if (!(r >= 0 && gnm_finite (r))) {
1526 g_printerr ("Range failure.\n");
1527 ok = FALSE;
1528 break;
1532 T = mean_target;
1533 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1534 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1535 g_printerr ("Mean failure.\n");
1536 ok = FALSE;
1539 T = var_target;
1540 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1541 if (!(var >= 0 && gnm_finite (var))) {
1542 /* That is a very simplistic test! */
1543 g_printerr ("Var failure.\n");
1544 ok = FALSE;
1547 T = skew_target;
1548 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1549 if (!gnm_finite (skew)) {
1550 /* That is a very simplistic test! */
1551 g_printerr ("Skew failure.\n");
1552 ok = FALSE;
1555 T = kurt_target;
1556 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1557 if (!(kurt >= -3 && gnm_finite (kurt))) {
1558 /* That is a very simplistic test! */
1559 g_printerr ("Kurt failure.\n");
1560 ok = FALSE;
1563 /* Fractile test */
1564 for (i = 1; i < nf; i++)
1565 fractiles[i] = qf (i / (double)nf, param_df1, param_df2, TRUE, FALSE);
1566 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1567 ok = FALSE;
1569 if (ok)
1570 g_printerr ("OK\n");
1571 else
1572 add_random_fail ("RANDFDIST");
1573 g_printerr ("\n");
1575 g_free (vals);
1578 static void
1579 test_random_randchisq (int N)
1581 gnm_float mean, var, skew, kurt;
1582 gnm_float *vals;
1583 gboolean ok;
1584 gnm_float param_df = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1585 gnm_float mean_target = param_df;
1586 gnm_float var_target = param_df * 2;
1587 gnm_float skew_target = gnm_sqrt (8 / param_df);
1588 gnm_float kurt_target = 12 / param_df;
1589 char *expr;
1590 gnm_float T;
1591 int i;
1592 gnm_float fractiles[10];
1593 const int nf = G_N_ELEMENTS (fractiles);
1595 expr = g_strdup_printf ("=RANDCHISQ(%.10" GNM_FORMAT_g ")", param_df);
1596 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1597 g_free (expr);
1599 ok = TRUE;
1600 for (i = 0; i < N; i++) {
1601 gnm_float r = vals[i];
1602 if (!(r >= 0 && gnm_finite (r))) {
1603 g_printerr ("Range failure.\n");
1604 ok = FALSE;
1605 break;
1609 T = mean_target;
1610 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1611 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1612 g_printerr ("Mean failure.\n");
1613 ok = FALSE;
1616 T = var_target;
1617 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1618 if (!(var >= 0 && gnm_finite (var))) {
1619 /* That is a very simplistic test! */
1620 g_printerr ("Var failure.\n");
1621 ok = FALSE;
1624 T = skew_target;
1625 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1626 if (!gnm_finite (skew)) {
1627 /* That is a very simplistic test! */
1628 g_printerr ("Skew failure.\n");
1629 ok = FALSE;
1632 T = kurt_target;
1633 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1634 if (!(kurt >= -3 && gnm_finite (kurt))) {
1635 /* That is a very simplistic test! */
1636 g_printerr ("Kurt failure.\n");
1637 ok = FALSE;
1640 /* Fractile test */
1641 for (i = 1; i < nf; i++)
1642 fractiles[i] = qchisq (i / (double)nf, param_df, TRUE, FALSE);
1643 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1644 ok = FALSE;
1646 if (ok)
1647 g_printerr ("OK\n");
1648 else
1649 add_random_fail ("RANDCHISQ");
1650 g_printerr ("\n");
1652 g_free (vals);
1655 static void
1656 test_random_randcauchy (int N)
1658 gnm_float mean, var, skew, kurt;
1659 gnm_float *vals;
1660 gboolean ok;
1661 gnm_float param_scale = 0.001 + gnm_pow (random_01 (), 4) * 1000;
1662 gnm_float mean_target = gnm_nan;
1663 gnm_float var_target = gnm_nan;
1664 gnm_float skew_target = gnm_nan;
1665 gnm_float kurt_target = gnm_nan;
1666 char *expr;
1667 gnm_float T;
1668 int i;
1669 gnm_float fractiles[10];
1670 const int nf = G_N_ELEMENTS (fractiles);
1673 * The distribution has no mean, no variance, no skew, and no kurtosis.
1674 * The support is all reals.
1677 expr = g_strdup_printf ("=RANDCAUCHY(%.10" GNM_FORMAT_g ")", param_scale);
1678 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1679 g_free (expr);
1681 ok = TRUE;
1682 for (i = 0; i < N; i++) {
1683 gnm_float r = vals[i];
1684 if (!(gnm_finite (r))) {
1685 g_printerr ("Range failure.\n");
1686 ok = FALSE;
1687 break;
1691 T = mean_target;
1692 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1693 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1694 g_printerr ("Mean failure.\n");
1695 ok = FALSE;
1698 T = var_target;
1699 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1700 if (!(var >= 0 && gnm_finite (var))) {
1701 /* That is a very simplistic test! */
1702 g_printerr ("Var failure.\n");
1703 ok = FALSE;
1706 T = skew_target;
1707 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1708 if (!gnm_finite (skew)) {
1709 /* That is a very simplistic test! */
1710 g_printerr ("Skew failure.\n");
1711 ok = FALSE;
1714 T = kurt_target;
1715 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1716 if (!(kurt >= -3 && gnm_finite (kurt))) {
1717 /* That is a very simplistic test! */
1718 g_printerr ("Kurt failure.\n");
1719 ok = FALSE;
1722 /* Fractile test */
1723 for (i = 1; i < nf; i++)
1724 fractiles[i] = qcauchy (i / (double)nf, 0.0, param_scale, TRUE, FALSE);
1725 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1726 ok = FALSE;
1728 if (ok)
1729 g_printerr ("OK\n");
1730 else
1731 add_random_fail ("RANDCAUCHY");
1732 g_printerr ("\n");
1734 g_free (vals);
1737 static void
1738 test_random_randbinom (int N)
1740 gnm_float mean, var, skew, kurt;
1741 gnm_float *vals;
1742 gboolean ok;
1743 gnm_float param_p = random_01 ();
1744 gnm_float param_trials = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
1745 gnm_float mean_target = param_trials * param_p;
1746 gnm_float var_target = mean_target * (1 - param_p);
1747 gnm_float skew_target = (1 - 2 * param_p) / gnm_sqrt (var_target);
1748 gnm_float kurt_target = (1 - 6 * param_p * (1 - param_p)) / var_target;
1749 char *expr;
1750 gnm_float T;
1751 int i;
1752 gnm_float fractiles[10], probs[10];
1753 const int nf = G_N_ELEMENTS (fractiles);
1755 expr = g_strdup_printf ("=RANDBINOM(%.10" GNM_FORMAT_g ",%.0" GNM_FORMAT_f ")", param_p, param_trials);
1756 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1757 g_free (expr);
1759 ok = TRUE;
1760 for (i = 0; i < N; i++) {
1761 gnm_float r = vals[i];
1762 if (!(r >= 0 && r <= param_trials && r == gnm_floor (r))) {
1763 g_printerr ("Range failure.\n");
1764 ok = FALSE;
1765 break;
1769 T = mean_target;
1770 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1771 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1772 g_printerr ("Mean failure.\n");
1773 ok = FALSE;
1776 T = var_target;
1777 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1778 if (!(var >= 0 && gnm_finite (var))) {
1779 /* That is a very simplistic test! */
1780 g_printerr ("Var failure.\n");
1781 ok = FALSE;
1784 T = skew_target;
1785 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1786 if (!gnm_finite (skew)) {
1787 /* That is a very simplistic test! */
1788 g_printerr ("Skew failure.\n");
1789 ok = FALSE;
1792 T = kurt_target;
1793 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1794 if (!(kurt >= -3 && gnm_finite (kurt))) {
1795 /* That is a very simplistic test! */
1796 g_printerr ("Kurt failure.\n");
1797 ok = FALSE;
1800 /* Fractile test */
1801 for (i = 1; i < nf; i++) {
1802 fractiles[i] = qbinom (i / (double)nf, param_trials, param_p, TRUE, FALSE);
1803 probs[i] = pbinom (fractiles[i], param_trials, param_p, TRUE, FALSE);
1805 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
1806 ok = FALSE;
1808 if (ok)
1809 g_printerr ("OK\n");
1810 else
1811 add_random_fail ("RANDBINOM");
1812 g_printerr ("\n");
1814 g_free (vals);
1817 static void
1818 test_random_randnegbinom (int N)
1820 gnm_float mean, var, skew, kurt;
1821 gnm_float *vals;
1822 gboolean ok;
1823 gnm_float param_p = random_01 ();
1824 gnm_float param_fails = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
1825 /* Warning: these differ from Wikipedia by swapping p and 1-p. */
1826 gnm_float mean_target = param_fails * (1 - param_p) / param_p;
1827 gnm_float var_target = mean_target / param_p;
1828 gnm_float skew_target = (2 - param_p) / gnm_sqrt (param_fails * (1 - param_p));
1829 gnm_float kurt_target = 6 / param_fails + 1 / var_target;
1830 char *expr;
1831 gnm_float T;
1832 int i;
1833 gnm_float fractiles[10], probs[10];
1834 const int nf = G_N_ELEMENTS (fractiles);
1836 expr = g_strdup_printf ("=RANDNEGBINOM(%.10" GNM_FORMAT_g ",%.0" GNM_FORMAT_f ")", param_p, param_fails);
1837 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1838 g_free (expr);
1840 ok = TRUE;
1841 for (i = 0; i < N; i++) {
1842 gnm_float r = vals[i];
1843 if (!(r >= 0 && gnm_finite (r) && r == gnm_floor (r))) {
1844 g_printerr ("Range failure.\n");
1845 ok = FALSE;
1846 break;
1850 T = mean_target;
1851 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1852 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1853 g_printerr ("Mean failure.\n");
1854 ok = FALSE;
1857 T = var_target;
1858 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1859 if (!(var >= 0 && gnm_finite (var))) {
1860 /* That is a very simplistic test! */
1861 g_printerr ("Var failure.\n");
1862 ok = FALSE;
1865 T = skew_target;
1866 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1867 if (!gnm_finite (skew)) {
1868 /* That is a very simplistic test! */
1869 g_printerr ("Skew failure.\n");
1870 ok = FALSE;
1873 T = kurt_target;
1874 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1875 if (!(kurt >= -3 && gnm_finite (kurt))) {
1876 /* That is a very simplistic test! */
1877 g_printerr ("Kurt failure.\n");
1878 ok = FALSE;
1881 /* Fractile test */
1882 for (i = 1; i < nf; i++) {
1883 fractiles[i] = qnbinom (i / (double)nf, param_fails, param_p, TRUE, FALSE);
1884 probs[i] = pnbinom (fractiles[i], param_fails, param_p, TRUE, FALSE);
1886 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
1887 ok = FALSE;
1889 if (ok)
1890 g_printerr ("OK\n");
1891 else
1892 add_random_fail ("RANDNEGBINOM");
1893 g_printerr ("\n");
1895 g_free (vals);
1898 static void
1899 test_random_randhyperg (int N)
1901 gnm_float mean, var, skew, kurt;
1902 gnm_float *vals;
1903 gboolean ok;
1904 gnm_float param_nr = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 4)));
1905 gnm_float param_nb = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 4)));
1906 gnm_float s = param_nr + param_nb;
1907 gnm_float param_n = gnm_floor (random_01 () * (s + 1));
1908 gnm_float mean_target = param_n * param_nr / s;
1909 gnm_float var_target = s > 1
1910 ? mean_target * (param_nb / s) * (s - param_n) / (s - 1)
1911 : 0;
1912 gnm_float skew_target = gnm_nan; /* Complicated */
1913 gnm_float kurt_target = gnm_nan; /* Complicated */
1914 char *expr;
1915 gnm_float T;
1916 int i;
1917 gnm_float fractiles[10], probs[10];
1918 const int nf = G_N_ELEMENTS (fractiles);
1920 expr = g_strdup_printf ("=RANDHYPERG(%.10" GNM_FORMAT_g ",%.0" GNM_FORMAT_f ",%.0" GNM_FORMAT_f ")", param_nr, param_nb, param_n);
1921 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1922 g_free (expr);
1924 ok = TRUE;
1925 for (i = 0; i < N; i++) {
1926 gnm_float r = vals[i];
1927 if (!(r >= 0 && r <= param_n && r == gnm_floor (r))) {
1928 g_printerr ("Range failure.\n");
1929 ok = FALSE;
1930 break;
1934 T = mean_target;
1935 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1936 if (gnm_finite (var_target) &&
1937 !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1938 g_printerr ("Mean failure.\n");
1939 ok = FALSE;
1942 T = var_target;
1943 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1944 if (!(var >= 0 && gnm_finite (var))) {
1945 /* That is a very simplistic test! */
1946 g_printerr ("Var failure.\n");
1947 ok = FALSE;
1950 T = skew_target;
1951 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1952 if (!gnm_finite (skew)) {
1953 /* That is a very simplistic test! */
1954 g_printerr ("Skew failure.\n");
1955 ok = FALSE;
1958 T = kurt_target;
1959 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1960 if (!(kurt >= -3 && gnm_finite (kurt))) {
1961 /* That is a very simplistic test! */
1962 g_printerr ("Kurt failure.\n");
1963 ok = FALSE;
1966 /* Fractile test */
1967 for (i = 1; i < nf; i++) {
1968 fractiles[i] = qhyper (i / (double)nf, param_nr, param_nb, param_n, TRUE, FALSE);
1969 probs[i] = phyper (fractiles[i], param_nr, param_nb, param_n, TRUE, FALSE);
1971 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
1972 ok = FALSE;
1974 if (ok)
1975 g_printerr ("OK\n");
1976 else
1977 add_random_fail ("RANDHYPERG");
1978 g_printerr ("\n");
1980 g_free (vals);
1983 static void
1984 test_random_randbetween (int N)
1986 gnm_float mean, var, skew, kurt;
1987 gnm_float *vals;
1988 gboolean ok;
1989 gnm_float lsign = (random_01 () > 0.75 ? 1 : -1);
1990 gnm_float param_l = lsign * gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
1991 gnm_float param_h = param_l + gnm_floor (1 / (0.0001 + gnm_pow (random_01 () / 2, 4)));
1992 gnm_float n = param_h - param_l + 1;
1993 gnm_float mean_target = (param_l + param_h) / 2;
1994 gnm_float var_target = (n * n - 1) / 12;
1995 gnm_float skew_target = 0;
1996 gnm_float kurt_target = (n * n + 1) / (n * n - 1) * -6 / 5;
1997 char *expr;
1998 gnm_float T;
1999 int i;
2001 expr = g_strdup_printf ("=RANDBETWEEN(%.0" GNM_FORMAT_f ",%.0" GNM_FORMAT_f ")", param_l, param_h);
2002 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2003 g_free (expr);
2005 ok = TRUE;
2006 for (i = 0; i < N; i++) {
2007 gnm_float r = vals[i];
2008 if (!(r >= param_l && r <= param_h && r == gnm_floor (r))) {
2009 g_printerr ("Range failure.\n");
2010 ok = FALSE;
2011 break;
2015 T = mean_target;
2016 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2017 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2018 g_printerr ("Mean failure.\n");
2019 ok = FALSE;
2022 T = var_target;
2023 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2024 if (!(var >= 0 && gnm_finite (var))) {
2025 /* That is a very simplistic test! */
2026 g_printerr ("Var failure.\n");
2027 ok = FALSE;
2030 T = skew_target;
2031 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2032 if (!gnm_finite (skew)) {
2033 /* That is a very simplistic test! */
2034 g_printerr ("Skew failure.\n");
2035 ok = FALSE;
2038 T = kurt_target;
2039 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2040 if (!(kurt >= -3 && gnm_finite (kurt))) {
2041 /* That is a very simplistic test! */
2042 g_printerr ("Kurt failure.\n");
2043 ok = FALSE;
2046 if (ok)
2047 g_printerr ("OK\n");
2048 else
2049 add_random_fail ("RANDBETWEEN");
2050 g_printerr ("\n");
2052 g_free (vals);
2055 static void
2056 test_random_randpoisson (int N)
2058 gnm_float mean, var, skew, kurt;
2059 gnm_float *vals;
2060 gboolean ok;
2061 gnm_float param_l = 1 / (0.0001 + gnm_pow (random_01 () / 2, 4));
2062 gnm_float mean_target = param_l;
2063 gnm_float var_target = param_l;
2064 gnm_float skew_target = 1 / gnm_sqrt (param_l);
2065 gnm_float kurt_target = 1 / param_l;
2066 char *expr;
2067 gnm_float T;
2068 int i;
2069 gnm_float fractiles[10], probs[10];
2070 const int nf = G_N_ELEMENTS (fractiles);
2072 expr = g_strdup_printf ("=RANDPOISSON(%.10" GNM_FORMAT_g ")", param_l);
2073 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2074 g_free (expr);
2076 ok = TRUE;
2077 for (i = 0; i < N; i++) {
2078 gnm_float r = vals[i];
2079 if (!(r >= 0 && gnm_finite (r) && r == gnm_floor (r))) {
2080 g_printerr ("Range failure.\n");
2081 ok = FALSE;
2082 break;
2086 T = mean_target;
2087 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2088 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2089 g_printerr ("Mean failure.\n");
2090 ok = FALSE;
2093 T = var_target;
2094 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2095 if (!(var >= 0 && gnm_finite (var))) {
2096 /* That is a very simplistic test! */
2097 g_printerr ("Var failure.\n");
2098 ok = FALSE;
2101 T = skew_target;
2102 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2103 if (!gnm_finite (skew)) {
2104 /* That is a very simplistic test! */
2105 g_printerr ("Skew failure.\n");
2106 ok = FALSE;
2109 T = kurt_target;
2110 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2111 if (!(kurt >= -3 && gnm_finite (kurt))) {
2112 /* That is a very simplistic test! */
2113 g_printerr ("Kurt failure.\n");
2114 ok = FALSE;
2117 /* Fractile test */
2118 for (i = 1; i < nf; i++) {
2119 fractiles[i] = qpois (i / (double)nf, param_l, TRUE, FALSE);
2120 probs[i] = ppois (fractiles[i], param_l, TRUE, FALSE);
2122 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
2123 ok = FALSE;
2125 if (ok)
2126 g_printerr ("OK\n");
2127 else
2128 add_random_fail ("RANDPOISSON");
2129 g_printerr ("\n");
2131 g_free (vals);
2135 * Note: this geometric distribution is the only with support {0,1,2,...}
2137 static void
2138 test_random_randgeom (int N)
2140 gnm_float mean, var, skew, kurt;
2141 gnm_float *vals;
2142 gboolean ok;
2143 gnm_float param_p = random_01 ();
2144 gnm_float mean_target = (1 - param_p) / param_p;
2145 gnm_float var_target = (1 - param_p) / (param_p * param_p);
2146 gnm_float skew_target = (2 - param_p) / gnm_sqrt (1 - param_p);
2147 gnm_float kurt_target = 6 + (param_p * param_p) / (1 - param_p);
2148 char *expr;
2149 gnm_float T;
2150 int i;
2151 gnm_float fractiles[10], probs[10];
2152 const int nf = G_N_ELEMENTS (fractiles);
2154 expr = g_strdup_printf ("=RANDGEOM(%.10" GNM_FORMAT_g ")", param_p);
2155 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2156 g_free (expr);
2158 ok = TRUE;
2159 for (i = 0; i < N; i++) {
2160 gnm_float r = vals[i];
2161 if (!(r >= 0 && gnm_finite (r) && r == gnm_floor (r))) {
2162 g_printerr ("Range failure.\n");
2163 ok = FALSE;
2164 break;
2168 T = mean_target;
2169 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2170 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2171 g_printerr ("Mean failure.\n");
2172 ok = FALSE;
2175 T = var_target;
2176 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2177 if (!(var >= 0 && gnm_finite (var))) {
2178 /* That is a very simplistic test! */
2179 g_printerr ("Var failure.\n");
2180 ok = FALSE;
2183 T = skew_target;
2184 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2185 if (!gnm_finite (skew)) {
2186 /* That is a very simplistic test! */
2187 g_printerr ("Skew failure.\n");
2188 ok = FALSE;
2191 T = kurt_target;
2192 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2193 if (!(kurt >= -3 && gnm_finite (kurt))) {
2194 /* That is a very simplistic test! */
2195 g_printerr ("Kurt failure.\n");
2196 ok = FALSE;
2199 /* Fractile test */
2200 for (i = 1; i < nf; i++) {
2201 fractiles[i] = qgeom (i / (double)nf, param_p, TRUE, FALSE);
2202 probs[i] = pgeom (fractiles[i], param_p, TRUE, FALSE);
2204 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
2205 ok = FALSE;
2207 if (ok)
2208 g_printerr ("OK\n");
2209 else
2210 add_random_fail ("RANDGEOM");
2211 g_printerr ("\n");
2213 g_free (vals);
2216 static void
2217 test_random_randlog (int N)
2219 gnm_float mean, var, skew, kurt;
2220 gnm_float *vals;
2221 gboolean ok;
2222 gnm_float param_p = random_01 ();
2223 gnm_float p = param_p;
2224 gnm_float l1mp = gnm_log1p (-p);
2225 gnm_float mean_target = -p / (1 - p) / l1mp;
2226 gnm_float var_target = -(p * (p + l1mp)) / gnm_pow ((1 - p) * l1mp, 2);
2227 /* See http://mathworld.wolfram.com/Log-SeriesDistribution.html */
2228 gnm_float skew_target =
2229 -l1mp *
2230 (2 * p * p + 3 * p * l1mp + (1 + p) * l1mp * l1mp) /
2231 (l1mp * (p + l1mp) * gnm_sqrt (-p * (p + l1mp)));
2232 gnm_float kurt_target =
2233 -(6 * p * p * p +
2234 12 * p * p * l1mp +
2235 p * (7 + 4 * p) * l1mp * l1mp +
2236 (1 + 4 * p + p * p) * l1mp * l1mp * l1mp) /
2237 (p * gnm_pow (p + l1mp, 2));
2238 char *expr;
2239 gnm_float T;
2240 int i;
2242 expr = g_strdup_printf ("=RANDLOG(%.10" GNM_FORMAT_g ")", param_p);
2243 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2244 g_free (expr);
2246 ok = TRUE;
2247 for (i = 0; i < N; i++) {
2248 gnm_float r = vals[i];
2249 if (!(r >= 1 && gnm_finite (r) && r == gnm_floor (r))) {
2250 g_printerr ("Range failure.\n");
2251 ok = FALSE;
2252 break;
2256 T = mean_target;
2257 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2258 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2259 g_printerr ("Mean failure.\n");
2260 ok = FALSE;
2263 T = var_target;
2264 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2265 if (!(var >= 0 && gnm_finite (var))) {
2266 /* That is a very simplistic test! */
2267 g_printerr ("Var failure.\n");
2268 ok = FALSE;
2271 T = skew_target;
2272 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2273 if (!gnm_finite (skew)) {
2274 /* That is a very simplistic test! */
2275 g_printerr ("Skew failure.\n");
2276 ok = FALSE;
2279 T = kurt_target;
2280 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2281 if (!(kurt >= -3 && gnm_finite (kurt))) {
2282 /* That is a very simplistic test! */
2283 g_printerr ("Kurt failure.\n");
2284 ok = FALSE;
2287 if (ok)
2288 g_printerr ("OK\n");
2289 else
2290 add_random_fail ("RANDLOG");
2291 g_printerr ("\n");
2293 g_free (vals);
2296 static void
2297 test_random_randweibull (int N)
2299 gnm_float mean, var, skew, kurt;
2300 gnm_float *vals;
2301 gboolean ok;
2302 gnm_float shape = 1 / (0.0001 + gnm_pow (random_01 () / 2, 2));
2303 gnm_float scale = 2 * random_01 ();
2304 gnm_float mean_target = scale * gnm_gamma (1 + 1 / shape);
2305 gnm_float var_target = scale * scale *
2306 (gnm_gamma (1 + 2 / shape) -
2307 gnm_pow (gnm_gamma (1 + 1 / shape), 2));
2308 /* See https://en.wikipedia.org/wiki/Weibull_distribution */
2309 gnm_float skew_target =
2310 (gnm_gamma (1 + 3 / shape) * gnm_pow (scale, 3) -
2311 3 * mean_target * var_target -
2312 gnm_pow (mean_target, 3)) /
2313 gnm_pow (var_target, 1.5);
2314 gnm_float kurt_target = gnm_nan; /* Complicated */
2315 char *expr;
2316 gnm_float T;
2317 int i;
2318 gnm_float fractiles[10];
2319 const int nf = G_N_ELEMENTS (fractiles);
2321 expr = g_strdup_printf ("=RANDWEIBULL(%.10" GNM_FORMAT_f ",%.10" GNM_FORMAT_f ")", scale, shape);
2322 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2323 g_free (expr);
2325 ok = TRUE;
2326 for (i = 0; i < N; i++) {
2327 gnm_float r = vals[i];
2328 if (!(r >= 0 && gnm_finite (r))) {
2329 g_printerr ("Range failure.\n");
2330 ok = FALSE;
2331 break;
2335 T = mean_target;
2336 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2337 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2338 g_printerr ("Mean failure.\n");
2339 ok = FALSE;
2342 T = var_target;
2343 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2344 if (!(var >= 0 && gnm_finite (var))) {
2345 /* That is a very simplistic test! */
2346 g_printerr ("Var failure.\n");
2347 ok = FALSE;
2350 T = skew_target;
2351 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2352 if (!gnm_finite (skew)) {
2353 /* That is a very simplistic test! */
2354 g_printerr ("Skew failure.\n");
2355 ok = FALSE;
2358 T = kurt_target;
2359 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2360 if (!(kurt >= -3 && gnm_finite (kurt))) {
2361 /* That is a very simplistic test! */
2362 g_printerr ("Kurt failure.\n");
2363 ok = FALSE;
2366 /* Fractile test */
2367 for (i = 1; i < nf; i++)
2368 fractiles[i] = qweibull (i / (double)nf, shape, scale, TRUE, FALSE);
2369 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
2370 ok = FALSE;
2372 if (ok)
2373 g_printerr ("OK\n");
2374 else
2375 add_random_fail ("RANDWEIBULL");
2376 g_printerr ("\n");
2378 g_free (vals);
2381 static void
2382 test_random_randlognorm (int N)
2384 gnm_float mean, var, skew, kurt;
2385 gnm_float *vals;
2386 gboolean ok;
2387 gnm_float lm = (random_01() - 0.5) / (0.1 + gnm_pow (random_01 () / 2, 2));
2388 gnm_float ls = 1 / (1 + gnm_pow (random_01 () / 2, 2));
2389 gnm_float mean_target = gnm_exp (lm + ls * ls / 2);
2390 gnm_float var_target = gnm_expm1 (ls * ls) * (mean_target * mean_target);
2391 /* See https://en.wikipedia.org/wiki/Log-normal_distribution */
2392 gnm_float skew_target = (gnm_exp (ls * ls) + 2) *
2393 gnm_sqrt (gnm_expm1 (ls * ls));
2394 gnm_float kurt_target = gnm_nan; /* Complicated */
2395 char *expr;
2396 gnm_float T;
2397 int i;
2398 gnm_float fractiles[10];
2399 const int nf = G_N_ELEMENTS (fractiles);
2401 expr = g_strdup_printf ("=RANDLOGNORM(%.10" GNM_FORMAT_f ",%.10" GNM_FORMAT_f ")", lm, ls);
2402 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2403 g_free (expr);
2405 ok = TRUE;
2406 for (i = 0; i < N; i++) {
2407 gnm_float r = vals[i];
2408 if (!(r >= 0 && r <= gnm_pinf)) {
2409 g_printerr ("Range failure.\n");
2410 ok = FALSE;
2411 break;
2415 T = mean_target;
2416 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2417 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2418 g_printerr ("Mean failure.\n");
2419 ok = FALSE;
2422 T = var_target;
2423 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2424 if (!(var >= 0 && gnm_finite (var))) {
2425 /* That is a very simplistic test! */
2426 g_printerr ("Var failure.\n");
2427 ok = FALSE;
2430 T = skew_target;
2431 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2432 if (!gnm_finite (skew)) {
2433 /* That is a very simplistic test! */
2434 g_printerr ("Skew failure.\n");
2435 ok = FALSE;
2438 T = kurt_target;
2439 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2440 if (!(kurt >= -3 && gnm_finite (kurt))) {
2441 /* That is a very simplistic test! */
2442 g_printerr ("Kurt failure.\n");
2443 ok = FALSE;
2446 /* Fractile test */
2447 for (i = 1; i < nf; i++)
2448 fractiles[i] = qlnorm (i / (double)nf, lm, ls, TRUE, FALSE);
2449 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
2450 ok = FALSE;
2452 if (ok)
2453 g_printerr ("OK\n");
2454 else
2455 add_random_fail ("RANDLOGNORM");
2456 g_printerr ("\n");
2458 g_free (vals);
2461 static void
2462 test_random_randrayleigh (int N)
2464 gnm_float mean, var, skew, kurt;
2465 gnm_float *vals;
2466 gboolean ok;
2467 gnm_float ls = 1 / (1 + gnm_pow (random_01 () / 2, 2));
2468 gnm_float mean_target = ls * gnm_sqrt (M_PIgnum / 2);
2469 gnm_float var_target = (4 - M_PIgnum) / 2 * ls * ls;
2470 gnm_float skew_target = 2 * gnm_sqrt (M_PIgnum) * (M_PIgnum - 3) /
2471 gnm_pow (4 - M_PIgnum, 1.5);
2472 gnm_float kurt_target = gnm_nan; /* Complicated */
2473 char *expr;
2474 gnm_float T;
2475 int i;
2476 gnm_float fractiles[10];
2477 const int nf = G_N_ELEMENTS (fractiles);
2479 expr = g_strdup_printf ("=RANDRAYLEIGH(%.10" GNM_FORMAT_f ")", ls);
2480 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2481 g_free (expr);
2483 ok = TRUE;
2484 for (i = 0; i < N; i++) {
2485 gnm_float r = vals[i];
2486 if (!(r >= 0 && gnm_finite (r))) {
2487 g_printerr ("Range failure.\n");
2488 ok = FALSE;
2489 break;
2493 T = mean_target;
2494 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2495 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2496 g_printerr ("Mean failure.\n");
2497 ok = FALSE;
2500 T = var_target;
2501 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2502 if (!(var >= 0 && gnm_finite (var))) {
2503 /* That is a very simplistic test! */
2504 g_printerr ("Var failure.\n");
2505 ok = FALSE;
2508 T = skew_target;
2509 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2510 if (!gnm_finite (skew)) {
2511 /* That is a very simplistic test! */
2512 g_printerr ("Skew failure.\n");
2513 ok = FALSE;
2516 T = kurt_target;
2517 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2518 if (!(kurt >= -3 && gnm_finite (kurt))) {
2519 /* That is a very simplistic test! */
2520 g_printerr ("Kurt failure.\n");
2521 ok = FALSE;
2524 /* Fractile test */
2525 for (i = 1; i < nf; i++)
2526 fractiles[i] = qrayleigh (i / (double)nf, ls, TRUE, FALSE);
2527 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
2528 ok = FALSE;
2530 if (ok)
2531 g_printerr ("OK\n");
2532 else
2533 add_random_fail ("RANDRAYLEIGH");
2534 g_printerr ("\n");
2536 g_free (vals);
2539 static void
2540 test_random (void)
2542 const char *test_name = "test_random";
2543 const int N = sstest_fast ? 2000 : 20000;
2544 const int High_N = N * 10;
2545 const char *single = g_getenv ("SSTEST_RANDOM");
2547 mark_test_start (test_name);
2549 #define CHECK1(NAME,C) \
2550 do { if (!single || strcmp(single,#NAME) == 0) test_random_ ## NAME (C); } while (0)
2552 /* Continuous */
2553 CHECK1 (rand, N);
2554 CHECK1 (randuniform, N);
2555 CHECK1 (randbeta, N);
2556 CHECK1 (randcauchy, N);
2557 CHECK1 (randchisq, N);
2558 CHECK1 (randexp, N);
2559 CHECK1 (randfdist, N);
2560 CHECK1 (randgamma, N);
2561 CHECK1 (randlog, N);
2562 CHECK1 (randlognorm, N);
2563 CHECK1 (randnorm, High_N);
2564 CHECK1 (randsnorm, High_N);
2565 CHECK1 (randtdist, N);
2566 CHECK1 (randweibull, N);
2567 CHECK1 (randrayleigh, N);
2568 #if 0
2569 CHECK1 (randexppow, N);
2570 CHECK1 (randgumbel, N);
2571 CHECK1 (randlandau, N);
2572 CHECK1 (randlaplace, N);
2573 CHECK1 (randlevy, N);
2574 CHECK1 (randlogistic, N);
2575 CHECK1 (randnormtail, N);
2576 CHECK1 (randpareto, N);
2577 CHECK1 (randrayleightail, N);
2578 CHECK1 (randstdist, N);
2579 #endif
2581 /* Discrete */
2582 CHECK1 (randbernoulli, N);
2583 CHECK1 (randbetween, N);
2584 CHECK1 (randbinom, N);
2585 CHECK1 (randdiscrete, N);
2586 CHECK1 (randgeom, High_N);
2587 CHECK1 (randhyperg, High_N);
2588 CHECK1 (randnegbinom, High_N);
2589 CHECK1 (randpoisson, High_N);
2591 #undef CHECK1
2593 if (!single) {
2594 if (random_summary)
2595 g_printerr ("SUMMARY: FAIL for %s\n\n", random_summary);
2596 else
2597 g_printerr ("SUMMARY: OK\n\n");
2599 g_free (random_summary);
2600 random_summary = NULL;
2602 mark_test_end (test_name);
2605 /*-------------------------------------------------------------------------- */
2607 #define MAYBE_DO(name) if (strcmp (testname, "all") != 0 && strcmp (testname, (name)) != 0) { } else
2610 main (int argc, char const **argv)
2612 GOErrorInfo *plugin_errs;
2613 GOCmdContext *cc;
2614 GOptionContext *ocontext;
2615 GError *error = NULL;
2616 const char *testname;
2618 /* No code before here, we need to init threads */
2619 argv = gnm_pre_parse_init (argc, argv);
2621 ocontext = g_option_context_new (_("[testname]"));
2622 g_option_context_add_main_entries (ocontext, sstest_options, GETTEXT_PACKAGE);
2623 g_option_context_add_group (ocontext, gnm_get_option_group ());
2624 g_option_context_parse (ocontext, &argc, (gchar ***)&argv, &error);
2625 g_option_context_free (ocontext);
2627 if (error) {
2628 g_printerr (_("%s\nRun '%s --help' to see a full list of available command line options.\n"),
2629 error->message, g_get_prgname ());
2630 g_error_free (error);
2631 return 1;
2634 if (sstest_show_version) {
2635 g_printerr (_("version '%s'\ndatadir := '%s'\nlibdir := '%s'\n"),
2636 GNM_VERSION_FULL, gnm_sys_data_dir (), gnm_sys_lib_dir ());
2637 return 0;
2640 gnm_init ();
2642 cc = gnm_cmd_context_stderr_new ();
2643 gnm_plugins_init (GO_CMD_CONTEXT (cc));
2644 go_plugin_db_activate_plugin_list (
2645 go_plugins_get_available_plugins (), &plugin_errs);
2646 if (plugin_errs) {
2647 /* FIXME: What do we want to do here? */
2648 go_error_info_free (plugin_errs);
2651 testname = argv[1];
2652 if (!testname) testname = "all";
2654 /* ---------------------------------------- */
2656 MAYBE_DO ("test_insdel_rowcol_names") test_insdel_rowcol_names ();
2657 MAYBE_DO ("test_insert_delete") test_insert_delete ();
2658 MAYBE_DO ("test_func_help") test_func_help ();
2659 MAYBE_DO ("test_nonascii_numbers") test_nonascii_numbers ();
2660 MAYBE_DO ("test_random") test_random ();
2662 /* ---------------------------------------- */
2664 g_object_unref (cc);
2665 gnm_shutdown ();
2666 gnm_pre_parse_shutdown ();
2668 return 0;