/* nag_quasi_init_scrambled (g05ync) Example Program. * * Copyright 2008, Numerical Algorithms Group. * * Mark 9, 2009. */ #include #include #include #include #include #define QUAS(I, J) quas[(order == Nag_ColMajor)?(J*pdquas + I):(I*pdquas + J)] int main(void) { /*Integer scalar and array declarations */ Integer exit_status = 0; Integer liref, d, i, j, lstate, q_size; Integer *iref = 0, *state = 0; /* NAG structures */ Integer pdquas; NagError fail; /*Double scalar and array declarations */ double sum, tmp, vsbl; double *quas = 0; /* Number of dimensions */ Integer idim = 8; /* Set the sample size */ Integer n = 200; /* Skip the first 1000 variates */ Integer iskip = 1000; /* Use row major order */ Nag_OrderType order = Nag_RowMajor; /* Choose the base pseudo generator */ Nag_BaseRNG pgenid = Nag_Basic; Integer psubid = 0; /* Set the seed */ Integer seed[] = { 1762543 }; Integer lseed = 1; /* Choose the quasi generator */ Nag_QuasiRandom_Sequence genid = Nag_QuasiRandom_Sobol; /* Use Owen type scrambling */ Nag_QuasiRandom_Scrambling stype = Nag_OwenLike; /* Scramble the default number of digits */ Integer nsdigi = 0; /* Initialise the error structure */ INIT_FAIL(fail); printf("nag_quasi_init_scrambled (g05ync) Example Program Results\n"); /* Get the length of the state array */ lstate = -1; nag_rand_init_repeatable(pgenid, psubid, seed, lseed, state, &lstate, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n", fail.message); exit_status = 1; goto END; } pdquas = (order == Nag_RowMajor)?idim:n; q_size = (order == Nag_RowMajor)?pdquas * n:pdquas * idim; /* Calculate the size of the reference vector */ liref = (genid == Nag_QuasiRandom_Faure)?407:32 * idim + 7; /* Allocate arrays */ if (!(quas = NAG_ALLOC(q_size, double)) || !(iref = NAG_ALLOC(liref, Integer)) || !(state = NAG_ALLOC(lstate, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Initialise the psuedo-random generator used in the scrambling to a repeatable sequence */ nag_rand_init_repeatable(pgenid, psubid, seed, lseed, state, &lstate, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_rand_init_repeatable (g05kfc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Initialise the quasi-random sequence */ nag_quasi_init_scrambled(Nag_QuasiRandom_Sobol, stype, idim, iref, liref, iskip, nsdigi, state, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_quasi_init_scrambled (g05ync).\n%s\n", fail.message); exit_status = 1; goto END; } /* Generate n quasi-random variates */ nag_quasi_rand_uniform(order, n, quas, pdquas, iref, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_quasi_rand_uniform (g05ymc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Estimate integral by evaluating function at each variate and summing */ sum = 0.0e0; for (i = 0; i < n; i++) { tmp = 1.0e0; for (d = 0; d < idim; d++) tmp *= fabs(4.0e0 * QUAS(i, d) - 2.0e0); sum += tmp; } /* Convert sum to mean value*/ vsbl = sum / (double) n; /* Print the estimated value of the integral */ printf("Value of integral = %8.4f\n\n", vsbl); /* Display the first 10 variates used */ printf("First 10 variates\n"); for (i = 0; i < 10; i++) { printf(" %3ld", i + 1); for (j = 0; j < idim; j++) printf("%8.4f%s", QUAS(i, j), ((j+1)%20)?" ":"\n"); if (idim%20) printf("\n"); } END: if (quas) NAG_FREE(quas); if (iref) NAG_FREE(iref); if (state) NAG_FREE(state); return exit_status; }