From 326e1549ca91f5649663c855d1a6978aaef77e90 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Wed, 10 Dec 2014 14:09:37 +0100 Subject: [PATCH 1/9] a large restructure of the file structure to accomodate the new c files for montecarlo and tests --- tardis/cmontecarlo.c | 948 ----------------- tardis/model.py | 2 +- tardis/montecarlo/__init__.py | 1 + tardis/{ => montecarlo}/montecarlo.pyx | 0 tardis/montecarlo/setup_package.py | 16 + tardis/{ => montecarlo/src}/cmontecarlo.h | 0 tardis/{ => montecarlo/src}/randomkit/LICENSE | 0 .../src}/randomkit/randomkit.h | 0 .../{ => montecarlo/src}/randomkit/rk_isaac.h | 0 tardis/{ => montecarlo/src}/randomkit/rk_mt.h | 0 .../src}/randomkit/rk_primitive.h | 0 .../{ => montecarlo/src}/randomkit/rk_sobol.h | 0 tardis/parallel.py | 22 - tardis/randomkit/rk_isaac.c | 258 ----- tardis/randomkit/rk_mt.c | 342 ------ tardis/randomkit/rk_primitive.c | 520 --------- tardis/randomkit/rk_sobol.c | 991 ------------------ tardis/setup_package.py | 13 - 18 files changed, 18 insertions(+), 3095 deletions(-) delete mode 100644 tardis/cmontecarlo.c create mode 100644 tardis/montecarlo/__init__.py rename tardis/{ => montecarlo}/montecarlo.pyx (100%) create mode 100644 tardis/montecarlo/setup_package.py rename tardis/{ => montecarlo/src}/cmontecarlo.h (100%) rename tardis/{ => montecarlo/src}/randomkit/LICENSE (100%) rename tardis/{ => montecarlo/src}/randomkit/randomkit.h (100%) rename tardis/{ => montecarlo/src}/randomkit/rk_isaac.h (100%) rename tardis/{ => montecarlo/src}/randomkit/rk_mt.h (100%) rename tardis/{ => montecarlo/src}/randomkit/rk_primitive.h (100%) rename tardis/{ => montecarlo/src}/randomkit/rk_sobol.h (100%) delete mode 100644 tardis/parallel.py delete mode 100644 tardis/randomkit/rk_isaac.c delete mode 100644 tardis/randomkit/rk_mt.c delete mode 100644 tardis/randomkit/rk_primitive.c delete mode 100644 tardis/randomkit/rk_sobol.c delete mode 100644 tardis/setup_package.py diff --git a/tardis/cmontecarlo.c b/tardis/cmontecarlo.c deleted file mode 100644 index bc2229a50e1..00000000000 --- a/tardis/cmontecarlo.c +++ /dev/null @@ -1,948 +0,0 @@ -#include "cmontecarlo.h" - -rk_state mt_state; - -void -initialize_random_kit (unsigned long seed) -{ - rk_seed (seed, &mt_state); -} - -INLINE tardis_error_t -line_search (double *nu, double nu_insert, int64_t number_of_lines, - int64_t * result) -{ - tardis_error_t ret_val = TARDIS_ERROR_OK; - int64_t imin, imax; - imin = 0; - imax = number_of_lines - 1; - if (nu_insert > nu[imin]) - { - *result = imin; - } - else if (nu_insert < nu[imax]) - { - *result = imax + 1; - } - else - { - ret_val = reverse_binary_search (nu, nu_insert, imin, imax, result); - *result = *result + 1; - } - return ret_val; -} - -inline tardis_error_t -reverse_binary_search (double *x, double x_insert, - int64_t imin, int64_t imax, int64_t * result) -{ - /* - Have in mind that *x points to a reverse sorted array. - That is large values will have small indices and small ones - will have large indices. - */ - tardis_error_t ret_val = TARDIS_ERROR_OK; - if (x_insert > x[imin] || x_insert < x[imax]) - { - ret_val = TARDIS_ERROR_BOUNDS_ERROR; - } - else - { - int imid = (imin + imax) / 2; - while (imax - imin > 2) - { - if (x[imid] < x_insert) - { - imax = imid + 1; - } - else - { - imin = imid; - } - imid = (imin + imax) / 2; - } - if (imax - imin == 2 && x_insert < x[imin + 1]) - { - *result = imin + 1; - } - else - { - *result = imin; - } - } - return ret_val; -} - -inline tardis_error_t -binary_search (double *x, double x_insert, int64_t imin, - int64_t imax, int64_t * result) -{ - /* - Have in mind that *x points to a sorted array. - Like [1,2,3,4,5,...] - */ - int imid; - tardis_error_t ret_val = TARDIS_ERROR_OK; - if (x_insert < x[imin] || x_insert > x[imax]) - { - ret_val = TARDIS_ERROR_BOUNDS_ERROR; - } - else - { - while (imax >= imin) - { - imid = (imin + imax) / 2; - if (x[imid] == x_insert) - { - *result = imid; - break; - } - else if (x[imid] < x_insert) - { - imin = imid + 1; - } - else - { - imax = imid - 1; - } - } - if (imax - imid == 2 && x_insert < x[imin + 1]) - { - *result = imin; - } - else - { - *result = imin; - } - } - return ret_val; -} - -INLINE double -rpacket_doppler_factor (rpacket_t * packet, storage_model_t * storage) -{ - return 1.0 - - rpacket_get_mu (packet) * rpacket_get_r (packet) * - storage->inverse_time_explosion * INVERSE_C; -} - -INLINE double -compute_distance2boundary (rpacket_t * packet, storage_model_t * storage) -{ - double r = rpacket_get_r (packet); - double mu = rpacket_get_mu (packet); - double r_outer = storage->r_outer[rpacket_get_current_shell_id (packet)]; - double r_inner = storage->r_inner[rpacket_get_current_shell_id (packet)]; - double d_outer = - sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); - double d_inner; - double result; - if (rpacket_get_recently_crossed_boundary (packet) == 1) - { - rpacket_set_next_shell_id (packet, 1); - return d_outer; - } - else - { - double check = r_inner * r_inner + (r * r * (mu * mu - 1.0)); - if (check < 0.0) - { - rpacket_set_next_shell_id (packet, 1); - return d_outer; - } - else - { - d_inner = mu < 0.0 ? -r * mu - sqrt (check) : MISS_DISTANCE; - } - } - if (d_inner < d_outer) - { - rpacket_set_next_shell_id (packet, -1); - return d_inner; - } - else - { - rpacket_set_next_shell_id (packet, 1); - return d_outer; - } -} - -INLINE tardis_error_t -compute_distance2line (rpacket_t * packet, storage_model_t * storage, - double *result) -{ - tardis_error_t ret_val = TARDIS_ERROR_OK; - if (rpacket_get_last_line (packet)) - { - *result = MISS_DISTANCE; - } - else - { - double r = rpacket_get_r (packet); - double mu = rpacket_get_mu (packet); - double nu = rpacket_get_nu (packet); - double nu_line = rpacket_get_nu_line (packet); - double t_exp = storage->time_explosion; - double inverse_t_exp = storage->inverse_time_explosion; - int64_t cur_zone_id = rpacket_get_current_shell_id (packet); - double comov_nu, doppler_factor; - doppler_factor = 1.0 - mu * r * inverse_t_exp * INVERSE_C; - comov_nu = nu * doppler_factor; - if (comov_nu < nu_line) - { - if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) - { - fprintf (stderr, "last_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) - 1]); - fprintf (stderr, "Last line in line list reached!"); - } - else if (rpacket_get_next_line_id (packet) == 0) - { - fprintf (stderr, "First line in line list!"); - fprintf (stderr, "next_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) + 1]); - } - else - { - fprintf (stderr, "last_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) - 1]); - fprintf (stderr, "next_line = %f\n", - storage-> - line_list_nu[rpacket_get_next_line_id (packet) + 1]); - } - fprintf (stderr, "ERROR: Comoving nu less than nu_line!\n"); - fprintf (stderr, "comov_nu = %f\n", comov_nu); - fprintf (stderr, "nu_line = %f\n", nu_line); - fprintf (stderr, "(comov_nu - nu_line) / nu_line = %f\n", - (comov_nu - nu_line) / nu_line); - fprintf (stderr, "r = %f\n", r); - fprintf (stderr, "mu = %f\n", mu); - fprintf (stderr, "nu = %f\n", nu); - fprintf (stderr, "doppler_factor = %f\n", doppler_factor); - fprintf (stderr, "cur_zone_id = %d\n", cur_zone_id); - ret_val = TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE; - } - else - { - *result = ((comov_nu - nu_line) / nu) * C * t_exp; - } - } - return ret_val; -} - -INLINE double -compute_distance2electron (rpacket_t * packet, storage_model_t * storage) -{ - if (rpacket_get_virtual_packet (packet) > 0) - { - return MISS_DISTANCE; - } - double inverse_ne = - storage-> - inverse_electron_densities[rpacket_get_current_shell_id (packet)] * - storage->inverse_sigma_thomson; - return rpacket_get_tau_event (packet) * inverse_ne; -} - -INLINE int64_t -macro_atom (rpacket_t * packet, storage_model_t * storage) -{ - int emit = 0, i = 0; - double p, event_random; - int activate_level = - storage->line2macro_level_upper[rpacket_get_next_line_id (packet) - 1]; - while (emit != -1) - { - event_random = rk_double (&mt_state); - i = storage->macro_block_references[activate_level] - 1; - p = 0.0; - do - { - p += - storage-> - transition_probabilities[rpacket_get_current_shell_id (packet) * - storage->transition_probabilities_nd + - (++i)]; - } - while (p <= event_random); - emit = storage->transition_type[i]; - activate_level = storage->destination_level_id[i]; - } - return storage->transition_line_id[i]; -} - -INLINE double -move_packet (rpacket_t * packet, storage_model_t * storage, double distance) -{ - double new_r, doppler_factor, comov_energy, comov_nu; - doppler_factor = rpacket_doppler_factor (packet, storage); - if (distance > 0.0) - { - double r = rpacket_get_r (packet); - new_r = - sqrt (r * r + distance * distance + - 2.0 * r * distance * rpacket_get_mu (packet)); - rpacket_set_mu (packet, - (rpacket_get_mu (packet) * r + distance) / new_r); - rpacket_set_r (packet, new_r); - if (rpacket_get_virtual_packet (packet) <= 0) - { - comov_energy = rpacket_get_energy (packet) * doppler_factor; - comov_nu = rpacket_get_nu (packet) * doppler_factor; - storage->js[rpacket_get_current_shell_id (packet)] += - comov_energy * distance; - storage->nubars[rpacket_get_current_shell_id (packet)] += - comov_energy * distance * comov_nu; - } - } - return doppler_factor; -} - -INLINE void -increment_j_blue_estimator (rpacket_t * packet, storage_model_t * storage, - double d_line, int64_t j_blue_idx) -{ - double comov_energy, r_interaction, mu_interaction, doppler_factor; - double r = rpacket_get_r (packet); - r_interaction = - sqrt (r * r + d_line * d_line + - 2.0 * r * d_line * rpacket_get_mu (packet)); - mu_interaction = (rpacket_get_mu (packet) * r + d_line) / r_interaction; - doppler_factor = 1.0 - mu_interaction * r_interaction * - storage->inverse_time_explosion * INVERSE_C; - comov_energy = rpacket_get_energy (packet) * doppler_factor; - storage->line_lists_j_blues[j_blue_idx] += - comov_energy / rpacket_get_nu (packet); -} - -int64_t -montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_mode) -{ - int64_t i; - rpacket_t virt_packet; - double mu_bin; - double mu_min; - double doppler_factor_ratio; - double weight; - int64_t virt_id_nu; - int64_t reabsorbed; - if (virtual_mode == 0) - { - reabsorbed = montecarlo_one_packet_loop (storage, packet, 0); - } - else - { - for (i = 0; i < rpacket_get_virtual_packet_flag (packet); i++) - { - memcpy ((void *) &virt_packet, (void *) packet, sizeof (rpacket_t)); - if (virt_packet.r > storage->r_inner[0]) - { - mu_min = - -1.0 * sqrt (1.0 - - (storage->r_inner[0] / virt_packet.r) * - (storage->r_inner[0] / virt_packet.r)); - } - else - { - mu_min = 0.0; - } - mu_bin = (1.0 - mu_min) / rpacket_get_virtual_packet_flag (packet); - virt_packet.mu = mu_min + (i + rk_double (&mt_state)) * mu_bin; - switch (virtual_mode) - { - case -2: - weight = 1.0 / rpacket_get_virtual_packet_flag (packet); - break; - case -1: - weight = - 2.0 * virt_packet.mu / - rpacket_get_virtual_packet_flag (packet); - break; - case 1: - weight = - (1.0 - - mu_min) / 2.0 / rpacket_get_virtual_packet_flag (packet); - break; - default: - fprintf (stderr, "Something has gone horribly wrong!\n"); - } - doppler_factor_ratio = - rpacket_doppler_factor (packet, storage) / - rpacket_doppler_factor (&virt_packet, storage); - virt_packet.energy = - rpacket_get_energy (packet) * doppler_factor_ratio; - virt_packet.nu = rpacket_get_nu (packet) * doppler_factor_ratio; - reabsorbed = montecarlo_one_packet_loop (storage, &virt_packet, 1); - if ((virt_packet.nu < storage->spectrum_end_nu) && - (virt_packet.nu > storage->spectrum_start_nu)) - { - virt_id_nu = - floor ((virt_packet.nu - - storage->spectrum_start_nu) / - storage->spectrum_delta_nu); - storage->spectrum_virt_nu[virt_id_nu] += - virt_packet.energy * weight; - } - } - } - return reabsorbed; -} - -void -move_packet_across_shell_boundary (rpacket_t * packet, - storage_model_t * storage, double distance) -{ - double comov_energy, doppler_factor, comov_nu, inverse_doppler_factor; - move_packet (packet, storage, distance); - if (rpacket_get_virtual_packet (packet) > 0) - { - double delta_tau_event = distance * - storage->electron_densities[rpacket_get_current_shell_id (packet)] * - storage->sigma_thomson; - rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) + - delta_tau_event); - } - else - { - rpacket_reset_tau_event (packet); - } - if ((rpacket_get_current_shell_id (packet) < storage->no_of_shells - 1 - && rpacket_get_next_shell_id (packet) == 1) - || (rpacket_get_current_shell_id (packet) > 0 - && rpacket_get_next_shell_id (packet) == -1)) - { - rpacket_set_current_shell_id (packet, - rpacket_get_current_shell_id (packet) + - rpacket_get_next_shell_id (packet)); - rpacket_set_recently_crossed_boundary (packet, - rpacket_get_next_shell_id - (packet)); - } - else if (rpacket_get_next_shell_id (packet) == 1) - { - rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); - } - else if ((storage->reflective_inner_boundary == 0) || - (rk_double (&mt_state) > storage->inner_boundary_albedo)) - { - rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); - } - else - { - doppler_factor = rpacket_doppler_factor (packet, storage); - comov_nu = rpacket_get_nu (packet) * doppler_factor; - comov_energy = rpacket_get_energy (packet) * doppler_factor; - rpacket_set_mu (packet, rk_double (&mt_state)); - inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage); - rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - rpacket_set_recently_crossed_boundary (packet, 1); - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - montecarlo_one_packet (storage, packet, -2); - } - } -} - -void -montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, - double distance) -{ - double comov_energy, doppler_factor, comov_nu, inverse_doppler_factor; - doppler_factor = move_packet (packet, storage, distance); - comov_nu = rpacket_get_nu (packet) * doppler_factor; - comov_energy = rpacket_get_energy (packet) * doppler_factor; - rpacket_set_mu (packet, 2.0 * rk_double (&mt_state) - 1.0); - inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage); - rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - rpacket_reset_tau_event (packet); - rpacket_set_recently_crossed_boundary (packet, 0); - storage->last_interaction_type[storage->current_packet_id] = 1; - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - montecarlo_one_packet (storage, packet, 1); - } -} - -void -montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, - double distance) -{ - double comov_energy = 0.0; - int64_t emission_line_id = 0; - double old_doppler_factor = 0.0; - double inverse_doppler_factor = 0.0; - double tau_line = 0.0; - double tau_electron = 0.0; - double tau_combined = 0.0; - bool virtual_close_line = false; - int64_t j_blue_idx = -1; - if (rpacket_get_virtual_packet (packet) == 0) - { - j_blue_idx = - rpacket_get_current_shell_id (packet) * - storage->line_lists_j_blues_nd + rpacket_get_next_line_id (packet); - increment_j_blue_estimator (packet, storage, distance, j_blue_idx); - } - tau_line = - storage->line_lists_tau_sobolevs[rpacket_get_current_shell_id (packet) * - storage->line_lists_tau_sobolevs_nd + - rpacket_get_next_line_id (packet)]; - tau_electron = - storage->sigma_thomson * - storage->electron_densities[rpacket_get_current_shell_id (packet)] * - distance; - tau_combined = tau_line + tau_electron; - rpacket_set_next_line_id (packet, rpacket_get_next_line_id (packet) + 1); - if (rpacket_get_next_line_id (packet) == storage->no_of_lines) - { - rpacket_set_last_line (packet, true); - } - if (rpacket_get_virtual_packet (packet) > 0) - { - rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) + tau_line); - } - else if (rpacket_get_tau_event (packet) < tau_combined) - { - old_doppler_factor = move_packet (packet, storage, distance); - rpacket_set_mu (packet, 2.0 * rk_double (&mt_state) - 1.0); - inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage); - comov_energy = rpacket_get_energy (packet) * old_doppler_factor; - rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); - storage->last_line_interaction_in_id[storage->current_packet_id] = - rpacket_get_next_line_id (packet) - 1; - storage->last_line_interaction_shell_id[storage->current_packet_id] = - rpacket_get_current_shell_id (packet); - storage->last_interaction_type[storage->current_packet_id] = 2; - if (storage->line_interaction_id == 0) - { - emission_line_id = rpacket_get_next_line_id (packet) - 1; - } - else if (storage->line_interaction_id >= 1) - { - emission_line_id = macro_atom (packet, storage); - } - storage->last_line_interaction_out_id[storage->current_packet_id] = - emission_line_id; - rpacket_set_nu (packet, - storage->line_list_nu[emission_line_id] * - inverse_doppler_factor); - rpacket_set_nu_line (packet, storage->line_list_nu[emission_line_id]); - rpacket_set_next_line_id (packet, emission_line_id + 1); - rpacket_reset_tau_event (packet); - rpacket_set_recently_crossed_boundary (packet, 0); - if (rpacket_get_virtual_packet_flag (packet) > 0) - { - virtual_close_line = false; - if (!rpacket_get_last_line (packet) && - fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - - rpacket_get_nu_line (packet)) / - rpacket_get_nu_line (packet) < 1e-7) - { - virtual_close_line = true; - } - // QUESTIONABLE!!! - bool old_close_line = rpacket_get_close_line (packet); - rpacket_set_close_line (packet, virtual_close_line); - montecarlo_one_packet (storage, packet, 1); - rpacket_set_close_line (packet, old_close_line); - virtual_close_line = false; - } - } - else - { - rpacket_set_tau_event (packet, - rpacket_get_tau_event (packet) - tau_line); - } - if (!rpacket_get_last_line (packet) && - fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - - rpacket_get_nu_line (packet)) / rpacket_get_nu_line (packet) < - 1e-7) - { - rpacket_set_close_line (packet, true); - } -} - -INLINE void -montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage) -{ - // Check if the last line was the same nu as the current line. - if (rpacket_get_close_line (packet)) - { - // If so set the distance to the line to 0.0 - rpacket_set_d_line (packet, 0.0); - // Reset close_line. - rpacket_set_close_line (packet, false); - } - else - { - rpacket_set_d_boundary (packet, - compute_distance2boundary (packet, storage)); - double d_line; - compute_distance2line (packet, storage, &d_line); - rpacket_set_d_line (packet, d_line); - rpacket_set_d_electron (packet, - compute_distance2electron (packet, storage)); - } -} - -INLINE montecarlo_event_handler_t -get_event_handler (rpacket_t * packet, storage_model_t * storage, - double *distance) -{ - double d_boundary, d_electron, d_line; - montecarlo_compute_distances (packet, storage); - d_boundary = rpacket_get_d_boundary (packet); - d_electron = rpacket_get_d_electron (packet); - d_line = rpacket_get_d_line (packet); - montecarlo_event_handler_t handler; - if (d_line <= d_boundary && d_line <= d_electron) - { - *distance = d_line; - handler = &montecarlo_line_scatter; - } - else if (d_boundary <= d_electron) - { - *distance = d_boundary; - handler = &move_packet_across_shell_boundary; - } - else - { - *distance = d_electron; - handler = &montecarlo_thomson_scatter; - } - return handler; -} - -int64_t -montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet, - int64_t virtual_packet) -{ - rpacket_set_tau_event (packet, 0.0); - rpacket_set_nu_line (packet, 0.0); - rpacket_set_virtual_packet (packet, virtual_packet); - rpacket_set_status (packet, TARDIS_PACKET_STATUS_IN_PROCESS); - // Initializing tau_event if it's a real packet. - if (virtual_packet == 0) - { - rpacket_reset_tau_event (packet); - } - // For a virtual packet tau_event is the sum of all the tau's that the packet passes. - while (rpacket_get_status (packet) == TARDIS_PACKET_STATUS_IN_PROCESS) - { - // Check if we are at the end of line list. - if (!rpacket_get_last_line (packet)) - { - rpacket_set_nu_line (packet, - storage-> - line_list_nu[rpacket_get_next_line_id - (packet)]); - } - double distance; - get_event_handler (packet, storage, &distance) (packet, storage, - distance); - if (virtual_packet > 0 && rpacket_get_tau_event (packet) > 10.0) - { - rpacket_set_tau_event (packet, 100.0); - rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); - } - } - if (virtual_packet > 0) - { - rpacket_set_energy (packet, - rpacket_get_energy (packet) * exp (-1.0 * - rpacket_get_tau_event - (packet))); - } - return rpacket_get_status (packet) == - TARDIS_PACKET_STATUS_REABSORBED ? 1 : 0; -} - -tardis_error_t -rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index, - int virtual_packet_flag) -{ - double nu_line; - double current_r; - double current_mu; - double current_nu; - double comov_current_nu; - double current_energy; - int64_t current_line_id; - int current_shell_id; - bool last_line; - bool close_line; - int recently_crossed_boundary; - tardis_error_t ret_val = TARDIS_ERROR_OK; - current_nu = storage->packet_nus[packet_index]; - current_energy = storage->packet_energies[packet_index]; - current_mu = storage->packet_mus[packet_index]; - comov_current_nu = current_nu; - current_shell_id = 0; - current_r = storage->r_inner[0]; - current_nu = - current_nu / (1 - - (current_mu * current_r * storage->inverse_time_explosion * - INVERSE_C)); - current_energy = - current_energy / (1 - - (current_mu * current_r * - storage->inverse_time_explosion * INVERSE_C)); - if ((ret_val = - line_search (storage->line_list_nu, comov_current_nu, - storage->no_of_lines, - ¤t_line_id)) != TARDIS_ERROR_OK) - { - return ret_val; - } - last_line = (current_line_id == storage->no_of_lines); - recently_crossed_boundary = true; - rpacket_set_nu (packet, current_nu); - rpacket_set_mu (packet, current_mu); - rpacket_set_energy (packet, current_energy); - rpacket_set_r (packet, current_r); - rpacket_set_current_shell_id (packet, current_shell_id); - rpacket_set_next_line_id (packet, current_line_id); - rpacket_set_last_line (packet, last_line); - rpacket_set_close_line (packet, false); - rpacket_set_recently_crossed_boundary (packet, recently_crossed_boundary); - rpacket_set_virtual_packet_flag (packet, virtual_packet_flag); - return ret_val; -} - -/* - Getter and setter methods. -*/ - -INLINE double -rpacket_get_nu (rpacket_t * packet) -{ - return packet->nu; -} - -INLINE void -rpacket_set_nu (rpacket_t * packet, double nu) -{ - packet->nu = nu; -} - -INLINE double -rpacket_get_mu (rpacket_t * packet) -{ - return packet->mu; -} - -INLINE void -rpacket_set_mu (rpacket_t * packet, double mu) -{ - packet->mu = mu; -} - -INLINE double -rpacket_get_energy (rpacket_t * packet) -{ - return packet->energy; -} - -INLINE void -rpacket_set_energy (rpacket_t * packet, double energy) -{ - packet->energy = energy; -} - -INLINE double -rpacket_get_r (rpacket_t * packet) -{ - return packet->r; -} - -INLINE void -rpacket_set_r (rpacket_t * packet, double r) -{ - packet->r = r; -} - -INLINE double -rpacket_get_tau_event (rpacket_t * packet) -{ - return packet->tau_event; -} - -INLINE void -rpacket_set_tau_event (rpacket_t * packet, double tau_event) -{ - packet->tau_event = tau_event; -} - -INLINE double -rpacket_get_nu_line (rpacket_t * packet) -{ - return packet->nu_line; -} - -INLINE void -rpacket_set_nu_line (rpacket_t * packet, double nu_line) -{ - packet->nu_line = nu_line; -} - -INLINE unsigned int -rpacket_get_current_shell_id (rpacket_t * packet) -{ - return packet->current_shell_id; -} - -INLINE void -rpacket_set_current_shell_id (rpacket_t * packet, - unsigned int current_shell_id) -{ - packet->current_shell_id = current_shell_id; -} - -INLINE unsigned int -rpacket_get_next_line_id (rpacket_t * packet) -{ - return packet->next_line_id; -} - -INLINE void -rpacket_set_next_line_id (rpacket_t * packet, unsigned int next_line_id) -{ - packet->next_line_id = next_line_id; -} - -INLINE bool -rpacket_get_last_line (rpacket_t * packet) -{ - return packet->last_line; -} - -INLINE void -rpacket_set_last_line (rpacket_t * packet, bool last_line) -{ - packet->last_line = last_line; -} - -INLINE bool -rpacket_get_close_line (rpacket_t * packet) -{ - return packet->close_line; -} - -INLINE void -rpacket_set_close_line (rpacket_t * packet, bool close_line) -{ - packet->close_line = close_line; -} - -INLINE int -rpacket_get_recently_crossed_boundary (rpacket_t * packet) -{ - return packet->recently_crossed_boundary; -} - -INLINE void -rpacket_set_recently_crossed_boundary (rpacket_t * packet, - int recently_crossed_boundary) -{ - packet->recently_crossed_boundary = recently_crossed_boundary; -} - -INLINE int -rpacket_get_virtual_packet_flag (rpacket_t * packet) -{ - return packet->virtual_packet_flag; -} - -INLINE void -rpacket_set_virtual_packet_flag (rpacket_t * packet, int virtual_packet_flag) -{ - packet->virtual_packet_flag = virtual_packet_flag; -} - -INLINE int -rpacket_get_virtual_packet (rpacket_t * packet) -{ - return packet->virtual_packet; -} - -INLINE void -rpacket_set_virtual_packet (rpacket_t * packet, int virtual_packet) -{ - packet->virtual_packet = virtual_packet; -} - -INLINE double -rpacket_get_d_boundary (rpacket_t * packet) -{ - return packet->d_boundary; -} - -INLINE void -rpacket_set_d_boundary (rpacket_t * packet, double d_boundary) -{ - packet->d_boundary = d_boundary; -} - -INLINE double -rpacket_get_d_electron (rpacket_t * packet) -{ - return packet->d_electron; -} - -INLINE void -rpacket_set_d_electron (rpacket_t * packet, double d_electron) -{ - packet->d_electron = d_electron; -} - -INLINE double -rpacket_get_d_line (rpacket_t * packet) -{ - return packet->d_line; -} - -INLINE void -rpacket_set_d_line (rpacket_t * packet, double d_line) -{ - packet->d_line = d_line; -} - -INLINE int -rpacket_get_next_shell_id (rpacket_t * packet) -{ - return packet->next_shell_id; -} - -INLINE void -rpacket_set_next_shell_id (rpacket_t * packet, int next_shell_id) -{ - packet->next_shell_id = next_shell_id; -} - -INLINE rpacket_status_t -rpacket_get_status (rpacket_t * packet) -{ - return packet->status; -} - -INLINE void -rpacket_set_status (rpacket_t * packet, rpacket_status_t status) -{ - packet->status = status; -} - -/* Other accessor methods. */ - -INLINE void -rpacket_reset_tau_event (rpacket_t * packet) -{ - rpacket_set_tau_event (packet, -log (rk_double (&mt_state))); -} diff --git a/tardis/model.py b/tardis/model.py index 753d21cc0a7..fdd7665dd5c 100644 --- a/tardis/model.py +++ b/tardis/model.py @@ -10,7 +10,7 @@ import scipy.special from tardis import packet_source, plasma_array -import montecarlo +from tardis.montecarlo import montecarlo from util import intensity_black_body diff --git a/tardis/montecarlo/__init__.py b/tardis/montecarlo/__init__.py new file mode 100644 index 00000000000..8b137891791 --- /dev/null +++ b/tardis/montecarlo/__init__.py @@ -0,0 +1 @@ + diff --git a/tardis/montecarlo.pyx b/tardis/montecarlo/montecarlo.pyx similarity index 100% rename from tardis/montecarlo.pyx rename to tardis/montecarlo/montecarlo.pyx diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py new file mode 100644 index 00000000000..93fccc684c6 --- /dev/null +++ b/tardis/montecarlo/setup_package.py @@ -0,0 +1,16 @@ +#setting the right include +from setuptools import Extension +import numpy as np +import os + +randomkit_files = ['rk_isaac.c', 'rk_mt.c', 'rk_primitive.c','rk_sobol.c'] + +def get_extensions(): + return [Extension('tardis.montecarlo.montecarlo', + ['tardis/montecarlo/montecarlo.pyx', + 'tardis/montecarlo/src/cmontecarlo.c'] + + [os.path.join('tardis/montecarlo/src/randomkit', fname) + for fname in randomkit_files], + include_dirs=['tardis/montecarlo/src', + 'tardis/montecarlo/src/randomkit', + np.get_include()])] diff --git a/tardis/cmontecarlo.h b/tardis/montecarlo/src/cmontecarlo.h similarity index 100% rename from tardis/cmontecarlo.h rename to tardis/montecarlo/src/cmontecarlo.h diff --git a/tardis/randomkit/LICENSE b/tardis/montecarlo/src/randomkit/LICENSE similarity index 100% rename from tardis/randomkit/LICENSE rename to tardis/montecarlo/src/randomkit/LICENSE diff --git a/tardis/randomkit/randomkit.h b/tardis/montecarlo/src/randomkit/randomkit.h similarity index 100% rename from tardis/randomkit/randomkit.h rename to tardis/montecarlo/src/randomkit/randomkit.h diff --git a/tardis/randomkit/rk_isaac.h b/tardis/montecarlo/src/randomkit/rk_isaac.h similarity index 100% rename from tardis/randomkit/rk_isaac.h rename to tardis/montecarlo/src/randomkit/rk_isaac.h diff --git a/tardis/randomkit/rk_mt.h b/tardis/montecarlo/src/randomkit/rk_mt.h similarity index 100% rename from tardis/randomkit/rk_mt.h rename to tardis/montecarlo/src/randomkit/rk_mt.h diff --git a/tardis/randomkit/rk_primitive.h b/tardis/montecarlo/src/randomkit/rk_primitive.h similarity index 100% rename from tardis/randomkit/rk_primitive.h rename to tardis/montecarlo/src/randomkit/rk_primitive.h diff --git a/tardis/randomkit/rk_sobol.h b/tardis/montecarlo/src/randomkit/rk_sobol.h similarity index 100% rename from tardis/randomkit/rk_sobol.h rename to tardis/montecarlo/src/randomkit/rk_sobol.h diff --git a/tardis/parallel.py b/tardis/parallel.py deleted file mode 100644 index f6f9fd3e603..00000000000 --- a/tardis/parallel.py +++ /dev/null @@ -1,22 +0,0 @@ -#parallel launcher scripts -import tardis.simulation -from IPython.parallel import require - -@require(tardis.simulation) -def simulation_launcher(config_dict): - return tardis.simulation.run_multizone(config_dict, globals()['amodel']) - - -def parallel_macro_multizone(lbv, default_dict, **kwargs): - #lbv = client.load_balanced_view() - - config_dicts = [] - - for i in xrange(len(kwargs.items()[0])): - tmp_config_dict = default_dict.copy() - for key in kwargs: - tmp_config_dict[key] = kwargs[key][i] - config_dicts.append(tmp_config_dict) - - specs = lbv.map_async(simulation_launcher, config_dicts) - return specs \ No newline at end of file diff --git a/tardis/randomkit/rk_isaac.c b/tardis/randomkit/rk_isaac.c deleted file mode 100644 index 9f52fafae4c..00000000000 --- a/tardis/randomkit/rk_isaac.c +++ /dev/null @@ -1,258 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * ISAAC RNG by By Bob Jenkins. Based on Bob Jenkins public domain code. - * - * Original algorithm for the implementation of rk_interval function from - * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by - * Magnus Jonsson. - * - * Constants used in the rk_double implementation by Isaku Wada. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_isaac.c,v 1.4 2006/02/20 19:13:37 js Exp $"; - -#include "rk_isaac.h" -#include -#include -#include - -/* use the contents of randrsl[0..RK_ISAAC_STATE_LEN-1] as the seed. */ -static void isaac_init(rk_isaac_state *rg); - -void rk_isaac_seed(unsigned long seed, rk_isaac_state *state) -{ - rk_knuth_fill(seed, state->randrsl, RK_ISAAC_STATE_LEN); - isaac_init(state); -} - -rk_error rk_isaac_randomseed(rk_isaac_state *state) -{ - if(rk_devfill(state->randrsl, sizeof(state->randrsl), 1) == RK_NOERR) - { - isaac_init(state); - return RK_NOERR; - } - - rk_isaac_seed(rk_seedfromsystem(), state); - - return RK_ENODEV; -} - -#define ind(mm,x) ((mm)[(x>>2)&(RK_ISAAC_STATE_LEN-1)]) -#define rngstep(mix,a,b,mm,m,m2,r,x) \ -{ \ - x = *m; \ - a = ((a^(mix)) + *(m2++)) & 0xffffffff; \ - *(m++) = y = (ind(mm,x) + a + b) & 0xffffffff; \ - *(r++) = b = (ind(mm,y>>RK_ISAAC_STATE_POW) + x) & 0xffffffff; \ -} - -/* Call rand(rk_isaac_state *r) to retrieve a single 32-bit random value */ -unsigned long rk_isaac_random(rk_isaac_state *state) -{ - if (!state->randcnt--) - { - register unsigned long a,b,x,y,*m,*mm,*m2,*r,*mend; - mm=state->randmem; r=state->randrsl; - a = state->randa; b = (state->randb + (++state->randc)) & 0xffffffff; - for (m = mm, mend = m2 = m+(RK_ISAAC_STATE_LEN/2); m>6 , a, b, mm, m, m2, r, x); - rngstep( a<<2 , a, b, mm, m, m2, r, x); - rngstep( a>>16, a, b, mm, m, m2, r, x); - } - for (m2 = mm; m2>6 , a, b, mm, m, m2, r, x); - rngstep( a<<2 , a, b, mm, m, m2, r, x); - rngstep( a>>16, a, b, mm, m, m2, r, x); - } - state->randb = b; state->randa = a; - state->randcnt=RK_ISAAC_STATE_LEN-1; - } - return state->randrsl[state->randcnt] & 0xFFFFFFFF; -} - -#define mix(a,b,c,d,e,f,g,h) \ -{ \ - a^=b<<11; d+=a; b+=c; \ - b^=(c & 0xFFFFFFFF)>>2; e+=b; c+=d; \ - c^=d<<8; f+=c; d+=e; \ - d^=(e & 0xFFFFFFFF)>>16; g+=d; e+=f; \ - e^=f<<10; h+=e; f+=g; \ - f^=(g & 0xFFFFFFFF)>>4; a+=f; g+=h; \ - g^=h<<8; b+=g; h+=a; \ - h^=(a & 0xFFFFFFFF)>>9; c+=h; a+=b; \ -} - -/* if (flag==1), then use the contents of randrsl[] to initialize mm[]. */ -void isaac_init(rk_isaac_state *rk_isaac_state) -{ - int i; - unsigned long a,b,c,d,e,f,g,h; - unsigned long *m,*r; - rk_isaac_state->randa = rk_isaac_state->randb = rk_isaac_state->randc = 0; - m=rk_isaac_state->randmem; - r=rk_isaac_state->randrsl; - a=b=c=d=e=f=g=h=0x9e3779b9; /* the golden ratio */ - - for (i=0; i<4; ++i) /* scramble it */ - mix(a,b,c,d,e,f,g,h); - - /* initialize using the contents of r[] as the seed */ - for (i=0; irandcnt=0; - rk_isaac_state->has_gauss=0; -} - -long rk_isaac_long(rk_isaac_state *state) -{ - return rk_isaac_ulong(state) >> 1; -} - -unsigned long rk_isaac_ulong(rk_isaac_state *state) -{ -#if ULONG_MAX <= 0xffffffffUL - return rk_isaac_random(state); -#else - /* Assumes 64 bits */ - return (rk_isaac_random(state) << 32) | (rk_isaac_random(state)); -#endif -} - -unsigned long rk_isaac_interval(unsigned long max, rk_isaac_state *state) -{ - unsigned long mask = max, value; - - if (max == 0) return 0; - - /* Smallest bit mask >= max */ - mask |= mask >> 1; - mask |= mask >> 2; - mask |= mask >> 4; - mask |= mask >> 8; - mask |= mask >> 16; -#if ULONG_MAX > 0xffffffffUL - mask |= mask >> 32; -#endif - - /* Search a random value in [0..mask] <= max */ - while ((value = (rk_isaac_ulong(state) & mask)) > max); - - return value; -} - -double rk_isaac_double(rk_isaac_state *state) -{ - /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ - long a = rk_isaac_random(state) >> 5, b = rk_isaac_random(state) >> 6; - return (a * 67108864.0 + b) / 9007199254740992.0; -} - -void rk_isaac_copy(rk_isaac_state *copy, rk_isaac_state *orig) -{ - memcpy(copy, orig, sizeof(rk_isaac_state)); -} - -double rk_isaac_gauss(rk_isaac_state *state) -{ - if (state->has_gauss) - { - state->has_gauss = 0; - return state->gauss; - } - else - { - double f, x1, x2, r2; - do - { - x1 = 2.0*rk_isaac_double(state) - 1.0; - x2 = 2.0*rk_isaac_double(state) - 1.0; - r2 = x1*x1 + x2*x2; - } - while (r2 >= 1.0 || r2 == 0.0); - - f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ - state->has_gauss = 1; - state->gauss = f*x1; /* Keep for next call */ - return f*x2; - } -} - -void rk_isaac_fill(void *buffer, size_t size, rk_isaac_state *state) -{ - unsigned long r; - unsigned char *buf = buffer; - rk_isaac_state tempstate; - - if (size > 0 && state == NULL) - { - rk_isaac_randomseed(&tempstate); - state = &tempstate; - } - - for (; size >= 4; size -= 4) - { - r = rk_isaac_random(state); - *(buf++) = r & 0xFF; - *(buf++) = (r >> 8) & 0xFF; - *(buf++) = (r >> 16) & 0xFF; - *(buf++) = (r >> 24) & 0xFF; - } - - if (!size) return; - - r = rk_isaac_random(state); - - for (; size; r >>= 8, size --) - *(buf++) = (unsigned char)(r & 0xFF); -} - -void rk_seed_isaac(rk_isaac_state *i_state, rk_state *state) -{ - rk_isaac_fill(state->key, sizeof(state->key), i_state); - state->pos = RK_STATE_LEN; - state->has_gauss = 0; -} diff --git a/tardis/randomkit/rk_mt.c b/tardis/randomkit/rk_mt.c deleted file mode 100644 index 6ef13513989..00000000000 --- a/tardis/randomkit/rk_mt.c +++ /dev/null @@ -1,342 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * The rk_random and rk_seed functions algorithms and the original design of - * the Mersenne Twister RNG: - * - * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, - * All rights reserved. - * - * Original algorithm for the implementation of rk_interval function from - * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by - * Magnus Jonsson. - * - * Constants used in the rk_double implementation by Isaku Wada. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_mt.c,v 1.6 2006/02/19 13:48:34 js Exp $"; - -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef _WIN32 -/* Windows */ -#include -#ifndef RK_NO_WINCRYPT -/* Windows crypto */ -#ifndef _WIN32_WINNT -#define _WIN32_WINNT 0x0400 -#endif -#include -#include -#endif -#else -/* Unix */ -#include -#include -#endif - -#include "randomkit.h" - -#ifndef RK_DEV_URANDOM -#define RK_DEV_URANDOM "/dev/urandom" -#endif - -#ifndef RK_DEV_RANDOM -#define RK_DEV_RANDOM "/dev/random" -#endif - -char *rk_strerror[RK_ERR_MAX] = -{ - "no error", - "random device unvavailable" -}; - -/* static functions */ -static unsigned long rk_hash(unsigned long key); - -void rk_knuth_fill(unsigned long seed, unsigned long *key, unsigned long len) -{ - int pos; - seed &= 0xffffffffUL; - - /* Knuth's PRNG as used in the Mersenne Twister reference implementation */ - for (pos=0; pos> 30)) + pos + 1) & 0xffffffffUL; - } -} - -void rk_seed(unsigned long seed, rk_state *state) -{ - rk_knuth_fill(seed, state->key, RK_STATE_LEN); - state->pos = RK_STATE_LEN; - state->has_gauss = 0; -} - -/* Thomas Wang 32 bits integer hash function */ -unsigned long rk_hash(unsigned long key) -{ - key += ~(key << 15); - key ^= (key >> 10); - key += (key << 3); - key ^= (key >> 6); - key += ~(key << 11); - key ^= (key >> 16); - return key; -} - -unsigned long rk_seedfromsystem() -{ -#ifndef _WIN32 - struct timeval tv; -#else - struct _timeb tv; -#endif - -#ifndef _WIN32 - gettimeofday(&tv, NULL); - return rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec) - ^ rk_hash(clock()); -#else - _ftime(&tv); - return rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()); -#endif -} - -rk_error rk_randomseed(rk_state *state) -{ - if(rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) - { - state->key[0] |= 0x80000000UL; /* ensures non-zero key */ - state->pos = RK_STATE_LEN; - state->has_gauss = 0; - return RK_NOERR; - } - - rk_seed(rk_seedfromsystem(), state); - - return RK_ENODEV; -} - -/* Magic Mersenne Twister constants */ -#define N 624 -#define M 397 -#define MATRIX_A 0x9908b0dfUL -#define UPPER_MASK 0x80000000UL -#define LOWER_MASK 0x7fffffffUL - -/* Slightly optimised reference implementation of the Mersenne Twister */ -unsigned long rk_random(rk_state *state) -{ - unsigned long y; - - if (state->pos == RK_STATE_LEN) - { - int i; - - for (i=0;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); - state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - } - for (;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); - state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - } - y = (state->key[N-1] & UPPER_MASK) | (state->key[0] & LOWER_MASK); - state->key[N-1] = state->key[M-1] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); - - state->pos = 0; - } - - y = state->key[state->pos++]; - - /* Tempering */ - y ^= (y >> 11); - y ^= (y << 7) & 0x9d2c5680UL; - y ^= (y << 15) & 0xefc60000UL; - y ^= (y >> 18); - - return y; -} - -long rk_long(rk_state *state) -{ - return rk_ulong(state) >> 1; -} - -unsigned long rk_ulong(rk_state *state) -{ -#if ULONG_MAX <= 0xffffffffUL - return rk_random(state); -#else - /* Assumes 64 bits */ - return (rk_random(state) << 32) | (rk_random(state)); -#endif -} - -unsigned long rk_interval(unsigned long max, rk_state *state) -{ - unsigned long mask = max, value; - - if (max == 0) return 0; - - /* Smallest bit mask >= max */ - mask |= mask >> 1; - mask |= mask >> 2; - mask |= mask >> 4; - mask |= mask >> 8; - mask |= mask >> 16; -#if ULONG_MAX > 0xffffffffUL - mask |= mask >> 32; -#endif - - /* Search a random value in [0..mask] <= max */ - while ((value = (rk_ulong(state) & mask)) > max); - - return value; -} - -double rk_double(rk_state *state) -{ - /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ - long a = rk_random(state) >> 5, b = rk_random(state) >> 6; - return (a * 67108864.0 + b) / 9007199254740992.0; -} - -void rk_copy(rk_state *copy, rk_state *orig) -{ - memcpy(copy, orig, sizeof(rk_state)); -} - -void rk_fill(void *buffer, size_t size, rk_state *state) -{ - unsigned long r; - unsigned char *buf = buffer; - rk_state tempstate; - - if (size > 0 && state == NULL) - { - rk_randomseed(&tempstate); - state = &tempstate; - } - - for (; size >= 4; size -= 4) - { - r = rk_random(state); - *(buf++) = r & 0xFF; - *(buf++) = (r >> 8) & 0xFF; - *(buf++) = (r >> 16) & 0xFF; - *(buf++) = (r >> 24) & 0xFF; - } - - if (!size) return; - - r = rk_random(state); - - for (; size; r >>= 8, size --) - *(buf++) = (unsigned char)(r & 0xFF); -} - -rk_error rk_devfill(void *buffer, size_t size, int strong) -{ -#ifndef _WIN32 - FILE *rfile; - int done; - - if (strong) - rfile = fopen(RK_DEV_RANDOM, "rb"); - else - rfile = fopen(RK_DEV_URANDOM, "rb"); - if (rfile == NULL) - return RK_ENODEV; - done = fread(buffer, size, 1, rfile); - fclose(rfile); - if (done) - return RK_NOERR; -#else - -#ifndef RK_NO_WINCRYPT - HCRYPTPROV hCryptProv; - BOOL done; - - if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL, - CRYPT_VERIFYCONTEXT) || !hCryptProv) - return RK_ENODEV; - done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer); - CryptReleaseContext(hCryptProv, 0); - if (done) - return RK_NOERR; -#endif - -#endif - - return RK_ENODEV; -} - -rk_error rk_altfill(void *buffer, size_t size, int strong, rk_state *state) -{ - rk_error err; - - err = rk_devfill(buffer, size, strong); - if (err) - rk_fill(buffer, size, state); - - return err; -} - -double rk_gauss(rk_state *state) -{ - if (state->has_gauss) - { - state->has_gauss = 0; - return state->gauss; - } - else - { - double f, x1, x2, r2; - do - { - x1 = 2.0*rk_double(state) - 1.0; - x2 = 2.0*rk_double(state) - 1.0; - r2 = x1*x1 + x2*x2; - } - while (r2 >= 1.0 || r2 == 0.0); - - f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ - state->has_gauss = 1; - state->gauss = f*x1; /* Keep for next call */ - return f*x2; - } -} diff --git a/tardis/randomkit/rk_primitive.c b/tardis/randomkit/rk_primitive.c deleted file mode 100644 index 25e70c6c7f5..00000000000 --- a/tardis/randomkit/rk_primitive.c +++ /dev/null @@ -1,520 +0,0 @@ -/* Random kit 1.6 */ -/* Primitivity test for binary polynomials of low degree */ - -/* - * Copyright (c) 2005-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Methodology inspired by scott duplichan's ppsearch code. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_primitive.c,v 1.7 2006/02/19 13:48:34 js Exp $"; - -#include -#include -#include "rk_primitive.h" - -#ifndef LONG_BIT -#if ULONG_MAX <= 0xffffffffUL -#define LONG_BIT 32 -#else -#define LONG_BIT 64 -#endif -#endif - -static unsigned long modmul(unsigned long poly1, unsigned long poly2, - unsigned long modulo, unsigned long mask); -static unsigned long modpow(unsigned long polynomial, unsigned long power, - unsigned long modulo, int degree); - -/* - * For all powers i of two up to 64, list all the number of the form - * (2^i-1)/p for p a prime factor of 2^i-1. - */ -static const unsigned long divisors[][12]={ -/* 2^0-1 */ - {1UL, - 0UL}, -/* 2^1-1 */ - {1UL, - 0UL}, -/* 2^2-1 */ - {1UL, - 0UL}, -/* 2^3-1 */ - {1UL, - 0UL}, -/* 2^4-1 */ - {5UL, - 3UL, - 0UL}, -/* 2^5-1 */ - {1UL, - 0UL}, -/* 2^6-1 */ - {21UL, - 9UL, - 0UL}, -/* 2^7-1 */ - {1UL, - 0UL}, -/* 2^8-1 */ - {85UL, - 51UL, - 15UL, - 0UL}, -/* 2^9-1 */ - {73UL, - 7UL, - 0UL}, -/* 2^10-1 */ - {341UL, - 93UL, - 33UL, - 0UL}, -/* 2^11-1 */ - {89UL, - 23UL, - 0UL}, -/* 2^12-1 */ - {1365UL, - 819UL, - 585UL, - 315UL, - 0UL}, -/* 2^13-1 */ - {1UL, - 0UL}, -/* 2^14-1 */ - {5461UL, - 381UL, - 129UL, - 0UL}, -/* 2^15-1 */ - {4681UL, - 1057UL, - 217UL, - 0UL}, -/* 2^16-1 */ - {21845UL, - 13107UL, - 3855UL, - 255UL, - 0UL}, -/* 2^17-1 */ - {1UL, - 0UL}, -/* 2^18-1 */ - {87381UL, - 37449UL, - 13797UL, - 3591UL, - 0UL}, -/* 2^19-1 */ - {1UL, - 0UL}, -/* 2^20-1 */ - {349525UL, - 209715UL, - 95325UL, - 33825UL, - 25575UL, - 0UL}, -/* 2^21-1 */ - {299593UL, - 16513UL, - 6223UL, - 0UL}, -/* 2^22-1 */ - {1398101UL, - 182361UL, - 47127UL, - 6141UL, - 0UL}, -/* 2^23-1 */ - {178481UL, - 47UL, - 0UL}, -/* 2^24-1 */ - {5592405UL, - 3355443UL, - 2396745UL, - 1290555UL, - 986895UL, - 69615UL, - 0UL}, -/* 2^25-1 */ - {1082401UL, - 55831UL, - 18631UL, - 0UL}, -/* 2^26-1 */ - {22369621UL, - 24573UL, - 8193UL, - 0UL}, -/* 2^27-1 */ - {19173961UL, - 1838599UL, - 511UL, - 0UL}, -/* 2^28-1 */ - {89478485UL, - 53687091UL, - 9256395UL, - 6242685UL, - 2375535UL, - 2113665UL, - 0UL}, -/* 2^29-1 */ - {2304167UL, - 486737UL, - 256999UL, - 0UL}, -/* 2^30-1 */ - {357913941UL, - 153391689UL, - 97612893UL, - 34636833UL, - 7110873UL, - 3243933UL, - 0UL}, -/* 2^31-1 */ - {1UL, - 0UL}, -/* 2^32-1 */ - {1431655765UL, - 858993459UL, - 252645135UL, - 16711935UL, - 65535UL, - 0UL}, -#if LONG_BIT > 32 -/* 2^33-1 */ - {1227133513UL, - 373475417UL, - 96516119UL, - 14329UL, - 0UL}, -/* 2^34-1 */ - {5726623061UL, - 393213UL, - 131073UL, - 0UL}, -/* 2^35-1 */ - {1108378657UL, - 483939977UL, - 270549121UL, - 279527UL, - 0UL}, -/* 2^36-1 */ - {22906492245UL, - 13743895347UL, - 9817068105UL, - 5286113595UL, - 3616814565UL, - 1857283155UL, - 941362695UL, - 630453915UL, - 0UL}, -/* 2^37-1 */ - {616318177UL, - 223UL, - 0UL}, -/* 2^38-1 */ - {91625968981UL, - 1572861UL, - 524289UL, - 0UL}, -/* 2^39-1 */ - {78536544841UL, - 6958934353UL, - 67117057UL, - 4529623UL, - 0UL}, -/* 2^40-1 */ - {366503875925UL, - 219902325555UL, - 99955602525UL, - 64677154575UL, - 35468117025UL, - 26817356775UL, - 17825775UL, - 0UL}, -/* 2^41-1 */ - {164511353UL, - 13367UL, - 0UL}, -/* 2^42-1 */ - {1466015503701UL, - 628292358729UL, - 102280151421UL, - 34630287489UL, - 13050583119UL, - 811597437UL, - 0UL}, -/* 2^43-1 */ - {20408568497UL, - 905040953UL, - 4188889UL, - 0UL}, -/* 2^44-1 */ - {5864062014805UL, - 3518437208883UL, - 764877654105UL, - 197665011735UL, - 44312811195UL, - 25757227005UL, - 8325691455UL, - 0UL}, -/* 2^45-1 */ - {5026338869833UL, - 1134979744801UL, - 481977699847UL, - 233009086681UL, - 55759702201UL, - 1509346321UL, - 0UL}, -/* 2^46-1 */ - {23456248059221UL, - 1497207322929UL, - 394264623UL, - 25165821UL, - 0UL}, -/* 2^47-1 */ - {59862819377UL, - 31184907679UL, - 10610063UL, - 0UL}, -/* 2^48-1 */ - {93824992236885UL, - 56294995342131UL, - 40210710958665UL, - 21651921285435UL, - 16557351571215UL, - 2901803883615UL, - 1167945961455UL, - 1095233372415UL, - 418239192735UL, - 0UL}, -/* 2^49-1 */ - {4432676798593UL, - 127UL, - 0UL}, -/* 2^50-1 */ - {375299968947541UL, - 102354536985693UL, - 36319351833633UL, - 4485656999373UL, - 1873377548823UL, - 625152641223UL, - 277931351973UL, - 0UL}, -/* 2^51-1 */ - {321685687669321UL, - 21862134113449UL, - 1050769861729UL, - 202518195313UL, - 17180000257UL, - 0UL}, -/* 2^52-1 */ - {1501199875790165UL, - 900719925474099UL, - 84973577874915UL, - 28685347945035UL, - 2792064245115UL, - 1649066139645UL, - 549822930945UL, - 0UL}, -/* 2^53-1 */ - {1416003655831UL, - 129728784761UL, - 441650591UL, - 0UL}, -/* 2^54-1 */ - {6004799503160661UL, - 2573485501354569UL, - 948126237341157UL, - 246772582321671UL, - 206561081853UL, - 68585259519UL, - 0UL}, -/* 2^55-1 */ - {1566469435607129UL, - 1162219258676257UL, - 404817944033303UL, - 40895342813807UL, - 11290754314937UL, - 178394823847UL, - 0UL}, -/* 2^56-1 */ - {24019198012642645UL, - 14411518807585587UL, - 4238682002231055UL, - 2484744621997515UL, - 1675758000882045UL, - 637677823344495UL, - 567382630219905UL, - 4563402735UL, - 0UL}, -/* 2^57-1 */ - {20587884010836553UL, - 4451159405623UL, - 274878431233UL, - 118823881393UL, - 0UL}, -/* 2^58-1 */ - {96076792050570581UL, - 4885260612740877UL, - 1237040240994471UL, - 261314937580881UL, - 137975287770087UL, - 95026151247UL, - 0UL}, -/* 2^59-1 */ - {3203431780337UL, - 179951UL, - 0UL}, -/* 2^60-1 */ - {384307168202282325UL, - 230584300921369395UL, - 164703072086692425UL, - 104811045873349725UL, - 88686269585142075UL, - 37191016277640225UL, - 28120036697727975UL, - 18900352534538475UL, - 7635241752363225UL, - 3483146539597725UL, - 872764197279975UL, - 0UL}, -/* 2^61-1 */ - {1UL, - 0UL}, -/* 2^62-1 */ - {1537228672809129301UL, - 6442450941UL, - 2147483649UL, - 0UL}, -/* 2^63-1 */ - {1317624576693539401UL, - 126347562148695559UL, - 72624976668147841UL, - 27369056489183311UL, - 99457304386111UL, - 14197294936951UL, - 0UL}, -/* 2^64-1 */ - {6148914691236517205UL, - 3689348814741910323UL, - 1085102592571150095UL, - 71777214294589695UL, - 28778071877862015UL, - 281470681808895UL, - 2753074036095UL, - 0UL}, -#if LONG_BIT > 64 -#error Factorization of numbers up to 2^LONG_BIT required -#endif -#endif -}; - -/* - * Modular multiply for two binary polynomial - * mask is 1UL << the degree of the modulus. - */ -unsigned long modmul(unsigned long poly1, unsigned long poly2, - unsigned long modulo, unsigned long mask) -{ - unsigned long result = 0; - - for (; poly1; poly1 >>= 1) - { - if (poly1 & 1) - result ^= poly2; - - poly2 <<= 1; - if (poly2 & mask) - poly2 ^= modulo; - } - return result; -} - -/* - * Modular exponentiation for a binary polynomial - * degree is the degree of the modulus. - */ -unsigned long modpow(unsigned long polynomial, unsigned long power, - unsigned long modulo, int degree) -{ - unsigned long result = 1, mask = 1UL << degree; - - for (; power; power >>= 1) - { - if (power & 1) - result = modmul(result, polynomial, modulo, mask); - polynomial = modmul(polynomial, polynomial, modulo, mask); - } - return result; -} - -/* - * Test the primitivity of a polynomial - */ -int rk_isprimitive(unsigned long polynomial) -{ - unsigned long pelement = 2, temp = polynomial >> 1; - int k, degree = 0, weight = 1; - - /* Special case for polynomials of degree < 2 */ - if (polynomial < 4) - return (polynomial == 3) || (polynomial == 1); - - /* A binary primitive polynomial has a constant term */ - if (!(polynomial & 1)) - return 0; - - /* - * A binary primitive polynomial of degree > 1 has an odd number of terms. - * temp ^= temp >> 16; temp ^= temp >> 8; ... would be sligthly faster. - * Compute the degree at the same time. - */ - for (; temp; degree++, temp >>= 1) - weight += temp & 1; - if (!(weight & 1)) - return 0; - - /* - * Check if the period is 2^degree-1. - * Sufficient if 2^degree-1 is prime. - */ - if (modpow(pelement, 1UL << degree, polynomial, degree) != pelement) - return 0; - - if (divisors[degree][0] != 1) - /* Primitivity test */ - for (k = 0; divisors[degree][k]; k++) - if (modpow(pelement, divisors[degree][k], polynomial, degree) == 1) - return 0; - - return 1; -} diff --git a/tardis/randomkit/rk_sobol.c b/tardis/randomkit/rk_sobol.c deleted file mode 100644 index c6a016bb033..00000000000 --- a/tardis/randomkit/rk_sobol.c +++ /dev/null @@ -1,991 +0,0 @@ -/* Random kit 1.6 */ - -/* - * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) - * - * Original algorithm from Numerical Recipes, 2nd edition, by Press et al. - * The inverse normal cdf formulas are from Peter J. Acklam. - * The initialization directions were found in Ferdinando Ametrano's - * implementation in QuantLib. - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject to - * the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS - * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. - * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY - * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, - * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE - * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. - */ - -static char const rcsid[] = - "@(#) $Jeannot: rk_sobol.c,v 1.9 2006/02/19 13:48:34 js Exp $"; - -#include -#include -#include -#include "rk_sobol.h" -#include "rk_mt.h" -#include "rk_primitive.h" - -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -#ifndef M_SQRT1_2 -#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ -#endif -#define RK_SOBOL_M_SQRT2PI 2.506628274631000502415 /* sqrt(2*pi) */ - -#ifndef LONG_BIT -#if ULONG_MAX <= 0xffffffffUL -#define LONG_BIT 32 -#else -#define LONG_BIT 64 -#endif -#endif - -char *rk_sobol_strerror[] = -{ - "no error", - "invalid dimension", - "too many numbers generated", - "not enough memory" -}; - -static double inverse_normal(double p); - -/* - * Sobol/Levitan coefficients of the free direction integers as given - * by Bratley, P., Fox, B.L. (1988) - */ - -const unsigned long rk_sobol_SLdirections[] = { - 1, - 1, 1, - 1, 3, 7, - 1, 1, 5, - 1, 3, 1, 1, - 1, 1, 3, 7, - 1, 3, 3, 9, 9, - 1, 3, 7, 13, 3, - 1, 1, 5, 11, 27, - 1, 3, 5, 1, 15, - 1, 1, 7, 3, 29, - 1, 3, 7, 7, 21, - 1, 1, 1, 9, 23, 37, - 1, 3, 3, 5, 19, 33, - 1, 1, 3, 13, 11, 7, - 1, 1, 7, 13, 25, 5, - 1, 3, 5, 11, 7, 11, - 1, 1, 1, 3, 13, 39, - 1, 3, 1, 15, 17, 63, 13, - 1, 1, 5, 5, 1, 27, 33, - 1, 3, 3, 3, 25, 17, 115, - 1, 1, 3, 15, 29, 15, 41, - 1, 3, 1, 7, 3, 23, 79, - 1, 3, 7, 9, 31, 29, 17, - 1, 1, 5, 13, 11, 3, 29, - 1, 3, 1, 9, 5, 21, 119, - 1, 1, 3, 1, 23, 13, 75, - 1, 3, 3, 11, 27, 31, 73, - 1, 1, 7, 7, 19, 25, 105, - 1, 3, 5, 5, 21, 9, 7, - 1, 1, 1, 15, 5, 49, 59, - 1, 1, 1, 1, 1, 33, 65, - 1, 3, 5, 15, 17, 19, 21, - 1, 1, 7, 11, 13, 29, 3, - 1, 3, 7, 5, 7, 11, 113, - 1, 1, 5, 3, 15, 19, 61, - 1, 3, 1, 1, 9, 27, 89, 7, - 1, 1, 3, 7, 31, 15, 45, 23, - 1, 3, 3, 9, 9, 25, 107, 39, - 0 -}; - -/* - * Lemieux coefficients of the free direction integers as given - * in QuantLib by Christiane Lemieux, private communication, September 2004 - */ - -const unsigned long rk_sobol_Ldirections[] = { - 1, - 1, 1, - 1, 3, 7, - 1, 1, 5, - 1, 3, 1, 1, - 1, 1, 3, 7, - 1, 3, 3, 9, 9, - 1, 3, 7, 13, 3, - 1, 1, 5, 11, 27, - 1, 3, 5, 1, 15, - 1, 1, 7, 3, 29, - 1, 3, 7, 7, 21, - 1, 1, 1, 9, 23, 37, - 1, 3, 3, 5, 19, 33, - 1, 1, 3, 13, 11, 7, - 1, 1, 7, 13, 25, 5, - 1, 3, 5, 11, 7, 11, - 1, 1, 1, 3, 13, 39, - 1, 3, 1, 15, 17, 63, 13, - 1, 1, 5, 5, 1, 27, 33, - 1, 3, 3, 3, 25, 17, 115, - 1, 1, 3, 15, 29, 15, 41, - 1, 3, 1, 7, 3, 23, 79, - 1, 3, 7, 9, 31, 29, 17, - 1, 1, 5, 13, 11, 3, 29, - 1, 3, 1, 9, 5, 21, 119, - 1, 1, 3, 1, 23, 13, 75, - 1, 3, 3, 11, 27, 31, 73, - 1, 1, 7, 7, 19, 25, 105, - 1, 3, 5, 5, 21, 9, 7, - 1, 1, 1, 15, 5, 49, 59, - 1, 1, 1, 1, 1, 33, 65, - 1, 3, 5, 15, 17, 19, 21, - 1, 1, 7, 11, 13, 29, 3, - 1, 3, 7, 5, 7, 11, 113, - 1, 1, 5, 3, 15, 19, 61, - 1, 3, 1, 1, 9, 27, 89, 7, - 1, 1, 3, 7, 31, 15, 45, 23, - 1, 3, 3, 9, 9, 25, 107, 39, - 1, 1, 3, 13, 7, 35, 61, 91, - 1, 1, 7, 11, 5, 35, 55, 75, - 1, 3, 5, 5, 11, 23, 29, 139, - 1, 1, 1, 7, 11, 15, 17, 81, - 1, 1, 7, 9, 5, 57, 79, 103, - 1, 1, 7, 13, 19, 5, 5, 185, - 1, 3, 1, 3, 13, 57, 97, 131, - 1, 1, 5, 5, 21, 25, 125, 197, - 1, 3, 3, 9, 31, 11, 103, 201, - 1, 1, 5, 3, 7, 25, 51, 121, - 1, 3, 7, 15, 19, 53, 73, 189, - 1, 1, 1, 15, 19, 55, 27, 183, - 1, 1, 7, 13, 3, 29, 109, 69, - 1, 1, 5, 15, 15, 23, 15, 1, 57, - 1, 3, 1, 3, 23, 55, 43, 143, 397, - 1, 1, 3, 11, 29, 9, 35, 131, 411, - 1, 3, 1, 7, 27, 39, 103, 199, 277, - 1, 3, 7, 3, 19, 55, 127, 67, 449, - 1, 3, 7, 3, 5, 29, 45, 85, 3, - 1, 3, 5, 5, 13, 23, 75, 245, 453, - 1, 3, 1, 15, 21, 47, 3, 77, 165, - 1, 1, 7, 9, 15, 5, 117, 73, 473, - 1, 3, 1, 9, 1, 21, 13, 173, 313, - 1, 1, 7, 3, 11, 45, 63, 77, 49, - 1, 1, 1, 1, 1, 25, 123, 39, 259, - 1, 1, 1, 5, 23, 11, 59, 11, 203, - 1, 3, 3, 15, 21, 1, 73, 71, 421, - 1, 1, 5, 11, 15, 31, 115, 95, 217, - 1, 1, 3, 3, 7, 53, 37, 43, 439, - 1, 1, 1, 1, 27, 53, 69, 159, 321, - 1, 1, 5, 15, 29, 17, 19, 43, 449, - 1, 1, 3, 9, 1, 55, 121, 205, 255, - 1, 1, 3, 11, 9, 47, 107, 11, 417, - 1, 1, 1, 5, 17, 25, 21, 83, 95, - 1, 3, 5, 13, 31, 25, 61, 157, 407, - 1, 1, 7, 9, 25, 33, 41, 35, 17, - 1, 3, 7, 15, 13, 39, 61, 187, 461, - 1, 3, 7, 13, 5, 57, 23, 177, 435, - 1, 1, 3, 15, 11, 27, 115, 5, 337, - 1, 3, 7, 3, 15, 63, 61, 171, 339, - 1, 3, 3, 13, 15, 61, 59, 47, 1, - 1, 1, 5, 15, 13, 5, 39, 83, 329, - 1, 1, 5, 5, 5, 27, 25, 39, 301, - 1, 1, 5, 11, 31, 41, 35, 233, 27, - 1, 3, 5, 15, 7, 37, 119, 171, 419, - 1, 3, 5, 5, 3, 29, 21, 189, 417, - 1, 1, 1, 1, 21, 41, 117, 119, 351, - 1, 1, 3, 1, 7, 27, 87, 19, 213, - 1, 1, 1, 1, 17, 7, 97, 217, 477, - 1, 1, 7, 1, 29, 61, 103, 231, 269, - 1, 1, 7, 13, 9, 27, 107, 207, 311, - 1, 1, 7, 5, 25, 21, 107, 179, 423, - 1, 3, 5, 11, 7, 1, 17, 245, 281, - 1, 3, 5, 9, 1, 5, 53, 59, 125, - 1, 1, 7, 1, 31, 57, 71, 245, 125, - 1, 1, 7, 5, 5, 57, 53, 253, 441, - 1, 3, 1, 13, 19, 35, 119, 235, 381, - 1, 3, 1, 7, 19, 59, 115, 33, 361, - 1, 1, 3, 5, 13, 1, 49, 143, 501, - 1, 1, 3, 5, 1, 63, 101, 85, 189, - 1, 1, 5, 11, 27, 63, 13, 131, 5, - 1, 1, 5, 7, 15, 45, 75, 59, 455, 585, - 1, 3, 1, 3, 7, 7, 111, 23, 119, 959, - 1, 3, 3, 9, 11, 41, 109, 163, 161, 879, - 1, 3, 5, 1, 21, 41, 121, 183, 315, 219, - 1, 1, 3, 9, 15, 3, 9, 223, 441, 929, - 1, 1, 7, 9, 3, 5, 93, 57, 253, 457, - 1, 1, 7, 13, 15, 29, 83, 21, 35, 45, - 1, 1, 3, 7, 13, 61, 119, 219, 85, 505, - 1, 1, 3, 3, 17, 13, 35, 197, 291, 109, - 1, 1, 3, 3, 5, 1, 113, 103, 217, 253, - 1, 1, 7, 1, 15, 39, 63, 223, 17, 9, - 1, 3, 7, 1, 17, 29, 67, 103, 495, 383, - 1, 3, 3, 15, 31, 59, 75, 165, 51, 913, - 1, 3, 7, 9, 5, 27, 79, 219, 233, 37, - 1, 3, 5, 15, 1, 11, 15, 211, 417, 811, - 1, 3, 5, 3, 29, 27, 39, 137, 407, 231, - 1, 1, 3, 5, 29, 43, 125, 135, 109, 67, - 1, 1, 1, 5, 11, 39, 107, 159, 323, 381, - 1, 1, 1, 1, 9, 11, 33, 55, 169, 253, - 1, 3, 5, 5, 11, 53, 63, 101, 251, 897, - 1, 3, 7, 1, 25, 15, 83, 119, 53, 157, - 1, 3, 5, 13, 5, 5, 3, 195, 111, 451, - 1, 3, 1, 15, 11, 1, 19, 11, 307, 777, - 1, 3, 7, 11, 5, 5, 17, 231, 345, 981, - 1, 1, 3, 3, 1, 33, 83, 201, 57, 475, - 1, 3, 7, 7, 17, 13, 35, 175, 499, 809, - 1, 1, 5, 3, 3, 17, 103, 119, 499, 865, - 1, 1, 1, 11, 27, 25, 37, 121, 401, 11, - 1, 1, 1, 11, 9, 25, 25, 241, 403, 3, - 1, 1, 1, 1, 11, 1, 39, 163, 231, 573, - 1, 1, 1, 13, 13, 21, 75, 185, 99, 545, - 1, 1, 1, 15, 3, 63, 69, 11, 173, 315, - 1, 3, 5, 15, 11, 3, 95, 49, 123, 765, - 1, 1, 1, 15, 3, 63, 77, 31, 425, 711, - 1, 1, 7, 15, 1, 37, 119, 145, 489, 583, - 1, 3, 5, 15, 3, 49, 117, 211, 165, 323, - 1, 3, 7, 1, 27, 63, 77, 201, 225, 803, - 1, 1, 1, 11, 23, 35, 67, 21, 469, 357, - 1, 1, 7, 7, 9, 7, 25, 237, 237, 571, - 1, 1, 3, 15, 29, 5, 107, 109, 241, 47, - 1, 3, 5, 11, 27, 63, 29, 13, 203, 675, - 1, 1, 3, 9, 9, 11, 103, 179, 449, 263, - 1, 3, 5, 11, 29, 63, 53, 151, 259, 223, - 1, 1, 3, 7, 9, 25, 5, 197, 237, 163, - 1, 3, 7, 13, 5, 57, 67, 193, 147, 241, - 1, 1, 5, 15, 15, 33, 17, 67, 161, 341, - 1, 1, 3, 13, 17, 43, 21, 197, 441, 985, - 1, 3, 1, 5, 15, 33, 33, 193, 305, 829, - 1, 1, 1, 13, 19, 27, 71, 187, 477, 239, - 1, 1, 1, 9, 9, 17, 41, 177, 229, 983, - 1, 3, 5, 9, 15, 45, 97, 205, 43, 767, - 1, 1, 1, 9, 31, 31, 77, 159, 395, 809, - 1, 3, 3, 3, 29, 19, 73, 123, 165, 307, - 1, 3, 1, 7, 5, 11, 77, 227, 355, 403, - 1, 3, 5, 5, 25, 31, 1, 215, 451, 195, - 1, 3, 7, 15, 29, 37, 101, 241, 17, 633, - 1, 1, 5, 1, 11, 3, 107, 137, 489, 5, - 1, 1, 1, 7, 19, 19, 75, 85, 471, 355, - 1, 1, 3, 3, 9, 13, 113, 167, 13, 27, - 1, 3, 5, 11, 21, 3, 89, 205, 377, 307, - 1, 1, 1, 9, 31, 61, 65, 9, 391, 141, 867, - 1, 1, 1, 9, 19, 19, 61, 227, 241, 55, 161, - 1, 1, 1, 11, 1, 19, 7, 233, 463, 171, 1941, - 1, 1, 5, 7, 25, 13, 103, 75, 19, 1021, 1063, - 1, 1, 1, 15, 17, 17, 79, 63, 391, 403, 1221, - 1, 3, 3, 11, 29, 25, 29, 107, 335, 475, 963, - 1, 3, 5, 1, 31, 33, 49, 43, 155, 9, 1285, - 1, 1, 5, 5, 15, 47, 39, 161, 357, 863, 1039, - 1, 3, 7, 15, 1, 39, 47, 109, 427, 393, 1103, - 1, 1, 1, 9, 9, 29, 121, 233, 157, 99, 701, - 1, 1, 1, 7, 1, 29, 75, 121, 439, 109, 993, - 1, 1, 1, 9, 5, 1, 39, 59, 89, 157, 1865, - 1, 1, 5, 1, 3, 37, 89, 93, 143, 533, 175, - 1, 1, 3, 5, 7, 33, 35, 173, 159, 135, 241, - 1, 1, 1, 15, 17, 37, 79, 131, 43, 891, 229, - 1, 1, 1, 1, 1, 35, 121, 177, 397, 1017, 583, - 1, 1, 3, 15, 31, 21, 43, 67, 467, 923, 1473, - 1, 1, 1, 7, 1, 33, 77, 111, 125, 771, 1975, - 1, 3, 7, 13, 1, 51, 113, 139, 245, 573, 503, - 1, 3, 1, 9, 21, 49, 15, 157, 49, 483, 291, - 1, 1, 1, 1, 29, 35, 17, 65, 403, 485, 1603, - 1, 1, 1, 7, 19, 1, 37, 129, 203, 321, 1809, - 1, 3, 7, 15, 15, 9, 5, 77, 29, 485, 581, - 1, 1, 3, 5, 15, 49, 97, 105, 309, 875, 1581, - 1, 3, 5, 1, 5, 19, 63, 35, 165, 399, 1489, - 1, 3, 5, 3, 23, 5, 79, 137, 115, 599, 1127, - 1, 1, 7, 5, 3, 61, 27, 177, 257, 91, 841, - 1, 1, 3, 5, 9, 31, 91, 209, 409, 661, 159, - 1, 3, 1, 15, 23, 39, 23, 195, 245, 203, 947, - 1, 1, 3, 1, 15, 59, 67, 95, 155, 461, 147, - 1, 3, 7, 5, 23, 25, 87, 11, 51, 449, 1631, - 1, 1, 1, 1, 17, 57, 7, 197, 409, 609, 135, - 1, 1, 1, 9, 1, 61, 115, 113, 495, 895, 1595, - 1, 3, 7, 15, 9, 47, 121, 211, 379, 985, 1755, - 1, 3, 1, 3, 7, 57, 27, 231, 339, 325, 1023, - 1, 1, 1, 1, 19, 63, 63, 239, 31, 643, 373, - 1, 3, 1, 11, 19, 9, 7, 171, 21, 691, 215, - 1, 1, 5, 13, 11, 57, 39, 211, 241, 893, 555, - 1, 1, 7, 5, 29, 21, 45, 59, 509, 223, 491, - 1, 1, 7, 9, 15, 61, 97, 75, 127, 779, 839, - 1, 1, 7, 15, 17, 33, 75, 237, 191, 925, 681, - 1, 3, 5, 7, 27, 57, 123, 111, 101, 371, 1129, - 1, 3, 5, 5, 29, 45, 59, 127, 229, 967, 2027, - 1, 1, 1, 1, 17, 7, 23, 199, 241, 455, 135, - 1, 1, 7, 15, 27, 29, 105, 171, 337, 503, 1817, - 1, 1, 3, 7, 21, 35, 61, 71, 405, 647, 2045, - 1, 1, 1, 1, 1, 15, 65, 167, 501, 79, 737, - 1, 1, 5, 1, 3, 49, 27, 189, 341, 615, 1287, - 1, 1, 1, 9, 1, 7, 31, 159, 503, 327, 1613, - 1, 3, 3, 3, 3, 23, 99, 115, 323, 997, 987, - 1, 1, 1, 9, 19, 33, 93, 247, 509, 453, 891, - 1, 1, 3, 1, 13, 19, 35, 153, 161, 633, 445, - 1, 3, 5, 15, 31, 5, 87, 197, 183, 783, 1823, - 1, 1, 7, 5, 19, 63, 69, 221, 129, 231, 1195, - 1, 1, 5, 5, 13, 23, 19, 231, 245, 917, 379, - 1, 3, 1, 15, 19, 43, 27, 223, 171, 413, 125, - 1, 1, 1, 9, 1, 59, 21, 15, 509, 207, 589, - 1, 3, 5, 3, 19, 31, 113, 19, 23, 733, 499, - 1, 1, 7, 1, 19, 51, 101, 165, 47, 925, 1093, - 1, 3, 3, 9, 15, 21, 43, 243, 237, 461, 1361, - 1, 1, 1, 9, 17, 15, 75, 75, 113, 715, 1419, - 1, 1, 7, 13, 17, 1, 99, 15, 347, 721, 1405, - 1, 1, 7, 15, 7, 27, 23, 183, 39, 59, 571, - 1, 3, 5, 9, 7, 43, 35, 165, 463, 567, 859, - 1, 3, 3, 11, 15, 19, 17, 129, 311, 343, 15, - 1, 1, 1, 15, 31, 59, 63, 39, 347, 359, 105, - 1, 1, 1, 15, 5, 43, 87, 241, 109, 61, 685, - 1, 1, 7, 7, 9, 39, 121, 127, 369, 579, 853, - 1, 1, 1, 1, 17, 15, 15, 95, 325, 627, 299, - 1, 1, 3, 13, 31, 53, 85, 111, 289, 811, 1635, - 1, 3, 7, 1, 19, 29, 75, 185, 153, 573, 653, - 1, 3, 7, 1, 29, 31, 55, 91, 249, 247, 1015, - 1, 3, 5, 7, 1, 49, 113, 139, 257, 127, 307, - 1, 3, 5, 9, 15, 15, 123, 105, 105, 225, 1893, - 1, 3, 3, 1, 15, 5, 105, 249, 73, 709, 1557, - 1, 1, 1, 9, 17, 31, 113, 73, 65, 701, 1439, - 1, 3, 5, 15, 13, 21, 117, 131, 243, 859, 323, - 1, 1, 1, 9, 19, 15, 69, 149, 89, 681, 515, - 1, 1, 1, 5, 29, 13, 21, 97, 301, 27, 967, - 1, 1, 3, 3, 15, 45, 107, 227, 495, 769, 1935, - 1, 1, 1, 11, 5, 27, 41, 173, 261, 703, 1349, - 1, 3, 3, 3, 11, 35, 97, 43, 501, 563, 1331, - 1, 1, 1, 7, 1, 17, 87, 17, 429, 245, 1941, - 1, 1, 7, 15, 29, 13, 1, 175, 425, 233, 797, - 1, 1, 3, 11, 21, 57, 49, 49, 163, 685, 701, - 1, 3, 3, 7, 11, 45, 107, 111, 379, 703, 1403, - 1, 1, 7, 3, 21, 7, 117, 49, 469, 37, 775, - 1, 1, 5, 15, 31, 63, 101, 77, 507, 489, 1955, - 1, 3, 3, 11, 19, 21, 101, 255, 203, 673, 665, - 1, 3, 3, 15, 17, 47, 125, 187, 271, 899, 2003, - 1, 1, 7, 7, 1, 35, 13, 235, 5, 337, 905, - 1, 3, 1, 15, 1, 43, 1, 27, 37, 695, 1429, - 1, 3, 1, 11, 21, 27, 93, 161, 299, 665, 495, - 1, 3, 3, 15, 3, 1, 81, 111, 105, 547, 897, - 1, 3, 5, 1, 3, 53, 97, 253, 401, 827, 1467, - 1, 1, 1, 5, 19, 59, 105, 125, 271, 351, 719, - 1, 3, 5, 13, 7, 11, 91, 41, 441, 759, 1827, - 1, 3, 7, 11, 29, 61, 61, 23, 307, 863, 363, - 1, 1, 7, 1, 15, 35, 29, 133, 415, 473, 1737, - 1, 1, 1, 13, 7, 33, 35, 225, 117, 681, 1545, - 1, 1, 1, 3, 5, 41, 83, 247, 13, 373, 1091, - 1, 3, 1, 13, 25, 61, 71, 217, 233, 313, 547, - 1, 3, 1, 7, 3, 29, 3, 49, 93, 465, 15, - 1, 1, 1, 9, 17, 61, 99, 163, 129, 485, 1087, - 1, 1, 1, 9, 9, 33, 31, 163, 145, 649, 253, - 1, 1, 1, 1, 17, 63, 43, 235, 287, 111, 567, - 1, 3, 5, 13, 29, 7, 11, 69, 153, 127, 449, - 1, 1, 5, 9, 11, 21, 15, 189, 431, 493, 1219, - 1, 1, 1, 15, 19, 5, 47, 91, 399, 293, 1743, - 1, 3, 3, 11, 29, 53, 53, 225, 409, 303, 333, - 1, 1, 1, 15, 31, 31, 21, 81, 147, 287, 1753, - 1, 3, 5, 5, 5, 63, 35, 125, 41, 687, 1793, - 1, 1, 1, 9, 19, 59, 107, 219, 455, 971, 297, - 1, 1, 3, 5, 3, 51, 121, 31, 245, 105, 1311, - 1, 3, 1, 5, 5, 57, 75, 107, 161, 431, 1693, - 1, 3, 1, 3, 19, 53, 27, 31, 191, 565, 1015, - 1, 3, 5, 13, 9, 41, 35, 249, 287, 49, 123, - 1, 1, 5, 7, 27, 17, 21, 3, 151, 885, 1165, - 1, 1, 7, 1, 15, 17, 65, 139, 427, 339, 1171, - 1, 1, 1, 5, 23, 5, 9, 89, 321, 907, 391, - 1, 1, 7, 9, 15, 1, 77, 71, 87, 701, 917, - 1, 1, 7, 1, 17, 37, 115, 127, 469, 779, 1543, - 1, 3, 7, 3, 5, 61, 15, 37, 301, 951, 1437, - 1, 1, 1, 13, 9, 51, 127, 145, 229, 55, 1567, - 1, 3, 7, 15, 19, 47, 53, 153, 295, 47, 1337, - 1, 3, 3, 5, 11, 31, 29, 133, 327, 287, 507, - 1, 1, 7, 7, 25, 31, 37, 199, 25, 927, 1317, - 1, 1, 7, 9, 3, 39, 127, 167, 345, 467, 759, - 1, 1, 1, 1, 31, 21, 15, 101, 293, 787, 1025, - 1, 1, 5, 3, 11, 41, 105, 109, 149, 837, 1813, - 1, 1, 3, 5, 29, 13, 19, 97, 309, 901, 753, - 1, 1, 7, 1, 19, 17, 31, 39, 173, 361, 1177, - 1, 3, 3, 3, 3, 41, 81, 7, 341, 491, 43, - 1, 1, 7, 7, 31, 35, 29, 77, 11, 335, 1275, - 1, 3, 3, 15, 17, 45, 19, 63, 151, 849, 129, - 1, 1, 7, 5, 7, 13, 47, 73, 79, 31, 499, - 1, 3, 1, 11, 1, 41, 59, 151, 247, 115, 1295, - 1, 1, 1, 9, 31, 37, 73, 23, 295, 483, 179, - 1, 3, 1, 15, 13, 63, 81, 27, 169, 825, 2037, - 1, 3, 5, 15, 7, 11, 73, 1, 451, 101, 2039, - 1, 3, 5, 3, 13, 53, 31, 137, 173, 319, 1521, - 1, 3, 1, 3, 29, 1, 73, 227, 377, 337, 1189, - 1, 3, 3, 13, 27, 9, 31, 101, 229, 165, 1983, - 1, 3, 1, 13, 13, 19, 19, 111, 319, 421, 223, - 1, 1, 7, 15, 25, 37, 61, 55, 359, 255, 1955, - 1, 1, 5, 13, 17, 43, 49, 215, 383, 915, 51, - 1, 1, 3, 1, 3, 7, 13, 119, 155, 585, 967, - 1, 3, 1, 13, 1, 63, 125, 21, 103, 287, 457, - 1, 1, 7, 1, 31, 17, 125, 137, 345, 379, 1925, - 1, 1, 3, 5, 5, 25, 119, 153, 455, 271, 2023, - 1, 1, 7, 9, 9, 37, 115, 47, 5, 255, 917, - 1, 3, 5, 3, 31, 21, 75, 203, 489, 593, 1, - 1, 3, 7, 15, 19, 63, 123, 153, 135, 977, 1875, - 1, 1, 1, 1, 5, 59, 31, 25, 127, 209, 745, - 1, 1, 1, 1, 19, 45, 67, 159, 301, 199, 535, - 1, 1, 7, 1, 31, 17, 19, 225, 369, 125, 421, - 1, 3, 3, 11, 7, 59, 115, 197, 459, 469, 1055, - 1, 3, 1, 3, 27, 45, 35, 131, 349, 101, 411, - 1, 3, 7, 11, 9, 3, 67, 145, 299, 253, 1339, - 1, 3, 3, 11, 9, 37, 123, 229, 273, 269, 515, - 1, 3, 7, 15, 11, 25, 75, 5, 367, 217, 951, - 1, 1, 3, 7, 9, 23, 63, 237, 385, 159, 1273, - 1, 1, 5, 11, 23, 5, 55, 193, 109, 865, 663, - 1, 1, 7, 15, 1, 57, 17, 141, 51, 217, 1259, - 1, 1, 3, 3, 15, 7, 89, 233, 71, 329, 203, - 1, 3, 7, 11, 11, 1, 19, 155, 89, 437, 573, - 1, 3, 1, 9, 27, 61, 47, 109, 161, 913, 1681, - 1, 1, 7, 15, 1, 33, 19, 15, 23, 913, 989, - 1, 3, 1, 1, 25, 39, 119, 193, 13, 571, 157, - 1, 1, 7, 13, 9, 55, 59, 147, 361, 935, 515, - 1, 1, 1, 9, 7, 59, 67, 117, 71, 855, 1493, - 1, 3, 1, 3, 13, 19, 57, 141, 305, 275, 1079, - 1, 1, 1, 9, 17, 61, 33, 7, 43, 931, 781, - 1, 1, 3, 1, 11, 17, 21, 97, 295, 277, 1721, - 1, 3, 1, 13, 15, 43, 11, 241, 147, 391, 1641, - 1, 1, 1, 1, 1, 19, 37, 21, 255, 263, 1571, - 1, 1, 3, 3, 23, 59, 89, 17, 475, 303, 757, 543, - 1, 3, 3, 9, 11, 55, 35, 159, 139, 203, 1531, 1825, - 1, 1, 5, 3, 17, 53, 51, 241, 269, 949, 1373, 325, - 1, 3, 7, 7, 5, 29, 91, 149, 239, 193, 1951, 2675, - 1, 3, 5, 1, 27, 33, 69, 11, 51, 371, 833, 2685, - 1, 1, 1, 15, 1, 17, 35, 57, 171, 1007, 449, 367, - 1, 1, 1, 7, 25, 61, 73, 219, 379, 53, 589, 4065, - 1, 3, 5, 13, 21, 29, 45, 19, 163, 169, 147, 597, - 1, 1, 5, 11, 21, 27, 7, 17, 237, 591, 255, 1235, - 1, 1, 7, 7, 17, 41, 69, 237, 397, 173, 1229, 2341, - 1, 1, 3, 1, 1, 33, 125, 47, 11, 783, 1323, 2469, - 1, 3, 1, 11, 3, 39, 35, 133, 153, 55, 1171, 3165, - 1, 1, 5, 11, 27, 23, 103, 245, 375, 753, 477, 2165, - 1, 3, 1, 15, 15, 49, 127, 223, 387, 771, 1719, 1465, - 1, 1, 1, 9, 11, 9, 17, 185, 239, 899, 1273, 3961, - 1, 1, 3, 13, 11, 51, 73, 81, 389, 647, 1767, 1215, - 1, 3, 5, 15, 19, 9, 69, 35, 349, 977, 1603, 1435, - 1, 1, 1, 1, 19, 59, 123, 37, 41, 961, 181, 1275, - 1, 1, 1, 1, 31, 29, 37, 71, 205, 947, 115, 3017, - 1, 1, 7, 15, 5, 37, 101, 169, 221, 245, 687, 195, - 1, 1, 1, 1, 19, 9, 125, 157, 119, 283, 1721, 743, - 1, 1, 7, 3, 1, 7, 61, 71, 119, 257, 1227, 2893, - 1, 3, 3, 3, 25, 41, 25, 225, 31, 57, 925, 2139, - 0 -}; - - -/* - * coefficients of the free direction integers as given in - * "Monte Carlo Methods in Finance", by Peter J�ckel, section 8.3 - */ - -const unsigned long rk_sobol_Jdirections[] = { - 1, - 1, 1, - 1, 3, 7, - 1, 1, 5, - 1, 3, 1, 1, - 1, 1, 3, 7, - 1, 3, 3, 9, 9, - 1, 3, 7, 7, 21, - 1, 1, 5, 11, 27, - 1, 1, 7, 3, 29, - 1, 3, 7, 13, 3, - 1, 3, 5, 1, 15, - 1, 1, 1, 9, 23, 37, - 1, 1, 3, 13, 11, 7, - 1, 3, 3, 5, 19, 33, - 1, 1, 7, 13, 25, 5, - 1, 1, 1, 3, 13, 39, - 1, 3, 5, 11, 7, 11, - 1, 3, 1, 7, 3, 23, 79, - 1, 3, 1, 15, 17, 63, 13, - 1, 3, 3, 3, 25, 17, 115, - 1, 3, 7, 9, 31, 29, 17, - 1, 1, 3, 15, 29, 15, 41, - 1, 3, 1, 9, 5, 21, 119, - 1, 1, 5, 5, 1, 27, 33, - 1, 1, 3, 1, 23, 13, 75, - 1, 1, 7, 7, 19, 25, 105, - 1, 3, 5, 5, 21, 9, 7, - 1, 1, 1, 15, 5, 49, 59, - 1, 3, 5, 15, 17, 19, 21, - 1, 1, 7, 11, 13, 29, 3, - 0 -}; - -/* - * 0 terminated list of primitive polynomials to speed up initialization - * All polynomials up to degree 13 (ie. 1111 polynomials) - */ -static const unsigned long rk_sobol_primitive_polynomials[] = { - 0x1UL, 0x3UL, 0x7UL, 0xBUL, 0xDUL, 0x13UL, 0x19UL, 0x25UL, 0x29UL, - 0x2FUL, 0x37UL, 0x3BUL, 0x3DUL, 0x43UL, 0x5BUL, 0x61UL, 0x67UL, 0x6DUL, - 0x73UL, 0x83UL, 0x89UL, 0x8FUL, 0x91UL, 0x9DUL, 0xA7UL, 0xABUL, 0xB9UL, - 0xBFUL, 0xC1UL, 0xCBUL, 0xD3UL, 0xD5UL, 0xE5UL, 0xEFUL, 0xF1UL, 0xF7UL, - 0xFDUL, 0x11DUL, 0x12BUL, 0x12DUL, 0x14DUL, 0x15FUL, 0x163UL, 0x165UL, - 0x169UL, 0x171UL, 0x187UL, 0x18DUL, 0x1A9UL, 0x1C3UL, 0x1CFUL, 0x1E7UL, - 0x1F5UL, 0x211UL, 0x21BUL, 0x221UL, 0x22DUL, 0x233UL, 0x259UL, 0x25FUL, - 0x269UL, 0x26FUL, 0x277UL, 0x27DUL, 0x287UL, 0x295UL, 0x2A3UL, 0x2A5UL, - 0x2AFUL, 0x2B7UL, 0x2BDUL, 0x2CFUL, 0x2D1UL, 0x2DBUL, 0x2F5UL, 0x2F9UL, - 0x313UL, 0x315UL, 0x31FUL, 0x323UL, 0x331UL, 0x33BUL, 0x34FUL, 0x35BUL, - 0x361UL, 0x36BUL, 0x36DUL, 0x373UL, 0x37FUL, 0x385UL, 0x38FUL, 0x3B5UL, - 0x3B9UL, 0x3C7UL, 0x3CBUL, 0x3CDUL, 0x3D5UL, 0x3D9UL, 0x3E3UL, 0x3E9UL, - 0x3FBUL, 0x409UL, 0x41BUL, 0x427UL, 0x42DUL, 0x465UL, 0x46FUL, 0x481UL, - 0x48BUL, 0x4C5UL, 0x4D7UL, 0x4E7UL, 0x4F3UL, 0x4FFUL, 0x50DUL, 0x519UL, - 0x523UL, 0x531UL, 0x53DUL, 0x543UL, 0x557UL, 0x56BUL, 0x585UL, 0x58FUL, - 0x597UL, 0x5A1UL, 0x5C7UL, 0x5E5UL, 0x5F7UL, 0x5FBUL, 0x613UL, 0x615UL, - 0x625UL, 0x637UL, 0x643UL, 0x64FUL, 0x65BUL, 0x679UL, 0x67FUL, 0x689UL, - 0x6B5UL, 0x6C1UL, 0x6D3UL, 0x6DFUL, 0x6FDUL, 0x717UL, 0x71DUL, 0x721UL, - 0x739UL, 0x747UL, 0x74DUL, 0x755UL, 0x759UL, 0x763UL, 0x77DUL, 0x78DUL, - 0x793UL, 0x7B1UL, 0x7DBUL, 0x7F3UL, 0x7F9UL, 0x805UL, 0x817UL, 0x82BUL, - 0x82DUL, 0x847UL, 0x863UL, 0x865UL, 0x871UL, 0x87BUL, 0x88DUL, 0x895UL, - 0x89FUL, 0x8A9UL, 0x8B1UL, 0x8CFUL, 0x8D1UL, 0x8E1UL, 0x8E7UL, 0x8EBUL, - 0x8F5UL, 0x90DUL, 0x913UL, 0x925UL, 0x929UL, 0x93BUL, 0x93DUL, 0x945UL, - 0x949UL, 0x951UL, 0x95BUL, 0x973UL, 0x975UL, 0x97FUL, 0x983UL, 0x98FUL, - 0x9ABUL, 0x9ADUL, 0x9B9UL, 0x9C7UL, 0x9D9UL, 0x9E5UL, 0x9F7UL, 0xA01UL, - 0xA07UL, 0xA13UL, 0xA15UL, 0xA29UL, 0xA49UL, 0xA61UL, 0xA6DUL, 0xA79UL, - 0xA7FUL, 0xA85UL, 0xA91UL, 0xA9DUL, 0xAA7UL, 0xAABUL, 0xAB3UL, 0xAB5UL, - 0xAD5UL, 0xADFUL, 0xAE9UL, 0xAEFUL, 0xAF1UL, 0xAFBUL, 0xB03UL, 0xB09UL, - 0xB11UL, 0xB33UL, 0xB3FUL, 0xB41UL, 0xB4BUL, 0xB59UL, 0xB5FUL, 0xB65UL, - 0xB6FUL, 0xB7DUL, 0xB87UL, 0xB8BUL, 0xB93UL, 0xB95UL, 0xBAFUL, 0xBB7UL, - 0xBBDUL, 0xBC9UL, 0xBDBUL, 0xBDDUL, 0xBE7UL, 0xBEDUL, 0xC0BUL, 0xC0DUL, - 0xC19UL, 0xC1FUL, 0xC57UL, 0xC61UL, 0xC6BUL, 0xC73UL, 0xC85UL, 0xC89UL, - 0xC97UL, 0xC9BUL, 0xC9DUL, 0xCB3UL, 0xCBFUL, 0xCC7UL, 0xCCDUL, 0xCD3UL, - 0xCD5UL, 0xCE3UL, 0xCE9UL, 0xCF7UL, 0xD03UL, 0xD0FUL, 0xD1DUL, 0xD27UL, - 0xD2DUL, 0xD41UL, 0xD47UL, 0xD55UL, 0xD59UL, 0xD63UL, 0xD6FUL, 0xD71UL, - 0xD93UL, 0xD9FUL, 0xDA9UL, 0xDBBUL, 0xDBDUL, 0xDC9UL, 0xDD7UL, 0xDDBUL, - 0xDE1UL, 0xDE7UL, 0xDF5UL, 0xE05UL, 0xE1DUL, 0xE21UL, 0xE27UL, 0xE2BUL, - 0xE33UL, 0xE39UL, 0xE47UL, 0xE4BUL, 0xE55UL, 0xE5FUL, 0xE71UL, 0xE7BUL, - 0xE7DUL, 0xE81UL, 0xE93UL, 0xE9FUL, 0xEA3UL, 0xEBBUL, 0xECFUL, 0xEDDUL, - 0xEF3UL, 0xEF9UL, 0xF0BUL, 0xF19UL, 0xF31UL, 0xF37UL, 0xF5DUL, 0xF6BUL, - 0xF6DUL, 0xF75UL, 0xF83UL, 0xF91UL, 0xF97UL, 0xF9BUL, 0xFA7UL, 0xFADUL, - 0xFB5UL, 0xFCDUL, 0xFD3UL, 0xFE5UL, 0xFE9UL, 0x1053UL, 0x1069UL, - 0x107BUL, 0x107DUL, 0x1099UL, 0x10D1UL, 0x10EBUL, 0x1107UL, 0x111FUL, - 0x1123UL, 0x113BUL, 0x114FUL, 0x1157UL, 0x1161UL, 0x116BUL, 0x1185UL, - 0x11B3UL, 0x11D9UL, 0x11DFUL, 0x120DUL, 0x1237UL, 0x123DUL, 0x1267UL, - 0x1273UL, 0x127FUL, 0x12B9UL, 0x12C1UL, 0x12CBUL, 0x130FUL, 0x131DUL, - 0x1321UL, 0x1339UL, 0x133FUL, 0x134DUL, 0x1371UL, 0x1399UL, 0x13A3UL, - 0x13A9UL, 0x1407UL, 0x1431UL, 0x1437UL, 0x144FUL, 0x145DUL, 0x1467UL, - 0x1475UL, 0x14A7UL, 0x14ADUL, 0x14D3UL, 0x150FUL, 0x151DUL, 0x154DUL, - 0x1593UL, 0x15C5UL, 0x15D7UL, 0x15DDUL, 0x15EBUL, 0x1609UL, 0x1647UL, - 0x1655UL, 0x1659UL, 0x16A5UL, 0x16BDUL, 0x1715UL, 0x1719UL, 0x1743UL, - 0x1745UL, 0x1775UL, 0x1789UL, 0x17ADUL, 0x17B3UL, 0x17BFUL, 0x17C1UL, - 0x1857UL, 0x185DUL, 0x1891UL, 0x1897UL, 0x18B9UL, 0x18EFUL, 0x191BUL, - 0x1935UL, 0x1941UL, 0x1965UL, 0x197BUL, 0x198BUL, 0x19B1UL, 0x19BDUL, - 0x19C9UL, 0x19CFUL, 0x19E7UL, 0x1A1BUL, 0x1A2BUL, 0x1A33UL, 0x1A69UL, - 0x1A8BUL, 0x1AD1UL, 0x1AE1UL, 0x1AF5UL, 0x1B0BUL, 0x1B13UL, 0x1B1FUL, - 0x1B57UL, 0x1B91UL, 0x1BA7UL, 0x1BBFUL, 0x1BC1UL, 0x1BD3UL, 0x1C05UL, - 0x1C11UL, 0x1C17UL, 0x1C27UL, 0x1C4DUL, 0x1C87UL, 0x1C9FUL, 0x1CA5UL, - 0x1CBBUL, 0x1CC5UL, 0x1CC9UL, 0x1CCFUL, 0x1CF3UL, 0x1D07UL, 0x1D23UL, - 0x1D43UL, 0x1D51UL, 0x1D5BUL, 0x1D75UL, 0x1D85UL, 0x1D89UL, 0x1E15UL, - 0x1E19UL, 0x1E2FUL, 0x1E45UL, 0x1E51UL, 0x1E67UL, 0x1E73UL, 0x1E8FUL, - 0x1EE3UL, 0x1F11UL, 0x1F1BUL, 0x1F27UL, 0x1F71UL, 0x1F99UL, 0x1FBBUL, - 0x1FBDUL, 0x1FC9UL, 0x201BUL, 0x2027UL, 0x2035UL, 0x2053UL, 0x2065UL, - 0x206FUL, 0x208BUL, 0x208DUL, 0x209FUL, 0x20A5UL, 0x20AFUL, 0x20BBUL, - 0x20BDUL, 0x20C3UL, 0x20C9UL, 0x20E1UL, 0x20F3UL, 0x210DUL, 0x2115UL, - 0x2129UL, 0x212FUL, 0x213BUL, 0x2143UL, 0x2167UL, 0x216BUL, 0x2179UL, - 0x2189UL, 0x2197UL, 0x219DUL, 0x21BFUL, 0x21C1UL, 0x21C7UL, 0x21CDUL, - 0x21DFUL, 0x21E3UL, 0x21F1UL, 0x21FBUL, 0x2219UL, 0x2225UL, 0x2237UL, - 0x223DUL, 0x2243UL, 0x225BUL, 0x225DUL, 0x2279UL, 0x227FUL, 0x2289UL, - 0x2297UL, 0x229BUL, 0x22B3UL, 0x22BFUL, 0x22CDUL, 0x22EFUL, 0x22F7UL, - 0x22FBUL, 0x2305UL, 0x2327UL, 0x232BUL, 0x2347UL, 0x2355UL, 0x2359UL, - 0x236FUL, 0x2371UL, 0x237DUL, 0x2387UL, 0x238DUL, 0x2395UL, 0x23A3UL, - 0x23A9UL, 0x23B1UL, 0x23B7UL, 0x23BBUL, 0x23E1UL, 0x23EDUL, 0x23F9UL, - 0x240BUL, 0x2413UL, 0x241FUL, 0x2425UL, 0x2429UL, 0x243DUL, 0x2451UL, - 0x2457UL, 0x2461UL, 0x246DUL, 0x247FUL, 0x2483UL, 0x249BUL, 0x249DUL, - 0x24B5UL, 0x24BFUL, 0x24C1UL, 0x24C7UL, 0x24CBUL, 0x24E3UL, 0x2509UL, - 0x2517UL, 0x251DUL, 0x2521UL, 0x252DUL, 0x2539UL, 0x2553UL, 0x2555UL, - 0x2563UL, 0x2571UL, 0x2577UL, 0x2587UL, 0x258BUL, 0x2595UL, 0x2599UL, - 0x259FUL, 0x25AFUL, 0x25BDUL, 0x25C5UL, 0x25CFUL, 0x25D7UL, 0x25EBUL, - 0x2603UL, 0x2605UL, 0x2611UL, 0x262DUL, 0x263FUL, 0x264BUL, 0x2653UL, - 0x2659UL, 0x2669UL, 0x2677UL, 0x267BUL, 0x2687UL, 0x2693UL, 0x2699UL, - 0x26B1UL, 0x26B7UL, 0x26BDUL, 0x26C3UL, 0x26EBUL, 0x26F5UL, 0x2713UL, - 0x2729UL, 0x273BUL, 0x274FUL, 0x2757UL, 0x275DUL, 0x276BUL, 0x2773UL, - 0x2779UL, 0x2783UL, 0x2791UL, 0x27A1UL, 0x27B9UL, 0x27C7UL, 0x27CBUL, - 0x27DFUL, 0x27EFUL, 0x27F1UL, 0x2807UL, 0x2819UL, 0x281FUL, 0x2823UL, - 0x2831UL, 0x283BUL, 0x283DUL, 0x2845UL, 0x2867UL, 0x2875UL, 0x2885UL, - 0x28ABUL, 0x28ADUL, 0x28BFUL, 0x28CDUL, 0x28D5UL, 0x28DFUL, 0x28E3UL, - 0x28E9UL, 0x28FBUL, 0x2909UL, 0x290FUL, 0x2911UL, 0x291BUL, 0x292BUL, - 0x2935UL, 0x293FUL, 0x2941UL, 0x294BUL, 0x2955UL, 0x2977UL, 0x297DUL, - 0x2981UL, 0x2993UL, 0x299FUL, 0x29AFUL, 0x29B7UL, 0x29BDUL, 0x29C3UL, - 0x29D7UL, 0x29F3UL, 0x29F5UL, 0x2A03UL, 0x2A0FUL, 0x2A1DUL, 0x2A21UL, - 0x2A33UL, 0x2A35UL, 0x2A4DUL, 0x2A69UL, 0x2A6FUL, 0x2A71UL, 0x2A7BUL, - 0x2A7DUL, 0x2AA5UL, 0x2AA9UL, 0x2AB1UL, 0x2AC5UL, 0x2AD7UL, 0x2ADBUL, - 0x2AEBUL, 0x2AF3UL, 0x2B01UL, 0x2B15UL, 0x2B23UL, 0x2B25UL, 0x2B2FUL, - 0x2B37UL, 0x2B43UL, 0x2B49UL, 0x2B6DUL, 0x2B7FUL, 0x2B85UL, 0x2B97UL, - 0x2B9BUL, 0x2BADUL, 0x2BB3UL, 0x2BD9UL, 0x2BE5UL, 0x2BFDUL, 0x2C0FUL, - 0x2C21UL, 0x2C2BUL, 0x2C2DUL, 0x2C3FUL, 0x2C41UL, 0x2C4DUL, 0x2C71UL, - 0x2C8BUL, 0x2C8DUL, 0x2C95UL, 0x2CA3UL, 0x2CAFUL, 0x2CBDUL, 0x2CC5UL, - 0x2CD1UL, 0x2CD7UL, 0x2CE1UL, 0x2CE7UL, 0x2CEBUL, 0x2D0DUL, 0x2D19UL, - 0x2D29UL, 0x2D2FUL, 0x2D37UL, 0x2D3BUL, 0x2D45UL, 0x2D5BUL, 0x2D67UL, - 0x2D75UL, 0x2D89UL, 0x2D8FUL, 0x2DA7UL, 0x2DABUL, 0x2DB5UL, 0x2DE3UL, - 0x2DF1UL, 0x2DFDUL, 0x2E07UL, 0x2E13UL, 0x2E15UL, 0x2E29UL, 0x2E49UL, - 0x2E4FUL, 0x2E5BUL, 0x2E5DUL, 0x2E61UL, 0x2E6BUL, 0x2E8FUL, 0x2E91UL, - 0x2E97UL, 0x2E9DUL, 0x2EABUL, 0x2EB3UL, 0x2EB9UL, 0x2EDFUL, 0x2EFBUL, - 0x2EFDUL, 0x2F05UL, 0x2F09UL, 0x2F11UL, 0x2F17UL, 0x2F3FUL, 0x2F41UL, - 0x2F4BUL, 0x2F4DUL, 0x2F59UL, 0x2F5FUL, 0x2F65UL, 0x2F69UL, 0x2F95UL, - 0x2FA5UL, 0x2FAFUL, 0x2FB1UL, 0x2FCFUL, 0x2FDDUL, 0x2FE7UL, 0x2FEDUL, - 0x2FF5UL, 0x2FFFUL, 0x3007UL, 0x3015UL, 0x3019UL, 0x302FUL, 0x3049UL, - 0x304FUL, 0x3067UL, 0x3079UL, 0x307FUL, 0x3091UL, 0x30A1UL, 0x30B5UL, - 0x30BFUL, 0x30C1UL, 0x30D3UL, 0x30D9UL, 0x30E5UL, 0x30EFUL, 0x3105UL, - 0x310FUL, 0x3135UL, 0x3147UL, 0x314DUL, 0x315FUL, 0x3163UL, 0x3171UL, - 0x317BUL, 0x31A3UL, 0x31A9UL, 0x31B7UL, 0x31C5UL, 0x31C9UL, 0x31DBUL, - 0x31E1UL, 0x31EBUL, 0x31EDUL, 0x31F3UL, 0x31FFUL, 0x3209UL, 0x320FUL, - 0x321DUL, 0x3227UL, 0x3239UL, 0x324BUL, 0x3253UL, 0x3259UL, 0x3265UL, - 0x3281UL, 0x3293UL, 0x3299UL, 0x329FUL, 0x32A9UL, 0x32B7UL, 0x32BBUL, - 0x32C3UL, 0x32D7UL, 0x32DBUL, 0x32E7UL, 0x3307UL, 0x3315UL, 0x332FUL, - 0x3351UL, 0x335DUL, 0x3375UL, 0x3397UL, 0x339BUL, 0x33ABUL, 0x33B9UL, - 0x33C1UL, 0x33C7UL, 0x33D5UL, 0x33E3UL, 0x33E5UL, 0x33F7UL, 0x33FBUL, - 0x3409UL, 0x341BUL, 0x3427UL, 0x3441UL, 0x344DUL, 0x345FUL, 0x3469UL, - 0x3477UL, 0x347BUL, 0x3487UL, 0x3493UL, 0x3499UL, 0x34A5UL, 0x34BDUL, - 0x34C9UL, 0x34DBUL, 0x34E7UL, 0x34F9UL, 0x350DUL, 0x351FUL, 0x3525UL, - 0x3531UL, 0x3537UL, 0x3545UL, 0x354FUL, 0x355DUL, 0x356DUL, 0x3573UL, - 0x357FUL, 0x359DUL, 0x35A1UL, 0x35B9UL, 0x35CDUL, 0x35D5UL, 0x35D9UL, - 0x35E3UL, 0x35E9UL, 0x35EFUL, 0x3601UL, 0x360BUL, 0x361FUL, 0x3625UL, - 0x362FUL, 0x363BUL, 0x3649UL, 0x3651UL, 0x365BUL, 0x3673UL, 0x3675UL, - 0x3691UL, 0x369BUL, 0x369DUL, 0x36ADUL, 0x36CBUL, 0x36D3UL, 0x36D5UL, - 0x36E3UL, 0x36EFUL, 0x3705UL, 0x370FUL, 0x371BUL, 0x3721UL, 0x372DUL, - 0x3739UL, 0x3741UL, 0x3747UL, 0x3753UL, 0x3771UL, 0x3777UL, 0x378BUL, - 0x3795UL, 0x3799UL, 0x37A3UL, 0x37C5UL, 0x37CFUL, 0x37D1UL, 0x37D7UL, - 0x37DDUL, 0x37E1UL, 0x37F3UL, 0x3803UL, 0x3805UL, 0x3817UL, 0x381DUL, - 0x3827UL, 0x3833UL, 0x384BUL, 0x3859UL, 0x3869UL, 0x3871UL, 0x38A3UL, - 0x38B1UL, 0x38BBUL, 0x38C9UL, 0x38CFUL, 0x38E1UL, 0x38F3UL, 0x38F9UL, - 0x3901UL, 0x3907UL, 0x390BUL, 0x3913UL, 0x3931UL, 0x394FUL, 0x3967UL, - 0x396DUL, 0x3983UL, 0x3985UL, 0x3997UL, 0x39A1UL, 0x39A7UL, 0x39ADUL, - 0x39CBUL, 0x39CDUL, 0x39D3UL, 0x39EFUL, 0x39F7UL, 0x39FDUL, 0x3A07UL, - 0x3A29UL, 0x3A2FUL, 0x3A3DUL, 0x3A51UL, 0x3A5DUL, 0x3A61UL, 0x3A67UL, - 0x3A73UL, 0x3A75UL, 0x3A89UL, 0x3AB9UL, 0x3ABFUL, 0x3ACDUL, 0x3AD3UL, - 0x3AD5UL, 0x3ADFUL, 0x3AE5UL, 0x3AE9UL, 0x3AFBUL, 0x3B11UL, 0x3B2BUL, - 0x3B2DUL, 0x3B35UL, 0x3B3FUL, 0x3B53UL, 0x3B59UL, 0x3B63UL, 0x3B65UL, - 0x3B6FUL, 0x3B71UL, 0x3B77UL, 0x3B8BUL, 0x3B99UL, 0x3BA5UL, 0x3BA9UL, - 0x3BB7UL, 0x3BBBUL, 0x3BD1UL, 0x3BE7UL, 0x3BF3UL, 0x3BFFUL, 0x3C0DUL, - 0x3C13UL, 0x3C15UL, 0x3C1FUL, 0x3C23UL, 0x3C25UL, 0x3C3BUL, 0x3C4FUL, - 0x3C5DUL, 0x3C6DUL, 0x3C83UL, 0x3C8FUL, 0x3C9DUL, 0x3CA7UL, 0x3CABUL, - 0x3CB9UL, 0x3CC7UL, 0x3CE9UL, 0x3CFBUL, 0x3CFDUL, 0x3D03UL, 0x3D17UL, - 0x3D1BUL, 0x3D21UL, 0x3D2DUL, 0x3D33UL, 0x3D35UL, 0x3D41UL, 0x3D4DUL, - 0x3D65UL, 0x3D69UL, 0x3D7DUL, 0x3D81UL, 0x3D95UL, 0x3DB1UL, 0x3DB7UL, - 0x3DC3UL, 0x3DD1UL, 0x3DDBUL, 0x3DE7UL, 0x3DEBUL, 0x3DF9UL, 0x3E05UL, - 0x3E09UL, 0x3E0FUL, 0x3E1BUL, 0x3E2BUL, 0x3E3FUL, 0x3E41UL, 0x3E53UL, - 0x3E65UL, 0x3E69UL, 0x3E8BUL, 0x3EA3UL, 0x3EBDUL, 0x3EC5UL, 0x3ED7UL, - 0x3EDDUL, 0x3EE1UL, 0x3EF9UL, 0x3F0DUL, 0x3F19UL, 0x3F1FUL, 0x3F25UL, - 0x3F37UL, 0x3F3DUL, 0x3F43UL, 0x3F45UL, 0x3F49UL, 0x3F51UL, 0x3F57UL, - 0x3F61UL, 0x3F83UL, 0x3F89UL, 0x3F91UL, 0x3FABUL, 0x3FB5UL, 0x3FE3UL, - 0x3FF7UL, 0x3FFDUL, - 0UL -}; - -rk_sobol_error rk_sobol_init(size_t dimension, rk_sobol_state *s, - rk_state *rs_dir, const unsigned long *directions, - const unsigned long *polynomials) -{ - rk_state rs_dir_temp; - int j, l, degree = 0, last_degree = 0, ooord = 0; - size_t k, cdir = 0, cpol = 0; - unsigned long polynomial = 1, rev = 0, last = 0; - - if (dimension == 0) - return RK_SOBOL_EINVAL; - - if (polynomials == NULL) - polynomials = rk_sobol_primitive_polynomials; - - /* Allocate the structure */ - s->direction = NULL; s->numerator = NULL; - s->direction = malloc(sizeof(*(s->direction))*dimension*LONG_BIT); - s->numerator = malloc(sizeof(*(s->numerator))*dimension); - if (!s->direction | !s->numerator) - { - if (!s->direction) free(s->direction); - if (!s->numerator) free(s->numerator); - return RK_SOBOL_ENOMEM; - } - - /* Initialize directions */ - /* Degree 0 */ - for (j = degree; j < LONG_BIT; j++) - s->direction[j*dimension] = 1UL << (LONG_BIT-j-1); - - /* Skip unused first polynomial */ - if (polynomials[cpol]) - cpol++; - - /* Degree >0 */ - for (k = 1; k < dimension; k++) - { - unsigned long temp; - - /* Find a new primitive polynomial */ - if (polynomials[cpol]) - polynomial = polynomials[cpol++]; - else if (rev) - { - /* We are generating polynomials out of order: - use the reverse of the previous polynomial */ - last = polynomial; - polynomial = rev; - rev = 0; - } - else - { - if (last) - { - polynomial = last; - last = 0; - } - - /* Find a new primitive polynomial */ - while(1) - { - if (polynomial == ULONG_MAX) - { - /* Not enough polynomials */ - free(s->direction); - free(s->numerator); - return RK_SOBOL_EINVAL; - } - - polynomial += 2; - - if (ooord) - { - unsigned long copy = polynomial; - /* We are generating polynomials out of order: - check if the reverse was already checked */ - for (rev = 0; copy; copy >>= 1) - rev = (rev << 1) | (copy & 1); - if (ooord && rev < polynomial) - continue; - } - - if (rk_isprimitive(polynomial)) - break; - } - - if (rev == polynomial) - /* We are generating polynomials out of order: - the reverse is not different, discard it */ - rev = 0; - } - - /* Compute the degree */ - for (temp = polynomial >> 1, degree = 0; temp; degree++, temp >>= 1); - - for (j=0; jdirection[j*dimension+k] = m << (LONG_BIT-j-1); - } - - /* Scaled recursion for directions */ - for (j = degree; j < LONG_BIT; j++) - { - unsigned long effdir = s->direction[(j-degree)*dimension+k], - ptemp = polynomial >> 1; - effdir ^= (effdir >> degree); - for (l = degree-1; l >= 1; l--, ptemp >>= 1) - if (ptemp & 1) - effdir ^= s->direction[(j-l)*dimension+k]; - s->direction[j*dimension+k] = effdir; - } - - /* Can we generate polynomials out of order ? */ - if (!ooord && polynomials[cpol] == 0 && degree > last_degree - && (directions == NULL || directions[cdir] == 0)) - ooord = 0; - else - last_degree = degree; - } - - /* Initialize numerator */ - for (k=0; knumerator[k] = 0; - - s->dimension = dimension; - s->gcount = 0; - s->count = 0; - return RK_SOBOL_OK; -} - -void rk_sobol_reinit(rk_sobol_state *s) -{ - size_t k; - - /* Initialize numerator */ - for (k=0; kdimension; k++) - s->numerator[k] = 0; - - s->count = 0; - s->gcount = 0; -} - -void rk_sobol_randomshift(rk_sobol_state *s, rk_state *rs_num) -{ - rk_state rs_num_temp; - size_t k; - - if (rs_num == NULL) - { - rs_num = &rs_num_temp; - rk_randomseed(rs_num); - } - - /* Initialize numerator */ - for (k=0; kdimension; k++) - s->numerator[k] = rk_ulong(rs_num); -} - -rk_sobol_error rk_sobol_copy(rk_sobol_state *copy, rk_sobol_state *orig) -{ - size_t k; - - /* Allocate the structure */ - copy->direction = NULL; copy->numerator = NULL; - copy->direction = malloc(sizeof(*(copy->direction))*orig->dimension*LONG_BIT); - copy->numerator = malloc(sizeof(*(copy->numerator))*orig->dimension); - if (!copy->direction | !copy->numerator) - { - if (!copy->direction) free(copy->direction); - if (!copy->numerator) free(copy->numerator); - return RK_SOBOL_ENOMEM; - } - - /* Initialize numerator */ - for (k=0; kdimension; k++) - copy->numerator[k] = orig->numerator[k]; - for (k=0; k<(orig->dimension*LONG_BIT); k++) - copy->direction[k] = orig->direction[k]; - - copy->count = orig->count; - copy->gcount = orig->gcount; - copy->dimension = orig->dimension; - - return RK_SOBOL_OK; -} - -rk_sobol_error rk_sobol_double(rk_sobol_state *s, double *x) -{ - int j; - size_t k; - unsigned long im; - const double inverse_denominator=1.0/(ULONG_MAX+1.0); - - if (s->count == ULONG_MAX) - j = 0; - else - for (im = s->count, j=0; im & 1; j++, im >>= 1); - s->count++; - - for (k=0; kdimension; k++) - { - s->numerator[k] ^= s->direction[j*s->dimension+k]; - x[k] = s->numerator[k]*inverse_denominator; - } - - if ((s->gcount++) == ULONG_MAX) return RK_SOBOL_EXHAUST; - return RK_SOBOL_OK; -} - -void rk_sobol_setcount(rk_sobol_state *s, unsigned long count) -{ - s->count = count; -} - -void rk_sobol_free(rk_sobol_state *s) -{ - free(s->direction); - free(s->numerator); -} - -double inverse_normal(double p) -{ - double q, t, x; - - /* Peter J. Acklam constants for the rational approximation */ - const double a[6] = - { - -3.969683028665376e+01, 2.209460984245205e+02, - -2.759285104469687e+02, 1.383577518672690e+02, - -3.066479806614716e+01, 2.506628277459239e+00 - }; - const double b[5] = - { - -5.447609879822406e+01, 1.615858368580409e+02, - -1.556989798598866e+02, 6.680131188771972e+01, - -1.328068155288572e+01 - }; - const double c[6] = - { - -7.784894002430293e-03, -3.223964580411365e-01, - -2.400758277161838e+00, -2.549732539343734e+00, - 4.374664141464968e+00, 2.938163982698783e+00 - }; - const double d[4] = - { - 7.784695709041462e-03, 3.224671290700398e-01, - 2.445134137142996e+00, 3.754408661907416e+00 - }; - - if (p <= 0) - return -HUGE_VAL; - else if (p >= 1) - return HUGE_VAL; - - q = p<0.5 ? p : 1-p; - if (q > 0.02425) - { - /* Rational approximation for central region */ - x = q-0.5; - t = x*x; - x = x*(((((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4])*t+a[5]) - /(((((b[0]*t+b[1])*t+b[2])*t+b[3])*t+b[4])*t+1); - } - else - { - /* Rational approximation for tail region */ - t = sqrt(-2*log(q)); - x = (((((c[0]*t+c[1])*t+c[2])*t+c[3])*t+c[4])*t+c[5]) - /((((d[0]*t+d[1])*t+d[2])*t+d[3])*t+1); - } - - /* If we have erfc, improve the precision */ -#ifndef WIN32 - /* Halley's rational method */ - t = (erfc(-x*M_SQRT1_2)/2 - q) * RK_SOBOL_M_SQRT2PI * exp(x*x/2); - x -= t/(1 + x*t/2); -#endif - - return p>0.5 ? -x : x; -} - -rk_sobol_error rk_sobol_gauss(rk_sobol_state *s, double *x) -{ - size_t k; - rk_sobol_error rc = rk_sobol_double(s, x); - - for (k=0; kdimension; k++) - x[k] = inverse_normal(x[k]); - - return rc; -} diff --git a/tardis/setup_package.py b/tardis/setup_package.py deleted file mode 100644 index defe4c24468..00000000000 --- a/tardis/setup_package.py +++ /dev/null @@ -1,13 +0,0 @@ -#setting the right include -from setuptools import Extension -import numpy as np - -randomkit_files = ['tardis/randomkit/rk_isaac.c', 'tardis/randomkit/rk_mt.c', - 'tardis/randomkit/rk_primitive.c', - 'tardis/randomkit/rk_sobol.c'] - -def get_extensions(): - return [Extension('tardis.montecarlo', - ['tardis/montecarlo.pyx', 'tardis/cmontecarlo.c'] + - randomkit_files, - include_dirs=['tardis/randomkit', np.get_include()])] From 599538d611d85b91d0ea30e8902f853dd6a187f3 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Wed, 10 Dec 2014 15:21:23 +0100 Subject: [PATCH 2/9] adding .c and .h files in the restructure montecarlo to setup_package.py --- tardis/montecarlo/setup_package.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py index 93fccc684c6..fff8c4cb8f9 100644 --- a/tardis/montecarlo/setup_package.py +++ b/tardis/montecarlo/setup_package.py @@ -14,3 +14,7 @@ def get_extensions(): include_dirs=['tardis/montecarlo/src', 'tardis/montecarlo/src/randomkit', np.get_include()])] +def get_package_data(): + return {'tardis.montecarlo.montecarlo': ['src/*.c', 'src/*.h', + 'src/randomkit/*.c', + 'src/randomkit/*.h']} \ No newline at end of file From 1b479042adc0a5b521273f0a0591f234fd840003 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Wed, 10 Dec 2014 15:33:34 +0100 Subject: [PATCH 3/9] changed to a relative path for the c files --- tardis/montecarlo/setup_package.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py index fff8c4cb8f9..639c7404882 100644 --- a/tardis/montecarlo/setup_package.py +++ b/tardis/montecarlo/setup_package.py @@ -8,13 +8,9 @@ def get_extensions(): return [Extension('tardis.montecarlo.montecarlo', ['tardis/montecarlo/montecarlo.pyx', - 'tardis/montecarlo/src/cmontecarlo.c'] + - [os.path.join('tardis/montecarlo/src/randomkit', fname) + 'src/cmontecarlo.c'] + + [os.path.join('src/randomkit', fname) for fname in randomkit_files], include_dirs=['tardis/montecarlo/src', 'tardis/montecarlo/src/randomkit', np.get_include()])] -def get_package_data(): - return {'tardis.montecarlo.montecarlo': ['src/*.c', 'src/*.h', - 'src/randomkit/*.c', - 'src/randomkit/*.h']} \ No newline at end of file From 98d4d57326b217a2774d6c71ae3e350e8d1817f9 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Wed, 10 Dec 2014 15:54:01 +0100 Subject: [PATCH 4/9] redid the c discovery --- tardis/montecarlo/setup_package.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/tardis/montecarlo/setup_package.py b/tardis/montecarlo/setup_package.py index 639c7404882..f3b6dab0fc2 100644 --- a/tardis/montecarlo/setup_package.py +++ b/tardis/montecarlo/setup_package.py @@ -3,14 +3,17 @@ import numpy as np import os -randomkit_files = ['rk_isaac.c', 'rk_mt.c', 'rk_primitive.c','rk_sobol.c'] +from glob import glob + def get_extensions(): - return [Extension('tardis.montecarlo.montecarlo', - ['tardis/montecarlo/montecarlo.pyx', - 'src/cmontecarlo.c'] + - [os.path.join('src/randomkit', fname) - for fname in randomkit_files], + sources = ['tardis/montecarlo/montecarlo.pyx'] + sources += [os.path.relpath(fname) for fname in glob( + os.path.join(os.path.dirname(__file__), 'src', '*.c'))] + sources += [os.path.relpath(fname) for fname in glob( + os.path.join(os.path.dirname(__file__), 'src/randomkit', '*.c'))] + + return [Extension('tardis.montecarlo.montecarlo', sources, include_dirs=['tardis/montecarlo/src', 'tardis/montecarlo/src/randomkit', np.get_include()])] From 9e8d82d4383d2748a8912124a8a74a1354831094 Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Thu, 11 Dec 2014 14:25:35 +0100 Subject: [PATCH 5/9] added cmontecarlo.c back --- tardis/montecarlo/montecarlo.pyx | 2 +- tardis/montecarlo/src/cmontecarlo.c | 948 ++++++++++++++++++++++++++++ 2 files changed, 949 insertions(+), 1 deletion(-) create mode 100644 tardis/montecarlo/src/cmontecarlo.c diff --git a/tardis/montecarlo/montecarlo.pyx b/tardis/montecarlo/montecarlo.pyx index aea165e1266..8e8fbfe7362 100644 --- a/tardis/montecarlo/montecarlo.pyx +++ b/tardis/montecarlo/montecarlo.pyx @@ -15,7 +15,7 @@ np.import_array() ctypedef np.int64_t int_type_t -cdef extern from "cmontecarlo.h": +cdef extern from "src/cmontecarlo.h": ctypedef enum rpacket_status_t: TARDIS_PACKET_STATUS_IN_PROCESS = 0 TARDIS_PACKET_STATUS_EMITTED = 1 diff --git a/tardis/montecarlo/src/cmontecarlo.c b/tardis/montecarlo/src/cmontecarlo.c new file mode 100644 index 00000000000..bc2229a50e1 --- /dev/null +++ b/tardis/montecarlo/src/cmontecarlo.c @@ -0,0 +1,948 @@ +#include "cmontecarlo.h" + +rk_state mt_state; + +void +initialize_random_kit (unsigned long seed) +{ + rk_seed (seed, &mt_state); +} + +INLINE tardis_error_t +line_search (double *nu, double nu_insert, int64_t number_of_lines, + int64_t * result) +{ + tardis_error_t ret_val = TARDIS_ERROR_OK; + int64_t imin, imax; + imin = 0; + imax = number_of_lines - 1; + if (nu_insert > nu[imin]) + { + *result = imin; + } + else if (nu_insert < nu[imax]) + { + *result = imax + 1; + } + else + { + ret_val = reverse_binary_search (nu, nu_insert, imin, imax, result); + *result = *result + 1; + } + return ret_val; +} + +inline tardis_error_t +reverse_binary_search (double *x, double x_insert, + int64_t imin, int64_t imax, int64_t * result) +{ + /* + Have in mind that *x points to a reverse sorted array. + That is large values will have small indices and small ones + will have large indices. + */ + tardis_error_t ret_val = TARDIS_ERROR_OK; + if (x_insert > x[imin] || x_insert < x[imax]) + { + ret_val = TARDIS_ERROR_BOUNDS_ERROR; + } + else + { + int imid = (imin + imax) / 2; + while (imax - imin > 2) + { + if (x[imid] < x_insert) + { + imax = imid + 1; + } + else + { + imin = imid; + } + imid = (imin + imax) / 2; + } + if (imax - imin == 2 && x_insert < x[imin + 1]) + { + *result = imin + 1; + } + else + { + *result = imin; + } + } + return ret_val; +} + +inline tardis_error_t +binary_search (double *x, double x_insert, int64_t imin, + int64_t imax, int64_t * result) +{ + /* + Have in mind that *x points to a sorted array. + Like [1,2,3,4,5,...] + */ + int imid; + tardis_error_t ret_val = TARDIS_ERROR_OK; + if (x_insert < x[imin] || x_insert > x[imax]) + { + ret_val = TARDIS_ERROR_BOUNDS_ERROR; + } + else + { + while (imax >= imin) + { + imid = (imin + imax) / 2; + if (x[imid] == x_insert) + { + *result = imid; + break; + } + else if (x[imid] < x_insert) + { + imin = imid + 1; + } + else + { + imax = imid - 1; + } + } + if (imax - imid == 2 && x_insert < x[imin + 1]) + { + *result = imin; + } + else + { + *result = imin; + } + } + return ret_val; +} + +INLINE double +rpacket_doppler_factor (rpacket_t * packet, storage_model_t * storage) +{ + return 1.0 - + rpacket_get_mu (packet) * rpacket_get_r (packet) * + storage->inverse_time_explosion * INVERSE_C; +} + +INLINE double +compute_distance2boundary (rpacket_t * packet, storage_model_t * storage) +{ + double r = rpacket_get_r (packet); + double mu = rpacket_get_mu (packet); + double r_outer = storage->r_outer[rpacket_get_current_shell_id (packet)]; + double r_inner = storage->r_inner[rpacket_get_current_shell_id (packet)]; + double d_outer = + sqrt (r_outer * r_outer + ((mu * mu - 1.0) * r * r)) - (r * mu); + double d_inner; + double result; + if (rpacket_get_recently_crossed_boundary (packet) == 1) + { + rpacket_set_next_shell_id (packet, 1); + return d_outer; + } + else + { + double check = r_inner * r_inner + (r * r * (mu * mu - 1.0)); + if (check < 0.0) + { + rpacket_set_next_shell_id (packet, 1); + return d_outer; + } + else + { + d_inner = mu < 0.0 ? -r * mu - sqrt (check) : MISS_DISTANCE; + } + } + if (d_inner < d_outer) + { + rpacket_set_next_shell_id (packet, -1); + return d_inner; + } + else + { + rpacket_set_next_shell_id (packet, 1); + return d_outer; + } +} + +INLINE tardis_error_t +compute_distance2line (rpacket_t * packet, storage_model_t * storage, + double *result) +{ + tardis_error_t ret_val = TARDIS_ERROR_OK; + if (rpacket_get_last_line (packet)) + { + *result = MISS_DISTANCE; + } + else + { + double r = rpacket_get_r (packet); + double mu = rpacket_get_mu (packet); + double nu = rpacket_get_nu (packet); + double nu_line = rpacket_get_nu_line (packet); + double t_exp = storage->time_explosion; + double inverse_t_exp = storage->inverse_time_explosion; + int64_t cur_zone_id = rpacket_get_current_shell_id (packet); + double comov_nu, doppler_factor; + doppler_factor = 1.0 - mu * r * inverse_t_exp * INVERSE_C; + comov_nu = nu * doppler_factor; + if (comov_nu < nu_line) + { + if (rpacket_get_next_line_id (packet) == storage->no_of_lines - 1) + { + fprintf (stderr, "last_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) - 1]); + fprintf (stderr, "Last line in line list reached!"); + } + else if (rpacket_get_next_line_id (packet) == 0) + { + fprintf (stderr, "First line in line list!"); + fprintf (stderr, "next_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) + 1]); + } + else + { + fprintf (stderr, "last_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) - 1]); + fprintf (stderr, "next_line = %f\n", + storage-> + line_list_nu[rpacket_get_next_line_id (packet) + 1]); + } + fprintf (stderr, "ERROR: Comoving nu less than nu_line!\n"); + fprintf (stderr, "comov_nu = %f\n", comov_nu); + fprintf (stderr, "nu_line = %f\n", nu_line); + fprintf (stderr, "(comov_nu - nu_line) / nu_line = %f\n", + (comov_nu - nu_line) / nu_line); + fprintf (stderr, "r = %f\n", r); + fprintf (stderr, "mu = %f\n", mu); + fprintf (stderr, "nu = %f\n", nu); + fprintf (stderr, "doppler_factor = %f\n", doppler_factor); + fprintf (stderr, "cur_zone_id = %d\n", cur_zone_id); + ret_val = TARDIS_ERROR_COMOV_NU_LESS_THAN_NU_LINE; + } + else + { + *result = ((comov_nu - nu_line) / nu) * C * t_exp; + } + } + return ret_val; +} + +INLINE double +compute_distance2electron (rpacket_t * packet, storage_model_t * storage) +{ + if (rpacket_get_virtual_packet (packet) > 0) + { + return MISS_DISTANCE; + } + double inverse_ne = + storage-> + inverse_electron_densities[rpacket_get_current_shell_id (packet)] * + storage->inverse_sigma_thomson; + return rpacket_get_tau_event (packet) * inverse_ne; +} + +INLINE int64_t +macro_atom (rpacket_t * packet, storage_model_t * storage) +{ + int emit = 0, i = 0; + double p, event_random; + int activate_level = + storage->line2macro_level_upper[rpacket_get_next_line_id (packet) - 1]; + while (emit != -1) + { + event_random = rk_double (&mt_state); + i = storage->macro_block_references[activate_level] - 1; + p = 0.0; + do + { + p += + storage-> + transition_probabilities[rpacket_get_current_shell_id (packet) * + storage->transition_probabilities_nd + + (++i)]; + } + while (p <= event_random); + emit = storage->transition_type[i]; + activate_level = storage->destination_level_id[i]; + } + return storage->transition_line_id[i]; +} + +INLINE double +move_packet (rpacket_t * packet, storage_model_t * storage, double distance) +{ + double new_r, doppler_factor, comov_energy, comov_nu; + doppler_factor = rpacket_doppler_factor (packet, storage); + if (distance > 0.0) + { + double r = rpacket_get_r (packet); + new_r = + sqrt (r * r + distance * distance + + 2.0 * r * distance * rpacket_get_mu (packet)); + rpacket_set_mu (packet, + (rpacket_get_mu (packet) * r + distance) / new_r); + rpacket_set_r (packet, new_r); + if (rpacket_get_virtual_packet (packet) <= 0) + { + comov_energy = rpacket_get_energy (packet) * doppler_factor; + comov_nu = rpacket_get_nu (packet) * doppler_factor; + storage->js[rpacket_get_current_shell_id (packet)] += + comov_energy * distance; + storage->nubars[rpacket_get_current_shell_id (packet)] += + comov_energy * distance * comov_nu; + } + } + return doppler_factor; +} + +INLINE void +increment_j_blue_estimator (rpacket_t * packet, storage_model_t * storage, + double d_line, int64_t j_blue_idx) +{ + double comov_energy, r_interaction, mu_interaction, doppler_factor; + double r = rpacket_get_r (packet); + r_interaction = + sqrt (r * r + d_line * d_line + + 2.0 * r * d_line * rpacket_get_mu (packet)); + mu_interaction = (rpacket_get_mu (packet) * r + d_line) / r_interaction; + doppler_factor = 1.0 - mu_interaction * r_interaction * + storage->inverse_time_explosion * INVERSE_C; + comov_energy = rpacket_get_energy (packet) * doppler_factor; + storage->line_lists_j_blues[j_blue_idx] += + comov_energy / rpacket_get_nu (packet); +} + +int64_t +montecarlo_one_packet (storage_model_t * storage, rpacket_t * packet, + int64_t virtual_mode) +{ + int64_t i; + rpacket_t virt_packet; + double mu_bin; + double mu_min; + double doppler_factor_ratio; + double weight; + int64_t virt_id_nu; + int64_t reabsorbed; + if (virtual_mode == 0) + { + reabsorbed = montecarlo_one_packet_loop (storage, packet, 0); + } + else + { + for (i = 0; i < rpacket_get_virtual_packet_flag (packet); i++) + { + memcpy ((void *) &virt_packet, (void *) packet, sizeof (rpacket_t)); + if (virt_packet.r > storage->r_inner[0]) + { + mu_min = + -1.0 * sqrt (1.0 - + (storage->r_inner[0] / virt_packet.r) * + (storage->r_inner[0] / virt_packet.r)); + } + else + { + mu_min = 0.0; + } + mu_bin = (1.0 - mu_min) / rpacket_get_virtual_packet_flag (packet); + virt_packet.mu = mu_min + (i + rk_double (&mt_state)) * mu_bin; + switch (virtual_mode) + { + case -2: + weight = 1.0 / rpacket_get_virtual_packet_flag (packet); + break; + case -1: + weight = + 2.0 * virt_packet.mu / + rpacket_get_virtual_packet_flag (packet); + break; + case 1: + weight = + (1.0 - + mu_min) / 2.0 / rpacket_get_virtual_packet_flag (packet); + break; + default: + fprintf (stderr, "Something has gone horribly wrong!\n"); + } + doppler_factor_ratio = + rpacket_doppler_factor (packet, storage) / + rpacket_doppler_factor (&virt_packet, storage); + virt_packet.energy = + rpacket_get_energy (packet) * doppler_factor_ratio; + virt_packet.nu = rpacket_get_nu (packet) * doppler_factor_ratio; + reabsorbed = montecarlo_one_packet_loop (storage, &virt_packet, 1); + if ((virt_packet.nu < storage->spectrum_end_nu) && + (virt_packet.nu > storage->spectrum_start_nu)) + { + virt_id_nu = + floor ((virt_packet.nu - + storage->spectrum_start_nu) / + storage->spectrum_delta_nu); + storage->spectrum_virt_nu[virt_id_nu] += + virt_packet.energy * weight; + } + } + } + return reabsorbed; +} + +void +move_packet_across_shell_boundary (rpacket_t * packet, + storage_model_t * storage, double distance) +{ + double comov_energy, doppler_factor, comov_nu, inverse_doppler_factor; + move_packet (packet, storage, distance); + if (rpacket_get_virtual_packet (packet) > 0) + { + double delta_tau_event = distance * + storage->electron_densities[rpacket_get_current_shell_id (packet)] * + storage->sigma_thomson; + rpacket_set_tau_event (packet, + rpacket_get_tau_event (packet) + + delta_tau_event); + } + else + { + rpacket_reset_tau_event (packet); + } + if ((rpacket_get_current_shell_id (packet) < storage->no_of_shells - 1 + && rpacket_get_next_shell_id (packet) == 1) + || (rpacket_get_current_shell_id (packet) > 0 + && rpacket_get_next_shell_id (packet) == -1)) + { + rpacket_set_current_shell_id (packet, + rpacket_get_current_shell_id (packet) + + rpacket_get_next_shell_id (packet)); + rpacket_set_recently_crossed_boundary (packet, + rpacket_get_next_shell_id + (packet)); + } + else if (rpacket_get_next_shell_id (packet) == 1) + { + rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); + } + else if ((storage->reflective_inner_boundary == 0) || + (rk_double (&mt_state) > storage->inner_boundary_albedo)) + { + rpacket_set_status (packet, TARDIS_PACKET_STATUS_REABSORBED); + } + else + { + doppler_factor = rpacket_doppler_factor (packet, storage); + comov_nu = rpacket_get_nu (packet) * doppler_factor; + comov_energy = rpacket_get_energy (packet) * doppler_factor; + rpacket_set_mu (packet, rk_double (&mt_state)); + inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage); + rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + rpacket_set_recently_crossed_boundary (packet, 1); + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + montecarlo_one_packet (storage, packet, -2); + } + } +} + +void +montecarlo_thomson_scatter (rpacket_t * packet, storage_model_t * storage, + double distance) +{ + double comov_energy, doppler_factor, comov_nu, inverse_doppler_factor; + doppler_factor = move_packet (packet, storage, distance); + comov_nu = rpacket_get_nu (packet) * doppler_factor; + comov_energy = rpacket_get_energy (packet) * doppler_factor; + rpacket_set_mu (packet, 2.0 * rk_double (&mt_state) - 1.0); + inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage); + rpacket_set_nu (packet, comov_nu * inverse_doppler_factor); + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + rpacket_reset_tau_event (packet); + rpacket_set_recently_crossed_boundary (packet, 0); + storage->last_interaction_type[storage->current_packet_id] = 1; + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + montecarlo_one_packet (storage, packet, 1); + } +} + +void +montecarlo_line_scatter (rpacket_t * packet, storage_model_t * storage, + double distance) +{ + double comov_energy = 0.0; + int64_t emission_line_id = 0; + double old_doppler_factor = 0.0; + double inverse_doppler_factor = 0.0; + double tau_line = 0.0; + double tau_electron = 0.0; + double tau_combined = 0.0; + bool virtual_close_line = false; + int64_t j_blue_idx = -1; + if (rpacket_get_virtual_packet (packet) == 0) + { + j_blue_idx = + rpacket_get_current_shell_id (packet) * + storage->line_lists_j_blues_nd + rpacket_get_next_line_id (packet); + increment_j_blue_estimator (packet, storage, distance, j_blue_idx); + } + tau_line = + storage->line_lists_tau_sobolevs[rpacket_get_current_shell_id (packet) * + storage->line_lists_tau_sobolevs_nd + + rpacket_get_next_line_id (packet)]; + tau_electron = + storage->sigma_thomson * + storage->electron_densities[rpacket_get_current_shell_id (packet)] * + distance; + tau_combined = tau_line + tau_electron; + rpacket_set_next_line_id (packet, rpacket_get_next_line_id (packet) + 1); + if (rpacket_get_next_line_id (packet) == storage->no_of_lines) + { + rpacket_set_last_line (packet, true); + } + if (rpacket_get_virtual_packet (packet) > 0) + { + rpacket_set_tau_event (packet, + rpacket_get_tau_event (packet) + tau_line); + } + else if (rpacket_get_tau_event (packet) < tau_combined) + { + old_doppler_factor = move_packet (packet, storage, distance); + rpacket_set_mu (packet, 2.0 * rk_double (&mt_state) - 1.0); + inverse_doppler_factor = 1.0 / rpacket_doppler_factor (packet, storage); + comov_energy = rpacket_get_energy (packet) * old_doppler_factor; + rpacket_set_energy (packet, comov_energy * inverse_doppler_factor); + storage->last_line_interaction_in_id[storage->current_packet_id] = + rpacket_get_next_line_id (packet) - 1; + storage->last_line_interaction_shell_id[storage->current_packet_id] = + rpacket_get_current_shell_id (packet); + storage->last_interaction_type[storage->current_packet_id] = 2; + if (storage->line_interaction_id == 0) + { + emission_line_id = rpacket_get_next_line_id (packet) - 1; + } + else if (storage->line_interaction_id >= 1) + { + emission_line_id = macro_atom (packet, storage); + } + storage->last_line_interaction_out_id[storage->current_packet_id] = + emission_line_id; + rpacket_set_nu (packet, + storage->line_list_nu[emission_line_id] * + inverse_doppler_factor); + rpacket_set_nu_line (packet, storage->line_list_nu[emission_line_id]); + rpacket_set_next_line_id (packet, emission_line_id + 1); + rpacket_reset_tau_event (packet); + rpacket_set_recently_crossed_boundary (packet, 0); + if (rpacket_get_virtual_packet_flag (packet) > 0) + { + virtual_close_line = false; + if (!rpacket_get_last_line (packet) && + fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - + rpacket_get_nu_line (packet)) / + rpacket_get_nu_line (packet) < 1e-7) + { + virtual_close_line = true; + } + // QUESTIONABLE!!! + bool old_close_line = rpacket_get_close_line (packet); + rpacket_set_close_line (packet, virtual_close_line); + montecarlo_one_packet (storage, packet, 1); + rpacket_set_close_line (packet, old_close_line); + virtual_close_line = false; + } + } + else + { + rpacket_set_tau_event (packet, + rpacket_get_tau_event (packet) - tau_line); + } + if (!rpacket_get_last_line (packet) && + fabs (storage->line_list_nu[rpacket_get_next_line_id (packet)] - + rpacket_get_nu_line (packet)) / rpacket_get_nu_line (packet) < + 1e-7) + { + rpacket_set_close_line (packet, true); + } +} + +INLINE void +montecarlo_compute_distances (rpacket_t * packet, storage_model_t * storage) +{ + // Check if the last line was the same nu as the current line. + if (rpacket_get_close_line (packet)) + { + // If so set the distance to the line to 0.0 + rpacket_set_d_line (packet, 0.0); + // Reset close_line. + rpacket_set_close_line (packet, false); + } + else + { + rpacket_set_d_boundary (packet, + compute_distance2boundary (packet, storage)); + double d_line; + compute_distance2line (packet, storage, &d_line); + rpacket_set_d_line (packet, d_line); + rpacket_set_d_electron (packet, + compute_distance2electron (packet, storage)); + } +} + +INLINE montecarlo_event_handler_t +get_event_handler (rpacket_t * packet, storage_model_t * storage, + double *distance) +{ + double d_boundary, d_electron, d_line; + montecarlo_compute_distances (packet, storage); + d_boundary = rpacket_get_d_boundary (packet); + d_electron = rpacket_get_d_electron (packet); + d_line = rpacket_get_d_line (packet); + montecarlo_event_handler_t handler; + if (d_line <= d_boundary && d_line <= d_electron) + { + *distance = d_line; + handler = &montecarlo_line_scatter; + } + else if (d_boundary <= d_electron) + { + *distance = d_boundary; + handler = &move_packet_across_shell_boundary; + } + else + { + *distance = d_electron; + handler = &montecarlo_thomson_scatter; + } + return handler; +} + +int64_t +montecarlo_one_packet_loop (storage_model_t * storage, rpacket_t * packet, + int64_t virtual_packet) +{ + rpacket_set_tau_event (packet, 0.0); + rpacket_set_nu_line (packet, 0.0); + rpacket_set_virtual_packet (packet, virtual_packet); + rpacket_set_status (packet, TARDIS_PACKET_STATUS_IN_PROCESS); + // Initializing tau_event if it's a real packet. + if (virtual_packet == 0) + { + rpacket_reset_tau_event (packet); + } + // For a virtual packet tau_event is the sum of all the tau's that the packet passes. + while (rpacket_get_status (packet) == TARDIS_PACKET_STATUS_IN_PROCESS) + { + // Check if we are at the end of line list. + if (!rpacket_get_last_line (packet)) + { + rpacket_set_nu_line (packet, + storage-> + line_list_nu[rpacket_get_next_line_id + (packet)]); + } + double distance; + get_event_handler (packet, storage, &distance) (packet, storage, + distance); + if (virtual_packet > 0 && rpacket_get_tau_event (packet) > 10.0) + { + rpacket_set_tau_event (packet, 100.0); + rpacket_set_status (packet, TARDIS_PACKET_STATUS_EMITTED); + } + } + if (virtual_packet > 0) + { + rpacket_set_energy (packet, + rpacket_get_energy (packet) * exp (-1.0 * + rpacket_get_tau_event + (packet))); + } + return rpacket_get_status (packet) == + TARDIS_PACKET_STATUS_REABSORBED ? 1 : 0; +} + +tardis_error_t +rpacket_init (rpacket_t * packet, storage_model_t * storage, int packet_index, + int virtual_packet_flag) +{ + double nu_line; + double current_r; + double current_mu; + double current_nu; + double comov_current_nu; + double current_energy; + int64_t current_line_id; + int current_shell_id; + bool last_line; + bool close_line; + int recently_crossed_boundary; + tardis_error_t ret_val = TARDIS_ERROR_OK; + current_nu = storage->packet_nus[packet_index]; + current_energy = storage->packet_energies[packet_index]; + current_mu = storage->packet_mus[packet_index]; + comov_current_nu = current_nu; + current_shell_id = 0; + current_r = storage->r_inner[0]; + current_nu = + current_nu / (1 - + (current_mu * current_r * storage->inverse_time_explosion * + INVERSE_C)); + current_energy = + current_energy / (1 - + (current_mu * current_r * + storage->inverse_time_explosion * INVERSE_C)); + if ((ret_val = + line_search (storage->line_list_nu, comov_current_nu, + storage->no_of_lines, + ¤t_line_id)) != TARDIS_ERROR_OK) + { + return ret_val; + } + last_line = (current_line_id == storage->no_of_lines); + recently_crossed_boundary = true; + rpacket_set_nu (packet, current_nu); + rpacket_set_mu (packet, current_mu); + rpacket_set_energy (packet, current_energy); + rpacket_set_r (packet, current_r); + rpacket_set_current_shell_id (packet, current_shell_id); + rpacket_set_next_line_id (packet, current_line_id); + rpacket_set_last_line (packet, last_line); + rpacket_set_close_line (packet, false); + rpacket_set_recently_crossed_boundary (packet, recently_crossed_boundary); + rpacket_set_virtual_packet_flag (packet, virtual_packet_flag); + return ret_val; +} + +/* + Getter and setter methods. +*/ + +INLINE double +rpacket_get_nu (rpacket_t * packet) +{ + return packet->nu; +} + +INLINE void +rpacket_set_nu (rpacket_t * packet, double nu) +{ + packet->nu = nu; +} + +INLINE double +rpacket_get_mu (rpacket_t * packet) +{ + return packet->mu; +} + +INLINE void +rpacket_set_mu (rpacket_t * packet, double mu) +{ + packet->mu = mu; +} + +INLINE double +rpacket_get_energy (rpacket_t * packet) +{ + return packet->energy; +} + +INLINE void +rpacket_set_energy (rpacket_t * packet, double energy) +{ + packet->energy = energy; +} + +INLINE double +rpacket_get_r (rpacket_t * packet) +{ + return packet->r; +} + +INLINE void +rpacket_set_r (rpacket_t * packet, double r) +{ + packet->r = r; +} + +INLINE double +rpacket_get_tau_event (rpacket_t * packet) +{ + return packet->tau_event; +} + +INLINE void +rpacket_set_tau_event (rpacket_t * packet, double tau_event) +{ + packet->tau_event = tau_event; +} + +INLINE double +rpacket_get_nu_line (rpacket_t * packet) +{ + return packet->nu_line; +} + +INLINE void +rpacket_set_nu_line (rpacket_t * packet, double nu_line) +{ + packet->nu_line = nu_line; +} + +INLINE unsigned int +rpacket_get_current_shell_id (rpacket_t * packet) +{ + return packet->current_shell_id; +} + +INLINE void +rpacket_set_current_shell_id (rpacket_t * packet, + unsigned int current_shell_id) +{ + packet->current_shell_id = current_shell_id; +} + +INLINE unsigned int +rpacket_get_next_line_id (rpacket_t * packet) +{ + return packet->next_line_id; +} + +INLINE void +rpacket_set_next_line_id (rpacket_t * packet, unsigned int next_line_id) +{ + packet->next_line_id = next_line_id; +} + +INLINE bool +rpacket_get_last_line (rpacket_t * packet) +{ + return packet->last_line; +} + +INLINE void +rpacket_set_last_line (rpacket_t * packet, bool last_line) +{ + packet->last_line = last_line; +} + +INLINE bool +rpacket_get_close_line (rpacket_t * packet) +{ + return packet->close_line; +} + +INLINE void +rpacket_set_close_line (rpacket_t * packet, bool close_line) +{ + packet->close_line = close_line; +} + +INLINE int +rpacket_get_recently_crossed_boundary (rpacket_t * packet) +{ + return packet->recently_crossed_boundary; +} + +INLINE void +rpacket_set_recently_crossed_boundary (rpacket_t * packet, + int recently_crossed_boundary) +{ + packet->recently_crossed_boundary = recently_crossed_boundary; +} + +INLINE int +rpacket_get_virtual_packet_flag (rpacket_t * packet) +{ + return packet->virtual_packet_flag; +} + +INLINE void +rpacket_set_virtual_packet_flag (rpacket_t * packet, int virtual_packet_flag) +{ + packet->virtual_packet_flag = virtual_packet_flag; +} + +INLINE int +rpacket_get_virtual_packet (rpacket_t * packet) +{ + return packet->virtual_packet; +} + +INLINE void +rpacket_set_virtual_packet (rpacket_t * packet, int virtual_packet) +{ + packet->virtual_packet = virtual_packet; +} + +INLINE double +rpacket_get_d_boundary (rpacket_t * packet) +{ + return packet->d_boundary; +} + +INLINE void +rpacket_set_d_boundary (rpacket_t * packet, double d_boundary) +{ + packet->d_boundary = d_boundary; +} + +INLINE double +rpacket_get_d_electron (rpacket_t * packet) +{ + return packet->d_electron; +} + +INLINE void +rpacket_set_d_electron (rpacket_t * packet, double d_electron) +{ + packet->d_electron = d_electron; +} + +INLINE double +rpacket_get_d_line (rpacket_t * packet) +{ + return packet->d_line; +} + +INLINE void +rpacket_set_d_line (rpacket_t * packet, double d_line) +{ + packet->d_line = d_line; +} + +INLINE int +rpacket_get_next_shell_id (rpacket_t * packet) +{ + return packet->next_shell_id; +} + +INLINE void +rpacket_set_next_shell_id (rpacket_t * packet, int next_shell_id) +{ + packet->next_shell_id = next_shell_id; +} + +INLINE rpacket_status_t +rpacket_get_status (rpacket_t * packet) +{ + return packet->status; +} + +INLINE void +rpacket_set_status (rpacket_t * packet, rpacket_status_t status) +{ + packet->status = status; +} + +/* Other accessor methods. */ + +INLINE void +rpacket_reset_tau_event (rpacket_t * packet) +{ + rpacket_set_tau_event (packet, -log (rk_double (&mt_state))); +} From 378f258ee33857b19ec51692d08578c0516283cd Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Thu, 11 Dec 2014 15:32:26 +0100 Subject: [PATCH 6/9] added right path for randomkit.h --- tardis/montecarlo/src/cmontecarlo.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tardis/montecarlo/src/cmontecarlo.h b/tardis/montecarlo/src/cmontecarlo.h index 3d67e92bde5..bbcf2707d9a 100644 --- a/tardis/montecarlo/src/cmontecarlo.h +++ b/tardis/montecarlo/src/cmontecarlo.h @@ -7,7 +7,7 @@ #include #include #include -#include "randomkit.h" +#include "randomkit/randomkit.h" #ifdef __clang__ #define INLINE extern inline From 83ac255ebb72059e7257b5642caf7b5fe17b4fed Mon Sep 17 00:00:00 2001 From: Wolfgang Kerzendorf Date: Thu, 11 Dec 2014 16:05:31 +0100 Subject: [PATCH 7/9] readded the c files for the randomkit --- tardis/montecarlo/src/randomkit/rk_isaac.c | 258 +++++ tardis/montecarlo/src/randomkit/rk_mt.c | 342 ++++++ .../montecarlo/src/randomkit/rk_primitive.c | 520 +++++++++ tardis/montecarlo/src/randomkit/rk_sobol.c | 991 ++++++++++++++++++ 4 files changed, 2111 insertions(+) create mode 100644 tardis/montecarlo/src/randomkit/rk_isaac.c create mode 100644 tardis/montecarlo/src/randomkit/rk_mt.c create mode 100644 tardis/montecarlo/src/randomkit/rk_primitive.c create mode 100644 tardis/montecarlo/src/randomkit/rk_sobol.c diff --git a/tardis/montecarlo/src/randomkit/rk_isaac.c b/tardis/montecarlo/src/randomkit/rk_isaac.c new file mode 100644 index 00000000000..9f52fafae4c --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_isaac.c @@ -0,0 +1,258 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * ISAAC RNG by By Bob Jenkins. Based on Bob Jenkins public domain code. + * + * Original algorithm for the implementation of rk_interval function from + * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by + * Magnus Jonsson. + * + * Constants used in the rk_double implementation by Isaku Wada. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_isaac.c,v 1.4 2006/02/20 19:13:37 js Exp $"; + +#include "rk_isaac.h" +#include +#include +#include + +/* use the contents of randrsl[0..RK_ISAAC_STATE_LEN-1] as the seed. */ +static void isaac_init(rk_isaac_state *rg); + +void rk_isaac_seed(unsigned long seed, rk_isaac_state *state) +{ + rk_knuth_fill(seed, state->randrsl, RK_ISAAC_STATE_LEN); + isaac_init(state); +} + +rk_error rk_isaac_randomseed(rk_isaac_state *state) +{ + if(rk_devfill(state->randrsl, sizeof(state->randrsl), 1) == RK_NOERR) + { + isaac_init(state); + return RK_NOERR; + } + + rk_isaac_seed(rk_seedfromsystem(), state); + + return RK_ENODEV; +} + +#define ind(mm,x) ((mm)[(x>>2)&(RK_ISAAC_STATE_LEN-1)]) +#define rngstep(mix,a,b,mm,m,m2,r,x) \ +{ \ + x = *m; \ + a = ((a^(mix)) + *(m2++)) & 0xffffffff; \ + *(m++) = y = (ind(mm,x) + a + b) & 0xffffffff; \ + *(r++) = b = (ind(mm,y>>RK_ISAAC_STATE_POW) + x) & 0xffffffff; \ +} + +/* Call rand(rk_isaac_state *r) to retrieve a single 32-bit random value */ +unsigned long rk_isaac_random(rk_isaac_state *state) +{ + if (!state->randcnt--) + { + register unsigned long a,b,x,y,*m,*mm,*m2,*r,*mend; + mm=state->randmem; r=state->randrsl; + a = state->randa; b = (state->randb + (++state->randc)) & 0xffffffff; + for (m = mm, mend = m2 = m+(RK_ISAAC_STATE_LEN/2); m>6 , a, b, mm, m, m2, r, x); + rngstep( a<<2 , a, b, mm, m, m2, r, x); + rngstep( a>>16, a, b, mm, m, m2, r, x); + } + for (m2 = mm; m2>6 , a, b, mm, m, m2, r, x); + rngstep( a<<2 , a, b, mm, m, m2, r, x); + rngstep( a>>16, a, b, mm, m, m2, r, x); + } + state->randb = b; state->randa = a; + state->randcnt=RK_ISAAC_STATE_LEN-1; + } + return state->randrsl[state->randcnt] & 0xFFFFFFFF; +} + +#define mix(a,b,c,d,e,f,g,h) \ +{ \ + a^=b<<11; d+=a; b+=c; \ + b^=(c & 0xFFFFFFFF)>>2; e+=b; c+=d; \ + c^=d<<8; f+=c; d+=e; \ + d^=(e & 0xFFFFFFFF)>>16; g+=d; e+=f; \ + e^=f<<10; h+=e; f+=g; \ + f^=(g & 0xFFFFFFFF)>>4; a+=f; g+=h; \ + g^=h<<8; b+=g; h+=a; \ + h^=(a & 0xFFFFFFFF)>>9; c+=h; a+=b; \ +} + +/* if (flag==1), then use the contents of randrsl[] to initialize mm[]. */ +void isaac_init(rk_isaac_state *rk_isaac_state) +{ + int i; + unsigned long a,b,c,d,e,f,g,h; + unsigned long *m,*r; + rk_isaac_state->randa = rk_isaac_state->randb = rk_isaac_state->randc = 0; + m=rk_isaac_state->randmem; + r=rk_isaac_state->randrsl; + a=b=c=d=e=f=g=h=0x9e3779b9; /* the golden ratio */ + + for (i=0; i<4; ++i) /* scramble it */ + mix(a,b,c,d,e,f,g,h); + + /* initialize using the contents of r[] as the seed */ + for (i=0; irandcnt=0; + rk_isaac_state->has_gauss=0; +} + +long rk_isaac_long(rk_isaac_state *state) +{ + return rk_isaac_ulong(state) >> 1; +} + +unsigned long rk_isaac_ulong(rk_isaac_state *state) +{ +#if ULONG_MAX <= 0xffffffffUL + return rk_isaac_random(state); +#else + /* Assumes 64 bits */ + return (rk_isaac_random(state) << 32) | (rk_isaac_random(state)); +#endif +} + +unsigned long rk_isaac_interval(unsigned long max, rk_isaac_state *state) +{ + unsigned long mask = max, value; + + if (max == 0) return 0; + + /* Smallest bit mask >= max */ + mask |= mask >> 1; + mask |= mask >> 2; + mask |= mask >> 4; + mask |= mask >> 8; + mask |= mask >> 16; +#if ULONG_MAX > 0xffffffffUL + mask |= mask >> 32; +#endif + + /* Search a random value in [0..mask] <= max */ + while ((value = (rk_isaac_ulong(state) & mask)) > max); + + return value; +} + +double rk_isaac_double(rk_isaac_state *state) +{ + /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ + long a = rk_isaac_random(state) >> 5, b = rk_isaac_random(state) >> 6; + return (a * 67108864.0 + b) / 9007199254740992.0; +} + +void rk_isaac_copy(rk_isaac_state *copy, rk_isaac_state *orig) +{ + memcpy(copy, orig, sizeof(rk_isaac_state)); +} + +double rk_isaac_gauss(rk_isaac_state *state) +{ + if (state->has_gauss) + { + state->has_gauss = 0; + return state->gauss; + } + else + { + double f, x1, x2, r2; + do + { + x1 = 2.0*rk_isaac_double(state) - 1.0; + x2 = 2.0*rk_isaac_double(state) - 1.0; + r2 = x1*x1 + x2*x2; + } + while (r2 >= 1.0 || r2 == 0.0); + + f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ + state->has_gauss = 1; + state->gauss = f*x1; /* Keep for next call */ + return f*x2; + } +} + +void rk_isaac_fill(void *buffer, size_t size, rk_isaac_state *state) +{ + unsigned long r; + unsigned char *buf = buffer; + rk_isaac_state tempstate; + + if (size > 0 && state == NULL) + { + rk_isaac_randomseed(&tempstate); + state = &tempstate; + } + + for (; size >= 4; size -= 4) + { + r = rk_isaac_random(state); + *(buf++) = r & 0xFF; + *(buf++) = (r >> 8) & 0xFF; + *(buf++) = (r >> 16) & 0xFF; + *(buf++) = (r >> 24) & 0xFF; + } + + if (!size) return; + + r = rk_isaac_random(state); + + for (; size; r >>= 8, size --) + *(buf++) = (unsigned char)(r & 0xFF); +} + +void rk_seed_isaac(rk_isaac_state *i_state, rk_state *state) +{ + rk_isaac_fill(state->key, sizeof(state->key), i_state); + state->pos = RK_STATE_LEN; + state->has_gauss = 0; +} diff --git a/tardis/montecarlo/src/randomkit/rk_mt.c b/tardis/montecarlo/src/randomkit/rk_mt.c new file mode 100644 index 00000000000..6ef13513989 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_mt.c @@ -0,0 +1,342 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2003-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * The rk_random and rk_seed functions algorithms and the original design of + * the Mersenne Twister RNG: + * + * Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + * All rights reserved. + * + * Original algorithm for the implementation of rk_interval function from + * Richard J. Wagner's implementation of the Mersenne Twister RNG, optimised by + * Magnus Jonsson. + * + * Constants used in the rk_double implementation by Isaku Wada. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_mt.c,v 1.6 2006/02/19 13:48:34 js Exp $"; + +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef _WIN32 +/* Windows */ +#include +#ifndef RK_NO_WINCRYPT +/* Windows crypto */ +#ifndef _WIN32_WINNT +#define _WIN32_WINNT 0x0400 +#endif +#include +#include +#endif +#else +/* Unix */ +#include +#include +#endif + +#include "randomkit.h" + +#ifndef RK_DEV_URANDOM +#define RK_DEV_URANDOM "/dev/urandom" +#endif + +#ifndef RK_DEV_RANDOM +#define RK_DEV_RANDOM "/dev/random" +#endif + +char *rk_strerror[RK_ERR_MAX] = +{ + "no error", + "random device unvavailable" +}; + +/* static functions */ +static unsigned long rk_hash(unsigned long key); + +void rk_knuth_fill(unsigned long seed, unsigned long *key, unsigned long len) +{ + int pos; + seed &= 0xffffffffUL; + + /* Knuth's PRNG as used in the Mersenne Twister reference implementation */ + for (pos=0; pos> 30)) + pos + 1) & 0xffffffffUL; + } +} + +void rk_seed(unsigned long seed, rk_state *state) +{ + rk_knuth_fill(seed, state->key, RK_STATE_LEN); + state->pos = RK_STATE_LEN; + state->has_gauss = 0; +} + +/* Thomas Wang 32 bits integer hash function */ +unsigned long rk_hash(unsigned long key) +{ + key += ~(key << 15); + key ^= (key >> 10); + key += (key << 3); + key ^= (key >> 6); + key += ~(key << 11); + key ^= (key >> 16); + return key; +} + +unsigned long rk_seedfromsystem() +{ +#ifndef _WIN32 + struct timeval tv; +#else + struct _timeb tv; +#endif + +#ifndef _WIN32 + gettimeofday(&tv, NULL); + return rk_hash(getpid()) ^ rk_hash(tv.tv_sec) ^ rk_hash(tv.tv_usec) + ^ rk_hash(clock()); +#else + _ftime(&tv); + return rk_hash(tv.time) ^ rk_hash(tv.millitm) ^ rk_hash(clock()); +#endif +} + +rk_error rk_randomseed(rk_state *state) +{ + if(rk_devfill(state->key, sizeof(state->key), 0) == RK_NOERR) + { + state->key[0] |= 0x80000000UL; /* ensures non-zero key */ + state->pos = RK_STATE_LEN; + state->has_gauss = 0; + return RK_NOERR; + } + + rk_seed(rk_seedfromsystem(), state); + + return RK_ENODEV; +} + +/* Magic Mersenne Twister constants */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL +#define UPPER_MASK 0x80000000UL +#define LOWER_MASK 0x7fffffffUL + +/* Slightly optimised reference implementation of the Mersenne Twister */ +unsigned long rk_random(rk_state *state) +{ + unsigned long y; + + if (state->pos == RK_STATE_LEN) + { + int i; + + for (i=0;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); + state->key[i] = state->key[i+M] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + } + for (;ikey[i] & UPPER_MASK) | (state->key[i+1] & LOWER_MASK); + state->key[i] = state->key[i+(M-N)] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + } + y = (state->key[N-1] & UPPER_MASK) | (state->key[0] & LOWER_MASK); + state->key[N-1] = state->key[M-1] ^ (y>>1) ^ (-(y & 1) & MATRIX_A); + + state->pos = 0; + } + + y = state->key[state->pos++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +long rk_long(rk_state *state) +{ + return rk_ulong(state) >> 1; +} + +unsigned long rk_ulong(rk_state *state) +{ +#if ULONG_MAX <= 0xffffffffUL + return rk_random(state); +#else + /* Assumes 64 bits */ + return (rk_random(state) << 32) | (rk_random(state)); +#endif +} + +unsigned long rk_interval(unsigned long max, rk_state *state) +{ + unsigned long mask = max, value; + + if (max == 0) return 0; + + /* Smallest bit mask >= max */ + mask |= mask >> 1; + mask |= mask >> 2; + mask |= mask >> 4; + mask |= mask >> 8; + mask |= mask >> 16; +#if ULONG_MAX > 0xffffffffUL + mask |= mask >> 32; +#endif + + /* Search a random value in [0..mask] <= max */ + while ((value = (rk_ulong(state) & mask)) > max); + + return value; +} + +double rk_double(rk_state *state) +{ + /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ + long a = rk_random(state) >> 5, b = rk_random(state) >> 6; + return (a * 67108864.0 + b) / 9007199254740992.0; +} + +void rk_copy(rk_state *copy, rk_state *orig) +{ + memcpy(copy, orig, sizeof(rk_state)); +} + +void rk_fill(void *buffer, size_t size, rk_state *state) +{ + unsigned long r; + unsigned char *buf = buffer; + rk_state tempstate; + + if (size > 0 && state == NULL) + { + rk_randomseed(&tempstate); + state = &tempstate; + } + + for (; size >= 4; size -= 4) + { + r = rk_random(state); + *(buf++) = r & 0xFF; + *(buf++) = (r >> 8) & 0xFF; + *(buf++) = (r >> 16) & 0xFF; + *(buf++) = (r >> 24) & 0xFF; + } + + if (!size) return; + + r = rk_random(state); + + for (; size; r >>= 8, size --) + *(buf++) = (unsigned char)(r & 0xFF); +} + +rk_error rk_devfill(void *buffer, size_t size, int strong) +{ +#ifndef _WIN32 + FILE *rfile; + int done; + + if (strong) + rfile = fopen(RK_DEV_RANDOM, "rb"); + else + rfile = fopen(RK_DEV_URANDOM, "rb"); + if (rfile == NULL) + return RK_ENODEV; + done = fread(buffer, size, 1, rfile); + fclose(rfile); + if (done) + return RK_NOERR; +#else + +#ifndef RK_NO_WINCRYPT + HCRYPTPROV hCryptProv; + BOOL done; + + if (!CryptAcquireContext(&hCryptProv, NULL, NULL, PROV_RSA_FULL, + CRYPT_VERIFYCONTEXT) || !hCryptProv) + return RK_ENODEV; + done = CryptGenRandom(hCryptProv, size, (unsigned char *)buffer); + CryptReleaseContext(hCryptProv, 0); + if (done) + return RK_NOERR; +#endif + +#endif + + return RK_ENODEV; +} + +rk_error rk_altfill(void *buffer, size_t size, int strong, rk_state *state) +{ + rk_error err; + + err = rk_devfill(buffer, size, strong); + if (err) + rk_fill(buffer, size, state); + + return err; +} + +double rk_gauss(rk_state *state) +{ + if (state->has_gauss) + { + state->has_gauss = 0; + return state->gauss; + } + else + { + double f, x1, x2, r2; + do + { + x1 = 2.0*rk_double(state) - 1.0; + x2 = 2.0*rk_double(state) - 1.0; + r2 = x1*x1 + x2*x2; + } + while (r2 >= 1.0 || r2 == 0.0); + + f = sqrt(-2.0*log(r2)/r2); /* Box-Muller transform */ + state->has_gauss = 1; + state->gauss = f*x1; /* Keep for next call */ + return f*x2; + } +} diff --git a/tardis/montecarlo/src/randomkit/rk_primitive.c b/tardis/montecarlo/src/randomkit/rk_primitive.c new file mode 100644 index 00000000000..25e70c6c7f5 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_primitive.c @@ -0,0 +1,520 @@ +/* Random kit 1.6 */ +/* Primitivity test for binary polynomials of low degree */ + +/* + * Copyright (c) 2005-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Methodology inspired by scott duplichan's ppsearch code. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_primitive.c,v 1.7 2006/02/19 13:48:34 js Exp $"; + +#include +#include +#include "rk_primitive.h" + +#ifndef LONG_BIT +#if ULONG_MAX <= 0xffffffffUL +#define LONG_BIT 32 +#else +#define LONG_BIT 64 +#endif +#endif + +static unsigned long modmul(unsigned long poly1, unsigned long poly2, + unsigned long modulo, unsigned long mask); +static unsigned long modpow(unsigned long polynomial, unsigned long power, + unsigned long modulo, int degree); + +/* + * For all powers i of two up to 64, list all the number of the form + * (2^i-1)/p for p a prime factor of 2^i-1. + */ +static const unsigned long divisors[][12]={ +/* 2^0-1 */ + {1UL, + 0UL}, +/* 2^1-1 */ + {1UL, + 0UL}, +/* 2^2-1 */ + {1UL, + 0UL}, +/* 2^3-1 */ + {1UL, + 0UL}, +/* 2^4-1 */ + {5UL, + 3UL, + 0UL}, +/* 2^5-1 */ + {1UL, + 0UL}, +/* 2^6-1 */ + {21UL, + 9UL, + 0UL}, +/* 2^7-1 */ + {1UL, + 0UL}, +/* 2^8-1 */ + {85UL, + 51UL, + 15UL, + 0UL}, +/* 2^9-1 */ + {73UL, + 7UL, + 0UL}, +/* 2^10-1 */ + {341UL, + 93UL, + 33UL, + 0UL}, +/* 2^11-1 */ + {89UL, + 23UL, + 0UL}, +/* 2^12-1 */ + {1365UL, + 819UL, + 585UL, + 315UL, + 0UL}, +/* 2^13-1 */ + {1UL, + 0UL}, +/* 2^14-1 */ + {5461UL, + 381UL, + 129UL, + 0UL}, +/* 2^15-1 */ + {4681UL, + 1057UL, + 217UL, + 0UL}, +/* 2^16-1 */ + {21845UL, + 13107UL, + 3855UL, + 255UL, + 0UL}, +/* 2^17-1 */ + {1UL, + 0UL}, +/* 2^18-1 */ + {87381UL, + 37449UL, + 13797UL, + 3591UL, + 0UL}, +/* 2^19-1 */ + {1UL, + 0UL}, +/* 2^20-1 */ + {349525UL, + 209715UL, + 95325UL, + 33825UL, + 25575UL, + 0UL}, +/* 2^21-1 */ + {299593UL, + 16513UL, + 6223UL, + 0UL}, +/* 2^22-1 */ + {1398101UL, + 182361UL, + 47127UL, + 6141UL, + 0UL}, +/* 2^23-1 */ + {178481UL, + 47UL, + 0UL}, +/* 2^24-1 */ + {5592405UL, + 3355443UL, + 2396745UL, + 1290555UL, + 986895UL, + 69615UL, + 0UL}, +/* 2^25-1 */ + {1082401UL, + 55831UL, + 18631UL, + 0UL}, +/* 2^26-1 */ + {22369621UL, + 24573UL, + 8193UL, + 0UL}, +/* 2^27-1 */ + {19173961UL, + 1838599UL, + 511UL, + 0UL}, +/* 2^28-1 */ + {89478485UL, + 53687091UL, + 9256395UL, + 6242685UL, + 2375535UL, + 2113665UL, + 0UL}, +/* 2^29-1 */ + {2304167UL, + 486737UL, + 256999UL, + 0UL}, +/* 2^30-1 */ + {357913941UL, + 153391689UL, + 97612893UL, + 34636833UL, + 7110873UL, + 3243933UL, + 0UL}, +/* 2^31-1 */ + {1UL, + 0UL}, +/* 2^32-1 */ + {1431655765UL, + 858993459UL, + 252645135UL, + 16711935UL, + 65535UL, + 0UL}, +#if LONG_BIT > 32 +/* 2^33-1 */ + {1227133513UL, + 373475417UL, + 96516119UL, + 14329UL, + 0UL}, +/* 2^34-1 */ + {5726623061UL, + 393213UL, + 131073UL, + 0UL}, +/* 2^35-1 */ + {1108378657UL, + 483939977UL, + 270549121UL, + 279527UL, + 0UL}, +/* 2^36-1 */ + {22906492245UL, + 13743895347UL, + 9817068105UL, + 5286113595UL, + 3616814565UL, + 1857283155UL, + 941362695UL, + 630453915UL, + 0UL}, +/* 2^37-1 */ + {616318177UL, + 223UL, + 0UL}, +/* 2^38-1 */ + {91625968981UL, + 1572861UL, + 524289UL, + 0UL}, +/* 2^39-1 */ + {78536544841UL, + 6958934353UL, + 67117057UL, + 4529623UL, + 0UL}, +/* 2^40-1 */ + {366503875925UL, + 219902325555UL, + 99955602525UL, + 64677154575UL, + 35468117025UL, + 26817356775UL, + 17825775UL, + 0UL}, +/* 2^41-1 */ + {164511353UL, + 13367UL, + 0UL}, +/* 2^42-1 */ + {1466015503701UL, + 628292358729UL, + 102280151421UL, + 34630287489UL, + 13050583119UL, + 811597437UL, + 0UL}, +/* 2^43-1 */ + {20408568497UL, + 905040953UL, + 4188889UL, + 0UL}, +/* 2^44-1 */ + {5864062014805UL, + 3518437208883UL, + 764877654105UL, + 197665011735UL, + 44312811195UL, + 25757227005UL, + 8325691455UL, + 0UL}, +/* 2^45-1 */ + {5026338869833UL, + 1134979744801UL, + 481977699847UL, + 233009086681UL, + 55759702201UL, + 1509346321UL, + 0UL}, +/* 2^46-1 */ + {23456248059221UL, + 1497207322929UL, + 394264623UL, + 25165821UL, + 0UL}, +/* 2^47-1 */ + {59862819377UL, + 31184907679UL, + 10610063UL, + 0UL}, +/* 2^48-1 */ + {93824992236885UL, + 56294995342131UL, + 40210710958665UL, + 21651921285435UL, + 16557351571215UL, + 2901803883615UL, + 1167945961455UL, + 1095233372415UL, + 418239192735UL, + 0UL}, +/* 2^49-1 */ + {4432676798593UL, + 127UL, + 0UL}, +/* 2^50-1 */ + {375299968947541UL, + 102354536985693UL, + 36319351833633UL, + 4485656999373UL, + 1873377548823UL, + 625152641223UL, + 277931351973UL, + 0UL}, +/* 2^51-1 */ + {321685687669321UL, + 21862134113449UL, + 1050769861729UL, + 202518195313UL, + 17180000257UL, + 0UL}, +/* 2^52-1 */ + {1501199875790165UL, + 900719925474099UL, + 84973577874915UL, + 28685347945035UL, + 2792064245115UL, + 1649066139645UL, + 549822930945UL, + 0UL}, +/* 2^53-1 */ + {1416003655831UL, + 129728784761UL, + 441650591UL, + 0UL}, +/* 2^54-1 */ + {6004799503160661UL, + 2573485501354569UL, + 948126237341157UL, + 246772582321671UL, + 206561081853UL, + 68585259519UL, + 0UL}, +/* 2^55-1 */ + {1566469435607129UL, + 1162219258676257UL, + 404817944033303UL, + 40895342813807UL, + 11290754314937UL, + 178394823847UL, + 0UL}, +/* 2^56-1 */ + {24019198012642645UL, + 14411518807585587UL, + 4238682002231055UL, + 2484744621997515UL, + 1675758000882045UL, + 637677823344495UL, + 567382630219905UL, + 4563402735UL, + 0UL}, +/* 2^57-1 */ + {20587884010836553UL, + 4451159405623UL, + 274878431233UL, + 118823881393UL, + 0UL}, +/* 2^58-1 */ + {96076792050570581UL, + 4885260612740877UL, + 1237040240994471UL, + 261314937580881UL, + 137975287770087UL, + 95026151247UL, + 0UL}, +/* 2^59-1 */ + {3203431780337UL, + 179951UL, + 0UL}, +/* 2^60-1 */ + {384307168202282325UL, + 230584300921369395UL, + 164703072086692425UL, + 104811045873349725UL, + 88686269585142075UL, + 37191016277640225UL, + 28120036697727975UL, + 18900352534538475UL, + 7635241752363225UL, + 3483146539597725UL, + 872764197279975UL, + 0UL}, +/* 2^61-1 */ + {1UL, + 0UL}, +/* 2^62-1 */ + {1537228672809129301UL, + 6442450941UL, + 2147483649UL, + 0UL}, +/* 2^63-1 */ + {1317624576693539401UL, + 126347562148695559UL, + 72624976668147841UL, + 27369056489183311UL, + 99457304386111UL, + 14197294936951UL, + 0UL}, +/* 2^64-1 */ + {6148914691236517205UL, + 3689348814741910323UL, + 1085102592571150095UL, + 71777214294589695UL, + 28778071877862015UL, + 281470681808895UL, + 2753074036095UL, + 0UL}, +#if LONG_BIT > 64 +#error Factorization of numbers up to 2^LONG_BIT required +#endif +#endif +}; + +/* + * Modular multiply for two binary polynomial + * mask is 1UL << the degree of the modulus. + */ +unsigned long modmul(unsigned long poly1, unsigned long poly2, + unsigned long modulo, unsigned long mask) +{ + unsigned long result = 0; + + for (; poly1; poly1 >>= 1) + { + if (poly1 & 1) + result ^= poly2; + + poly2 <<= 1; + if (poly2 & mask) + poly2 ^= modulo; + } + return result; +} + +/* + * Modular exponentiation for a binary polynomial + * degree is the degree of the modulus. + */ +unsigned long modpow(unsigned long polynomial, unsigned long power, + unsigned long modulo, int degree) +{ + unsigned long result = 1, mask = 1UL << degree; + + for (; power; power >>= 1) + { + if (power & 1) + result = modmul(result, polynomial, modulo, mask); + polynomial = modmul(polynomial, polynomial, modulo, mask); + } + return result; +} + +/* + * Test the primitivity of a polynomial + */ +int rk_isprimitive(unsigned long polynomial) +{ + unsigned long pelement = 2, temp = polynomial >> 1; + int k, degree = 0, weight = 1; + + /* Special case for polynomials of degree < 2 */ + if (polynomial < 4) + return (polynomial == 3) || (polynomial == 1); + + /* A binary primitive polynomial has a constant term */ + if (!(polynomial & 1)) + return 0; + + /* + * A binary primitive polynomial of degree > 1 has an odd number of terms. + * temp ^= temp >> 16; temp ^= temp >> 8; ... would be sligthly faster. + * Compute the degree at the same time. + */ + for (; temp; degree++, temp >>= 1) + weight += temp & 1; + if (!(weight & 1)) + return 0; + + /* + * Check if the period is 2^degree-1. + * Sufficient if 2^degree-1 is prime. + */ + if (modpow(pelement, 1UL << degree, polynomial, degree) != pelement) + return 0; + + if (divisors[degree][0] != 1) + /* Primitivity test */ + for (k = 0; divisors[degree][k]; k++) + if (modpow(pelement, divisors[degree][k], polynomial, degree) == 1) + return 0; + + return 1; +} diff --git a/tardis/montecarlo/src/randomkit/rk_sobol.c b/tardis/montecarlo/src/randomkit/rk_sobol.c new file mode 100644 index 00000000000..c6a016bb033 --- /dev/null +++ b/tardis/montecarlo/src/randomkit/rk_sobol.c @@ -0,0 +1,991 @@ +/* Random kit 1.6 */ + +/* + * Copyright (c) 2004-2006, Jean-Sebastien Roy (js@jeannot.org) + * + * Original algorithm from Numerical Recipes, 2nd edition, by Press et al. + * The inverse normal cdf formulas are from Peter J. Acklam. + * The initialization directions were found in Ferdinando Ametrano's + * implementation in QuantLib. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. + * IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY + * CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +static char const rcsid[] = + "@(#) $Jeannot: rk_sobol.c,v 1.9 2006/02/19 13:48:34 js Exp $"; + +#include +#include +#include +#include "rk_sobol.h" +#include "rk_mt.h" +#include "rk_primitive.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif +#ifndef M_SQRT1_2 +#define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */ +#endif +#define RK_SOBOL_M_SQRT2PI 2.506628274631000502415 /* sqrt(2*pi) */ + +#ifndef LONG_BIT +#if ULONG_MAX <= 0xffffffffUL +#define LONG_BIT 32 +#else +#define LONG_BIT 64 +#endif +#endif + +char *rk_sobol_strerror[] = +{ + "no error", + "invalid dimension", + "too many numbers generated", + "not enough memory" +}; + +static double inverse_normal(double p); + +/* + * Sobol/Levitan coefficients of the free direction integers as given + * by Bratley, P., Fox, B.L. (1988) + */ + +const unsigned long rk_sobol_SLdirections[] = { + 1, + 1, 1, + 1, 3, 7, + 1, 1, 5, + 1, 3, 1, 1, + 1, 1, 3, 7, + 1, 3, 3, 9, 9, + 1, 3, 7, 13, 3, + 1, 1, 5, 11, 27, + 1, 3, 5, 1, 15, + 1, 1, 7, 3, 29, + 1, 3, 7, 7, 21, + 1, 1, 1, 9, 23, 37, + 1, 3, 3, 5, 19, 33, + 1, 1, 3, 13, 11, 7, + 1, 1, 7, 13, 25, 5, + 1, 3, 5, 11, 7, 11, + 1, 1, 1, 3, 13, 39, + 1, 3, 1, 15, 17, 63, 13, + 1, 1, 5, 5, 1, 27, 33, + 1, 3, 3, 3, 25, 17, 115, + 1, 1, 3, 15, 29, 15, 41, + 1, 3, 1, 7, 3, 23, 79, + 1, 3, 7, 9, 31, 29, 17, + 1, 1, 5, 13, 11, 3, 29, + 1, 3, 1, 9, 5, 21, 119, + 1, 1, 3, 1, 23, 13, 75, + 1, 3, 3, 11, 27, 31, 73, + 1, 1, 7, 7, 19, 25, 105, + 1, 3, 5, 5, 21, 9, 7, + 1, 1, 1, 15, 5, 49, 59, + 1, 1, 1, 1, 1, 33, 65, + 1, 3, 5, 15, 17, 19, 21, + 1, 1, 7, 11, 13, 29, 3, + 1, 3, 7, 5, 7, 11, 113, + 1, 1, 5, 3, 15, 19, 61, + 1, 3, 1, 1, 9, 27, 89, 7, + 1, 1, 3, 7, 31, 15, 45, 23, + 1, 3, 3, 9, 9, 25, 107, 39, + 0 +}; + +/* + * Lemieux coefficients of the free direction integers as given + * in QuantLib by Christiane Lemieux, private communication, September 2004 + */ + +const unsigned long rk_sobol_Ldirections[] = { + 1, + 1, 1, + 1, 3, 7, + 1, 1, 5, + 1, 3, 1, 1, + 1, 1, 3, 7, + 1, 3, 3, 9, 9, + 1, 3, 7, 13, 3, + 1, 1, 5, 11, 27, + 1, 3, 5, 1, 15, + 1, 1, 7, 3, 29, + 1, 3, 7, 7, 21, + 1, 1, 1, 9, 23, 37, + 1, 3, 3, 5, 19, 33, + 1, 1, 3, 13, 11, 7, + 1, 1, 7, 13, 25, 5, + 1, 3, 5, 11, 7, 11, + 1, 1, 1, 3, 13, 39, + 1, 3, 1, 15, 17, 63, 13, + 1, 1, 5, 5, 1, 27, 33, + 1, 3, 3, 3, 25, 17, 115, + 1, 1, 3, 15, 29, 15, 41, + 1, 3, 1, 7, 3, 23, 79, + 1, 3, 7, 9, 31, 29, 17, + 1, 1, 5, 13, 11, 3, 29, + 1, 3, 1, 9, 5, 21, 119, + 1, 1, 3, 1, 23, 13, 75, + 1, 3, 3, 11, 27, 31, 73, + 1, 1, 7, 7, 19, 25, 105, + 1, 3, 5, 5, 21, 9, 7, + 1, 1, 1, 15, 5, 49, 59, + 1, 1, 1, 1, 1, 33, 65, + 1, 3, 5, 15, 17, 19, 21, + 1, 1, 7, 11, 13, 29, 3, + 1, 3, 7, 5, 7, 11, 113, + 1, 1, 5, 3, 15, 19, 61, + 1, 3, 1, 1, 9, 27, 89, 7, + 1, 1, 3, 7, 31, 15, 45, 23, + 1, 3, 3, 9, 9, 25, 107, 39, + 1, 1, 3, 13, 7, 35, 61, 91, + 1, 1, 7, 11, 5, 35, 55, 75, + 1, 3, 5, 5, 11, 23, 29, 139, + 1, 1, 1, 7, 11, 15, 17, 81, + 1, 1, 7, 9, 5, 57, 79, 103, + 1, 1, 7, 13, 19, 5, 5, 185, + 1, 3, 1, 3, 13, 57, 97, 131, + 1, 1, 5, 5, 21, 25, 125, 197, + 1, 3, 3, 9, 31, 11, 103, 201, + 1, 1, 5, 3, 7, 25, 51, 121, + 1, 3, 7, 15, 19, 53, 73, 189, + 1, 1, 1, 15, 19, 55, 27, 183, + 1, 1, 7, 13, 3, 29, 109, 69, + 1, 1, 5, 15, 15, 23, 15, 1, 57, + 1, 3, 1, 3, 23, 55, 43, 143, 397, + 1, 1, 3, 11, 29, 9, 35, 131, 411, + 1, 3, 1, 7, 27, 39, 103, 199, 277, + 1, 3, 7, 3, 19, 55, 127, 67, 449, + 1, 3, 7, 3, 5, 29, 45, 85, 3, + 1, 3, 5, 5, 13, 23, 75, 245, 453, + 1, 3, 1, 15, 21, 47, 3, 77, 165, + 1, 1, 7, 9, 15, 5, 117, 73, 473, + 1, 3, 1, 9, 1, 21, 13, 173, 313, + 1, 1, 7, 3, 11, 45, 63, 77, 49, + 1, 1, 1, 1, 1, 25, 123, 39, 259, + 1, 1, 1, 5, 23, 11, 59, 11, 203, + 1, 3, 3, 15, 21, 1, 73, 71, 421, + 1, 1, 5, 11, 15, 31, 115, 95, 217, + 1, 1, 3, 3, 7, 53, 37, 43, 439, + 1, 1, 1, 1, 27, 53, 69, 159, 321, + 1, 1, 5, 15, 29, 17, 19, 43, 449, + 1, 1, 3, 9, 1, 55, 121, 205, 255, + 1, 1, 3, 11, 9, 47, 107, 11, 417, + 1, 1, 1, 5, 17, 25, 21, 83, 95, + 1, 3, 5, 13, 31, 25, 61, 157, 407, + 1, 1, 7, 9, 25, 33, 41, 35, 17, + 1, 3, 7, 15, 13, 39, 61, 187, 461, + 1, 3, 7, 13, 5, 57, 23, 177, 435, + 1, 1, 3, 15, 11, 27, 115, 5, 337, + 1, 3, 7, 3, 15, 63, 61, 171, 339, + 1, 3, 3, 13, 15, 61, 59, 47, 1, + 1, 1, 5, 15, 13, 5, 39, 83, 329, + 1, 1, 5, 5, 5, 27, 25, 39, 301, + 1, 1, 5, 11, 31, 41, 35, 233, 27, + 1, 3, 5, 15, 7, 37, 119, 171, 419, + 1, 3, 5, 5, 3, 29, 21, 189, 417, + 1, 1, 1, 1, 21, 41, 117, 119, 351, + 1, 1, 3, 1, 7, 27, 87, 19, 213, + 1, 1, 1, 1, 17, 7, 97, 217, 477, + 1, 1, 7, 1, 29, 61, 103, 231, 269, + 1, 1, 7, 13, 9, 27, 107, 207, 311, + 1, 1, 7, 5, 25, 21, 107, 179, 423, + 1, 3, 5, 11, 7, 1, 17, 245, 281, + 1, 3, 5, 9, 1, 5, 53, 59, 125, + 1, 1, 7, 1, 31, 57, 71, 245, 125, + 1, 1, 7, 5, 5, 57, 53, 253, 441, + 1, 3, 1, 13, 19, 35, 119, 235, 381, + 1, 3, 1, 7, 19, 59, 115, 33, 361, + 1, 1, 3, 5, 13, 1, 49, 143, 501, + 1, 1, 3, 5, 1, 63, 101, 85, 189, + 1, 1, 5, 11, 27, 63, 13, 131, 5, + 1, 1, 5, 7, 15, 45, 75, 59, 455, 585, + 1, 3, 1, 3, 7, 7, 111, 23, 119, 959, + 1, 3, 3, 9, 11, 41, 109, 163, 161, 879, + 1, 3, 5, 1, 21, 41, 121, 183, 315, 219, + 1, 1, 3, 9, 15, 3, 9, 223, 441, 929, + 1, 1, 7, 9, 3, 5, 93, 57, 253, 457, + 1, 1, 7, 13, 15, 29, 83, 21, 35, 45, + 1, 1, 3, 7, 13, 61, 119, 219, 85, 505, + 1, 1, 3, 3, 17, 13, 35, 197, 291, 109, + 1, 1, 3, 3, 5, 1, 113, 103, 217, 253, + 1, 1, 7, 1, 15, 39, 63, 223, 17, 9, + 1, 3, 7, 1, 17, 29, 67, 103, 495, 383, + 1, 3, 3, 15, 31, 59, 75, 165, 51, 913, + 1, 3, 7, 9, 5, 27, 79, 219, 233, 37, + 1, 3, 5, 15, 1, 11, 15, 211, 417, 811, + 1, 3, 5, 3, 29, 27, 39, 137, 407, 231, + 1, 1, 3, 5, 29, 43, 125, 135, 109, 67, + 1, 1, 1, 5, 11, 39, 107, 159, 323, 381, + 1, 1, 1, 1, 9, 11, 33, 55, 169, 253, + 1, 3, 5, 5, 11, 53, 63, 101, 251, 897, + 1, 3, 7, 1, 25, 15, 83, 119, 53, 157, + 1, 3, 5, 13, 5, 5, 3, 195, 111, 451, + 1, 3, 1, 15, 11, 1, 19, 11, 307, 777, + 1, 3, 7, 11, 5, 5, 17, 231, 345, 981, + 1, 1, 3, 3, 1, 33, 83, 201, 57, 475, + 1, 3, 7, 7, 17, 13, 35, 175, 499, 809, + 1, 1, 5, 3, 3, 17, 103, 119, 499, 865, + 1, 1, 1, 11, 27, 25, 37, 121, 401, 11, + 1, 1, 1, 11, 9, 25, 25, 241, 403, 3, + 1, 1, 1, 1, 11, 1, 39, 163, 231, 573, + 1, 1, 1, 13, 13, 21, 75, 185, 99, 545, + 1, 1, 1, 15, 3, 63, 69, 11, 173, 315, + 1, 3, 5, 15, 11, 3, 95, 49, 123, 765, + 1, 1, 1, 15, 3, 63, 77, 31, 425, 711, + 1, 1, 7, 15, 1, 37, 119, 145, 489, 583, + 1, 3, 5, 15, 3, 49, 117, 211, 165, 323, + 1, 3, 7, 1, 27, 63, 77, 201, 225, 803, + 1, 1, 1, 11, 23, 35, 67, 21, 469, 357, + 1, 1, 7, 7, 9, 7, 25, 237, 237, 571, + 1, 1, 3, 15, 29, 5, 107, 109, 241, 47, + 1, 3, 5, 11, 27, 63, 29, 13, 203, 675, + 1, 1, 3, 9, 9, 11, 103, 179, 449, 263, + 1, 3, 5, 11, 29, 63, 53, 151, 259, 223, + 1, 1, 3, 7, 9, 25, 5, 197, 237, 163, + 1, 3, 7, 13, 5, 57, 67, 193, 147, 241, + 1, 1, 5, 15, 15, 33, 17, 67, 161, 341, + 1, 1, 3, 13, 17, 43, 21, 197, 441, 985, + 1, 3, 1, 5, 15, 33, 33, 193, 305, 829, + 1, 1, 1, 13, 19, 27, 71, 187, 477, 239, + 1, 1, 1, 9, 9, 17, 41, 177, 229, 983, + 1, 3, 5, 9, 15, 45, 97, 205, 43, 767, + 1, 1, 1, 9, 31, 31, 77, 159, 395, 809, + 1, 3, 3, 3, 29, 19, 73, 123, 165, 307, + 1, 3, 1, 7, 5, 11, 77, 227, 355, 403, + 1, 3, 5, 5, 25, 31, 1, 215, 451, 195, + 1, 3, 7, 15, 29, 37, 101, 241, 17, 633, + 1, 1, 5, 1, 11, 3, 107, 137, 489, 5, + 1, 1, 1, 7, 19, 19, 75, 85, 471, 355, + 1, 1, 3, 3, 9, 13, 113, 167, 13, 27, + 1, 3, 5, 11, 21, 3, 89, 205, 377, 307, + 1, 1, 1, 9, 31, 61, 65, 9, 391, 141, 867, + 1, 1, 1, 9, 19, 19, 61, 227, 241, 55, 161, + 1, 1, 1, 11, 1, 19, 7, 233, 463, 171, 1941, + 1, 1, 5, 7, 25, 13, 103, 75, 19, 1021, 1063, + 1, 1, 1, 15, 17, 17, 79, 63, 391, 403, 1221, + 1, 3, 3, 11, 29, 25, 29, 107, 335, 475, 963, + 1, 3, 5, 1, 31, 33, 49, 43, 155, 9, 1285, + 1, 1, 5, 5, 15, 47, 39, 161, 357, 863, 1039, + 1, 3, 7, 15, 1, 39, 47, 109, 427, 393, 1103, + 1, 1, 1, 9, 9, 29, 121, 233, 157, 99, 701, + 1, 1, 1, 7, 1, 29, 75, 121, 439, 109, 993, + 1, 1, 1, 9, 5, 1, 39, 59, 89, 157, 1865, + 1, 1, 5, 1, 3, 37, 89, 93, 143, 533, 175, + 1, 1, 3, 5, 7, 33, 35, 173, 159, 135, 241, + 1, 1, 1, 15, 17, 37, 79, 131, 43, 891, 229, + 1, 1, 1, 1, 1, 35, 121, 177, 397, 1017, 583, + 1, 1, 3, 15, 31, 21, 43, 67, 467, 923, 1473, + 1, 1, 1, 7, 1, 33, 77, 111, 125, 771, 1975, + 1, 3, 7, 13, 1, 51, 113, 139, 245, 573, 503, + 1, 3, 1, 9, 21, 49, 15, 157, 49, 483, 291, + 1, 1, 1, 1, 29, 35, 17, 65, 403, 485, 1603, + 1, 1, 1, 7, 19, 1, 37, 129, 203, 321, 1809, + 1, 3, 7, 15, 15, 9, 5, 77, 29, 485, 581, + 1, 1, 3, 5, 15, 49, 97, 105, 309, 875, 1581, + 1, 3, 5, 1, 5, 19, 63, 35, 165, 399, 1489, + 1, 3, 5, 3, 23, 5, 79, 137, 115, 599, 1127, + 1, 1, 7, 5, 3, 61, 27, 177, 257, 91, 841, + 1, 1, 3, 5, 9, 31, 91, 209, 409, 661, 159, + 1, 3, 1, 15, 23, 39, 23, 195, 245, 203, 947, + 1, 1, 3, 1, 15, 59, 67, 95, 155, 461, 147, + 1, 3, 7, 5, 23, 25, 87, 11, 51, 449, 1631, + 1, 1, 1, 1, 17, 57, 7, 197, 409, 609, 135, + 1, 1, 1, 9, 1, 61, 115, 113, 495, 895, 1595, + 1, 3, 7, 15, 9, 47, 121, 211, 379, 985, 1755, + 1, 3, 1, 3, 7, 57, 27, 231, 339, 325, 1023, + 1, 1, 1, 1, 19, 63, 63, 239, 31, 643, 373, + 1, 3, 1, 11, 19, 9, 7, 171, 21, 691, 215, + 1, 1, 5, 13, 11, 57, 39, 211, 241, 893, 555, + 1, 1, 7, 5, 29, 21, 45, 59, 509, 223, 491, + 1, 1, 7, 9, 15, 61, 97, 75, 127, 779, 839, + 1, 1, 7, 15, 17, 33, 75, 237, 191, 925, 681, + 1, 3, 5, 7, 27, 57, 123, 111, 101, 371, 1129, + 1, 3, 5, 5, 29, 45, 59, 127, 229, 967, 2027, + 1, 1, 1, 1, 17, 7, 23, 199, 241, 455, 135, + 1, 1, 7, 15, 27, 29, 105, 171, 337, 503, 1817, + 1, 1, 3, 7, 21, 35, 61, 71, 405, 647, 2045, + 1, 1, 1, 1, 1, 15, 65, 167, 501, 79, 737, + 1, 1, 5, 1, 3, 49, 27, 189, 341, 615, 1287, + 1, 1, 1, 9, 1, 7, 31, 159, 503, 327, 1613, + 1, 3, 3, 3, 3, 23, 99, 115, 323, 997, 987, + 1, 1, 1, 9, 19, 33, 93, 247, 509, 453, 891, + 1, 1, 3, 1, 13, 19, 35, 153, 161, 633, 445, + 1, 3, 5, 15, 31, 5, 87, 197, 183, 783, 1823, + 1, 1, 7, 5, 19, 63, 69, 221, 129, 231, 1195, + 1, 1, 5, 5, 13, 23, 19, 231, 245, 917, 379, + 1, 3, 1, 15, 19, 43, 27, 223, 171, 413, 125, + 1, 1, 1, 9, 1, 59, 21, 15, 509, 207, 589, + 1, 3, 5, 3, 19, 31, 113, 19, 23, 733, 499, + 1, 1, 7, 1, 19, 51, 101, 165, 47, 925, 1093, + 1, 3, 3, 9, 15, 21, 43, 243, 237, 461, 1361, + 1, 1, 1, 9, 17, 15, 75, 75, 113, 715, 1419, + 1, 1, 7, 13, 17, 1, 99, 15, 347, 721, 1405, + 1, 1, 7, 15, 7, 27, 23, 183, 39, 59, 571, + 1, 3, 5, 9, 7, 43, 35, 165, 463, 567, 859, + 1, 3, 3, 11, 15, 19, 17, 129, 311, 343, 15, + 1, 1, 1, 15, 31, 59, 63, 39, 347, 359, 105, + 1, 1, 1, 15, 5, 43, 87, 241, 109, 61, 685, + 1, 1, 7, 7, 9, 39, 121, 127, 369, 579, 853, + 1, 1, 1, 1, 17, 15, 15, 95, 325, 627, 299, + 1, 1, 3, 13, 31, 53, 85, 111, 289, 811, 1635, + 1, 3, 7, 1, 19, 29, 75, 185, 153, 573, 653, + 1, 3, 7, 1, 29, 31, 55, 91, 249, 247, 1015, + 1, 3, 5, 7, 1, 49, 113, 139, 257, 127, 307, + 1, 3, 5, 9, 15, 15, 123, 105, 105, 225, 1893, + 1, 3, 3, 1, 15, 5, 105, 249, 73, 709, 1557, + 1, 1, 1, 9, 17, 31, 113, 73, 65, 701, 1439, + 1, 3, 5, 15, 13, 21, 117, 131, 243, 859, 323, + 1, 1, 1, 9, 19, 15, 69, 149, 89, 681, 515, + 1, 1, 1, 5, 29, 13, 21, 97, 301, 27, 967, + 1, 1, 3, 3, 15, 45, 107, 227, 495, 769, 1935, + 1, 1, 1, 11, 5, 27, 41, 173, 261, 703, 1349, + 1, 3, 3, 3, 11, 35, 97, 43, 501, 563, 1331, + 1, 1, 1, 7, 1, 17, 87, 17, 429, 245, 1941, + 1, 1, 7, 15, 29, 13, 1, 175, 425, 233, 797, + 1, 1, 3, 11, 21, 57, 49, 49, 163, 685, 701, + 1, 3, 3, 7, 11, 45, 107, 111, 379, 703, 1403, + 1, 1, 7, 3, 21, 7, 117, 49, 469, 37, 775, + 1, 1, 5, 15, 31, 63, 101, 77, 507, 489, 1955, + 1, 3, 3, 11, 19, 21, 101, 255, 203, 673, 665, + 1, 3, 3, 15, 17, 47, 125, 187, 271, 899, 2003, + 1, 1, 7, 7, 1, 35, 13, 235, 5, 337, 905, + 1, 3, 1, 15, 1, 43, 1, 27, 37, 695, 1429, + 1, 3, 1, 11, 21, 27, 93, 161, 299, 665, 495, + 1, 3, 3, 15, 3, 1, 81, 111, 105, 547, 897, + 1, 3, 5, 1, 3, 53, 97, 253, 401, 827, 1467, + 1, 1, 1, 5, 19, 59, 105, 125, 271, 351, 719, + 1, 3, 5, 13, 7, 11, 91, 41, 441, 759, 1827, + 1, 3, 7, 11, 29, 61, 61, 23, 307, 863, 363, + 1, 1, 7, 1, 15, 35, 29, 133, 415, 473, 1737, + 1, 1, 1, 13, 7, 33, 35, 225, 117, 681, 1545, + 1, 1, 1, 3, 5, 41, 83, 247, 13, 373, 1091, + 1, 3, 1, 13, 25, 61, 71, 217, 233, 313, 547, + 1, 3, 1, 7, 3, 29, 3, 49, 93, 465, 15, + 1, 1, 1, 9, 17, 61, 99, 163, 129, 485, 1087, + 1, 1, 1, 9, 9, 33, 31, 163, 145, 649, 253, + 1, 1, 1, 1, 17, 63, 43, 235, 287, 111, 567, + 1, 3, 5, 13, 29, 7, 11, 69, 153, 127, 449, + 1, 1, 5, 9, 11, 21, 15, 189, 431, 493, 1219, + 1, 1, 1, 15, 19, 5, 47, 91, 399, 293, 1743, + 1, 3, 3, 11, 29, 53, 53, 225, 409, 303, 333, + 1, 1, 1, 15, 31, 31, 21, 81, 147, 287, 1753, + 1, 3, 5, 5, 5, 63, 35, 125, 41, 687, 1793, + 1, 1, 1, 9, 19, 59, 107, 219, 455, 971, 297, + 1, 1, 3, 5, 3, 51, 121, 31, 245, 105, 1311, + 1, 3, 1, 5, 5, 57, 75, 107, 161, 431, 1693, + 1, 3, 1, 3, 19, 53, 27, 31, 191, 565, 1015, + 1, 3, 5, 13, 9, 41, 35, 249, 287, 49, 123, + 1, 1, 5, 7, 27, 17, 21, 3, 151, 885, 1165, + 1, 1, 7, 1, 15, 17, 65, 139, 427, 339, 1171, + 1, 1, 1, 5, 23, 5, 9, 89, 321, 907, 391, + 1, 1, 7, 9, 15, 1, 77, 71, 87, 701, 917, + 1, 1, 7, 1, 17, 37, 115, 127, 469, 779, 1543, + 1, 3, 7, 3, 5, 61, 15, 37, 301, 951, 1437, + 1, 1, 1, 13, 9, 51, 127, 145, 229, 55, 1567, + 1, 3, 7, 15, 19, 47, 53, 153, 295, 47, 1337, + 1, 3, 3, 5, 11, 31, 29, 133, 327, 287, 507, + 1, 1, 7, 7, 25, 31, 37, 199, 25, 927, 1317, + 1, 1, 7, 9, 3, 39, 127, 167, 345, 467, 759, + 1, 1, 1, 1, 31, 21, 15, 101, 293, 787, 1025, + 1, 1, 5, 3, 11, 41, 105, 109, 149, 837, 1813, + 1, 1, 3, 5, 29, 13, 19, 97, 309, 901, 753, + 1, 1, 7, 1, 19, 17, 31, 39, 173, 361, 1177, + 1, 3, 3, 3, 3, 41, 81, 7, 341, 491, 43, + 1, 1, 7, 7, 31, 35, 29, 77, 11, 335, 1275, + 1, 3, 3, 15, 17, 45, 19, 63, 151, 849, 129, + 1, 1, 7, 5, 7, 13, 47, 73, 79, 31, 499, + 1, 3, 1, 11, 1, 41, 59, 151, 247, 115, 1295, + 1, 1, 1, 9, 31, 37, 73, 23, 295, 483, 179, + 1, 3, 1, 15, 13, 63, 81, 27, 169, 825, 2037, + 1, 3, 5, 15, 7, 11, 73, 1, 451, 101, 2039, + 1, 3, 5, 3, 13, 53, 31, 137, 173, 319, 1521, + 1, 3, 1, 3, 29, 1, 73, 227, 377, 337, 1189, + 1, 3, 3, 13, 27, 9, 31, 101, 229, 165, 1983, + 1, 3, 1, 13, 13, 19, 19, 111, 319, 421, 223, + 1, 1, 7, 15, 25, 37, 61, 55, 359, 255, 1955, + 1, 1, 5, 13, 17, 43, 49, 215, 383, 915, 51, + 1, 1, 3, 1, 3, 7, 13, 119, 155, 585, 967, + 1, 3, 1, 13, 1, 63, 125, 21, 103, 287, 457, + 1, 1, 7, 1, 31, 17, 125, 137, 345, 379, 1925, + 1, 1, 3, 5, 5, 25, 119, 153, 455, 271, 2023, + 1, 1, 7, 9, 9, 37, 115, 47, 5, 255, 917, + 1, 3, 5, 3, 31, 21, 75, 203, 489, 593, 1, + 1, 3, 7, 15, 19, 63, 123, 153, 135, 977, 1875, + 1, 1, 1, 1, 5, 59, 31, 25, 127, 209, 745, + 1, 1, 1, 1, 19, 45, 67, 159, 301, 199, 535, + 1, 1, 7, 1, 31, 17, 19, 225, 369, 125, 421, + 1, 3, 3, 11, 7, 59, 115, 197, 459, 469, 1055, + 1, 3, 1, 3, 27, 45, 35, 131, 349, 101, 411, + 1, 3, 7, 11, 9, 3, 67, 145, 299, 253, 1339, + 1, 3, 3, 11, 9, 37, 123, 229, 273, 269, 515, + 1, 3, 7, 15, 11, 25, 75, 5, 367, 217, 951, + 1, 1, 3, 7, 9, 23, 63, 237, 385, 159, 1273, + 1, 1, 5, 11, 23, 5, 55, 193, 109, 865, 663, + 1, 1, 7, 15, 1, 57, 17, 141, 51, 217, 1259, + 1, 1, 3, 3, 15, 7, 89, 233, 71, 329, 203, + 1, 3, 7, 11, 11, 1, 19, 155, 89, 437, 573, + 1, 3, 1, 9, 27, 61, 47, 109, 161, 913, 1681, + 1, 1, 7, 15, 1, 33, 19, 15, 23, 913, 989, + 1, 3, 1, 1, 25, 39, 119, 193, 13, 571, 157, + 1, 1, 7, 13, 9, 55, 59, 147, 361, 935, 515, + 1, 1, 1, 9, 7, 59, 67, 117, 71, 855, 1493, + 1, 3, 1, 3, 13, 19, 57, 141, 305, 275, 1079, + 1, 1, 1, 9, 17, 61, 33, 7, 43, 931, 781, + 1, 1, 3, 1, 11, 17, 21, 97, 295, 277, 1721, + 1, 3, 1, 13, 15, 43, 11, 241, 147, 391, 1641, + 1, 1, 1, 1, 1, 19, 37, 21, 255, 263, 1571, + 1, 1, 3, 3, 23, 59, 89, 17, 475, 303, 757, 543, + 1, 3, 3, 9, 11, 55, 35, 159, 139, 203, 1531, 1825, + 1, 1, 5, 3, 17, 53, 51, 241, 269, 949, 1373, 325, + 1, 3, 7, 7, 5, 29, 91, 149, 239, 193, 1951, 2675, + 1, 3, 5, 1, 27, 33, 69, 11, 51, 371, 833, 2685, + 1, 1, 1, 15, 1, 17, 35, 57, 171, 1007, 449, 367, + 1, 1, 1, 7, 25, 61, 73, 219, 379, 53, 589, 4065, + 1, 3, 5, 13, 21, 29, 45, 19, 163, 169, 147, 597, + 1, 1, 5, 11, 21, 27, 7, 17, 237, 591, 255, 1235, + 1, 1, 7, 7, 17, 41, 69, 237, 397, 173, 1229, 2341, + 1, 1, 3, 1, 1, 33, 125, 47, 11, 783, 1323, 2469, + 1, 3, 1, 11, 3, 39, 35, 133, 153, 55, 1171, 3165, + 1, 1, 5, 11, 27, 23, 103, 245, 375, 753, 477, 2165, + 1, 3, 1, 15, 15, 49, 127, 223, 387, 771, 1719, 1465, + 1, 1, 1, 9, 11, 9, 17, 185, 239, 899, 1273, 3961, + 1, 1, 3, 13, 11, 51, 73, 81, 389, 647, 1767, 1215, + 1, 3, 5, 15, 19, 9, 69, 35, 349, 977, 1603, 1435, + 1, 1, 1, 1, 19, 59, 123, 37, 41, 961, 181, 1275, + 1, 1, 1, 1, 31, 29, 37, 71, 205, 947, 115, 3017, + 1, 1, 7, 15, 5, 37, 101, 169, 221, 245, 687, 195, + 1, 1, 1, 1, 19, 9, 125, 157, 119, 283, 1721, 743, + 1, 1, 7, 3, 1, 7, 61, 71, 119, 257, 1227, 2893, + 1, 3, 3, 3, 25, 41, 25, 225, 31, 57, 925, 2139, + 0 +}; + + +/* + * coefficients of the free direction integers as given in + * "Monte Carlo Methods in Finance", by Peter J�ckel, section 8.3 + */ + +const unsigned long rk_sobol_Jdirections[] = { + 1, + 1, 1, + 1, 3, 7, + 1, 1, 5, + 1, 3, 1, 1, + 1, 1, 3, 7, + 1, 3, 3, 9, 9, + 1, 3, 7, 7, 21, + 1, 1, 5, 11, 27, + 1, 1, 7, 3, 29, + 1, 3, 7, 13, 3, + 1, 3, 5, 1, 15, + 1, 1, 1, 9, 23, 37, + 1, 1, 3, 13, 11, 7, + 1, 3, 3, 5, 19, 33, + 1, 1, 7, 13, 25, 5, + 1, 1, 1, 3, 13, 39, + 1, 3, 5, 11, 7, 11, + 1, 3, 1, 7, 3, 23, 79, + 1, 3, 1, 15, 17, 63, 13, + 1, 3, 3, 3, 25, 17, 115, + 1, 3, 7, 9, 31, 29, 17, + 1, 1, 3, 15, 29, 15, 41, + 1, 3, 1, 9, 5, 21, 119, + 1, 1, 5, 5, 1, 27, 33, + 1, 1, 3, 1, 23, 13, 75, + 1, 1, 7, 7, 19, 25, 105, + 1, 3, 5, 5, 21, 9, 7, + 1, 1, 1, 15, 5, 49, 59, + 1, 3, 5, 15, 17, 19, 21, + 1, 1, 7, 11, 13, 29, 3, + 0 +}; + +/* + * 0 terminated list of primitive polynomials to speed up initialization + * All polynomials up to degree 13 (ie. 1111 polynomials) + */ +static const unsigned long rk_sobol_primitive_polynomials[] = { + 0x1UL, 0x3UL, 0x7UL, 0xBUL, 0xDUL, 0x13UL, 0x19UL, 0x25UL, 0x29UL, + 0x2FUL, 0x37UL, 0x3BUL, 0x3DUL, 0x43UL, 0x5BUL, 0x61UL, 0x67UL, 0x6DUL, + 0x73UL, 0x83UL, 0x89UL, 0x8FUL, 0x91UL, 0x9DUL, 0xA7UL, 0xABUL, 0xB9UL, + 0xBFUL, 0xC1UL, 0xCBUL, 0xD3UL, 0xD5UL, 0xE5UL, 0xEFUL, 0xF1UL, 0xF7UL, + 0xFDUL, 0x11DUL, 0x12BUL, 0x12DUL, 0x14DUL, 0x15FUL, 0x163UL, 0x165UL, + 0x169UL, 0x171UL, 0x187UL, 0x18DUL, 0x1A9UL, 0x1C3UL, 0x1CFUL, 0x1E7UL, + 0x1F5UL, 0x211UL, 0x21BUL, 0x221UL, 0x22DUL, 0x233UL, 0x259UL, 0x25FUL, + 0x269UL, 0x26FUL, 0x277UL, 0x27DUL, 0x287UL, 0x295UL, 0x2A3UL, 0x2A5UL, + 0x2AFUL, 0x2B7UL, 0x2BDUL, 0x2CFUL, 0x2D1UL, 0x2DBUL, 0x2F5UL, 0x2F9UL, + 0x313UL, 0x315UL, 0x31FUL, 0x323UL, 0x331UL, 0x33BUL, 0x34FUL, 0x35BUL, + 0x361UL, 0x36BUL, 0x36DUL, 0x373UL, 0x37FUL, 0x385UL, 0x38FUL, 0x3B5UL, + 0x3B9UL, 0x3C7UL, 0x3CBUL, 0x3CDUL, 0x3D5UL, 0x3D9UL, 0x3E3UL, 0x3E9UL, + 0x3FBUL, 0x409UL, 0x41BUL, 0x427UL, 0x42DUL, 0x465UL, 0x46FUL, 0x481UL, + 0x48BUL, 0x4C5UL, 0x4D7UL, 0x4E7UL, 0x4F3UL, 0x4FFUL, 0x50DUL, 0x519UL, + 0x523UL, 0x531UL, 0x53DUL, 0x543UL, 0x557UL, 0x56BUL, 0x585UL, 0x58FUL, + 0x597UL, 0x5A1UL, 0x5C7UL, 0x5E5UL, 0x5F7UL, 0x5FBUL, 0x613UL, 0x615UL, + 0x625UL, 0x637UL, 0x643UL, 0x64FUL, 0x65BUL, 0x679UL, 0x67FUL, 0x689UL, + 0x6B5UL, 0x6C1UL, 0x6D3UL, 0x6DFUL, 0x6FDUL, 0x717UL, 0x71DUL, 0x721UL, + 0x739UL, 0x747UL, 0x74DUL, 0x755UL, 0x759UL, 0x763UL, 0x77DUL, 0x78DUL, + 0x793UL, 0x7B1UL, 0x7DBUL, 0x7F3UL, 0x7F9UL, 0x805UL, 0x817UL, 0x82BUL, + 0x82DUL, 0x847UL, 0x863UL, 0x865UL, 0x871UL, 0x87BUL, 0x88DUL, 0x895UL, + 0x89FUL, 0x8A9UL, 0x8B1UL, 0x8CFUL, 0x8D1UL, 0x8E1UL, 0x8E7UL, 0x8EBUL, + 0x8F5UL, 0x90DUL, 0x913UL, 0x925UL, 0x929UL, 0x93BUL, 0x93DUL, 0x945UL, + 0x949UL, 0x951UL, 0x95BUL, 0x973UL, 0x975UL, 0x97FUL, 0x983UL, 0x98FUL, + 0x9ABUL, 0x9ADUL, 0x9B9UL, 0x9C7UL, 0x9D9UL, 0x9E5UL, 0x9F7UL, 0xA01UL, + 0xA07UL, 0xA13UL, 0xA15UL, 0xA29UL, 0xA49UL, 0xA61UL, 0xA6DUL, 0xA79UL, + 0xA7FUL, 0xA85UL, 0xA91UL, 0xA9DUL, 0xAA7UL, 0xAABUL, 0xAB3UL, 0xAB5UL, + 0xAD5UL, 0xADFUL, 0xAE9UL, 0xAEFUL, 0xAF1UL, 0xAFBUL, 0xB03UL, 0xB09UL, + 0xB11UL, 0xB33UL, 0xB3FUL, 0xB41UL, 0xB4BUL, 0xB59UL, 0xB5FUL, 0xB65UL, + 0xB6FUL, 0xB7DUL, 0xB87UL, 0xB8BUL, 0xB93UL, 0xB95UL, 0xBAFUL, 0xBB7UL, + 0xBBDUL, 0xBC9UL, 0xBDBUL, 0xBDDUL, 0xBE7UL, 0xBEDUL, 0xC0BUL, 0xC0DUL, + 0xC19UL, 0xC1FUL, 0xC57UL, 0xC61UL, 0xC6BUL, 0xC73UL, 0xC85UL, 0xC89UL, + 0xC97UL, 0xC9BUL, 0xC9DUL, 0xCB3UL, 0xCBFUL, 0xCC7UL, 0xCCDUL, 0xCD3UL, + 0xCD5UL, 0xCE3UL, 0xCE9UL, 0xCF7UL, 0xD03UL, 0xD0FUL, 0xD1DUL, 0xD27UL, + 0xD2DUL, 0xD41UL, 0xD47UL, 0xD55UL, 0xD59UL, 0xD63UL, 0xD6FUL, 0xD71UL, + 0xD93UL, 0xD9FUL, 0xDA9UL, 0xDBBUL, 0xDBDUL, 0xDC9UL, 0xDD7UL, 0xDDBUL, + 0xDE1UL, 0xDE7UL, 0xDF5UL, 0xE05UL, 0xE1DUL, 0xE21UL, 0xE27UL, 0xE2BUL, + 0xE33UL, 0xE39UL, 0xE47UL, 0xE4BUL, 0xE55UL, 0xE5FUL, 0xE71UL, 0xE7BUL, + 0xE7DUL, 0xE81UL, 0xE93UL, 0xE9FUL, 0xEA3UL, 0xEBBUL, 0xECFUL, 0xEDDUL, + 0xEF3UL, 0xEF9UL, 0xF0BUL, 0xF19UL, 0xF31UL, 0xF37UL, 0xF5DUL, 0xF6BUL, + 0xF6DUL, 0xF75UL, 0xF83UL, 0xF91UL, 0xF97UL, 0xF9BUL, 0xFA7UL, 0xFADUL, + 0xFB5UL, 0xFCDUL, 0xFD3UL, 0xFE5UL, 0xFE9UL, 0x1053UL, 0x1069UL, + 0x107BUL, 0x107DUL, 0x1099UL, 0x10D1UL, 0x10EBUL, 0x1107UL, 0x111FUL, + 0x1123UL, 0x113BUL, 0x114FUL, 0x1157UL, 0x1161UL, 0x116BUL, 0x1185UL, + 0x11B3UL, 0x11D9UL, 0x11DFUL, 0x120DUL, 0x1237UL, 0x123DUL, 0x1267UL, + 0x1273UL, 0x127FUL, 0x12B9UL, 0x12C1UL, 0x12CBUL, 0x130FUL, 0x131DUL, + 0x1321UL, 0x1339UL, 0x133FUL, 0x134DUL, 0x1371UL, 0x1399UL, 0x13A3UL, + 0x13A9UL, 0x1407UL, 0x1431UL, 0x1437UL, 0x144FUL, 0x145DUL, 0x1467UL, + 0x1475UL, 0x14A7UL, 0x14ADUL, 0x14D3UL, 0x150FUL, 0x151DUL, 0x154DUL, + 0x1593UL, 0x15C5UL, 0x15D7UL, 0x15DDUL, 0x15EBUL, 0x1609UL, 0x1647UL, + 0x1655UL, 0x1659UL, 0x16A5UL, 0x16BDUL, 0x1715UL, 0x1719UL, 0x1743UL, + 0x1745UL, 0x1775UL, 0x1789UL, 0x17ADUL, 0x17B3UL, 0x17BFUL, 0x17C1UL, + 0x1857UL, 0x185DUL, 0x1891UL, 0x1897UL, 0x18B9UL, 0x18EFUL, 0x191BUL, + 0x1935UL, 0x1941UL, 0x1965UL, 0x197BUL, 0x198BUL, 0x19B1UL, 0x19BDUL, + 0x19C9UL, 0x19CFUL, 0x19E7UL, 0x1A1BUL, 0x1A2BUL, 0x1A33UL, 0x1A69UL, + 0x1A8BUL, 0x1AD1UL, 0x1AE1UL, 0x1AF5UL, 0x1B0BUL, 0x1B13UL, 0x1B1FUL, + 0x1B57UL, 0x1B91UL, 0x1BA7UL, 0x1BBFUL, 0x1BC1UL, 0x1BD3UL, 0x1C05UL, + 0x1C11UL, 0x1C17UL, 0x1C27UL, 0x1C4DUL, 0x1C87UL, 0x1C9FUL, 0x1CA5UL, + 0x1CBBUL, 0x1CC5UL, 0x1CC9UL, 0x1CCFUL, 0x1CF3UL, 0x1D07UL, 0x1D23UL, + 0x1D43UL, 0x1D51UL, 0x1D5BUL, 0x1D75UL, 0x1D85UL, 0x1D89UL, 0x1E15UL, + 0x1E19UL, 0x1E2FUL, 0x1E45UL, 0x1E51UL, 0x1E67UL, 0x1E73UL, 0x1E8FUL, + 0x1EE3UL, 0x1F11UL, 0x1F1BUL, 0x1F27UL, 0x1F71UL, 0x1F99UL, 0x1FBBUL, + 0x1FBDUL, 0x1FC9UL, 0x201BUL, 0x2027UL, 0x2035UL, 0x2053UL, 0x2065UL, + 0x206FUL, 0x208BUL, 0x208DUL, 0x209FUL, 0x20A5UL, 0x20AFUL, 0x20BBUL, + 0x20BDUL, 0x20C3UL, 0x20C9UL, 0x20E1UL, 0x20F3UL, 0x210DUL, 0x2115UL, + 0x2129UL, 0x212FUL, 0x213BUL, 0x2143UL, 0x2167UL, 0x216BUL, 0x2179UL, + 0x2189UL, 0x2197UL, 0x219DUL, 0x21BFUL, 0x21C1UL, 0x21C7UL, 0x21CDUL, + 0x21DFUL, 0x21E3UL, 0x21F1UL, 0x21FBUL, 0x2219UL, 0x2225UL, 0x2237UL, + 0x223DUL, 0x2243UL, 0x225BUL, 0x225DUL, 0x2279UL, 0x227FUL, 0x2289UL, + 0x2297UL, 0x229BUL, 0x22B3UL, 0x22BFUL, 0x22CDUL, 0x22EFUL, 0x22F7UL, + 0x22FBUL, 0x2305UL, 0x2327UL, 0x232BUL, 0x2347UL, 0x2355UL, 0x2359UL, + 0x236FUL, 0x2371UL, 0x237DUL, 0x2387UL, 0x238DUL, 0x2395UL, 0x23A3UL, + 0x23A9UL, 0x23B1UL, 0x23B7UL, 0x23BBUL, 0x23E1UL, 0x23EDUL, 0x23F9UL, + 0x240BUL, 0x2413UL, 0x241FUL, 0x2425UL, 0x2429UL, 0x243DUL, 0x2451UL, + 0x2457UL, 0x2461UL, 0x246DUL, 0x247FUL, 0x2483UL, 0x249BUL, 0x249DUL, + 0x24B5UL, 0x24BFUL, 0x24C1UL, 0x24C7UL, 0x24CBUL, 0x24E3UL, 0x2509UL, + 0x2517UL, 0x251DUL, 0x2521UL, 0x252DUL, 0x2539UL, 0x2553UL, 0x2555UL, + 0x2563UL, 0x2571UL, 0x2577UL, 0x2587UL, 0x258BUL, 0x2595UL, 0x2599UL, + 0x259FUL, 0x25AFUL, 0x25BDUL, 0x25C5UL, 0x25CFUL, 0x25D7UL, 0x25EBUL, + 0x2603UL, 0x2605UL, 0x2611UL, 0x262DUL, 0x263FUL, 0x264BUL, 0x2653UL, + 0x2659UL, 0x2669UL, 0x2677UL, 0x267BUL, 0x2687UL, 0x2693UL, 0x2699UL, + 0x26B1UL, 0x26B7UL, 0x26BDUL, 0x26C3UL, 0x26EBUL, 0x26F5UL, 0x2713UL, + 0x2729UL, 0x273BUL, 0x274FUL, 0x2757UL, 0x275DUL, 0x276BUL, 0x2773UL, + 0x2779UL, 0x2783UL, 0x2791UL, 0x27A1UL, 0x27B9UL, 0x27C7UL, 0x27CBUL, + 0x27DFUL, 0x27EFUL, 0x27F1UL, 0x2807UL, 0x2819UL, 0x281FUL, 0x2823UL, + 0x2831UL, 0x283BUL, 0x283DUL, 0x2845UL, 0x2867UL, 0x2875UL, 0x2885UL, + 0x28ABUL, 0x28ADUL, 0x28BFUL, 0x28CDUL, 0x28D5UL, 0x28DFUL, 0x28E3UL, + 0x28E9UL, 0x28FBUL, 0x2909UL, 0x290FUL, 0x2911UL, 0x291BUL, 0x292BUL, + 0x2935UL, 0x293FUL, 0x2941UL, 0x294BUL, 0x2955UL, 0x2977UL, 0x297DUL, + 0x2981UL, 0x2993UL, 0x299FUL, 0x29AFUL, 0x29B7UL, 0x29BDUL, 0x29C3UL, + 0x29D7UL, 0x29F3UL, 0x29F5UL, 0x2A03UL, 0x2A0FUL, 0x2A1DUL, 0x2A21UL, + 0x2A33UL, 0x2A35UL, 0x2A4DUL, 0x2A69UL, 0x2A6FUL, 0x2A71UL, 0x2A7BUL, + 0x2A7DUL, 0x2AA5UL, 0x2AA9UL, 0x2AB1UL, 0x2AC5UL, 0x2AD7UL, 0x2ADBUL, + 0x2AEBUL, 0x2AF3UL, 0x2B01UL, 0x2B15UL, 0x2B23UL, 0x2B25UL, 0x2B2FUL, + 0x2B37UL, 0x2B43UL, 0x2B49UL, 0x2B6DUL, 0x2B7FUL, 0x2B85UL, 0x2B97UL, + 0x2B9BUL, 0x2BADUL, 0x2BB3UL, 0x2BD9UL, 0x2BE5UL, 0x2BFDUL, 0x2C0FUL, + 0x2C21UL, 0x2C2BUL, 0x2C2DUL, 0x2C3FUL, 0x2C41UL, 0x2C4DUL, 0x2C71UL, + 0x2C8BUL, 0x2C8DUL, 0x2C95UL, 0x2CA3UL, 0x2CAFUL, 0x2CBDUL, 0x2CC5UL, + 0x2CD1UL, 0x2CD7UL, 0x2CE1UL, 0x2CE7UL, 0x2CEBUL, 0x2D0DUL, 0x2D19UL, + 0x2D29UL, 0x2D2FUL, 0x2D37UL, 0x2D3BUL, 0x2D45UL, 0x2D5BUL, 0x2D67UL, + 0x2D75UL, 0x2D89UL, 0x2D8FUL, 0x2DA7UL, 0x2DABUL, 0x2DB5UL, 0x2DE3UL, + 0x2DF1UL, 0x2DFDUL, 0x2E07UL, 0x2E13UL, 0x2E15UL, 0x2E29UL, 0x2E49UL, + 0x2E4FUL, 0x2E5BUL, 0x2E5DUL, 0x2E61UL, 0x2E6BUL, 0x2E8FUL, 0x2E91UL, + 0x2E97UL, 0x2E9DUL, 0x2EABUL, 0x2EB3UL, 0x2EB9UL, 0x2EDFUL, 0x2EFBUL, + 0x2EFDUL, 0x2F05UL, 0x2F09UL, 0x2F11UL, 0x2F17UL, 0x2F3FUL, 0x2F41UL, + 0x2F4BUL, 0x2F4DUL, 0x2F59UL, 0x2F5FUL, 0x2F65UL, 0x2F69UL, 0x2F95UL, + 0x2FA5UL, 0x2FAFUL, 0x2FB1UL, 0x2FCFUL, 0x2FDDUL, 0x2FE7UL, 0x2FEDUL, + 0x2FF5UL, 0x2FFFUL, 0x3007UL, 0x3015UL, 0x3019UL, 0x302FUL, 0x3049UL, + 0x304FUL, 0x3067UL, 0x3079UL, 0x307FUL, 0x3091UL, 0x30A1UL, 0x30B5UL, + 0x30BFUL, 0x30C1UL, 0x30D3UL, 0x30D9UL, 0x30E5UL, 0x30EFUL, 0x3105UL, + 0x310FUL, 0x3135UL, 0x3147UL, 0x314DUL, 0x315FUL, 0x3163UL, 0x3171UL, + 0x317BUL, 0x31A3UL, 0x31A9UL, 0x31B7UL, 0x31C5UL, 0x31C9UL, 0x31DBUL, + 0x31E1UL, 0x31EBUL, 0x31EDUL, 0x31F3UL, 0x31FFUL, 0x3209UL, 0x320FUL, + 0x321DUL, 0x3227UL, 0x3239UL, 0x324BUL, 0x3253UL, 0x3259UL, 0x3265UL, + 0x3281UL, 0x3293UL, 0x3299UL, 0x329FUL, 0x32A9UL, 0x32B7UL, 0x32BBUL, + 0x32C3UL, 0x32D7UL, 0x32DBUL, 0x32E7UL, 0x3307UL, 0x3315UL, 0x332FUL, + 0x3351UL, 0x335DUL, 0x3375UL, 0x3397UL, 0x339BUL, 0x33ABUL, 0x33B9UL, + 0x33C1UL, 0x33C7UL, 0x33D5UL, 0x33E3UL, 0x33E5UL, 0x33F7UL, 0x33FBUL, + 0x3409UL, 0x341BUL, 0x3427UL, 0x3441UL, 0x344DUL, 0x345FUL, 0x3469UL, + 0x3477UL, 0x347BUL, 0x3487UL, 0x3493UL, 0x3499UL, 0x34A5UL, 0x34BDUL, + 0x34C9UL, 0x34DBUL, 0x34E7UL, 0x34F9UL, 0x350DUL, 0x351FUL, 0x3525UL, + 0x3531UL, 0x3537UL, 0x3545UL, 0x354FUL, 0x355DUL, 0x356DUL, 0x3573UL, + 0x357FUL, 0x359DUL, 0x35A1UL, 0x35B9UL, 0x35CDUL, 0x35D5UL, 0x35D9UL, + 0x35E3UL, 0x35E9UL, 0x35EFUL, 0x3601UL, 0x360BUL, 0x361FUL, 0x3625UL, + 0x362FUL, 0x363BUL, 0x3649UL, 0x3651UL, 0x365BUL, 0x3673UL, 0x3675UL, + 0x3691UL, 0x369BUL, 0x369DUL, 0x36ADUL, 0x36CBUL, 0x36D3UL, 0x36D5UL, + 0x36E3UL, 0x36EFUL, 0x3705UL, 0x370FUL, 0x371BUL, 0x3721UL, 0x372DUL, + 0x3739UL, 0x3741UL, 0x3747UL, 0x3753UL, 0x3771UL, 0x3777UL, 0x378BUL, + 0x3795UL, 0x3799UL, 0x37A3UL, 0x37C5UL, 0x37CFUL, 0x37D1UL, 0x37D7UL, + 0x37DDUL, 0x37E1UL, 0x37F3UL, 0x3803UL, 0x3805UL, 0x3817UL, 0x381DUL, + 0x3827UL, 0x3833UL, 0x384BUL, 0x3859UL, 0x3869UL, 0x3871UL, 0x38A3UL, + 0x38B1UL, 0x38BBUL, 0x38C9UL, 0x38CFUL, 0x38E1UL, 0x38F3UL, 0x38F9UL, + 0x3901UL, 0x3907UL, 0x390BUL, 0x3913UL, 0x3931UL, 0x394FUL, 0x3967UL, + 0x396DUL, 0x3983UL, 0x3985UL, 0x3997UL, 0x39A1UL, 0x39A7UL, 0x39ADUL, + 0x39CBUL, 0x39CDUL, 0x39D3UL, 0x39EFUL, 0x39F7UL, 0x39FDUL, 0x3A07UL, + 0x3A29UL, 0x3A2FUL, 0x3A3DUL, 0x3A51UL, 0x3A5DUL, 0x3A61UL, 0x3A67UL, + 0x3A73UL, 0x3A75UL, 0x3A89UL, 0x3AB9UL, 0x3ABFUL, 0x3ACDUL, 0x3AD3UL, + 0x3AD5UL, 0x3ADFUL, 0x3AE5UL, 0x3AE9UL, 0x3AFBUL, 0x3B11UL, 0x3B2BUL, + 0x3B2DUL, 0x3B35UL, 0x3B3FUL, 0x3B53UL, 0x3B59UL, 0x3B63UL, 0x3B65UL, + 0x3B6FUL, 0x3B71UL, 0x3B77UL, 0x3B8BUL, 0x3B99UL, 0x3BA5UL, 0x3BA9UL, + 0x3BB7UL, 0x3BBBUL, 0x3BD1UL, 0x3BE7UL, 0x3BF3UL, 0x3BFFUL, 0x3C0DUL, + 0x3C13UL, 0x3C15UL, 0x3C1FUL, 0x3C23UL, 0x3C25UL, 0x3C3BUL, 0x3C4FUL, + 0x3C5DUL, 0x3C6DUL, 0x3C83UL, 0x3C8FUL, 0x3C9DUL, 0x3CA7UL, 0x3CABUL, + 0x3CB9UL, 0x3CC7UL, 0x3CE9UL, 0x3CFBUL, 0x3CFDUL, 0x3D03UL, 0x3D17UL, + 0x3D1BUL, 0x3D21UL, 0x3D2DUL, 0x3D33UL, 0x3D35UL, 0x3D41UL, 0x3D4DUL, + 0x3D65UL, 0x3D69UL, 0x3D7DUL, 0x3D81UL, 0x3D95UL, 0x3DB1UL, 0x3DB7UL, + 0x3DC3UL, 0x3DD1UL, 0x3DDBUL, 0x3DE7UL, 0x3DEBUL, 0x3DF9UL, 0x3E05UL, + 0x3E09UL, 0x3E0FUL, 0x3E1BUL, 0x3E2BUL, 0x3E3FUL, 0x3E41UL, 0x3E53UL, + 0x3E65UL, 0x3E69UL, 0x3E8BUL, 0x3EA3UL, 0x3EBDUL, 0x3EC5UL, 0x3ED7UL, + 0x3EDDUL, 0x3EE1UL, 0x3EF9UL, 0x3F0DUL, 0x3F19UL, 0x3F1FUL, 0x3F25UL, + 0x3F37UL, 0x3F3DUL, 0x3F43UL, 0x3F45UL, 0x3F49UL, 0x3F51UL, 0x3F57UL, + 0x3F61UL, 0x3F83UL, 0x3F89UL, 0x3F91UL, 0x3FABUL, 0x3FB5UL, 0x3FE3UL, + 0x3FF7UL, 0x3FFDUL, + 0UL +}; + +rk_sobol_error rk_sobol_init(size_t dimension, rk_sobol_state *s, + rk_state *rs_dir, const unsigned long *directions, + const unsigned long *polynomials) +{ + rk_state rs_dir_temp; + int j, l, degree = 0, last_degree = 0, ooord = 0; + size_t k, cdir = 0, cpol = 0; + unsigned long polynomial = 1, rev = 0, last = 0; + + if (dimension == 0) + return RK_SOBOL_EINVAL; + + if (polynomials == NULL) + polynomials = rk_sobol_primitive_polynomials; + + /* Allocate the structure */ + s->direction = NULL; s->numerator = NULL; + s->direction = malloc(sizeof(*(s->direction))*dimension*LONG_BIT); + s->numerator = malloc(sizeof(*(s->numerator))*dimension); + if (!s->direction | !s->numerator) + { + if (!s->direction) free(s->direction); + if (!s->numerator) free(s->numerator); + return RK_SOBOL_ENOMEM; + } + + /* Initialize directions */ + /* Degree 0 */ + for (j = degree; j < LONG_BIT; j++) + s->direction[j*dimension] = 1UL << (LONG_BIT-j-1); + + /* Skip unused first polynomial */ + if (polynomials[cpol]) + cpol++; + + /* Degree >0 */ + for (k = 1; k < dimension; k++) + { + unsigned long temp; + + /* Find a new primitive polynomial */ + if (polynomials[cpol]) + polynomial = polynomials[cpol++]; + else if (rev) + { + /* We are generating polynomials out of order: + use the reverse of the previous polynomial */ + last = polynomial; + polynomial = rev; + rev = 0; + } + else + { + if (last) + { + polynomial = last; + last = 0; + } + + /* Find a new primitive polynomial */ + while(1) + { + if (polynomial == ULONG_MAX) + { + /* Not enough polynomials */ + free(s->direction); + free(s->numerator); + return RK_SOBOL_EINVAL; + } + + polynomial += 2; + + if (ooord) + { + unsigned long copy = polynomial; + /* We are generating polynomials out of order: + check if the reverse was already checked */ + for (rev = 0; copy; copy >>= 1) + rev = (rev << 1) | (copy & 1); + if (ooord && rev < polynomial) + continue; + } + + if (rk_isprimitive(polynomial)) + break; + } + + if (rev == polynomial) + /* We are generating polynomials out of order: + the reverse is not different, discard it */ + rev = 0; + } + + /* Compute the degree */ + for (temp = polynomial >> 1, degree = 0; temp; degree++, temp >>= 1); + + for (j=0; jdirection[j*dimension+k] = m << (LONG_BIT-j-1); + } + + /* Scaled recursion for directions */ + for (j = degree; j < LONG_BIT; j++) + { + unsigned long effdir = s->direction[(j-degree)*dimension+k], + ptemp = polynomial >> 1; + effdir ^= (effdir >> degree); + for (l = degree-1; l >= 1; l--, ptemp >>= 1) + if (ptemp & 1) + effdir ^= s->direction[(j-l)*dimension+k]; + s->direction[j*dimension+k] = effdir; + } + + /* Can we generate polynomials out of order ? */ + if (!ooord && polynomials[cpol] == 0 && degree > last_degree + && (directions == NULL || directions[cdir] == 0)) + ooord = 0; + else + last_degree = degree; + } + + /* Initialize numerator */ + for (k=0; knumerator[k] = 0; + + s->dimension = dimension; + s->gcount = 0; + s->count = 0; + return RK_SOBOL_OK; +} + +void rk_sobol_reinit(rk_sobol_state *s) +{ + size_t k; + + /* Initialize numerator */ + for (k=0; kdimension; k++) + s->numerator[k] = 0; + + s->count = 0; + s->gcount = 0; +} + +void rk_sobol_randomshift(rk_sobol_state *s, rk_state *rs_num) +{ + rk_state rs_num_temp; + size_t k; + + if (rs_num == NULL) + { + rs_num = &rs_num_temp; + rk_randomseed(rs_num); + } + + /* Initialize numerator */ + for (k=0; kdimension; k++) + s->numerator[k] = rk_ulong(rs_num); +} + +rk_sobol_error rk_sobol_copy(rk_sobol_state *copy, rk_sobol_state *orig) +{ + size_t k; + + /* Allocate the structure */ + copy->direction = NULL; copy->numerator = NULL; + copy->direction = malloc(sizeof(*(copy->direction))*orig->dimension*LONG_BIT); + copy->numerator = malloc(sizeof(*(copy->numerator))*orig->dimension); + if (!copy->direction | !copy->numerator) + { + if (!copy->direction) free(copy->direction); + if (!copy->numerator) free(copy->numerator); + return RK_SOBOL_ENOMEM; + } + + /* Initialize numerator */ + for (k=0; kdimension; k++) + copy->numerator[k] = orig->numerator[k]; + for (k=0; k<(orig->dimension*LONG_BIT); k++) + copy->direction[k] = orig->direction[k]; + + copy->count = orig->count; + copy->gcount = orig->gcount; + copy->dimension = orig->dimension; + + return RK_SOBOL_OK; +} + +rk_sobol_error rk_sobol_double(rk_sobol_state *s, double *x) +{ + int j; + size_t k; + unsigned long im; + const double inverse_denominator=1.0/(ULONG_MAX+1.0); + + if (s->count == ULONG_MAX) + j = 0; + else + for (im = s->count, j=0; im & 1; j++, im >>= 1); + s->count++; + + for (k=0; kdimension; k++) + { + s->numerator[k] ^= s->direction[j*s->dimension+k]; + x[k] = s->numerator[k]*inverse_denominator; + } + + if ((s->gcount++) == ULONG_MAX) return RK_SOBOL_EXHAUST; + return RK_SOBOL_OK; +} + +void rk_sobol_setcount(rk_sobol_state *s, unsigned long count) +{ + s->count = count; +} + +void rk_sobol_free(rk_sobol_state *s) +{ + free(s->direction); + free(s->numerator); +} + +double inverse_normal(double p) +{ + double q, t, x; + + /* Peter J. Acklam constants for the rational approximation */ + const double a[6] = + { + -3.969683028665376e+01, 2.209460984245205e+02, + -2.759285104469687e+02, 1.383577518672690e+02, + -3.066479806614716e+01, 2.506628277459239e+00 + }; + const double b[5] = + { + -5.447609879822406e+01, 1.615858368580409e+02, + -1.556989798598866e+02, 6.680131188771972e+01, + -1.328068155288572e+01 + }; + const double c[6] = + { + -7.784894002430293e-03, -3.223964580411365e-01, + -2.400758277161838e+00, -2.549732539343734e+00, + 4.374664141464968e+00, 2.938163982698783e+00 + }; + const double d[4] = + { + 7.784695709041462e-03, 3.224671290700398e-01, + 2.445134137142996e+00, 3.754408661907416e+00 + }; + + if (p <= 0) + return -HUGE_VAL; + else if (p >= 1) + return HUGE_VAL; + + q = p<0.5 ? p : 1-p; + if (q > 0.02425) + { + /* Rational approximation for central region */ + x = q-0.5; + t = x*x; + x = x*(((((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4])*t+a[5]) + /(((((b[0]*t+b[1])*t+b[2])*t+b[3])*t+b[4])*t+1); + } + else + { + /* Rational approximation for tail region */ + t = sqrt(-2*log(q)); + x = (((((c[0]*t+c[1])*t+c[2])*t+c[3])*t+c[4])*t+c[5]) + /((((d[0]*t+d[1])*t+d[2])*t+d[3])*t+1); + } + + /* If we have erfc, improve the precision */ +#ifndef WIN32 + /* Halley's rational method */ + t = (erfc(-x*M_SQRT1_2)/2 - q) * RK_SOBOL_M_SQRT2PI * exp(x*x/2); + x -= t/(1 + x*t/2); +#endif + + return p>0.5 ? -x : x; +} + +rk_sobol_error rk_sobol_gauss(rk_sobol_state *s, double *x) +{ + size_t k; + rk_sobol_error rc = rk_sobol_double(s, x); + + for (k=0; kdimension; k++) + x[k] = inverse_normal(x[k]); + + return rc; +} From 4bce6db08ca9199eca3293c6a03f4b2babc09841 Mon Sep 17 00:00:00 2001 From: Michi Date: Thu, 11 Dec 2014 16:09:44 +0100 Subject: [PATCH 8/9] include .c in montecarlo/src/ for the randomkit files --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 7b4979ae3b9..e0f489f9d44 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ __pycache__ # Ignore .c files by default to avoid including generated code. If you want to # add a non-generated .c extension, use `git add -f filename.c`. *.c +!tardis/montecarlo/src/.c # Other generated files */version.py From 6a14e242072e696cfa50c193ecedaa5d92d110eb Mon Sep 17 00:00:00 2001 From: Michi Date: Thu, 11 Dec 2014 16:14:47 +0100 Subject: [PATCH 9/9] add a * to .c --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index e0f489f9d44..e4d132bcfd7 100644 --- a/.gitignore +++ b/.gitignore @@ -8,7 +8,7 @@ __pycache__ # Ignore .c files by default to avoid including generated code. If you want to # add a non-generated .c extension, use `git add -f filename.c`. *.c -!tardis/montecarlo/src/.c +!tardis/montecarlo/src/*.c # Other generated files */version.py