?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
sstest: move handling of --ext-refs-file to here from main binary.
[gnumeric.git] / src / sstest.c
blob346fa8f5abe3b683c2c1175cb768dda0fb8219c5
1 /*
2 * sstest.c: Test code for Gnumeric
3 *
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"
30
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>
36
37 static gboolean sstest_show_version = FALSE;
38 static gboolean sstest_fast = FALSE;
39 static gchar *ext_refs_file = NULL;
40 static gchar *samples_file = NULL;
41
42 static GOptionEntry const sstest_options [] = {
43 {
44 "fast", 'f',
45 0, G_OPTION_ARG_NONE, &sstest_fast,
46 N_("Run fewer iterations"),
47 NULL
48 },
49
50 {
51 "ext-refs-file", 0,
52 0, G_OPTION_ARG_FILENAME, &ext_refs_file,
53 N_("Dumps web page for function help"),
54 N_("FILE")
55 },
56
57 {
58 "samples-file", 0,
59 0, G_OPTION_ARG_FILENAME, &samples_file,
60 N_("Dumps list of samples in function help"),
61 N_("FILE")
62 },
63
64 {
65 "version", 'V',
66 0, G_OPTION_ARG_NONE, &sstest_show_version,
67 N_("Display program version"),
68 NULL
69 },
70
71 { NULL }
72 };
73
74 static void
75 mark_test_start (const char *name)
76 {
77 g_printerr ("-----------------------------------------------------------------------------\nStart: %s\n-----------------------------------------------------------------------------\n\n", name);
78 }
79
80 static void
81 mark_test_end (const char *name)
82 {
83 g_printerr ("End: %s\n\n", name);
84 }
85
86 static void
87 cb_collect_names (G_GNUC_UNUSED const char *name, GnmNamedExpr *nexpr, GSList **names)
88 {
89 *names = g_slist_prepend (*names, nexpr);
90 }
91
92 static GnmCell *
93 fetch_cell (Sheet *sheet, const char *where)
94 {
95 GnmCellPos cp;
96 gboolean ok = cellpos_parse (where,
97 gnm_sheet_get_size (sheet),
98 &cp, TRUE) != NULL;
99 g_return_val_if_fail (ok, NULL);
100 return sheet_cell_fetch (sheet, cp.col, cp.row);
101 }
102
103 static void
104 set_cell (Sheet *sheet, const char *where, const char *what)
105 {
106 GnmCell *cell = fetch_cell (sheet, where);
107 if (cell)
108 gnm_cell_set_text (cell, what);
109 }
110
111 static void
112 dump_sheet (Sheet *sheet, const char *header)
113 {
114 GPtrArray *cells = sheet_cells (sheet, NULL);
115 unsigned ui;
116
117 if (header)
118 g_printerr ("# %s\n", header);
119 for (ui = 0; ui < cells->len; ui++) {
120 GnmCell *cell = g_ptr_array_index (cells, ui);
121 char *txt = gnm_cell_get_entered_text (cell);
122 g_printerr ("%s: %s\n",
123 cellpos_as_string (&cell->pos), txt);
124 g_free (txt);
125 }
126 g_ptr_array_free (cells, TRUE);
127 }
128
129
130 static void
131 dump_names (Workbook *wb)
132 {
133 GSList *l, *names = NULL;
134
135 workbook_foreach_name (wb, FALSE, (GHFunc)cb_collect_names, &names);
136 names = g_slist_sort (names, (GCompareFunc)expr_name_cmp_by_name);
137
138 g_printerr ("Dumping names...\n");
139 for (l = names; l; l = l->next) {
140 GnmNamedExpr *nexpr = l->data;
141 GnmConventionsOut out;
142
143 out.accum = g_string_new (NULL);
144 out.pp = &nexpr->pos;
145 out.convs = gnm_conventions_default;
146
147 g_string_append (out.accum, "Scope=");
148 if (out.pp->sheet)
149 g_string_append (out.accum, out.pp->sheet->name_quoted);
150 else
151 g_string_append (out.accum, "Global");
152
153 g_string_append (out.accum, " Name=");
154 go_strescape (out.accum, expr_name_name (nexpr));
155
156 g_string_append (out.accum, " Expr=");
157 gnm_expr_top_as_gstring (nexpr->texpr, &out);
158
159 g_printerr ("%s\n", out.accum->str);
160 g_string_free (out.accum, TRUE);
161 }
162 g_printerr ("Dumping names... Done\n");
163
164 g_slist_free (names);
165 }
166
167 static void
168 define_name (const char *name, const char *expr_txt, gpointer scope)
169 {
170 GnmParsePos pos;
171 GnmExprTop const *texpr;
172 GnmNamedExpr const *nexpr;
173 GnmConventions const *convs;
174
175 if (IS_SHEET (scope)) {
176 parse_pos_init_sheet (&pos, scope);
177 convs = sheet_get_conventions (pos.sheet);
178 } else {
179 parse_pos_init (&pos, WORKBOOK (scope), NULL, 0, 0);
180 convs = gnm_conventions_default;
181 }
182
183 texpr = gnm_expr_parse_str (expr_txt, &pos,
184 GNM_EXPR_PARSE_DEFAULT,
185 convs, NULL);
186 if (!texpr) {
187 g_printerr ("Failed to parse %s for name %s\n",
188 expr_txt, name);
189 return;
190 }
191
192 nexpr = expr_name_add (&pos, name, texpr, NULL, TRUE, NULL);
193 if (!nexpr)
194 g_printerr ("Failed to add name %s\n", name);
195 }
196
197 static void
198 test_insdel_rowcol_names (void)
199 {
200 Workbook *wb;
201 Sheet *sheet1,*sheet2;
202 const char *test_name = "test_insdel_rowcol_names";
203 GOUndo *undo;
204 int i;
205
206 mark_test_start (test_name);
207
208 wb = workbook_new ();
209 sheet1 = workbook_sheet_add (wb, -1,
210 GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
211 sheet2 = workbook_sheet_add (wb, -1,
212 GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
213
214 define_name ("Print_Area", "Sheet1!$A$1:$IV$65536", sheet1);
215 define_name ("Print_Area", "Sheet2!$A$1:$IV$65536", sheet2);
216
217 define_name ("NAMEGA1", "A1", wb);
218 define_name ("NAMEG2", "$A$14+Sheet1!$A$14+Sheet2!$A$14", wb);
219
220 define_name ("NAMEA1", "A1", sheet1);
221 define_name ("NAMEA2", "A2", sheet1);
222 define_name ("NAMEA1ABS", "$A$1", sheet1);
223 define_name ("NAMEA2ABS", "$A$2", sheet1);
224
225 dump_names (wb);
226
227 for (i = 3; i >= 0; i--) {
228 g_printerr ("About to insert before column %s on %s\n",
229 col_name (i), sheet1->name_unquoted);
230 sheet_insert_cols (sheet1, i, 12, &undo, NULL);
231 dump_names (wb);
232 g_printerr ("Undoing.\n");
233 go_undo_undo_with_data (undo, NULL);
234 g_object_unref (undo);
235 g_printerr ("Done.\n");
236 }
237
238 for (i = 3; i >= 0; i--) {
239 g_printerr ("About to insert before column %s on %s\n",
240 col_name (i), sheet2->name_unquoted);
241 sheet_insert_cols (sheet2, i, 12, &undo, NULL);
242 dump_names (wb);
243 g_printerr ("Undoing.\n");
244 go_undo_undo_with_data (undo, NULL);
245 g_object_unref (undo);
246 g_printerr ("Done.\n");
247 }
248
249 for (i = 3; i >= 0; i--) {
250 g_printerr ("About to delete column %s on %s\n",
251 col_name (i), sheet1->name_unquoted);
252 sheet_delete_cols (sheet1, i, 1, &undo, NULL);
253 dump_names (wb);
254 g_printerr ("Undoing.\n");
255 go_undo_undo_with_data (undo, NULL);
256 g_object_unref (undo);
257 g_printerr ("Done.\n");
258 }
259
260 g_object_unref (wb);
261
262 mark_test_end (test_name);
263 }
264
265 /*-------------------------------------------------------------------------- */
266
267 static void
268 test_insert_delete (void)
269 {
270 const char *test_name = "test_insert_delete";
271 Workbook *wb;
272 Sheet *sheet1;
273 int i;
274 GOUndo *u = NULL, *u1;
275
276 mark_test_start (test_name);
277
278 wb = workbook_new ();
279 sheet1 = workbook_sheet_add (wb, -1,
280 GNM_DEFAULT_COLS, GNM_DEFAULT_ROWS);
281 set_cell (sheet1, "B2", "=D4+1");
282 set_cell (sheet1, "D2", "=if(TRUE,B2,2)");
283
284 dump_sheet (sheet1, "Init");
285
286 for (i = 5; i >= 0; i--) {
287 g_printerr ("# About to insert column before %s\n",
288 col_name (i));
289 sheet_insert_cols (sheet1, i, 1, &u1, NULL);
290 u = go_undo_combine (u, u1);
291 dump_sheet (sheet1, NULL);
292 }
293
294 for (i = 5; i >= 0; i--) {
295 g_printerr ("# About to insert row before %s\n",
296 row_name (i));
297 sheet_insert_rows (sheet1, i, 1, &u1, NULL);
298 u = go_undo_combine (u, u1);
299 dump_sheet (sheet1, NULL);
300 }
301
302 go_undo_undo (u);
303 g_object_unref (u);
304 u = NULL;
305 dump_sheet (sheet1, "Undo the lot");
306
307 for (i = 5; i >= 0; i--) {
308 g_printerr ("# About to delete column %s\n",
309 col_name (i));
310 sheet_delete_cols (sheet1, i, 1, &u1, NULL);
311 u = go_undo_combine (u, u1);
312 dump_sheet (sheet1, NULL);
313 }
314
315 for (i = 5; i >= 0; i--) {
316 g_printerr ("# About to delete row %s\n",
317 row_name (i));
318 sheet_delete_rows (sheet1, i, 1, &u1, NULL);
319 u = go_undo_combine (u, u1);
320 dump_sheet (sheet1, NULL);
321 }
322
323 go_undo_undo (u);
324 g_object_unref (u);
325 u = NULL;
326 dump_sheet (sheet1, "Undo the lot");
327
328 g_object_unref (wb);
329
330 mark_test_end (test_name);
331 }
332
333 /*-------------------------------------------------------------------------- */
334
335 static void
336 test_func_help (void)
337 {
338 const char *test_name = "test_func_help";
339 int res;
340
341 mark_test_start (test_name);
342
343 res = gnm_func_sanity_check ();
344 g_printerr ("Result = %d\n", res);
345
346 mark_test_end (test_name);
347 }
348
349 /*-------------------------------------------------------------------------- */
350
351 static int
352 test_strtol_ok (const char *s, long l, size_t expected_len)
353 {
354 long l2;
355 char *end;
356 int save_errno;
357
358 l2 = gnm_utf8_strtol (s, &end);
359 save_errno = errno;
360
361 if (end != s + expected_len) {
362 g_printerr ("Unexpect conversion end of [%s]\n", s);
363 return 1;
364 }
365 if (l != l2) {
366 g_printerr ("Unexpect conversion result of [%s]\n", s);
367 return 1;
368 }
369 if (save_errno != 0) {
370 g_printerr ("Unexpect conversion errno of [%s]\n", s);
371 return 1;
372 }
373
374 return 0;
375 }
376
377 static int
378 test_strtol_noconv (const char *s)
379 {
380 long l;
381 char *end;
382 int save_errno;
383
384 l = gnm_utf8_strtol (s, &end);
385 save_errno = errno;
386
387 if (end != s) {
388 g_printerr ("Unexpect conversion end of [%s]\n", s);
389 return 1;
390 }
391 if (l != 0) {
392 g_printerr ("Unexpect conversion result of [%s]\n", s);
393 return 1;
394 }
395 if (save_errno != 0) {
396 g_printerr ("Unexpect conversion errno of [%s]\n", s);
397 return 1;
398 }
399
400 return 0;
401 }
402
403 static int
404 test_strtol_overflow (const char *s, gboolean pos)
405 {
406 long l;
407 char *end;
408 int save_errno;
409 size_t expected_len = strlen (s);
410
411 l = gnm_utf8_strtol (s, &end);
412 save_errno = errno;
413
414 if (end != s + expected_len) {
415 g_printerr ("Unexpect conversion end of [%s]\n", s);
416 return 1;
417 }
418 if (l != (pos ? LONG_MAX : LONG_MIN)) {
419 g_printerr ("Unexpect conversion result of [%s]\n", s);
420 return 1;
421 }
422 if (save_errno != ERANGE) {
423 g_printerr ("Unexpect conversion errno of [%s]\n", s);
424 return 1;
425 }
426
427 return 0;
428 }
429
430 static int
431 test_strtol_reverse (long l)
432 {
433 char buffer[4*sizeof(l) + 4];
434 int res = 0;
435
436 sprintf(buffer, "%ld", l);
437 res |= test_strtol_ok (buffer, l, strlen (buffer));
438
439 sprintf(buffer, " %ld", l);
440 res |= test_strtol_ok (buffer, l, strlen (buffer));
441
442 sprintf(buffer, "\xc2\xa0\n\t%ld", l);
443 res |= test_strtol_ok (buffer, l, strlen (buffer));
444
445 sprintf(buffer, " \t%ldx", l);
446 res |= test_strtol_ok (buffer, l, strlen (buffer) - 1);
447
448 return res;
449 }
450
451 static int
452 test_strtod_ok (const char *s, double d, size_t expected_len)
453 {
454 gnm_float d2;
455 char *end;
456 int save_errno;
457
458 d2 = gnm_utf8_strto (s, &end);
459 save_errno = errno;
460
461 if (end != s + expected_len) {
462 g_printerr ("Unexpect conversion end of [%s]\n", s);
463 return 1;
464 }
465 if (d != d2) {
466 g_printerr ("Unexpect conversion result of [%s]\n", s);
467 return 1;
468 }
469 if (save_errno != 0) {
470 g_printerr ("Unexpect conversion errno of [%s]\n", s);
471 return 1;
472 }
473
474 return 0;
475 }
476
477 static void
478 test_nonascii_numbers (void)
479 {
480 const char *test_name = "test_nonascii_numbers";
481 int res = 0;
482
483 mark_test_start (test_name);
484
485 res |= test_strtol_reverse (0);
486 res |= test_strtol_reverse (1);
487 res |= test_strtol_reverse (-1);
488 res |= test_strtol_reverse (LONG_MIN);
489 res |= test_strtol_reverse (LONG_MIN + 1);
490 res |= test_strtol_reverse (LONG_MAX - 1);
491
492 res |= test_strtol_ok ("\xef\xbc\x8d\xef\xbc\x91", -1, 6);
493 res |= test_strtol_ok ("\xc2\xa0+1", 1, 4);
494
495 res |= test_strtol_ok ("000000000000000000000000000000", 0, 30);
496
497 res |= test_strtol_noconv ("");
498 res |= test_strtol_noconv (" ");
499 res |= test_strtol_noconv (" +");
500 res |= test_strtol_noconv (" -");
501 res |= test_strtol_noconv (" .00");
502 res |= test_strtol_noconv (" e0");
503 res |= test_strtol_noconv ("--0");
504 res |= test_strtol_noconv ("+-0");
505 res |= test_strtol_noconv ("+ 0");
506 res |= test_strtol_noconv ("- 0");
507
508 {
509 char buffer[4 * sizeof (long) + 2];
510
511 sprintf (buffer, "-%lu", 1 + (unsigned long)LONG_MIN);
512 res |= test_strtol_overflow (buffer, FALSE);
513 sprintf (buffer, "-%lu", 10 + (unsigned long)LONG_MIN);
514 res |= test_strtol_overflow (buffer, FALSE);
515
516 sprintf (buffer, "%lu", 1 + (unsigned long)LONG_MAX);
517 res |= test_strtol_overflow (buffer, TRUE);
518 sprintf (buffer, "%lu", 10 + (unsigned long)LONG_MAX);
519 res |= test_strtol_overflow (buffer, TRUE);
520 }
521
522 /* -------------------- */
523
524 res |= test_strtod_ok ("0", 0, 1);
525 res |= test_strtod_ok ("1", 1, 1);
526 res |= test_strtod_ok ("-1", -1, 2);
527 res |= test_strtod_ok ("+1", 1, 2);
528 res |= test_strtod_ok (" +1", 1, 3);
529 res |= test_strtod_ok ("\xc2\xa0+1", 1, 4);
530 res |= test_strtod_ok ("\xc2\xa0+1x", 1, 4);
531 res |= test_strtod_ok ("\xc2\xa0+1e", 1, 4);
532 res |= test_strtod_ok ("\xc2\xa0+1e+", 1, 4);
533 res |= test_strtod_ok ("\xc2\xa0+1e+0", 1, 7);
534 res |= test_strtod_ok ("-1e1", -10, 4);
535 res |= test_strtod_ok ("100e-2", 1, 6);
536 res |= test_strtod_ok ("100e+2", 10000, 6);
537 res |= test_strtod_ok ("1x0p0", 1, 1);
538 res |= test_strtod_ok ("+inf", gnm_pinf, 4);
539 res |= test_strtod_ok ("-inf", gnm_ninf, 4);
540 res |= test_strtod_ok ("1.25", 1.25, 4);
541 res |= test_strtod_ok ("1.25e1", 12.5, 6);
542 res |= test_strtod_ok ("12.5e-1", 1.25, 7);
543
544 g_printerr ("Result = %d\n", res);
545
546 mark_test_end (test_name);
547 }
548
549 /*-------------------------------------------------------------------------- */
550
551 static char *random_summary = NULL;
552
553 static void
554 add_random_fail (const char *s)
555 {
556 if (random_summary) {
557 char *t = g_strconcat (random_summary, ", ", s, NULL);
558 g_free (random_summary);
559 random_summary = t;
560 } else
561 random_summary = g_strdup (s);
562 }
563
564 static void
565 define_cell (Sheet *sheet, int c, int r, const char *expr)
566 {
567 GnmCell *cell = sheet_cell_fetch (sheet, c, r);
568 sheet_cell_set_text (cell, expr, NULL);
569 }
570
571 #define GET_PROB(i_) ((i_) <= 0 ? 0 : ((i_) >= nf ? 1 : probs[(i_)]))
572
573 static gboolean
574 rand_fractile_test (gnm_float const *vals, int N, int nf,
575 gnm_float const *fractiles, gnm_float const *probs)
576 {
577 gnm_float f = 1.0 / nf;
578 int *fractilecount = g_new (int, nf + 1);
579 int *expected = g_new (int, nf + 1);
580 int i;
581 gboolean ok = TRUE;
582 gboolean debug = TRUE;
583
584 if (debug) {
585 g_printerr ("Bin upper limit:");
586 for (i = 1; i <= nf; i++) {
587 gnm_float U = (i == nf) ? gnm_pinf : fractiles[i];
588 g_printerr ("%s%" GNM_FORMAT_g,
589 (i == 1) ? " " : ", ",
590 U);
591 }
592 g_printerr (".\n");
593 }
594
595 if (debug && probs) {
596 g_printerr ("Cumulative probabilities:");
597 for (i = 1; i <= nf; i++)
598 g_printerr ("%s%.1" GNM_FORMAT_f "%%",
599 (i == 1) ? " " : ", ", 100 * GET_PROB (i));
600 g_printerr (".\n");
601 }
602
603 for (i = 1; i < nf - 1; i++) {
604 if (!(fractiles[i] <= fractiles[i + 1])) {
605 g_printerr ("Severe fractile ordering problem.\n");
606 return FALSE;
607 }
608
609 if (probs && !(probs[i] <= probs[i + 1])) {
610 g_printerr ("Severe cumulative probabilities ordering problem.\n");
611 return FALSE;
612 }
613 }
614 if (probs && (probs[1] < 0 || probs[nf - 1] > 1)) {
615 g_printerr ("Severe cumulative probabilities range problem.\n");
616 return FALSE;
617 }
618
619 for (i = 0; i <= nf; i++)
620 fractilecount[i] = 0;
621
622 for (i = 0; i < N; i++) {
623 gnm_float r = vals[i];
624 int j;
625 for (j = 1; j < nf; j++)
626 if (r <= fractiles[j])
627 break;
628 fractilecount[j]++;
629 }
630 g_printerr ("Fractile counts:");
631 for (i = 1; i <= nf; i++)
632 g_printerr ("%s%d", (i == 1) ? " " : ", ", fractilecount[i]);
633 g_printerr (".\n");
634
635 if (probs) {
636 g_printerr ("Expected counts:");
637 for (i = 1; i <= nf; i++) {
638 gnm_float p = GET_PROB (i) - GET_PROB (i-1);
639 expected[i] = gnm_floor (p * N + 0.5);
640 g_printerr ("%s%d", (i == 1) ? " " : ", ", expected[i]);
641 }
642 g_printerr (".\n");
643 } else {
644 gnm_float T = f * N;
645 g_printerr ("Expected count in each fractile: %.10" GNM_FORMAT_g "\n", T);
646 for (i = 0; i <= nf; i++)
647 expected[i] = T;
648 }
649
650 for (i = 1; i <= nf; i++) {
651 gnm_float T = expected[i];
652 if (!(gnm_abs (fractilecount[i] - T) <= 4 * gnm_sqrt (T))) {
653 g_printerr ("Fractile test failure for bin %d.\n", i);
654 ok = FALSE;
655 }
656 }
657
658 g_free (fractilecount);
659 g_free (expected);
660
661 return ok;
662 }
663
664 #undef GET_PROB
665
666 static gnm_float *
667 test_random_1 (int N, const char *expr,
668 gnm_float *mean, gnm_float *var,
669 gnm_float *skew, gnm_float *kurt)
670 {
671 Workbook *wb = workbook_new ();
672 Sheet *sheet;
673 gnm_float *res = g_new (gnm_float, N);
674 int i;
675 char *s;
676 int cols = 2, rows = N;
677
678 g_printerr ("Testing %s\n", expr);
679
680 gnm_sheet_suggest_size (&cols, &rows);
681 sheet = workbook_sheet_add (wb, -1, cols, rows);
682
683 for (i = 0; i < N; i++)
684 define_cell (sheet, 0, i, expr);
685
686 s = g_strdup_printf ("=average(a1:a%d)", N);
687 define_cell (sheet, 1, 0, s);
688 g_free (s);
689
690 s = g_strdup_printf ("=var(a1:a%d)", N);
691 define_cell (sheet, 1, 1, s);
692 g_free (s);
693
694 s = g_strdup_printf ("=skew(a1:a%d)", N);
695 define_cell (sheet, 1, 2, s);
696 g_free (s);
697
698 s = g_strdup_printf ("=kurt(a1:a%d)", N);
699 define_cell (sheet, 1, 3, s);
700 g_free (s);
701
702 /* Force recalc of all dirty cells even in manual mode. */
703 workbook_recalc (sheet->workbook);
704
705 for (i = 0; i < N; i++)
706 res[i] = value_get_as_float (sheet_cell_get (sheet, 0, i)->value);
707 *mean = value_get_as_float (sheet_cell_get (sheet, 1, 0)->value);
708 g_printerr ("Mean: %.10" GNM_FORMAT_g "\n", *mean);
709
710 *var = value_get_as_float (sheet_cell_get (sheet, 1, 1)->value);
711 g_printerr ("Var: %.10" GNM_FORMAT_g "\n", *var);
712
713 *skew = value_get_as_float (sheet_cell_get (sheet, 1, 2)->value);
714 g_printerr ("Skew: %.10" GNM_FORMAT_g "\n", *skew);
715
716 *kurt = value_get_as_float (sheet_cell_get (sheet, 1, 3)->value);
717 g_printerr ("Kurt: %.10" GNM_FORMAT_g "\n", *kurt);
718
719 g_object_unref (wb);
720 return res;
721 }
722
723 static gnm_float *
724 test_random_normality (int N, const char *expr,
725 gnm_float *mean, gnm_float *var,
726 gnm_float *adtest, gnm_float *cvmtest,
727 gnm_float *lkstest, gnm_float *sftest)
728 {
729 Workbook *wb = workbook_new ();
730 Sheet *sheet;
731 gnm_float *res = g_new (gnm_float, N);
732 int i;
733 char *s;
734 int cols = 2, rows = N;
735
736 g_printerr ("Testing %s\n", expr);
737
738 gnm_sheet_suggest_size (&cols, &rows);
739 sheet = workbook_sheet_add (wb, -1, cols, rows);
740
741 for (i = 0; i < N; i++)
742 define_cell (sheet, 0, i, expr);
743
744 s = g_strdup_printf ("=average(a1:a%d)", N);
745 define_cell (sheet, 1, 0, s);
746 g_free (s);
747
748 s = g_strdup_printf ("=var(a1:a%d)", N);
749 define_cell (sheet, 1, 1, s);
750 g_free (s);
751
752 s = g_strdup_printf ("=adtest(a1:a%d)", N);
753 define_cell (sheet, 1, 2, s);
754 g_free (s);
755
756 s = g_strdup_printf ("=cvmtest(a1:a%d)", N);
757 define_cell (sheet, 1, 3, s);
758 g_free (s);
759
760 s = g_strdup_printf ("=lkstest(a1:a%d)", N);
761 define_cell (sheet, 1, 4, s);
762 g_free (s);
763
764 s = g_strdup_printf ("=sftest(a1:a%d)", N > 5000 ? 5000 : N);
765 define_cell (sheet, 1, 5, s);
766 g_free (s);
767
768 /* Force recalc of all dirty cells even in manual mode. */
769 workbook_recalc (sheet->workbook);
770
771 for (i = 0; i < N; i++)
772 res[i] = value_get_as_float (sheet_cell_get (sheet, 0, i)->value);
773 *mean = value_get_as_float (sheet_cell_get (sheet, 1, 0)->value);
774 g_printerr ("Mean: %.10" GNM_FORMAT_g "\n", *mean);
775
776 *var = value_get_as_float (sheet_cell_get (sheet, 1, 1)->value);
777 g_printerr ("Var: %.10" GNM_FORMAT_g "\n", *var);
778
779 *adtest = value_get_as_float (sheet_cell_get (sheet, 1, 2)->value);
780 g_printerr ("ADTest: %.10" GNM_FORMAT_g "\n", *adtest);
781
782 *cvmtest = value_get_as_float (sheet_cell_get (sheet, 1, 3)->value);
783 g_printerr ("CVMTest: %.10" GNM_FORMAT_g "\n", *cvmtest);
784
785 *lkstest = value_get_as_float (sheet_cell_get (sheet, 1, 4)->value);
786 g_printerr ("LKSTest: %.10" GNM_FORMAT_g "\n", *lkstest);
787
788 *sftest = value_get_as_float (sheet_cell_get (sheet, 1, 5)->value);
789 g_printerr ("SFTest: %.10" GNM_FORMAT_g "\n", *sftest);
790
791 g_object_unref (wb);
792 return res;
793 }
794
795 static void
796 test_random_rand (int N)
797 {
798 gnm_float mean, var, skew, kurt;
799 gnm_float mean_target = 0.5;
800 gnm_float var_target = 1.0 / 12;
801 gnm_float skew_target = 0;
802 gnm_float kurt_target = -6.0 / 5;
803 gnm_float *vals;
804 int i;
805 gboolean ok;
806 gnm_float T;
807 gnm_float fractiles[10];
808 const int nf = G_N_ELEMENTS (fractiles);
809
810 vals = test_random_1 (N, "=RAND()", &mean, &var, &skew, &kurt);
811 ok = TRUE;
812 for (i = 0; i < N; i++) {
813 gnm_float r = vals[i];
814 if (!(r >= 0 && r < 1)) {
815 g_printerr ("Range failure.\n");
816 ok = FALSE;
817 break;
818 }
819 }
820
821 T = mean_target;
822 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
823 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
824 ok = FALSE;
825 }
826 T = var_target;
827 if (gnm_abs (var - T) > 0.01) {
828 g_printerr ("Var failure.\n");
829 ok = FALSE;
830 }
831 T = skew_target;
832 if (gnm_abs (skew - T) > 0.05) {
833 g_printerr ("Skew failure.\n");
834 ok = FALSE;
835 }
836 T = kurt_target;
837 if (gnm_abs (kurt - T) > 0.05) {
838 g_printerr ("Kurt failure.\n");
839 ok = FALSE;
840 }
841
842 /* Fractile test */
843 for (i = 1; i < nf; i++)
844 fractiles[i] = i / (double)nf;
845 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
846 ok = FALSE;
847
848 if (ok)
849 g_printerr ("OK\n");
850 else
851 add_random_fail ("RAND");
852 g_printerr ("\n");
853
854 g_free (vals);
855 }
856
857 static void
858 test_random_randuniform (int N)
859 {
860 gnm_float mean, var, skew, kurt;
861 gnm_float *vals;
862 gboolean ok;
863 gnm_float lsign = (random_01 () > 0.75 ? 1 : -1);
864 gnm_float param_l = lsign * gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
865 gnm_float param_h = param_l + gnm_floor (1 / (0.0001 + gnm_pow (random_01 () / 2, 4)));
866 gnm_float n = param_h - param_l;
867 gnm_float mean_target = (param_l + param_h) / 2;
868 gnm_float var_target = (n * n) / 12;
869 gnm_float skew_target = 0;
870 gnm_float kurt_target = -6 / 5.0;
871 char *expr;
872 gnm_float T;
873 int i;
874 gnm_float fractiles[10];
875 const int nf = G_N_ELEMENTS (fractiles);
876
877 expr = g_strdup_printf ("=RANDUNIFORM(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")", param_l, param_h);
878 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
879 g_free (expr);
880
881 ok = TRUE;
882 for (i = 0; i < N; i++) {
883 gnm_float r = vals[i];
884 if (!(r >= param_l && r < param_h)) {
885 g_printerr ("Range failure.\n");
886 ok = FALSE;
887 break;
888 }
889 }
890
891 T = mean_target;
892 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
893 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
894 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
895 ok = FALSE;
896 }
897
898 T = var_target;
899 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
900 if (!(var >= 0 && gnm_finite (var))) {
901 /* That is a very simplistic test! */
902 g_printerr ("Var failure.\n");
903 ok = FALSE;
904 }
905
906 T = skew_target;
907 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
908 if (!gnm_finite (skew)) {
909 /* That is a very simplistic test! */
910 g_printerr ("Skew failure.\n");
911 ok = FALSE;
912 }
913
914 T = kurt_target;
915 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
916 if (!(kurt >= -3 && gnm_finite (kurt))) {
917 /* That is a very simplistic test! */
918 g_printerr ("Kurt failure.\n");
919 ok = FALSE;
920 }
921
922 /* Fractile test */
923 for (i = 1; i < nf; i++)
924 fractiles[i] = param_l + n * i / (double)nf;
925 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
926 ok = FALSE;
927
928 if (ok)
929 g_printerr ("OK\n");
930 else
931 add_random_fail ("RANDUNIFORM");
932 g_printerr ("\n");
933
934 g_free (vals);
935 }
936
937 static void
938 test_random_randbernoulli (int N)
939 {
940 gnm_float p = 0.3;
941 gnm_float q = 1 - p;
942 gnm_float mean, var, skew, kurt;
943 gnm_float mean_target = p;
944 gnm_float var_target = p * (1 - p);
945 gnm_float skew_target = (q - p) / gnm_sqrt (p * q);
946 gnm_float kurt_target = (1 - 6 * p * q) / (p * q);
947 gnm_float *vals;
948 int i;
949 gboolean ok;
950 char *expr;
951 gnm_float T;
952
953 expr = g_strdup_printf ("=RANDBERNOULLI(%.10" GNM_FORMAT_g ")", p);
954 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
955 g_free (expr);
956
957 ok = TRUE;
958 for (i = 0; i < N; i++) {
959 gnm_float r = vals[i];
960 if (!(r == 0 || r == 1)) {
961 g_printerr ("Range failure.\n");
962 ok = FALSE;
963 break;
964 }
965 }
966
967 T = mean_target;
968 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
969 if (gnm_abs (mean - p) > 0.01) {
970 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
971 ok = FALSE;
972 }
973
974 T = var_target;
975 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
976 if (gnm_abs (var - T) > 0.01) {
977 g_printerr ("Var failure.\n");
978 ok = FALSE;
979 }
980
981 T = skew_target;
982 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
983 if (!(gnm_abs (skew - T) <= 0.10 * gnm_abs (T))) {
984 g_printerr ("Skew failure.\n");
985 ok = FALSE;
986 }
987
988 T = kurt_target;
989 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
990 if (!(gnm_abs (kurt - T) <= 0.15 * gnm_abs (T))) {
991 g_printerr ("Kurt failure.\n");
992 ok = FALSE;
993 }
994 if (ok)
995 g_printerr ("OK\n");
996 else
997 add_random_fail ("RANDBERNOULLI");
998 g_printerr ("\n");
999
1000 g_free (vals);
1001 }
1002
1003 static void
1004 test_random_randdiscrete (int N)
1005 {
1006 gnm_float mean, var, skew, kurt;
1007 gnm_float *vals;
1008 int i;
1009 gboolean ok;
1010 gnm_float mean_target = 13;
1011 gnm_float var_target = 156;
1012 gnm_float skew_target = 0.6748;
1013 gnm_float kurt_target = -0.9057;
1014 char *expr;
1015 gnm_float T;
1016
1017 expr = g_strdup_printf ("=RANDDISCRETE({0;1;4;9;16;25;36})");
1018 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1019 g_free (expr);
1020
1021 ok = TRUE;
1022 for (i = 0; i < N; i++) {
1023 gnm_float r = vals[i];
1024 if (!(r >= 0 && r <= 36 && gnm_sqrt (r) == gnm_floor (gnm_sqrt (r)))) {
1025 g_printerr ("Range failure.\n");
1026 ok = FALSE;
1027 break;
1028 }
1029 }
1030
1031 T = mean_target;
1032 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1033 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1034 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1035 ok = FALSE;
1036 }
1037
1038 T = var_target;
1039 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1040 if (!(var >= 0 && gnm_finite (var))) {
1041 /* That is a very simplistic test! */
1042 g_printerr ("Var failure.\n");
1043 ok = FALSE;
1044 }
1045
1046 T = skew_target;
1047 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1048 if (!gnm_finite (skew)) {
1049 /* That is a very simplistic test! */
1050 g_printerr ("Skew failure.\n");
1051 ok = FALSE;
1052 }
1053
1054 T = kurt_target;
1055 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1056 if (!(kurt >= -3 && gnm_finite (kurt))) {
1057 /* That is a very simplistic test! */
1058 g_printerr ("Kurt failure.\n");
1059 ok = FALSE;
1060 }
1061
1062 if (ok)
1063 g_printerr ("OK\n");
1064 else
1065 add_random_fail ("RANDDISCRETE");
1066 g_printerr ("\n");
1067
1068 g_free (vals);
1069 }
1070
1071 static void
1072 test_random_randnorm (int N)
1073 {
1074 gnm_float mean, var, adtest, cvmtest, lkstest, sftest;
1075 gnm_float mean_target = 0, var_target = 1;
1076 gnm_float *vals;
1077 gboolean ok;
1078 char *expr;
1079 gnm_float T;
1080 int i;
1081 gnm_float fractiles[10];
1082 const int nf = G_N_ELEMENTS (fractiles);
1083
1084 expr = g_strdup_printf ("=RANDNORM(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")",
1085 mean_target, var_target);
1086 vals = test_random_normality (N, expr, &mean, &var, &adtest, &cvmtest, &lkstest, &sftest);
1087 g_free (expr);
1088
1089 ok = TRUE;
1090 for (i = 0; i < N; i++) {
1091 gnm_float r = vals[i];
1092 if (!gnm_finite (r)) {
1093 g_printerr ("Range failure.\n");
1094 ok = FALSE;
1095 break;
1096 }
1097 }
1098
1099 T = mean_target;
1100 if (gnm_abs (mean - T) > 0.02) {
1101 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1102 ok = FALSE;
1103 }
1104 T = var_target;
1105 if (gnm_abs (var - T) > 0.02) {
1106 g_printerr ("Var failure.\n");
1107 ok = FALSE;
1108 }
1109
1110 /* Fractile test */
1111 for (i = 1; i < nf; i++)
1112 fractiles[i] = qnorm (i / (double)nf, mean_target, gnm_sqrt (var_target), TRUE, FALSE);
1113 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1114 ok = FALSE;
1115
1116 if (adtest < 0.01) {
1117 g_printerr ("Anderson Darling Test rejected [%.10" GNM_FORMAT_g "]\n", adtest);
1118 ok = FALSE;
1119 }
1120 if (cvmtest < 0.01) {
1121 g_printerr ("Cramér-von Mises Test rejected [%.10" GNM_FORMAT_g "]\n", cvmtest);
1122 ok = FALSE;
1123 }
1124 if (lkstest < 0.01) {
1125 g_printerr ("Lilliefors (Kolmogorov-Smirnov) Test rejected [%.10" GNM_FORMAT_g "]\n",
1126 lkstest);
1127 ok = FALSE;
1128 }
1129 if (sftest < 0.01) {
1130 g_printerr ("Shapiro-Francia Test rejected [%.10" GNM_FORMAT_g "]\n", sftest);
1131 ok = FALSE;
1132 }
1133
1134 if (ok)
1135 g_printerr ("OK\n");
1136 else
1137 add_random_fail ("RANDNORM");
1138 g_printerr ("\n");
1139
1140 g_free (vals);
1141 }
1142
1143 static void
1144 test_random_randsnorm (int N)
1145 {
1146 gnm_float mean, var, skew, kurt;
1147 gnm_float *vals;
1148 gboolean ok;
1149 gnm_float alpha = 5;
1150 gnm_float delta = alpha/gnm_sqrt(1+alpha*alpha);
1151 gnm_float mean_target = delta * gnm_sqrt (2/M_PIgnum);
1152 gnm_float var_target = 1-mean_target*mean_target;
1153 char *expr;
1154 gnm_float T;
1155 int i;
1156
1157 expr = g_strdup_printf ("=RANDSNORM(%.10" GNM_FORMAT_g ")", alpha);
1158 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1159 g_free (expr);
1160
1161 ok = TRUE;
1162 for (i = 0; i < N; i++) {
1163 gnm_float r = vals[i];
1164 if (!gnm_finite (r)) {
1165 g_printerr ("Range failure.\n");
1166 ok = FALSE;
1167 break;
1168 }
1169 }
1170
1171 T = mean_target;
1172 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1173 if (gnm_abs (mean - T) > 0.01) {
1174 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1175 ok = FALSE;
1176 }
1177
1178 T = var_target;
1179 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1180 if (gnm_abs (var - T) > 0.01) {
1181 g_printerr ("Var failure.\n");
1182 ok = FALSE;
1183 }
1184
1185 T = mean_target/gnm_sqrt(var_target);
1186 T = T*T*T*(4-M_PIgnum)/2;
1187 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1188 if (gnm_abs (skew - T) > 0.05) {
1189 g_printerr ("Skew failure.\n");
1190 ok = FALSE;
1191 }
1192
1193 T = 2*(M_PIgnum - 3)*mean_target*mean_target*mean_target*mean_target/(var_target*var_target);
1194 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1195 if (gnm_abs (kurt - T) > 0.15) {
1196 g_printerr ("Kurt failure.\n");
1197 ok = FALSE;
1198 }
1199
1200 if (ok)
1201 g_printerr ("OK\n");
1202 else
1203 add_random_fail ("RANDSNORM");
1204 g_printerr ("\n");
1205
1206 g_free (vals);
1207 }
1208
1209 static void
1210 test_random_randexp (int N)
1211 {
1212 gnm_float mean, var, skew, kurt;
1213 gnm_float *vals;
1214 gboolean ok;
1215 gnm_float param_l = 1 / (0.0001 + gnm_pow (random_01 () / 2, 4));
1216 gnm_float mean_target = param_l;
1217 gnm_float var_target = mean_target * mean_target;
1218 gnm_float skew_target = 2;
1219 gnm_float kurt_target = 6;
1220 char *expr;
1221 gnm_float T;
1222 int i;
1223 gnm_float fractiles[10];
1224 const int nf = G_N_ELEMENTS (fractiles);
1225
1226 expr = g_strdup_printf ("=RANDEXP(%.10" GNM_FORMAT_g ")", param_l);
1227 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1228 g_free (expr);
1229
1230 ok = TRUE;
1231 for (i = 0; i < N; i++) {
1232 gnm_float r = vals[i];
1233 if (!(r >= 0 && gnm_finite (r))) {
1234 g_printerr ("Range failure.\n");
1235 ok = FALSE;
1236 break;
1237 }
1238 }
1239
1240 T = mean_target;
1241 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1242 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1243 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1244 ok = FALSE;
1245 }
1246
1247 T = var_target;
1248 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1249 if (!(var >= 0 && gnm_finite (var))) {
1250 /* That is a very simplistic test! */
1251 g_printerr ("Var failure.\n");
1252 ok = FALSE;
1253 }
1254
1255 T = skew_target;
1256 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1257 if (!gnm_finite (skew)) {
1258 /* That is a very simplistic test! */
1259 g_printerr ("Skew failure.\n");
1260 ok = FALSE;
1261 }
1262
1263 T = kurt_target;
1264 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1265 if (!(kurt >= -3 && gnm_finite (kurt))) {
1266 /* That is a very simplistic test! */
1267 g_printerr ("Kurt failure.\n");
1268 ok = FALSE;
1269 }
1270
1271 /* Fractile test */
1272 for (i = 1; i < nf; i++)
1273 fractiles[i] = qexp (i / (double)nf, param_l, TRUE, FALSE);
1274 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1275 ok = FALSE;
1276
1277 if (ok)
1278 g_printerr ("OK\n");
1279 else
1280 add_random_fail ("RANDEXP");
1281 g_printerr ("\n");
1282
1283 g_free (vals);
1284 }
1285
1286 static void
1287 test_random_randgamma (int N)
1288 {
1289 gnm_float mean, var, skew, kurt;
1290 gnm_float *vals;
1291 gboolean ok;
1292 gnm_float param_shape = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 6)));
1293 gnm_float param_scale = 0.001 + gnm_pow (random_01 (), 4) * 1000;
1294 gnm_float mean_target = param_shape * param_scale;
1295 gnm_float var_target = mean_target * param_scale;
1296 gnm_float skew_target = 2 / gnm_sqrt (param_shape);
1297 gnm_float kurt_target = 6 / param_shape;
1298 char *expr;
1299 gnm_float T;
1300 int i;
1301 gnm_float fractiles[10];
1302 const int nf = G_N_ELEMENTS (fractiles);
1303
1304 expr = g_strdup_printf ("=RANDGAMMA(%.0" GNM_FORMAT_f ",%.10" GNM_FORMAT_g ")", param_shape, param_scale);
1305 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1306 g_free (expr);
1307
1308 ok = TRUE;
1309 for (i = 0; i < N; i++) {
1310 gnm_float r = vals[i];
1311 if (!(r > 0 && gnm_finite (r))) {
1312 g_printerr ("Range failure.\n");
1313 ok = FALSE;
1314 break;
1315 }
1316 }
1317
1318 T = mean_target;
1319 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1320 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1321 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1322 ok = FALSE;
1323 }
1324
1325 T = var_target;
1326 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1327 if (!(var >= 0 && gnm_finite (var))) {
1328 /* That is a very simplistic test! */
1329 g_printerr ("Var failure.\n");
1330 ok = FALSE;
1331 }
1332
1333 T = skew_target;
1334 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1335 if (!gnm_finite (skew)) {
1336 /* That is a very simplistic test! */
1337 g_printerr ("Skew failure.\n");
1338 ok = FALSE;
1339 }
1340
1341 T = kurt_target;
1342 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1343 if (!(kurt >= -3 && gnm_finite (kurt))) {
1344 /* That is a very simplistic test! */
1345 g_printerr ("Kurt failure.\n");
1346 ok = FALSE;
1347 }
1348
1349 /* Fractile test */
1350 for (i = 1; i < nf; i++)
1351 fractiles[i] = qgamma (i / (double)nf, param_shape, param_scale, TRUE, FALSE);
1352 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1353 ok = FALSE;
1354
1355 if (ok)
1356 g_printerr ("OK\n");
1357 else
1358 add_random_fail ("RANDGAMMA");
1359 g_printerr ("\n");
1360
1361 g_free (vals);
1362 }
1363
1364 static void
1365 test_random_randbeta (int N)
1366 {
1367 gnm_float mean, var, skew, kurt;
1368 gnm_float *vals;
1369 gboolean ok;
1370 gnm_float param_a = 1 / (0.0001 + gnm_pow (random_01 (), 6));
1371 gnm_float param_b = 1 / (0.0001 + gnm_pow (random_01 (), 6));
1372 gnm_float s = param_a + param_b;
1373 gnm_float mean_target = param_a / s;
1374 gnm_float var_target = mean_target * param_b / (s * (s + 1));
1375 gnm_float skew_target =
1376 (2 * (param_b - param_a) * gnm_sqrt (s + 1))/
1377 ((s + 2) * gnm_sqrt (param_a * param_b));
1378 gnm_float kurt_target = gnm_nan; /* Complicated */
1379 char *expr;
1380 gnm_float T;
1381 int i;
1382 gnm_float fractiles[10];
1383 const int nf = G_N_ELEMENTS (fractiles);
1384
1385 expr = g_strdup_printf ("=RANDBETA(%.10" GNM_FORMAT_g ",%.10" GNM_FORMAT_g ")", param_a, param_b);
1386 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1387 g_free (expr);
1388
1389 ok = TRUE;
1390 for (i = 0; i < N; i++) {
1391 gnm_float r = vals[i];
1392 if (!(r >= 0 && r <= 1)) {
1393 g_printerr ("Range failure.\n");
1394 ok = FALSE;
1395 break;
1396 }
1397 }
1398
1399 T = mean_target;
1400 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1401 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1402 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1403 ok = FALSE;
1404 }
1405
1406 T = var_target;
1407 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1408 if (!(var >= 0 && gnm_finite (var))) {
1409 /* That is a very simplistic test! */
1410 g_printerr ("Var failure.\n");
1411 ok = FALSE;
1412 }
1413
1414 T = skew_target;
1415 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1416 if (!gnm_finite (skew)) {
1417 /* That is a very simplistic test! */
1418 g_printerr ("Skew failure.\n");
1419 ok = FALSE;
1420 }
1421
1422 T = kurt_target;
1423 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1424 if (!(kurt >= -3 && gnm_finite (kurt))) {
1425 /* That is a very simplistic test! */
1426 g_printerr ("Kurt failure.\n");
1427 ok = FALSE;
1428 }
1429
1430 /* Fractile test */
1431 for (i = 1; i < nf; i++)
1432 fractiles[i] = qbeta (i / (double)nf, param_a, param_b, TRUE, FALSE);
1433 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1434 ok = FALSE;
1435
1436 if (ok)
1437 g_printerr ("OK\n");
1438 else
1439 add_random_fail ("RANDBETA");
1440 g_printerr ("\n");
1441
1442 g_free (vals);
1443 }
1444
1445 static void
1446 test_random_randtdist (int N)
1447 {
1448 gnm_float mean, var, skew, kurt;
1449 gnm_float *vals;
1450 gboolean ok;
1451 gnm_float param_df = 1 + gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1452 gnm_float mean_target = 0;
1453 gnm_float var_target = param_df > 2 ? param_df / (param_df - 2) : gnm_nan;
1454 gnm_float skew_target = param_df > 3 ? 0 : gnm_nan;
1455 gnm_float kurt_target = param_df > 4 ? 6 / (param_df - 4) : gnm_nan;
1456 char *expr;
1457 gnm_float T;
1458 int i;
1459 gnm_float fractiles[10];
1460 const int nf = G_N_ELEMENTS (fractiles);
1461
1462 expr = g_strdup_printf ("=RANDTDIST(%.0" GNM_FORMAT_f ")", param_df);
1463 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1464 g_free (expr);
1465
1466 ok = TRUE;
1467 for (i = 0; i < N; i++) {
1468 gnm_float r = vals[i];
1469 if (!(gnm_finite (r))) {
1470 g_printerr ("Range failure.\n");
1471 ok = FALSE;
1472 break;
1473 }
1474 }
1475
1476 T = mean_target;
1477 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1478 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1479 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1480 ok = FALSE;
1481 }
1482
1483 T = var_target;
1484 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1485 if (!(var >= 0 && gnm_finite (var))) {
1486 /* That is a very simplistic test! */
1487 g_printerr ("Var failure.\n");
1488 ok = FALSE;
1489 }
1490
1491 T = skew_target;
1492 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1493 if (!gnm_finite (skew)) {
1494 /* That is a very simplistic test! */
1495 g_printerr ("Skew failure.\n");
1496 ok = FALSE;
1497 }
1498
1499 T = kurt_target;
1500 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1501 if (!(kurt >= -3 && gnm_finite (kurt))) {
1502 /* That is a very simplistic test! */
1503 g_printerr ("Kurt failure.\n");
1504 ok = FALSE;
1505 }
1506
1507 /* Fractile test */
1508 for (i = 1; i < nf; i++)
1509 fractiles[i] = qt (i / (double)nf, param_df, TRUE, FALSE);
1510 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1511 ok = FALSE;
1512
1513 if (ok)
1514 g_printerr ("OK\n");
1515 else
1516 add_random_fail ("RANDTDIST");
1517 g_printerr ("\n");
1518
1519 g_free (vals);
1520 }
1521
1522 static void
1523 test_random_randfdist (int N)
1524 {
1525 gnm_float mean, var, skew, kurt;
1526 gnm_float *vals;
1527 gboolean ok;
1528 gnm_float param_df1 = 1 + gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1529 gnm_float param_df2 = 1 + gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1530 gnm_float mean_target = param_df2 > 2 ? param_df2 / (param_df2 - 2) : gnm_nan;
1531 gnm_float var_target = param_df2 > 4
1532 ? (2 * param_df2 * param_df2 * (param_df1 + param_df2 - 2) /
1533 (param_df1 * (param_df2 - 2) * (param_df2 - 2) * (param_df2 - 4)))
1534 : gnm_nan;
1535 gnm_float skew_target = gnm_nan; /* Complicated */
1536 gnm_float kurt_target = gnm_nan; /* Complicated */
1537 char *expr;
1538 gnm_float T;
1539 int i;
1540 gnm_float fractiles[10];
1541 const int nf = G_N_ELEMENTS (fractiles);
1542
1543 expr = g_strdup_printf ("=RANDFDIST(%.0" GNM_FORMAT_f ",%.0" GNM_FORMAT_f ")", param_df1, param_df2);
1544 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1545 g_free (expr);
1546
1547 ok = TRUE;
1548 for (i = 0; i < N; i++) {
1549 gnm_float r = vals[i];
1550 if (!(r >= 0 && gnm_finite (r))) {
1551 g_printerr ("Range failure.\n");
1552 ok = FALSE;
1553 break;
1554 }
1555 }
1556
1557 T = mean_target;
1558 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1559 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1560 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1561 ok = FALSE;
1562 }
1563
1564 T = var_target;
1565 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1566 if (!(var >= 0 && gnm_finite (var))) {
1567 /* That is a very simplistic test! */
1568 g_printerr ("Var failure.\n");
1569 ok = FALSE;
1570 }
1571
1572 T = skew_target;
1573 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1574 if (!gnm_finite (skew)) {
1575 /* That is a very simplistic test! */
1576 g_printerr ("Skew failure.\n");
1577 ok = FALSE;
1578 }
1579
1580 T = kurt_target;
1581 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1582 if (!(kurt >= -3 && gnm_finite (kurt))) {
1583 /* That is a very simplistic test! */
1584 g_printerr ("Kurt failure.\n");
1585 ok = FALSE;
1586 }
1587
1588 /* Fractile test */
1589 for (i = 1; i < nf; i++)
1590 fractiles[i] = qf (i / (double)nf, param_df1, param_df2, TRUE, FALSE);
1591 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1592 ok = FALSE;
1593
1594 if (ok)
1595 g_printerr ("OK\n");
1596 else
1597 add_random_fail ("RANDFDIST");
1598 g_printerr ("\n");
1599
1600 g_free (vals);
1601 }
1602
1603 static void
1604 test_random_randchisq (int N)
1605 {
1606 gnm_float mean, var, skew, kurt;
1607 gnm_float *vals;
1608 gboolean ok;
1609 gnm_float param_df = 1 + gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 6)));
1610 gnm_float mean_target = param_df;
1611 gnm_float var_target = param_df * 2;
1612 gnm_float skew_target = gnm_sqrt (8 / param_df);
1613 gnm_float kurt_target = 12 / param_df;
1614 char *expr;
1615 gnm_float T;
1616 int i;
1617 gnm_float fractiles[10];
1618 const int nf = G_N_ELEMENTS (fractiles);
1619
1620 expr = g_strdup_printf ("=RANDCHISQ(%.10" GNM_FORMAT_g ")", param_df);
1621 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1622 g_free (expr);
1623
1624 ok = TRUE;
1625 for (i = 0; i < N; i++) {
1626 gnm_float r = vals[i];
1627 if (!(r >= 0 && gnm_finite (r))) {
1628 g_printerr ("Range failure.\n");
1629 ok = FALSE;
1630 break;
1631 }
1632 }
1633
1634 T = mean_target;
1635 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1636 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1637 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1638 ok = FALSE;
1639 }
1640
1641 T = var_target;
1642 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1643 if (!(var >= 0 && gnm_finite (var))) {
1644 /* That is a very simplistic test! */
1645 g_printerr ("Var failure.\n");
1646 ok = FALSE;
1647 }
1648
1649 T = skew_target;
1650 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1651 if (!gnm_finite (skew)) {
1652 /* That is a very simplistic test! */
1653 g_printerr ("Skew failure.\n");
1654 ok = FALSE;
1655 }
1656
1657 T = kurt_target;
1658 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1659 if (!(kurt >= -3 && gnm_finite (kurt))) {
1660 /* That is a very simplistic test! */
1661 g_printerr ("Kurt failure.\n");
1662 ok = FALSE;
1663 }
1664
1665 /* Fractile test */
1666 for (i = 1; i < nf; i++)
1667 fractiles[i] = qchisq (i / (double)nf, param_df, TRUE, FALSE);
1668 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1669 ok = FALSE;
1670
1671 if (ok)
1672 g_printerr ("OK\n");
1673 else
1674 add_random_fail ("RANDCHISQ");
1675 g_printerr ("\n");
1676
1677 g_free (vals);
1678 }
1679
1680 static void
1681 test_random_randcauchy (int N)
1682 {
1683 gnm_float mean, var, skew, kurt;
1684 gnm_float *vals;
1685 gboolean ok;
1686 gnm_float param_scale = 0.001 + gnm_pow (random_01 (), 4) * 1000;
1687 gnm_float mean_target = gnm_nan;
1688 gnm_float var_target = gnm_nan;
1689 gnm_float skew_target = gnm_nan;
1690 gnm_float kurt_target = gnm_nan;
1691 char *expr;
1692 gnm_float T;
1693 int i;
1694 gnm_float fractiles[10];
1695 const int nf = G_N_ELEMENTS (fractiles);
1696
1697 /*
1698 * The distribution has no mean, no variance, no skew, and no kurtosis.
1699 * The support is all reals.
1700 */
1701
1702 expr = g_strdup_printf ("=RANDCAUCHY(%.10" GNM_FORMAT_g ")", param_scale);
1703 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1704 g_free (expr);
1705
1706 ok = TRUE;
1707 for (i = 0; i < N; i++) {
1708 gnm_float r = vals[i];
1709 if (!(gnm_finite (r))) {
1710 g_printerr ("Range failure.\n");
1711 ok = FALSE;
1712 break;
1713 }
1714 }
1715
1716 T = mean_target;
1717 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1718 if (gnm_finite (var_target) && !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1719 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1720 ok = FALSE;
1721 }
1722
1723 T = var_target;
1724 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1725 if (!(var >= 0 && gnm_finite (var))) {
1726 /* That is a very simplistic test! */
1727 g_printerr ("Var failure.\n");
1728 ok = FALSE;
1729 }
1730
1731 T = skew_target;
1732 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1733 if (!gnm_finite (skew)) {
1734 /* That is a very simplistic test! */
1735 g_printerr ("Skew failure.\n");
1736 ok = FALSE;
1737 }
1738
1739 T = kurt_target;
1740 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1741 if (!(kurt >= -3 && gnm_finite (kurt))) {
1742 /* That is a very simplistic test! */
1743 g_printerr ("Kurt failure.\n");
1744 ok = FALSE;
1745 }
1746
1747 /* Fractile test */
1748 for (i = 1; i < nf; i++)
1749 fractiles[i] = qcauchy (i / (double)nf, 0.0, param_scale, TRUE, FALSE);
1750 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
1751 ok = FALSE;
1752
1753 if (ok)
1754 g_printerr ("OK\n");
1755 else
1756 add_random_fail ("RANDCAUCHY");
1757 g_printerr ("\n");
1758
1759 g_free (vals);
1760 }
1761
1762 static void
1763 test_random_randbinom (int N)
1764 {
1765 gnm_float mean, var, skew, kurt;
1766 gnm_float *vals;
1767 gboolean ok;
1768 gnm_float param_p = random_01 ();
1769 gnm_float param_trials = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
1770 gnm_float mean_target = param_trials * param_p;
1771 gnm_float var_target = mean_target * (1 - param_p);
1772 gnm_float skew_target = (1 - 2 * param_p) / gnm_sqrt (var_target);
1773 gnm_float kurt_target = (1 - 6 * param_p * (1 - param_p)) / var_target;
1774 char *expr;
1775 gnm_float T;
1776 int i;
1777 gnm_float fractiles[10], probs[10];
1778 const int nf = G_N_ELEMENTS (fractiles);
1779
1780 expr = g_strdup_printf ("=RANDBINOM(%.10" GNM_FORMAT_g ",%.0" GNM_FORMAT_f ")", param_p, param_trials);
1781 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1782 g_free (expr);
1783
1784 ok = TRUE;
1785 for (i = 0; i < N; i++) {
1786 gnm_float r = vals[i];
1787 if (!(r >= 0 && r <= param_trials && r == gnm_floor (r))) {
1788 g_printerr ("Range failure.\n");
1789 ok = FALSE;
1790 break;
1791 }
1792 }
1793
1794 T = mean_target;
1795 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1796 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1797 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1798 ok = FALSE;
1799 }
1800
1801 T = var_target;
1802 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1803 if (!(var >= 0 && gnm_finite (var))) {
1804 /* That is a very simplistic test! */
1805 g_printerr ("Var failure.\n");
1806 ok = FALSE;
1807 }
1808
1809 T = skew_target;
1810 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1811 if (!gnm_finite (skew)) {
1812 /* That is a very simplistic test! */
1813 g_printerr ("Skew failure.\n");
1814 ok = FALSE;
1815 }
1816
1817 T = kurt_target;
1818 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1819 if (!(kurt >= -3 && gnm_finite (kurt))) {
1820 /* That is a very simplistic test! */
1821 g_printerr ("Kurt failure.\n");
1822 ok = FALSE;
1823 }
1824
1825 /* Fractile test */
1826 for (i = 1; i < nf; i++) {
1827 fractiles[i] = qbinom (i / (double)nf, param_trials, param_p, TRUE, FALSE);
1828 probs[i] = pbinom (fractiles[i], param_trials, param_p, TRUE, FALSE);
1829 }
1830 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
1831 ok = FALSE;
1832
1833 if (ok)
1834 g_printerr ("OK\n");
1835 else
1836 add_random_fail ("RANDBINOM");
1837 g_printerr ("\n");
1838
1839 g_free (vals);
1840 }
1841
1842 static void
1843 test_random_randnegbinom (int N)
1844 {
1845 gnm_float mean, var, skew, kurt;
1846 gnm_float *vals;
1847 gboolean ok;
1848 gnm_float param_p = random_01 ();
1849 gnm_float param_fails = gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
1850 /* Warning: these differ from Wikipedia by swapping p and 1-p. */
1851 gnm_float mean_target = param_fails * (1 - param_p) / param_p;
1852 gnm_float var_target = mean_target / param_p;
1853 gnm_float skew_target = (2 - param_p) / gnm_sqrt (param_fails * (1 - param_p));
1854 gnm_float kurt_target = 6 / param_fails + 1 / var_target;
1855 char *expr;
1856 gnm_float T;
1857 int i;
1858 gnm_float fractiles[10], probs[10];
1859 const int nf = G_N_ELEMENTS (fractiles);
1860
1861 expr = g_strdup_printf ("=RANDNEGBINOM(%.10" GNM_FORMAT_g ",%.0" GNM_FORMAT_f ")", param_p, param_fails);
1862 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1863 g_free (expr);
1864
1865 ok = TRUE;
1866 for (i = 0; i < N; i++) {
1867 gnm_float r = vals[i];
1868 if (!(r >= 0 && gnm_finite (r) && r == gnm_floor (r))) {
1869 g_printerr ("Range failure.\n");
1870 ok = FALSE;
1871 break;
1872 }
1873 }
1874
1875 T = mean_target;
1876 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1877 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1878 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1879 ok = FALSE;
1880 }
1881
1882 T = var_target;
1883 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1884 if (!(var >= 0 && gnm_finite (var))) {
1885 /* That is a very simplistic test! */
1886 g_printerr ("Var failure.\n");
1887 ok = FALSE;
1888 }
1889
1890 T = skew_target;
1891 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1892 if (!gnm_finite (skew)) {
1893 /* That is a very simplistic test! */
1894 g_printerr ("Skew failure.\n");
1895 ok = FALSE;
1896 }
1897
1898 T = kurt_target;
1899 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1900 if (!(kurt >= -3 && gnm_finite (kurt))) {
1901 /* That is a very simplistic test! */
1902 g_printerr ("Kurt failure.\n");
1903 ok = FALSE;
1904 }
1905
1906 /* Fractile test */
1907 for (i = 1; i < nf; i++) {
1908 fractiles[i] = qnbinom (i / (double)nf, param_fails, param_p, TRUE, FALSE);
1909 probs[i] = pnbinom (fractiles[i], param_fails, param_p, TRUE, FALSE);
1910 }
1911 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
1912 ok = FALSE;
1913
1914 if (ok)
1915 g_printerr ("OK\n");
1916 else
1917 add_random_fail ("RANDNEGBINOM");
1918 g_printerr ("\n");
1919
1920 g_free (vals);
1921 }
1922
1923 static void
1924 test_random_randhyperg (int N)
1925 {
1926 gnm_float mean, var, skew, kurt;
1927 gnm_float *vals;
1928 gboolean ok;
1929 gnm_float param_nr = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 4)));
1930 gnm_float param_nb = gnm_floor (1 / (0.01 + gnm_pow (random_01 (), 4)));
1931 gnm_float s = param_nr + param_nb;
1932 gnm_float param_n = gnm_floor (random_01 () * (s + 1));
1933 gnm_float mean_target = param_n * param_nr / s;
1934 gnm_float var_target = s > 1
1935 ? mean_target * (param_nb / s) * (s - param_n) / (s - 1)
1936 : 0;
1937 gnm_float skew_target = gnm_nan; /* Complicated */
1938 gnm_float kurt_target = gnm_nan; /* Complicated */
1939 char *expr;
1940 gnm_float T;
1941 int i;
1942 gnm_float fractiles[10], probs[10];
1943 const int nf = G_N_ELEMENTS (fractiles);
1944
1945 expr = g_strdup_printf ("=RANDHYPERG(%.10" GNM_FORMAT_g ",%.0" GNM_FORMAT_f ",%.0" GNM_FORMAT_f ")", param_nr, param_nb, param_n);
1946 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
1947 g_free (expr);
1948
1949 ok = TRUE;
1950 for (i = 0; i < N; i++) {
1951 gnm_float r = vals[i];
1952 if (!(r >= 0 && r <= param_n && r == gnm_floor (r))) {
1953 g_printerr ("Range failure.\n");
1954 ok = FALSE;
1955 break;
1956 }
1957 }
1958
1959 T = mean_target;
1960 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
1961 if (gnm_finite (var_target) &&
1962 !(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
1963 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
1964 ok = FALSE;
1965 }
1966
1967 T = var_target;
1968 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
1969 if (!(var >= 0 && gnm_finite (var))) {
1970 /* That is a very simplistic test! */
1971 g_printerr ("Var failure.\n");
1972 ok = FALSE;
1973 }
1974
1975 T = skew_target;
1976 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
1977 if (!gnm_finite (skew)) {
1978 /* That is a very simplistic test! */
1979 g_printerr ("Skew failure.\n");
1980 ok = FALSE;
1981 }
1982
1983 T = kurt_target;
1984 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
1985 if (!(kurt >= -3 && gnm_finite (kurt))) {
1986 /* That is a very simplistic test! */
1987 g_printerr ("Kurt failure.\n");
1988 ok = FALSE;
1989 }
1990
1991 /* Fractile test */
1992 for (i = 1; i < nf; i++) {
1993 fractiles[i] = qhyper (i / (double)nf, param_nr, param_nb, param_n, TRUE, FALSE);
1994 probs[i] = phyper (fractiles[i], param_nr, param_nb, param_n, TRUE, FALSE);
1995 }
1996 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
1997 ok = FALSE;
1998
1999 if (ok)
2000 g_printerr ("OK\n");
2001 else
2002 add_random_fail ("RANDHYPERG");
2003 g_printerr ("\n");
2004
2005 g_free (vals);
2006 }
2007
2008 static void
2009 test_random_randbetween (int N)
2010 {
2011 gnm_float mean, var, skew, kurt;
2012 gnm_float *vals;
2013 gboolean ok;
2014 gnm_float lsign = (random_01 () > 0.75 ? 1 : -1);
2015 gnm_float param_l = lsign * gnm_floor (1 / (0.0001 + gnm_pow (random_01 (), 4)));
2016 gnm_float param_h = param_l + gnm_floor (1 / (0.0001 + gnm_pow (random_01 () / 2, 4)));
2017 gnm_float n = param_h - param_l + 1;
2018 gnm_float mean_target = (param_l + param_h) / 2;
2019 gnm_float var_target = (n * n - 1) / 12;
2020 gnm_float skew_target = 0;
2021 gnm_float kurt_target = (n * n + 1) / (n * n - 1) * -6 / 5;
2022 char *expr;
2023 gnm_float T;
2024 int i;
2025
2026 expr = g_strdup_printf ("=RANDBETWEEN(%.0" GNM_FORMAT_f ",%.0" GNM_FORMAT_f ")", param_l, param_h);
2027 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2028 g_free (expr);
2029
2030 ok = TRUE;
2031 for (i = 0; i < N; i++) {
2032 gnm_float r = vals[i];
2033 if (!(r >= param_l && r <= param_h && r == gnm_floor (r))) {
2034 g_printerr ("Range failure.\n");
2035 ok = FALSE;
2036 break;
2037 }
2038 }
2039
2040 T = mean_target;
2041 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2042 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2043 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2044 ok = FALSE;
2045 }
2046
2047 T = var_target;
2048 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2049 if (!(var >= 0 && gnm_finite (var))) {
2050 /* That is a very simplistic test! */
2051 g_printerr ("Var failure.\n");
2052 ok = FALSE;
2053 }
2054
2055 T = skew_target;
2056 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2057 if (!gnm_finite (skew)) {
2058 /* That is a very simplistic test! */
2059 g_printerr ("Skew failure.\n");
2060 ok = FALSE;
2061 }
2062
2063 T = kurt_target;
2064 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2065 if (!(kurt >= -3 && gnm_finite (kurt))) {
2066 /* That is a very simplistic test! */
2067 g_printerr ("Kurt failure.\n");
2068 ok = FALSE;
2069 }
2070
2071 if (ok)
2072 g_printerr ("OK\n");
2073 else
2074 add_random_fail ("RANDBETWEEN");
2075 g_printerr ("\n");
2076
2077 g_free (vals);
2078 }
2079
2080 static void
2081 test_random_randpoisson (int N)
2082 {
2083 gnm_float mean, var, skew, kurt;
2084 gnm_float *vals;
2085 gboolean ok;
2086 gnm_float param_l = 1 / (0.0001 + gnm_pow (random_01 () / 2, 4));
2087 gnm_float mean_target = param_l;
2088 gnm_float var_target = param_l;
2089 gnm_float skew_target = 1 / gnm_sqrt (param_l);
2090 gnm_float kurt_target = 1 / param_l;
2091 char *expr;
2092 gnm_float T;
2093 int i;
2094 gnm_float fractiles[10], probs[10];
2095 const int nf = G_N_ELEMENTS (fractiles);
2096
2097 expr = g_strdup_printf ("=RANDPOISSON(%.10" GNM_FORMAT_g ")", param_l);
2098 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2099 g_free (expr);
2100
2101 ok = TRUE;
2102 for (i = 0; i < N; i++) {
2103 gnm_float r = vals[i];
2104 if (!(r >= 0 && gnm_finite (r) && r == gnm_floor (r))) {
2105 g_printerr ("Range failure.\n");
2106 ok = FALSE;
2107 break;
2108 }
2109 }
2110
2111 T = mean_target;
2112 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2113 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2114 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2115 ok = FALSE;
2116 }
2117
2118 T = var_target;
2119 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2120 if (!(var >= 0 && gnm_finite (var))) {
2121 /* That is a very simplistic test! */
2122 g_printerr ("Var failure.\n");
2123 ok = FALSE;
2124 }
2125
2126 T = skew_target;
2127 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2128 if (!gnm_finite (skew)) {
2129 /* That is a very simplistic test! */
2130 g_printerr ("Skew failure.\n");
2131 ok = FALSE;
2132 }
2133
2134 T = kurt_target;
2135 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2136 if (!(kurt >= -3 && gnm_finite (kurt))) {
2137 /* That is a very simplistic test! */
2138 g_printerr ("Kurt failure.\n");
2139 ok = FALSE;
2140 }
2141
2142 /* Fractile test */
2143 for (i = 1; i < nf; i++) {
2144 fractiles[i] = qpois (i / (double)nf, param_l, TRUE, FALSE);
2145 probs[i] = ppois (fractiles[i], param_l, TRUE, FALSE);
2146 }
2147 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
2148 ok = FALSE;
2149
2150 if (ok)
2151 g_printerr ("OK\n");
2152 else
2153 add_random_fail ("RANDPOISSON");
2154 g_printerr ("\n");
2155
2156 g_free (vals);
2157 }
2158
2159 /*
2160 * Note: this geometric distribution is the only with support {0,1,2,...}
2161 */
2162 static void
2163 test_random_randgeom (int N)
2164 {
2165 gnm_float mean, var, skew, kurt;
2166 gnm_float *vals;
2167 gboolean ok;
2168 gnm_float param_p = random_01 ();
2169 gnm_float mean_target = (1 - param_p) / param_p;
2170 gnm_float var_target = (1 - param_p) / (param_p * param_p);
2171 gnm_float skew_target = (2 - param_p) / gnm_sqrt (1 - param_p);
2172 gnm_float kurt_target = 6 + (param_p * param_p) / (1 - param_p);
2173 char *expr;
2174 gnm_float T;
2175 int i;
2176 gnm_float fractiles[10], probs[10];
2177 const int nf = G_N_ELEMENTS (fractiles);
2178
2179 expr = g_strdup_printf ("=RANDGEOM(%.10" GNM_FORMAT_g ")", param_p);
2180 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2181 g_free (expr);
2182
2183 ok = TRUE;
2184 for (i = 0; i < N; i++) {
2185 gnm_float r = vals[i];
2186 if (!(r >= 0 && gnm_finite (r) && r == gnm_floor (r))) {
2187 g_printerr ("Range failure.\n");
2188 ok = FALSE;
2189 break;
2190 }
2191 }
2192
2193 T = mean_target;
2194 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2195 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2196 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2197 ok = FALSE;
2198 }
2199
2200 T = var_target;
2201 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2202 if (!(var >= 0 && gnm_finite (var))) {
2203 /* That is a very simplistic test! */
2204 g_printerr ("Var failure.\n");
2205 ok = FALSE;
2206 }
2207
2208 T = skew_target;
2209 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2210 if (!gnm_finite (skew)) {
2211 /* That is a very simplistic test! */
2212 g_printerr ("Skew failure.\n");
2213 ok = FALSE;
2214 }
2215
2216 T = kurt_target;
2217 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2218 if (!(kurt >= -3 && gnm_finite (kurt))) {
2219 /* That is a very simplistic test! */
2220 g_printerr ("Kurt failure.\n");
2221 ok = FALSE;
2222 }
2223
2224 /* Fractile test */
2225 for (i = 1; i < nf; i++) {
2226 fractiles[i] = qgeom (i / (double)nf, param_p, TRUE, FALSE);
2227 probs[i] = pgeom (fractiles[i], param_p, TRUE, FALSE);
2228 }
2229 if (!rand_fractile_test (vals, N, nf, fractiles, probs))
2230 ok = FALSE;
2231
2232 if (ok)
2233 g_printerr ("OK\n");
2234 else
2235 add_random_fail ("RANDGEOM");
2236 g_printerr ("\n");
2237
2238 g_free (vals);
2239 }
2240
2241 static void
2242 test_random_randlog (int N)
2243 {
2244 gnm_float mean, var, skew, kurt;
2245 gnm_float *vals;
2246 gboolean ok;
2247 gnm_float param_p = random_01 ();
2248 gnm_float p = param_p;
2249 gnm_float l1mp = gnm_log1p (-p);
2250 gnm_float mean_target = -p / (1 - p) / l1mp;
2251 gnm_float var_target = -(p * (p + l1mp)) / gnm_pow ((1 - p) * l1mp, 2);
2252 /* See http://mathworld.wolfram.com/Log-SeriesDistribution.html */
2253 gnm_float skew_target =
2254 -l1mp *
2255 (2 * p * p + 3 * p * l1mp + (1 + p) * l1mp * l1mp) /
2256 (l1mp * (p + l1mp) * gnm_sqrt (-p * (p + l1mp)));
2257 gnm_float kurt_target =
2258 -(6 * p * p * p +
2259 12 * p * p * l1mp +
2260 p * (7 + 4 * p) * l1mp * l1mp +
2261 (1 + 4 * p + p * p) * l1mp * l1mp * l1mp) /
2262 (p * gnm_pow (p + l1mp, 2));
2263 char *expr;
2264 gnm_float T;
2265 int i;
2266
2267 expr = g_strdup_printf ("=RANDLOG(%.10" GNM_FORMAT_g ")", param_p);
2268 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2269 g_free (expr);
2270
2271 ok = TRUE;
2272 for (i = 0; i < N; i++) {
2273 gnm_float r = vals[i];
2274 if (!(r >= 1 && gnm_finite (r) && r == gnm_floor (r))) {
2275 g_printerr ("Range failure.\n");
2276 ok = FALSE;
2277 break;
2278 }
2279 }
2280
2281 T = mean_target;
2282 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2283 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2284 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2285 ok = FALSE;
2286 }
2287
2288 T = var_target;
2289 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2290 if (!(var >= 0 && gnm_finite (var))) {
2291 /* That is a very simplistic test! */
2292 g_printerr ("Var failure.\n");
2293 ok = FALSE;
2294 }
2295
2296 T = skew_target;
2297 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2298 if (!gnm_finite (skew)) {
2299 /* That is a very simplistic test! */
2300 g_printerr ("Skew failure.\n");
2301 ok = FALSE;
2302 }
2303
2304 T = kurt_target;
2305 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2306 if (!(kurt >= -3 && gnm_finite (kurt))) {
2307 /* That is a very simplistic test! */
2308 g_printerr ("Kurt failure.\n");
2309 ok = FALSE;
2310 }
2311
2312 if (ok)
2313 g_printerr ("OK\n");
2314 else
2315 add_random_fail ("RANDLOG");
2316 g_printerr ("\n");
2317
2318 g_free (vals);
2319 }
2320
2321 static void
2322 test_random_randweibull (int N)
2323 {
2324 gnm_float mean, var, skew, kurt;
2325 gnm_float *vals;
2326 gboolean ok;
2327 gnm_float shape = 1 / (0.0001 + gnm_pow (random_01 () / 2, 2));
2328 gnm_float scale = 2 * random_01 ();
2329 gnm_float mean_target = scale * gnm_gamma (1 + 1 / shape);
2330 gnm_float var_target = scale * scale *
2331 (gnm_gamma (1 + 2 / shape) -
2332 gnm_pow (gnm_gamma (1 + 1 / shape), 2));
2333 /* See https://en.wikipedia.org/wiki/Weibull_distribution */
2334 gnm_float skew_target =
2335 (gnm_gamma (1 + 3 / shape) * gnm_pow (scale, 3) -
2336 3 * mean_target * var_target -
2337 gnm_pow (mean_target, 3)) /
2338 gnm_pow (var_target, 1.5);
2339 gnm_float kurt_target = gnm_nan; /* Complicated */
2340 char *expr;
2341 gnm_float T;
2342 int i;
2343 gnm_float fractiles[10];
2344 const int nf = G_N_ELEMENTS (fractiles);
2345
2346 expr = g_strdup_printf ("=RANDWEIBULL(%.10" GNM_FORMAT_f ",%.10" GNM_FORMAT_f ")", scale, shape);
2347 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2348 g_free (expr);
2349
2350 ok = TRUE;
2351 for (i = 0; i < N; i++) {
2352 gnm_float r = vals[i];
2353 if (!(r >= 0 && gnm_finite (r))) {
2354 g_printerr ("Range failure.\n");
2355 ok = FALSE;
2356 break;
2357 }
2358 }
2359
2360 T = mean_target;
2361 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2362 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2363 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2364 ok = FALSE;
2365 }
2366
2367 T = var_target;
2368 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2369 if (!(var >= 0 && gnm_finite (var))) {
2370 /* That is a very simplistic test! */
2371 g_printerr ("Var failure.\n");
2372 ok = FALSE;
2373 }
2374
2375 T = skew_target;
2376 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2377 if (!gnm_finite (skew)) {
2378 /* That is a very simplistic test! */
2379 g_printerr ("Skew failure.\n");
2380 ok = FALSE;
2381 }
2382
2383 T = kurt_target;
2384 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2385 if (!(kurt >= -3 && gnm_finite (kurt))) {
2386 /* That is a very simplistic test! */
2387 g_printerr ("Kurt failure.\n");
2388 ok = FALSE;
2389 }
2390
2391 /* Fractile test */
2392 for (i = 1; i < nf; i++)
2393 fractiles[i] = qweibull (i / (double)nf, shape, scale, TRUE, FALSE);
2394 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
2395 ok = FALSE;
2396
2397 if (ok)
2398 g_printerr ("OK\n");
2399 else
2400 add_random_fail ("RANDWEIBULL");
2401 g_printerr ("\n");
2402
2403 g_free (vals);
2404 }
2405
2406 static void
2407 test_random_randlognorm (int N)
2408 {
2409 gnm_float mean, var, skew, kurt;
2410 gnm_float *vals;
2411 gboolean ok;
2412 gnm_float lm = (random_01() - 0.5) / (0.1 + gnm_pow (random_01 () / 2, 2));
2413 gnm_float ls = 1 / (1 + gnm_pow (random_01 () / 2, 2));
2414 gnm_float mean_target = gnm_exp (lm + ls * ls / 2);
2415 gnm_float var_target = gnm_expm1 (ls * ls) * (mean_target * mean_target);
2416 /* See https://en.wikipedia.org/wiki/Log-normal_distribution */
2417 gnm_float skew_target = (gnm_exp (ls * ls) + 2) *
2418 gnm_sqrt (gnm_expm1 (ls * ls));
2419 gnm_float kurt_target = gnm_nan; /* Complicated */
2420 char *expr;
2421 gnm_float T;
2422 int i;
2423 gnm_float fractiles[10];
2424 const int nf = G_N_ELEMENTS (fractiles);
2425
2426 expr = g_strdup_printf ("=RANDLOGNORM(%.10" GNM_FORMAT_f ",%.10" GNM_FORMAT_f ")", lm, ls);
2427 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2428 g_free (expr);
2429
2430 ok = TRUE;
2431 for (i = 0; i < N; i++) {
2432 gnm_float r = vals[i];
2433 if (!(r >= 0 && r <= gnm_pinf)) {
2434 g_printerr ("Range failure.\n");
2435 ok = FALSE;
2436 break;
2437 }
2438 }
2439
2440 T = mean_target;
2441 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2442 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2443 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2444 ok = FALSE;
2445 }
2446
2447 T = var_target;
2448 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2449 if (!(var >= 0 && gnm_finite (var))) {
2450 /* That is a very simplistic test! */
2451 g_printerr ("Var failure.\n");
2452 ok = FALSE;
2453 }
2454
2455 T = skew_target;
2456 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2457 if (!gnm_finite (skew)) {
2458 /* That is a very simplistic test! */
2459 g_printerr ("Skew failure.\n");
2460 ok = FALSE;
2461 }
2462
2463 T = kurt_target;
2464 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2465 if (!(kurt >= -3 && gnm_finite (kurt))) {
2466 /* That is a very simplistic test! */
2467 g_printerr ("Kurt failure.\n");
2468 ok = FALSE;
2469 }
2470
2471 /* Fractile test */
2472 for (i = 1; i < nf; i++)
2473 fractiles[i] = qlnorm (i / (double)nf, lm, ls, TRUE, FALSE);
2474 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
2475 ok = FALSE;
2476
2477 if (ok)
2478 g_printerr ("OK\n");
2479 else
2480 add_random_fail ("RANDLOGNORM");
2481 g_printerr ("\n");
2482
2483 g_free (vals);
2484 }
2485
2486 static void
2487 test_random_randrayleigh (int N)
2488 {
2489 gnm_float mean, var, skew, kurt;
2490 gnm_float *vals;
2491 gboolean ok;
2492 gnm_float ls = 1 / (1 + gnm_pow (random_01 () / 2, 2));
2493 gnm_float mean_target = ls * gnm_sqrt (M_PIgnum / 2);
2494 gnm_float var_target = (4 - M_PIgnum) / 2 * ls * ls;
2495 gnm_float skew_target = 2 * gnm_sqrt (M_PIgnum) * (M_PIgnum - 3) /
2496 gnm_pow (4 - M_PIgnum, 1.5);
2497 gnm_float kurt_target = gnm_nan; /* Complicated */
2498 char *expr;
2499 gnm_float T;
2500 int i;
2501 gnm_float fractiles[10];
2502 const int nf = G_N_ELEMENTS (fractiles);
2503
2504 expr = g_strdup_printf ("=RANDRAYLEIGH(%.10" GNM_FORMAT_f ")", ls);
2505 vals = test_random_1 (N, expr, &mean, &var, &skew, &kurt);
2506 g_free (expr);
2507
2508 ok = TRUE;
2509 for (i = 0; i < N; i++) {
2510 gnm_float r = vals[i];
2511 if (!(r >= 0 && gnm_finite (r))) {
2512 g_printerr ("Range failure.\n");
2513 ok = FALSE;
2514 break;
2515 }
2516 }
2517
2518 T = mean_target;
2519 g_printerr ("Expected mean: %.10" GNM_FORMAT_g "\n", T);
2520 if (!(gnm_abs (mean - T) <= 3 * gnm_sqrt (var_target / N))) {
2521 g_printerr ("Mean failure [%.1" GNM_FORMAT_f " stdev]\n", (mean - T) / gnm_sqrt (var_target / N));
2522 ok = FALSE;
2523 }
2524
2525 T = var_target;
2526 g_printerr ("Expected var: %.10" GNM_FORMAT_g "\n", T);
2527 if (!(var >= 0 && gnm_finite (var))) {
2528 /* That is a very simplistic test! */
2529 g_printerr ("Var failure.\n");
2530 ok = FALSE;
2531 }
2532
2533 T = skew_target;
2534 g_printerr ("Expected skew: %.10" GNM_FORMAT_g "\n", T);
2535 if (!gnm_finite (skew)) {
2536 /* That is a very simplistic test! */
2537 g_printerr ("Skew failure.\n");
2538 ok = FALSE;
2539 }
2540
2541 T = kurt_target;
2542 g_printerr ("Expected kurt: %.10" GNM_FORMAT_g "\n", T);
2543 if (!(kurt >= -3 && gnm_finite (kurt))) {
2544 /* That is a very simplistic test! */
2545 g_printerr ("Kurt failure.\n");
2546 ok = FALSE;
2547 }
2548
2549 /* Fractile test */
2550 for (i = 1; i < nf; i++)
2551 fractiles[i] = qrayleigh (i / (double)nf, ls, TRUE, FALSE);
2552 if (!rand_fractile_test (vals, N, nf, fractiles, NULL))
2553 ok = FALSE;
2554
2555 if (ok)
2556 g_printerr ("OK\n");
2557 else
2558 add_random_fail ("RANDRAYLEIGH");
2559 g_printerr ("\n");
2560
2561 g_free (vals);
2562 }
2563
2564 static void
2565 test_random (void)
2566 {
2567 const char *test_name = "test_random";
2568 const int N = sstest_fast ? 2000 : 20000;
2569 const int High_N = N * 10;
2570 const char *single = g_getenv ("SSTEST_RANDOM");
2571
2572 mark_test_start (test_name);
2573
2574 #define CHECK1(NAME,C) \
2575 do { if (!single || strcmp(single,#NAME) == 0) test_random_ ## NAME (C); } while (0)
2576
2577 /* Continuous */
2578 CHECK1 (rand, N);
2579 CHECK1 (randuniform, N);
2580 CHECK1 (randbeta, N);
2581 CHECK1 (randcauchy, N);
2582 CHECK1 (randchisq, N);
2583 CHECK1 (randexp, N);
2584 CHECK1 (randfdist, N);
2585 CHECK1 (randgamma, N);
2586 CHECK1 (randlog, N);
2587 CHECK1 (randlognorm, N);
2588 CHECK1 (randnorm, High_N);
2589 CHECK1 (randsnorm, High_N);
2590 CHECK1 (randtdist, N);
2591 CHECK1 (randweibull, N);
2592 CHECK1 (randrayleigh, N);
2593 #if 0
2594 CHECK1 (randexppow, N);
2595 CHECK1 (randgumbel, N);
2596 CHECK1 (randlandau, N);
2597 CHECK1 (randlaplace, N);
2598 CHECK1 (randlevy, N);
2599 CHECK1 (randlogistic, N);
2600 CHECK1 (randnormtail, N);
2601 CHECK1 (randpareto, N);
2602 CHECK1 (randrayleightail, N);
2603 CHECK1 (randstdist, N);
2604 #endif
2605
2606 /* Discrete */
2607 CHECK1 (randbernoulli, N);
2608 CHECK1 (randbetween, N);
2609 CHECK1 (randbinom, N);
2610 CHECK1 (randdiscrete, N);
2611 CHECK1 (randgeom, High_N);
2612 CHECK1 (randhyperg, High_N);
2613 CHECK1 (randnegbinom, High_N);
2614 CHECK1 (randpoisson, High_N);
2615
2616 #undef CHECK1
2617
2618 if (!single) {
2619 if (random_summary)
2620 g_printerr ("SUMMARY: FAIL for %s\n\n", random_summary);
2621 else
2622 g_printerr ("SUMMARY: OK\n\n");
2623 }
2624 g_free (random_summary);
2625 random_summary = NULL;
2626
2627 mark_test_end (test_name);
2628 }
2629
2630 /*-------------------------------------------------------------------------- */
2631
2632 #define MAYBE_DO(name) if (strcmp (testname, "all") != 0 && strcmp (testname, (name)) != 0) { } else
2633
2634 int
2635 main (int argc, char const **argv)
2636 {
2637 GOErrorInfo *plugin_errs;
2638 GOCmdContext *cc;
2639 GOptionContext *ocontext;
2640 GError *error = NULL;
2641 const char *testname;
2642
2643 /* No code before here, we need to init threads */
2644 argv = gnm_pre_parse_init (argc, argv);
2645
2646 ocontext = g_option_context_new (_("[testname]"));
2647 g_option_context_add_main_entries (ocontext, sstest_options, GETTEXT_PACKAGE);
2648 g_option_context_add_group (ocontext, gnm_get_option_group ());
2649 g_option_context_parse (ocontext, &argc, (gchar ***)&argv, &error);
2650 g_option_context_free (ocontext);
2651
2652 if (error) {
2653 g_printerr (_("%s\nRun '%s --help' to see a full list of available command line options.\n"),
2654 error->message, g_get_prgname ());
2655 g_error_free (error);
2656 return 1;
2657 }
2658
2659 if (sstest_show_version) {
2660 g_printerr (_("version '%s'\ndatadir := '%s'\nlibdir := '%s'\n"),
2661 GNM_VERSION_FULL, gnm_sys_data_dir (), gnm_sys_lib_dir ());
2662 return 0;
2663 }
2664
2665 gnm_init ();
2666
2667 cc = gnm_cmd_context_stderr_new ();
2668 gnm_plugins_init (GO_CMD_CONTEXT (cc));
2669 go_plugin_db_activate_plugin_list (
2670 go_plugins_get_available_plugins (), &plugin_errs);
2671 if (plugin_errs) {
2672 /* FIXME: What do we want to do here? */
2673 go_error_info_free (plugin_errs);
2674 }
2675
2676 if (ext_refs_file)
2677 return gnm_dump_func_defs (ext_refs_file, 4);
2678 if (samples_file)
2679 return gnm_dump_func_defs (samples_file, 5);
2680
2681 testname = argv[1];
2682 if (!testname) testname = "all";
2683
2684 /* ---------------------------------------- */
2685
2686 MAYBE_DO ("test_insdel_rowcol_names") test_insdel_rowcol_names ();
2687 MAYBE_DO ("test_insert_delete") test_insert_delete ();
2688 MAYBE_DO ("test_func_help") test_func_help ();
2689 MAYBE_DO ("test_nonascii_numbers") test_nonascii_numbers ();
2690 MAYBE_DO ("test_random") test_random ();
2691
2692 /* ---------------------------------------- */
2693
2694 g_object_unref (cc);
2695 gnm_shutdown ();
2696 gnm_pre_parse_shutdown ();
2697
2698 return 0;
2699 }
gnumeric mirror