| 853 | } |
| 854 | |
| 855 | double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) { |
| 856 | double s; |
| 857 | double U, V, W, Y, Z; |
| 858 | double result, mod; |
| 859 | int neg; |
| 860 | if (npy_isnan(kappa)) { |
| 861 | return NPY_NAN; |
| 862 | } |
| 863 | if (kappa < 1e-8) { |
| 864 | /* Use a uniform for very small values of kappa */ |
| 865 | return M_PI * (2 * next_double(bitgen_state) - 1); |
| 866 | } else { |
| 867 | /* with double precision rho is zero until 1.4e-8 */ |
| 868 | if (kappa < 1e-5) { |
| 869 | /* |
| 870 | * second order taylor expansion around kappa = 0 |
| 871 | * precise until relatively large kappas as second order is 0 |
| 872 | */ |
| 873 | s = (1. / kappa + kappa); |
| 874 | } else { |
| 875 | if (kappa <= 1e6) { |
| 876 | /* Path for 1e-5 <= kappa <= 1e6 */ |
| 877 | double r = 1 + sqrt(1 + 4 * kappa * kappa); |
| 878 | double rho = (r - sqrt(2 * r)) / (2 * kappa); |
| 879 | s = (1 + rho * rho) / (2 * rho); |
| 880 | } else { |
| 881 | /* Fallback to wrapped normal distribution for kappa > 1e6 */ |
| 882 | result = mu + sqrt(1. / kappa) * random_standard_normal(bitgen_state); |
| 883 | /* Ensure result is within bounds */ |
| 884 | if (result < -M_PI) { |
| 885 | result += 2*M_PI; |
| 886 | } |
| 887 | if (result > M_PI) { |
| 888 | result -= 2*M_PI; |
| 889 | } |
| 890 | return result; |
| 891 | } |
| 892 | } |
| 893 | |
| 894 | while (1) { |
| 895 | U = next_double(bitgen_state); |
| 896 | Z = cos(M_PI * U); |
| 897 | W = (1 + s * Z) / (s + Z); |
| 898 | Y = kappa * (s - W); |
| 899 | V = next_double(bitgen_state); |
| 900 | /* |
| 901 | * V==0.0 is ok here since Y >= 0 always leads |
| 902 | * to accept, while Y < 0 always rejects |
| 903 | */ |
| 904 | if ((Y * (2 - Y) - V >= 0) || (log(Y / V) + 1 - Y >= 0)) { |
| 905 | break; |
| 906 | } |
| 907 | } |
| 908 | |
| 909 | U = next_double(bitgen_state); |
| 910 | |
| 911 | result = acos(W); |
| 912 | if (U < 0.5) { |
nothing calls this directly
no test coverage detected