diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2014, Christopher Chalmers + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Christopher Chalmers nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,14 @@ +## pcg-random + +[![Build Status](https://travis-ci.org/cchalmers/pcg-random.svg)](https://travis-ci.org/cchalmers/pcg-random) + +Haskell bindings to the PCG random number generator http://www.pcg-random.org. + +> PCG is a family of simple fast space-efficient statistically good algorithms for random number generation with better-than-typical cryptographic security + +Implements the standard multiple stream generator as well as the fast, single and unique variants. + +The api is very similar to [mwc-random] but the pcg generator appears to be faster. There is also a pure interface via the [random] libray. + +[mwc-random]: https://hackage.haskell.org/package/mwc-random +[random]: http://hackage.haskell.org/package/random diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/c/LICENSE.txt b/c/LICENSE.txt new file mode 100644 --- /dev/null +++ b/c/LICENSE.txt @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "{}" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright {yyyy} {name of copyright owner} + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/c/pcg-advance-64.c b/c/pcg-advance-64.c new file mode 100644 --- /dev/null +++ b/c/pcg-advance-64.c @@ -0,0 +1,62 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the canonical C++ PCG implementation, which + * has many additional features and is preferable if you can use C++ in + * your project. + * + * Repetative C code is derived using C preprocessor metaprogramming + * techniques. + */ + +#include "pcg_variants.h" + +/* Multi-step advance functions (jump-ahead, jump-back) + * + * The method used here is based on Brown, "Random Number Generation + * with Arbitrary Stride,", Transactions of the American Nuclear + * Society (Nov. 1994). The algorithm is very similar to fast + * exponentiation. + * + * Even though delta is an unsigned integer, we can pass a + * signed integer to go backwards, it just goes "the long way round". + */ + +uint64_t pcg_advance_lcg_64(uint64_t state, uint64_t delta, uint64_t cur_mult, + uint64_t cur_plus) +{ + uint64_t acc_mult = 1u; + uint64_t acc_plus = 0u; + while (delta > 0) { + if (delta & 1) { + acc_mult *= cur_mult; + acc_plus = acc_plus * cur_mult + cur_plus; + } + cur_plus = (cur_mult + 1) * cur_plus; + cur_mult *= cur_mult; + delta /= 2; + } + return acc_mult * state + acc_plus; +} + diff --git a/c/pcg-global-64.c b/c/pcg-global-64.c new file mode 100644 --- /dev/null +++ b/c/pcg-global-64.c @@ -0,0 +1,59 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the canonical C++ PCG implementation, which + * has many additional features and is preferable if you can use C++ in + * your project. + * + * The contents of this file were mechanically derived from pcg_variants.h + * (every inline function defined there gets an exern declaration here). + */ + +#include "pcg_variants.h" + +#if PCG_HAS_128BIT_OPS + +static pcg64_random_t pcg64_global = PCG64_INITIALIZER; + +uint64_t pcg64_random() +{ + return pcg64_random_r(&pcg64_global); +} + +uint64_t pcg64_boundedrand(uint64_t bound) +{ + return pcg64_boundedrand_r(&pcg64_global, bound); +} + +void pcg64_srandom(pcg128_t seed, pcg128_t seq) +{ + pcg64_srandom_r(&pcg64_global, seed, seq); +} + +void pcg64_advance(pcg128_t delta) +{ + pcg64_advance_r(&pcg64_global, delta); +} + +#endif diff --git a/c/pcg-output-64.c b/c/pcg-output-64.c new file mode 100644 --- /dev/null +++ b/c/pcg-output-64.c @@ -0,0 +1,70 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the canonical C++ PCG implementation, which + * has many additional features and is preferable if you can use C++ in + * your project. + * + * The contents of this file were mechanically derived from pcg_variants.h + * (every inline function defined there gets an exern declaration here). + */ + +#include "pcg_variants.h" + +/* + * Rotate helper functions. + */ + +extern inline uint64_t pcg_rotr_64(uint64_t value, unsigned int rot); + +/* + * Output functions. These are the core of the PCG generation scheme. + */ + +// XSH RS + +#if PCG_HAS_128BIT_OPS +extern inline uint64_t pcg_output_xsh_rs_128_64(pcg128_t state); +#endif + +// XSH RR + +#if PCG_HAS_128BIT_OPS +extern inline uint64_t pcg_output_xsh_rr_128_64(pcg128_t state); +#endif + +// RXS M XS + +extern inline uint64_t pcg_output_rxs_m_xs_64_64(uint64_t state); + +// XSL RR (only defined for >= 64 bits) + +#if PCG_HAS_128BIT_OPS +extern inline uint64_t pcg_output_xsl_rr_128_64(pcg128_t state); +#endif + +// XSL RR RR (only defined for >= 64 bits) + +extern inline uint64_t pcg_output_xsl_rr_rr_64_64(uint64_t state); + diff --git a/c/pcg-rngs-64.c b/c/pcg-rngs-64.c new file mode 100644 --- /dev/null +++ b/c/pcg-rngs-64.c @@ -0,0 +1,232 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the canonical C++ PCG implementation, which + * has many additional features and is preferable if you can use C++ in + * your project. + * + * The contents of this file were mechanically derived from pcg_variants.h + * (every inline function defined there gets an exern declaration here). + */ + +#include "pcg_variants.h" + +/* Functions to advance the underlying LCG, one version for each size and + * each style. These functions are considered semi-private. There is rarely + * a good reason to call them directly. + */ + +extern inline void pcg_oneseq_64_step_r(struct pcg_state_64* rng); + +extern inline void pcg_oneseq_64_advance_r(struct pcg_state_64* rng, + uint64_t delta); + +extern inline void pcg_mcg_64_step_r(struct pcg_state_64* rng); + +extern inline void pcg_mcg_64_advance_r(struct pcg_state_64* rng, + uint64_t delta); + +extern inline void pcg_unique_64_step_r(struct pcg_state_64* rng); + +extern inline void pcg_unique_64_advance_r(struct pcg_state_64* rng, + uint64_t delta); + +extern inline void pcg_setseq_64_step_r(struct pcg_state_setseq_64* rng); + +extern inline void pcg_setseq_64_advance_r(struct pcg_state_setseq_64* rng, + uint64_t delta); + +/* Functions to seed the RNG state, one version for each size and each + * style. Unlike the step functions, regular users can and should call + * these functions. + */ + +extern inline void pcg_oneseq_64_srandom_r(struct pcg_state_64* rng, + uint64_t initstate); + +extern inline void pcg_mcg_64_srandom_r(struct pcg_state_64* rng, + uint64_t initstate); + +extern inline void pcg_unique_64_srandom_r(struct pcg_state_64* rng, + uint64_t initstate); + +extern inline void pcg_setseq_64_srandom_r(struct pcg_state_setseq_64* rng, + uint64_t initstate, + uint64_t initseq); + +/* Now, finally we create each of the individual generators. We provide + * a random_r function that provides a random number of the appropriate + * type (using the full range of the type) and a boundedrand_r version + * that provides + * + * Implementation notes for boundedrand_r: + * + * To avoid bias, we need to make the range of the RNG a multiple of + * bound, which we do by dropping output less than a threshold. + * Let's consider a 32-bit case... A naive scheme to calculate the + * threshold would be to do + * + * uint32_t threshold = 0x100000000ull % bound; + * + * but 64-bit div/mod is slower than 32-bit div/mod (especially on + * 32-bit platforms). In essence, we do + * + * uint32_t threshold = (0x100000000ull-bound) % bound; + * + * because this version will calculate the same modulus, but the LHS + * value is less than 2^32. + * + * (Note that using modulo is only wise for good RNGs, poorer RNGs + * such as raw LCGs do better using a technique based on division.) + * Empricical tests show that division is preferable to modulus for + * reducting the range of an RNG. It's faster, and sometimes it can + * even be statistically prefereable. + */ + +/* Generation functions for XSH RS */ + +extern inline uint32_t +pcg_oneseq_64_xsh_rs_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_oneseq_64_xsh_rs_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +extern inline uint32_t +pcg_unique_64_xsh_rs_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_unique_64_xsh_rs_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +extern inline uint32_t +pcg_setseq_64_xsh_rs_32_random_r(struct pcg_state_setseq_64* rng); + +extern inline uint32_t +pcg_setseq_64_xsh_rs_32_boundedrand_r(struct pcg_state_setseq_64* rng, + uint32_t bound); + +extern inline uint32_t pcg_mcg_64_xsh_rs_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_mcg_64_xsh_rs_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +/* Generation functions for XSH RR */ + +extern inline uint32_t +pcg_oneseq_64_xsh_rr_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_oneseq_64_xsh_rr_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +extern inline uint32_t +pcg_unique_64_xsh_rr_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_unique_64_xsh_rr_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +extern inline uint32_t +pcg_setseq_64_xsh_rr_32_random_r(struct pcg_state_setseq_64* rng); + +extern inline uint32_t +pcg_setseq_64_xsh_rr_32_boundedrand_r(struct pcg_state_setseq_64* rng, + uint32_t bound); + +extern inline uint32_t pcg_mcg_64_xsh_rr_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_mcg_64_xsh_rr_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +/* Generation functions for RXS M XS (no MCG versions because they + * don't make sense when you want to use the entire state) + */ + +extern inline uint64_t +pcg_oneseq_64_rxs_m_xs_64_random_r(struct pcg_state_64* rng); + +extern inline uint64_t +pcg_oneseq_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound); + +extern inline uint64_t +pcg_unique_64_rxs_m_xs_64_random_r(struct pcg_state_64* rng); + +extern inline uint64_t +pcg_unique_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound); + +extern inline uint64_t +pcg_setseq_64_rxs_m_xs_64_random_r(struct pcg_state_setseq_64* rng); + +extern inline uint64_t +pcg_setseq_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_setseq_64* rng, + uint64_t bound); + +/* Generation functions for XSL RR (only defined for "large" types) */ + +extern inline uint32_t +pcg_oneseq_64_xsl_rr_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_oneseq_64_xsl_rr_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +extern inline uint32_t +pcg_unique_64_xsl_rr_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_unique_64_xsl_rr_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +extern inline uint32_t +pcg_setseq_64_xsl_rr_32_random_r(struct pcg_state_setseq_64* rng); + +extern inline uint32_t +pcg_setseq_64_xsl_rr_32_boundedrand_r(struct pcg_state_setseq_64* rng, + uint32_t bound); + +extern inline uint32_t pcg_mcg_64_xsl_rr_32_random_r(struct pcg_state_64* rng); + +extern inline uint32_t +pcg_mcg_64_xsl_rr_32_boundedrand_r(struct pcg_state_64* rng, uint32_t bound); + +/* Generation functions for XSL RR RR (only defined for "large" types) */ + +extern inline uint64_t +pcg_oneseq_64_xsl_rr_rr_64_random_r(struct pcg_state_64* rng); + +extern inline uint64_t +pcg_oneseq_64_xsl_rr_rr_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound); + +extern inline uint64_t +pcg_unique_64_xsl_rr_rr_64_random_r(struct pcg_state_64* rng); + +extern inline uint64_t +pcg_unique_64_xsl_rr_rr_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound); + +extern inline uint64_t +pcg_setseq_64_xsl_rr_rr_64_random_r(struct pcg_state_setseq_64* rng); + +extern inline uint64_t +pcg_setseq_64_xsl_rr_rr_64_boundedrand_r(struct pcg_state_setseq_64* rng, + uint64_t bound); + diff --git a/c/pcg_variants.h b/c/pcg_variants.h new file mode 100644 --- /dev/null +++ b/c/pcg_variants.h @@ -0,0 +1,2211 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the canonical C++ PCG implementation, which + * has many additional features and is preferable if you can use C++ in + * your project. + * + * Much of the derivation was performed mechanically. In particular, the + * output functions were generated by compiling the C++ output functions + * into LLVM bitcode and then transforming that using the LLVM C backend + * (from https://github.com/draperlaboratory/llvm-cbe), and then + * postprocessing and hand editing the output. + * + * Much of the remaining code was generated by C-preprocessor metaprogramming. + */ + +#ifndef PCG_VARIANTS_H_INCLUDED +#define PCG_VARIANTS_H_INCLUDED 1 + +#include + +#if __SIZEOF_INT128__ + typedef __uint128_t pcg128_t; + #define PCG_128BIT_CONSTANT(high,low) \ + ((((pcg128_t)high) << 64) + low) + #define PCG_HAS_128BIT_OPS 1 +#endif + +#if __GNUC_GNU_INLINE__ && !defined(__cplusplus) + #error Nonstandard GNU inlining semanatics. Compile with -std=c99 or better. + // We could instead use macros PCG_INLINE and PCG_EXTERN_INLINE + // but better to just reject ancient C code. +#endif + +#if __cplusplus +extern "C" { +#endif + +/* + * Rotate helper functions. + */ + +inline uint8_t pcg_rotr_8(uint8_t value, unsigned int rot) +{ +/* Unfortunately, clang is kinda pathetic when it comes to properly + * recognizing idiomatic rotate code, so for clang we actually provide + * assembler directives (enabled with PCG_USE_INLINE_ASM). Boo, hiss. + */ +#if PCG_USE_INLINE_ASM && __clang__ && (__x86_64__ || __i386__) + asm ("rorb %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +#else + return (value >> rot) | (value << ((- rot) & 7)); +#endif +} + +inline uint16_t pcg_rotr_16(uint16_t value, unsigned int rot) +{ +#if PCG_USE_INLINE_ASM && __clang__ && (__x86_64__ || __i386__) + asm ("rorw %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +#else + return (value >> rot) | (value << ((- rot) & 15)); +#endif +} + +inline uint32_t pcg_rotr_32(uint32_t value, unsigned int rot) +{ +#if PCG_USE_INLINE_ASM && __clang__ && (__x86_64__ || __i386__) + asm ("rorl %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +#else + return (value >> rot) | (value << ((- rot) & 31)); +#endif +} + +inline uint64_t pcg_rotr_64(uint64_t value, unsigned int rot) +{ +#if 0 && PCG_USE_INLINE_ASM && __clang__ && __x86_64__ + // For whatever reason, clang actually *does* generator rotq by + // itself, so we don't need this code. + asm ("rorq %%cl, %0" : "=r" (value) : "0" (value), "c" (rot)); + return value; +#else + return (value >> rot) | (value << ((- rot) & 63)); +#endif +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_rotr_128(pcg128_t value, unsigned int rot) +{ + return (value >> rot) | (value << ((- rot) & 127)); +} +#endif + +/* + * Output functions. These are the core of the PCG generation scheme. + */ + +// XSH RS + +inline uint8_t pcg_output_xsh_rs_16_8(uint16_t state) +{ + return (uint8_t)(((state >> 7u) ^ state) >> ((state >> 14u) + 3u)); +} + +inline uint16_t pcg_output_xsh_rs_32_16(uint32_t state) +{ + return (uint16_t)(((state >> 11u) ^ state) >> ((state >> 30u) + 11u)); +} + +inline uint32_t pcg_output_xsh_rs_64_32(uint64_t state) +{ + + return (uint32_t)(((state >> 22u) ^ state) >> ((state >> 61u) + 22u)); +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_output_xsh_rs_128_64(pcg128_t state) +{ + return (uint64_t)(((state >> 43u) ^ state) >> ((state >> 124u) + 45u)); +} +#endif + +// XSH RR + +inline uint8_t pcg_output_xsh_rr_16_8(uint16_t state) +{ + return pcg_rotr_8(((state >> 5u) ^ state) >> 5u, state >> 13u); +} + +inline uint16_t pcg_output_xsh_rr_32_16(uint32_t state) +{ + return pcg_rotr_16(((state >> 10u) ^ state) >> 12u, state >> 28u); +} + +inline uint32_t pcg_output_xsh_rr_64_32(uint64_t state) +{ + return pcg_rotr_32(((state >> 18u) ^ state) >> 27u, state >> 59u); +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_output_xsh_rr_128_64(pcg128_t state) +{ + return pcg_rotr_64(((state >> 29u) ^ state) >> 58u, state >> 122u); +} +#endif + +// RXS M XS + +inline uint8_t pcg_output_rxs_m_xs_8_8(uint8_t state) +{ + uint8_t word = ((state >> ((state >> 6u) + 2u)) ^ state) * 217u; + return (word >> 6u) ^ word; +} + +inline uint16_t pcg_output_rxs_m_xs_16_16(uint16_t state) +{ + uint16_t word = ((state >> ((state >> 13u) + 3u)) ^ state) * 62169u; + return (word >> 11u) ^ word; +} + +inline uint32_t pcg_output_rxs_m_xs_32_32(uint32_t state) +{ + uint32_t word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u; + return (word >> 22u) ^ word; +} + +inline uint64_t pcg_output_rxs_m_xs_64_64(uint64_t state) +{ + uint64_t word = ((state >> ((state >> 59u) + 5u)) ^ state) + * 12605985483714917081ull; + return (word >> 43u) ^ word; +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_output_rxs_m_xs_128_128(pcg128_t state) +{ + pcg128_t word = ((state >> ((state >> 122u) + 6u)) ^ state) + * (PCG_128BIT_CONSTANT(17766728186571221404ULL, + 12605985483714917081ULL)); + // 327738287884841127335028083622016905945 + return (word >> 86u) ^ word; +} +#endif + +// XSL RR (only defined for >= 64 bits) + +inline uint32_t pcg_output_xsl_rr_64_32(uint64_t state) +{ + return pcg_rotr_32(((uint32_t)(state >> 32u)) ^ (uint32_t)state, + state >> 59u); +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_output_xsl_rr_128_64(pcg128_t state) +{ + return pcg_rotr_64(((uint64_t)(state >> 64u)) ^ (uint64_t)state, + state >> 122u); +} +#endif + +// XSL RR RR (only defined for >= 64 bits) + +inline uint64_t pcg_output_xsl_rr_rr_64_64(uint64_t state) +{ + uint32_t rot1 = (uint32_t)(state >> 59u); + uint32_t high = (uint32_t)(state >> 32u); + uint32_t low = (uint32_t)state; + uint32_t xored = high ^ low; + uint32_t newlow = pcg_rotr_32(xored, rot1); + uint32_t newhigh = pcg_rotr_32(high, newlow & 31u); + return (((uint64_t)newhigh) << 32u) | newlow; +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_output_xsl_rr_rr_128_128(pcg128_t state) +{ + uint32_t rot1 = (uint32_t)(state >> 122u); + uint64_t high = (uint64_t)(state >> 64u); + uint64_t low = (uint64_t)state; + uint64_t xored = high ^ low; + uint64_t newlow = pcg_rotr_64(xored, rot1); + uint64_t newhigh = pcg_rotr_64(high, newlow & 63u); + return (((pcg128_t)newhigh) << 64u) | newlow; +} +#endif + +#define PCG_DEFAULT_MULTIPLIER_8 141U +#define PCG_DEFAULT_MULTIPLIER_16 25385U +#define PCG_DEFAULT_MULTIPLIER_32 747796405U +#define PCG_DEFAULT_MULTIPLIER_64 6364136223846793005ULL + +#define PCG_DEFAULT_INCREMENT_8 77U +#define PCG_DEFAULT_INCREMENT_16 47989U +#define PCG_DEFAULT_INCREMENT_32 2891336453U +#define PCG_DEFAULT_INCREMENT_64 1442695040888963407ULL + +#if PCG_HAS_128BIT_OPS +#define PCG_DEFAULT_MULTIPLIER_128 \ + PCG_128BIT_CONSTANT(2549297995355413924ULL,4865540595714422341ULL) +#define PCG_DEFAULT_INCREMENT_128 \ + PCG_128BIT_CONSTANT(6364136223846793005ULL,1442695040888963407ULL) +#endif + +/* + * Static initialization constants (if you can't call srandom for some + * bizarre reason). + */ + +#if PCG_HAS_128BIT_OPS +#define PCG_STATE_ONESEQ_8_INITIALIZER { 0xd7U } +#define PCG_STATE_ONESEQ_16_INITIALIZER { 0x20dfU } +#define PCG_STATE_ONESEQ_32_INITIALIZER { 0x46b56677U } +#define PCG_STATE_ONESEQ_64_INITIALIZER { 0x4d595df4d0f33173ULL } +#define PCG_STATE_ONESEQ_128_INITIALIZER \ + { PCG_128BIT_CONSTANT(0xb8dc10e158a92392ULL, 0x98046df007ec0a53ULL) } +#endif + +#if PCG_HAS_128BIT_OPS +#define PCG_STATE_UNIQUE_8_INITIALIZER PCG_STATE_ONESEQ_8_INITIALIZER +#define PCG_STATE_UNIQUE_16_INITIALIZER PCG_STATE_ONESEQ_16_INITIALIZER +#define PCG_STATE_UNIQUE_32_INITIALIZER PCG_STATE_ONESEQ_32_INITIALIZER +#define PCG_STATE_UNIQUE_64_INITIALIZER PCG_STATE_ONESEQ_64_INITIALIZER +#define PCG_STATE_UNIQUE_128_INITIALIZER PCG_STATE_ONESEQ_128_INITIALIZER +#endif + +#if PCG_HAS_128BIT_OPS +#define PCG_STATE_MCG_8_INITIALIZER { 0xe5U } +#define PCG_STATE_MCG_16_INITIALIZER { 0xa5e5U } +#define PCG_STATE_MCG_32_INITIALIZER { 0xd15ea5e5U } +#define PCG_STATE_MCG_64_INITIALIZER { 0xcafef00dd15ea5e5ULL } +#define PCG_STATE_MCG_128_INITIALIZER \ + { PCG_128BIT_CONSTANT(0x0000000000000000ULL, 0xcafef00dd15ea5e5ULL) } +#endif + +#if PCG_HAS_128BIT_OPS +#define PCG_STATE_SETSEQ_8_INITIALIZER { 0x9bU, 0xdbU } +#define PCG_STATE_SETSEQ_16_INITIALIZER { 0xe39bU, 0x5bdbU } +#define PCG_STATE_SETSEQ_32_INITIALIZER { 0xec02d89bU, 0x94b95bdbU } +#define PCG_STATE_SETSEQ_64_INITIALIZER \ + { 0x853c49e6748fea9bULL, 0xda3e39cb94b95bdbULL } +#define PCG_STATE_SETSEQ_128_INITIALIZER \ + { PCG_128BIT_CONSTANT(0x979c9a98d8462005ULL, 0x7d3e9cb6cfe0549bULL), \ + PCG_128BIT_CONSTANT(0x0000000000000001ULL, 0xda3e39cb94b95bdbULL) } +#endif + +/* Representations for the oneseq, mcg, and unique variants */ + +struct pcg_state_8 { + uint8_t state; +}; + +struct pcg_state_16 { + uint16_t state; +}; + +struct pcg_state_32 { + uint32_t state; +}; + +struct pcg_state_64 { + uint64_t state; +}; + +#if PCG_HAS_128BIT_OPS +struct pcg_state_128 { + pcg128_t state; +}; +#endif + +/* Representations setseq variants */ + +struct pcg_state_setseq_8 { + uint8_t state; + uint8_t inc; +}; + +struct pcg_state_setseq_16 { + uint16_t state; + uint16_t inc; +}; + +struct pcg_state_setseq_32 { + uint32_t state; + uint32_t inc; +}; + +struct pcg_state_setseq_64 { + uint64_t state; + uint64_t inc; +}; + +#if PCG_HAS_128BIT_OPS +struct pcg_state_setseq_128 { + pcg128_t state; + pcg128_t inc; +}; +#endif + +/* Multi-step advance functions (jump-ahead, jump-back) */ + +extern uint8_t pcg_advance_lcg_8(uint8_t state, uint8_t delta, uint8_t cur_mult, + uint8_t cur_plus); +extern uint16_t pcg_advance_lcg_16(uint16_t state, uint16_t delta, + uint16_t cur_mult, uint16_t cur_plus); +extern uint32_t pcg_advance_lcg_32(uint32_t state, uint32_t delta, + uint32_t cur_mult, uint32_t cur_plus); +extern uint64_t pcg_advance_lcg_64(uint64_t state, uint64_t delta, + uint64_t cur_mult, uint64_t cur_plus); + +#if PCG_HAS_128BIT_OPS +extern pcg128_t pcg_advance_lcg_128(pcg128_t state, pcg128_t delta, + pcg128_t cur_mult, pcg128_t cur_plus); +#endif + +/* Functions to advance the underlying LCG, one version for each size and + * each style. These functions are considered semi-private. There is rarely + * a good reason to call them directly. + */ + +inline void pcg_oneseq_8_step_r(struct pcg_state_8* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_8 + + PCG_DEFAULT_INCREMENT_8; +} + +inline void pcg_oneseq_8_advance_r(struct pcg_state_8* rng, uint8_t delta) +{ + rng->state = pcg_advance_lcg_8(rng->state, delta, PCG_DEFAULT_MULTIPLIER_8, + PCG_DEFAULT_INCREMENT_8); +} + +inline void pcg_mcg_8_step_r(struct pcg_state_8* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_8; +} + +inline void pcg_mcg_8_advance_r(struct pcg_state_8* rng, uint8_t delta) +{ + rng->state + = pcg_advance_lcg_8(rng->state, delta, PCG_DEFAULT_MULTIPLIER_8, 0u); +} + +inline void pcg_unique_8_step_r(struct pcg_state_8* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_8 + + (uint8_t)(((intptr_t)rng) | 1u); +} + +inline void pcg_unique_8_advance_r(struct pcg_state_8* rng, uint8_t delta) +{ + rng->state = pcg_advance_lcg_8(rng->state, delta, PCG_DEFAULT_MULTIPLIER_8, + (uint8_t)(((intptr_t)rng) | 1u)); +} + +inline void pcg_setseq_8_step_r(struct pcg_state_setseq_8* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_8 + rng->inc; +} + +inline void pcg_setseq_8_advance_r(struct pcg_state_setseq_8* rng, + uint8_t delta) +{ + rng->state = pcg_advance_lcg_8(rng->state, delta, PCG_DEFAULT_MULTIPLIER_8, + rng->inc); +} + +inline void pcg_oneseq_16_step_r(struct pcg_state_16* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_16 + + PCG_DEFAULT_INCREMENT_16; +} + +inline void pcg_oneseq_16_advance_r(struct pcg_state_16* rng, uint16_t delta) +{ + rng->state = pcg_advance_lcg_16( + rng->state, delta, PCG_DEFAULT_MULTIPLIER_16, PCG_DEFAULT_INCREMENT_16); +} + +inline void pcg_mcg_16_step_r(struct pcg_state_16* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_16; +} + +inline void pcg_mcg_16_advance_r(struct pcg_state_16* rng, uint16_t delta) +{ + rng->state + = pcg_advance_lcg_16(rng->state, delta, PCG_DEFAULT_MULTIPLIER_16, 0u); +} + +inline void pcg_unique_16_step_r(struct pcg_state_16* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_16 + + (uint16_t)(((intptr_t)rng) | 1u); +} + +inline void pcg_unique_16_advance_r(struct pcg_state_16* rng, uint16_t delta) +{ + rng->state + = pcg_advance_lcg_16(rng->state, delta, PCG_DEFAULT_MULTIPLIER_16, + (uint16_t)(((intptr_t)rng) | 1u)); +} + +inline void pcg_setseq_16_step_r(struct pcg_state_setseq_16* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_16 + rng->inc; +} + +inline void pcg_setseq_16_advance_r(struct pcg_state_setseq_16* rng, + uint16_t delta) +{ + rng->state = pcg_advance_lcg_16(rng->state, delta, + PCG_DEFAULT_MULTIPLIER_16, rng->inc); +} + +inline void pcg_oneseq_32_step_r(struct pcg_state_32* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_32 + + PCG_DEFAULT_INCREMENT_32; +} + +inline void pcg_oneseq_32_advance_r(struct pcg_state_32* rng, uint32_t delta) +{ + rng->state = pcg_advance_lcg_32( + rng->state, delta, PCG_DEFAULT_MULTIPLIER_32, PCG_DEFAULT_INCREMENT_32); +} + +inline void pcg_mcg_32_step_r(struct pcg_state_32* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_32; +} + +inline void pcg_mcg_32_advance_r(struct pcg_state_32* rng, uint32_t delta) +{ + rng->state + = pcg_advance_lcg_32(rng->state, delta, PCG_DEFAULT_MULTIPLIER_32, 0u); +} + +inline void pcg_unique_32_step_r(struct pcg_state_32* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_32 + + (uint32_t)(((intptr_t)rng) | 1u); +} + +inline void pcg_unique_32_advance_r(struct pcg_state_32* rng, uint32_t delta) +{ + rng->state + = pcg_advance_lcg_32(rng->state, delta, PCG_DEFAULT_MULTIPLIER_32, + (uint32_t)(((intptr_t)rng) | 1u)); +} + +inline void pcg_setseq_32_step_r(struct pcg_state_setseq_32* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_32 + rng->inc; +} + +inline void pcg_setseq_32_advance_r(struct pcg_state_setseq_32* rng, + uint32_t delta) +{ + rng->state = pcg_advance_lcg_32(rng->state, delta, + PCG_DEFAULT_MULTIPLIER_32, rng->inc); +} + +inline void pcg_oneseq_64_step_r(struct pcg_state_64* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_64 + + PCG_DEFAULT_INCREMENT_64; +} + +inline void pcg_oneseq_64_advance_r(struct pcg_state_64* rng, uint64_t delta) +{ + rng->state = pcg_advance_lcg_64( + rng->state, delta, PCG_DEFAULT_MULTIPLIER_64, PCG_DEFAULT_INCREMENT_64); +} + +inline void pcg_mcg_64_step_r(struct pcg_state_64* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_64; +} + +inline void pcg_mcg_64_advance_r(struct pcg_state_64* rng, uint64_t delta) +{ + rng->state + = pcg_advance_lcg_64(rng->state, delta, PCG_DEFAULT_MULTIPLIER_64, 0u); +} + +inline void pcg_unique_64_step_r(struct pcg_state_64* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_64 + + (uint64_t)(((intptr_t)rng) | 1u); +} + +inline void pcg_unique_64_advance_r(struct pcg_state_64* rng, uint64_t delta) +{ + rng->state + = pcg_advance_lcg_64(rng->state, delta, PCG_DEFAULT_MULTIPLIER_64, + (uint64_t)(((intptr_t)rng) | 1u)); +} + +inline void pcg_setseq_64_step_r(struct pcg_state_setseq_64* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_64 + rng->inc; +} + +inline void pcg_setseq_64_advance_r(struct pcg_state_setseq_64* rng, + uint64_t delta) +{ + rng->state = pcg_advance_lcg_64(rng->state, delta, + PCG_DEFAULT_MULTIPLIER_64, rng->inc); +} + +#if PCG_HAS_128BIT_OPS +inline void pcg_oneseq_128_step_r(struct pcg_state_128* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_128 + + PCG_DEFAULT_INCREMENT_128; +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_oneseq_128_advance_r(struct pcg_state_128* rng, pcg128_t delta) +{ + rng->state + = pcg_advance_lcg_128(rng->state, delta, PCG_DEFAULT_MULTIPLIER_128, + PCG_DEFAULT_INCREMENT_128); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_mcg_128_step_r(struct pcg_state_128* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_128; +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_mcg_128_advance_r(struct pcg_state_128* rng, pcg128_t delta) +{ + rng->state = pcg_advance_lcg_128(rng->state, delta, + PCG_DEFAULT_MULTIPLIER_128, 0u); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_unique_128_step_r(struct pcg_state_128* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_128 + + (pcg128_t)(((intptr_t)rng) | 1u); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_unique_128_advance_r(struct pcg_state_128* rng, pcg128_t delta) +{ + rng->state + = pcg_advance_lcg_128(rng->state, delta, PCG_DEFAULT_MULTIPLIER_128, + (pcg128_t)(((intptr_t)rng) | 1u)); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_setseq_128_step_r(struct pcg_state_setseq_128* rng) +{ + rng->state = rng->state * PCG_DEFAULT_MULTIPLIER_128 + rng->inc; +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_setseq_128_advance_r(struct pcg_state_setseq_128* rng, + pcg128_t delta) +{ + rng->state = pcg_advance_lcg_128(rng->state, delta, + PCG_DEFAULT_MULTIPLIER_128, rng->inc); +} +#endif + +/* Functions to seed the RNG state, one version for each size and each + * style. Unlike the step functions, regular users can and should call + * these functions. + */ + +inline void pcg_oneseq_8_srandom_r(struct pcg_state_8* rng, uint8_t initstate) +{ + rng->state = 0U; + pcg_oneseq_8_step_r(rng); + rng->state += initstate; + pcg_oneseq_8_step_r(rng); +} + +inline void pcg_mcg_8_srandom_r(struct pcg_state_8* rng, uint8_t initstate) +{ + rng->state = initstate | 1u; +} + +inline void pcg_unique_8_srandom_r(struct pcg_state_8* rng, uint8_t initstate) +{ + rng->state = 0U; + pcg_unique_8_step_r(rng); + rng->state += initstate; + pcg_unique_8_step_r(rng); +} + +inline void pcg_setseq_8_srandom_r(struct pcg_state_setseq_8* rng, + uint8_t initstate, uint8_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg_setseq_8_step_r(rng); + rng->state += initstate; + pcg_setseq_8_step_r(rng); +} + +inline void pcg_oneseq_16_srandom_r(struct pcg_state_16* rng, + uint16_t initstate) +{ + rng->state = 0U; + pcg_oneseq_16_step_r(rng); + rng->state += initstate; + pcg_oneseq_16_step_r(rng); +} + +inline void pcg_mcg_16_srandom_r(struct pcg_state_16* rng, uint16_t initstate) +{ + rng->state = initstate | 1u; +} + +inline void pcg_unique_16_srandom_r(struct pcg_state_16* rng, + uint16_t initstate) +{ + rng->state = 0U; + pcg_unique_16_step_r(rng); + rng->state += initstate; + pcg_unique_16_step_r(rng); +} + +inline void pcg_setseq_16_srandom_r(struct pcg_state_setseq_16* rng, + uint16_t initstate, uint16_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg_setseq_16_step_r(rng); + rng->state += initstate; + pcg_setseq_16_step_r(rng); +} + +inline void pcg_oneseq_32_srandom_r(struct pcg_state_32* rng, + uint32_t initstate) +{ + rng->state = 0U; + pcg_oneseq_32_step_r(rng); + rng->state += initstate; + pcg_oneseq_32_step_r(rng); +} + +inline void pcg_mcg_32_srandom_r(struct pcg_state_32* rng, uint32_t initstate) +{ + rng->state = initstate | 1u; +} + +inline void pcg_unique_32_srandom_r(struct pcg_state_32* rng, + uint32_t initstate) +{ + rng->state = 0U; + pcg_unique_32_step_r(rng); + rng->state += initstate; + pcg_unique_32_step_r(rng); +} + +inline void pcg_setseq_32_srandom_r(struct pcg_state_setseq_32* rng, + uint32_t initstate, uint32_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg_setseq_32_step_r(rng); + rng->state += initstate; + pcg_setseq_32_step_r(rng); +} + +inline void pcg_oneseq_64_srandom_r(struct pcg_state_64* rng, + uint64_t initstate) +{ + rng->state = 0U; + pcg_oneseq_64_step_r(rng); + rng->state += initstate; + pcg_oneseq_64_step_r(rng); +} + +inline void pcg_mcg_64_srandom_r(struct pcg_state_64* rng, uint64_t initstate) +{ + rng->state = initstate | 1u; +} + +inline void pcg_unique_64_srandom_r(struct pcg_state_64* rng, + uint64_t initstate) +{ + rng->state = 0U; + pcg_unique_64_step_r(rng); + rng->state += initstate; + pcg_unique_64_step_r(rng); +} + +inline void pcg_setseq_64_srandom_r(struct pcg_state_setseq_64* rng, + uint64_t initstate, uint64_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg_setseq_64_step_r(rng); + rng->state += initstate; + pcg_setseq_64_step_r(rng); +} + +#if PCG_HAS_128BIT_OPS +inline void pcg_oneseq_128_srandom_r(struct pcg_state_128* rng, + pcg128_t initstate) +{ + rng->state = 0U; + pcg_oneseq_128_step_r(rng); + rng->state += initstate; + pcg_oneseq_128_step_r(rng); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_mcg_128_srandom_r(struct pcg_state_128* rng, pcg128_t initstate) +{ + rng->state = initstate | 1u; +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_unique_128_srandom_r(struct pcg_state_128* rng, + pcg128_t initstate) +{ + rng->state = 0U; + pcg_unique_128_step_r(rng); + rng->state += initstate; + pcg_unique_128_step_r(rng); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline void pcg_setseq_128_srandom_r(struct pcg_state_setseq_128* rng, + pcg128_t initstate, pcg128_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg_setseq_128_step_r(rng); + rng->state += initstate; + pcg_setseq_128_step_r(rng); +} +#endif + +/* Now, finally we create each of the individual generators. We provide + * a random_r function that provides a random number of the appropriate + * type (using the full range of the type) and a boundedrand_r version + * that provides + * + * Implementation notes for boundedrand_r: + * + * To avoid bias, we need to make the range of the RNG a multiple of + * bound, which we do by dropping output less than a threshold. + * Let's consider a 32-bit case... A naive scheme to calculate the + * threshold would be to do + * + * uint32_t threshold = 0x100000000ull % bound; + * + * but 64-bit div/mod is slower than 32-bit div/mod (especially on + * 32-bit platforms). In essence, we do + * + * uint32_t threshold = (0x100000000ull-bound) % bound; + * + * because this version will calculate the same modulus, but the LHS + * value is less than 2^32. + * + * (Note that using modulo is only wise for good RNGs, poorer RNGs + * such as raw LCGs do better using a technique based on division.) + * Empricical tests show that division is preferable to modulus for + * reducting the range of an RNG. It's faster, and sometimes it can + * even be statistically prefereable. + */ + +/* Generation functions for XSH RS */ + +inline uint8_t pcg_oneseq_16_xsh_rs_8_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_oneseq_16_step_r(rng); + return pcg_output_xsh_rs_16_8(oldstate); +} + +inline uint8_t pcg_oneseq_16_xsh_rs_8_boundedrand_r(struct pcg_state_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_oneseq_16_xsh_rs_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_oneseq_32_xsh_rs_16_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_oneseq_32_step_r(rng); + return pcg_output_xsh_rs_32_16(oldstate); +} + +inline uint16_t pcg_oneseq_32_xsh_rs_16_boundedrand_r(struct pcg_state_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_oneseq_32_xsh_rs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_oneseq_64_xsh_rs_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_oneseq_64_step_r(rng); + return pcg_output_xsh_rs_64_32(oldstate); +} + +inline uint32_t pcg_oneseq_64_xsh_rs_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_oneseq_64_xsh_rs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_oneseq_128_xsh_rs_64_random_r(struct pcg_state_128* rng) +{ + pcg_oneseq_128_step_r(rng); + return pcg_output_xsh_rs_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_oneseq_128_xsh_rs_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_oneseq_128_xsh_rs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_unique_16_xsh_rs_8_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_unique_16_step_r(rng); + return pcg_output_xsh_rs_16_8(oldstate); +} + +inline uint8_t pcg_unique_16_xsh_rs_8_boundedrand_r(struct pcg_state_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_unique_16_xsh_rs_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_unique_32_xsh_rs_16_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_unique_32_step_r(rng); + return pcg_output_xsh_rs_32_16(oldstate); +} + +inline uint16_t pcg_unique_32_xsh_rs_16_boundedrand_r(struct pcg_state_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_unique_32_xsh_rs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_unique_64_xsh_rs_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_unique_64_step_r(rng); + return pcg_output_xsh_rs_64_32(oldstate); +} + +inline uint32_t pcg_unique_64_xsh_rs_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_unique_64_xsh_rs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_unique_128_xsh_rs_64_random_r(struct pcg_state_128* rng) +{ + pcg_unique_128_step_r(rng); + return pcg_output_xsh_rs_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_unique_128_xsh_rs_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_unique_128_xsh_rs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_setseq_16_xsh_rs_8_random_r(struct pcg_state_setseq_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_setseq_16_step_r(rng); + return pcg_output_xsh_rs_16_8(oldstate); +} + +inline uint8_t +pcg_setseq_16_xsh_rs_8_boundedrand_r(struct pcg_state_setseq_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_setseq_16_xsh_rs_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t +pcg_setseq_32_xsh_rs_16_random_r(struct pcg_state_setseq_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_setseq_32_step_r(rng); + return pcg_output_xsh_rs_32_16(oldstate); +} + +inline uint16_t +pcg_setseq_32_xsh_rs_16_boundedrand_r(struct pcg_state_setseq_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_setseq_32_xsh_rs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t +pcg_setseq_64_xsh_rs_32_random_r(struct pcg_state_setseq_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_setseq_64_step_r(rng); + return pcg_output_xsh_rs_64_32(oldstate); +} + +inline uint32_t +pcg_setseq_64_xsh_rs_32_boundedrand_r(struct pcg_state_setseq_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_setseq_64_xsh_rs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_setseq_128_xsh_rs_64_random_r(struct pcg_state_setseq_128* rng) +{ + pcg_setseq_128_step_r(rng); + return pcg_output_xsh_rs_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_setseq_128_xsh_rs_64_boundedrand_r(struct pcg_state_setseq_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_setseq_128_xsh_rs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_mcg_16_xsh_rs_8_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_mcg_16_step_r(rng); + return pcg_output_xsh_rs_16_8(oldstate); +} + +inline uint8_t pcg_mcg_16_xsh_rs_8_boundedrand_r(struct pcg_state_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_mcg_16_xsh_rs_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_mcg_32_xsh_rs_16_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_mcg_32_step_r(rng); + return pcg_output_xsh_rs_32_16(oldstate); +} + +inline uint16_t pcg_mcg_32_xsh_rs_16_boundedrand_r(struct pcg_state_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_mcg_32_xsh_rs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_mcg_64_xsh_rs_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_mcg_64_step_r(rng); + return pcg_output_xsh_rs_64_32(oldstate); +} + +inline uint32_t pcg_mcg_64_xsh_rs_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_mcg_64_xsh_rs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_mcg_128_xsh_rs_64_random_r(struct pcg_state_128* rng) +{ + pcg_mcg_128_step_r(rng); + return pcg_output_xsh_rs_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_mcg_128_xsh_rs_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_mcg_128_xsh_rs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +/* Generation functions for XSH RR */ + +inline uint8_t pcg_oneseq_16_xsh_rr_8_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_oneseq_16_step_r(rng); + return pcg_output_xsh_rr_16_8(oldstate); +} + +inline uint8_t pcg_oneseq_16_xsh_rr_8_boundedrand_r(struct pcg_state_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_oneseq_16_xsh_rr_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_oneseq_32_xsh_rr_16_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_oneseq_32_step_r(rng); + return pcg_output_xsh_rr_32_16(oldstate); +} + +inline uint16_t pcg_oneseq_32_xsh_rr_16_boundedrand_r(struct pcg_state_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_oneseq_32_xsh_rr_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_oneseq_64_xsh_rr_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_oneseq_64_step_r(rng); + return pcg_output_xsh_rr_64_32(oldstate); +} + +inline uint32_t pcg_oneseq_64_xsh_rr_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_oneseq_64_xsh_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_oneseq_128_xsh_rr_64_random_r(struct pcg_state_128* rng) +{ + pcg_oneseq_128_step_r(rng); + return pcg_output_xsh_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_oneseq_128_xsh_rr_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_oneseq_128_xsh_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_unique_16_xsh_rr_8_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_unique_16_step_r(rng); + return pcg_output_xsh_rr_16_8(oldstate); +} + +inline uint8_t pcg_unique_16_xsh_rr_8_boundedrand_r(struct pcg_state_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_unique_16_xsh_rr_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_unique_32_xsh_rr_16_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_unique_32_step_r(rng); + return pcg_output_xsh_rr_32_16(oldstate); +} + +inline uint16_t pcg_unique_32_xsh_rr_16_boundedrand_r(struct pcg_state_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_unique_32_xsh_rr_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_unique_64_xsh_rr_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_unique_64_step_r(rng); + return pcg_output_xsh_rr_64_32(oldstate); +} + +inline uint32_t pcg_unique_64_xsh_rr_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_unique_64_xsh_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_unique_128_xsh_rr_64_random_r(struct pcg_state_128* rng) +{ + pcg_unique_128_step_r(rng); + return pcg_output_xsh_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_unique_128_xsh_rr_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_unique_128_xsh_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_setseq_16_xsh_rr_8_random_r(struct pcg_state_setseq_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_setseq_16_step_r(rng); + return pcg_output_xsh_rr_16_8(oldstate); +} + +inline uint8_t +pcg_setseq_16_xsh_rr_8_boundedrand_r(struct pcg_state_setseq_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_setseq_16_xsh_rr_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t +pcg_setseq_32_xsh_rr_16_random_r(struct pcg_state_setseq_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_setseq_32_step_r(rng); + return pcg_output_xsh_rr_32_16(oldstate); +} + +inline uint16_t +pcg_setseq_32_xsh_rr_16_boundedrand_r(struct pcg_state_setseq_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_setseq_32_xsh_rr_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t +pcg_setseq_64_xsh_rr_32_random_r(struct pcg_state_setseq_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_setseq_64_step_r(rng); + return pcg_output_xsh_rr_64_32(oldstate); +} + +inline uint32_t +pcg_setseq_64_xsh_rr_32_boundedrand_r(struct pcg_state_setseq_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_setseq_64_xsh_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_setseq_128_xsh_rr_64_random_r(struct pcg_state_setseq_128* rng) +{ + pcg_setseq_128_step_r(rng); + return pcg_output_xsh_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_setseq_128_xsh_rr_64_boundedrand_r(struct pcg_state_setseq_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_setseq_128_xsh_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_mcg_16_xsh_rr_8_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_mcg_16_step_r(rng); + return pcg_output_xsh_rr_16_8(oldstate); +} + +inline uint8_t pcg_mcg_16_xsh_rr_8_boundedrand_r(struct pcg_state_16* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_mcg_16_xsh_rr_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_mcg_32_xsh_rr_16_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_mcg_32_step_r(rng); + return pcg_output_xsh_rr_32_16(oldstate); +} + +inline uint16_t pcg_mcg_32_xsh_rr_16_boundedrand_r(struct pcg_state_32* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_mcg_32_xsh_rr_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_mcg_64_xsh_rr_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_mcg_64_step_r(rng); + return pcg_output_xsh_rr_64_32(oldstate); +} + +inline uint32_t pcg_mcg_64_xsh_rr_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_mcg_64_xsh_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_mcg_128_xsh_rr_64_random_r(struct pcg_state_128* rng) +{ + pcg_mcg_128_step_r(rng); + return pcg_output_xsh_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_mcg_128_xsh_rr_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_mcg_128_xsh_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +/* Generation functions for RXS M XS (no MCG versions because they + * don't make sense when you want to use the entire state) + */ + +inline uint8_t pcg_oneseq_8_rxs_m_xs_8_random_r(struct pcg_state_8* rng) +{ + uint8_t oldstate = rng->state; + pcg_oneseq_8_step_r(rng); + return pcg_output_rxs_m_xs_8_8(oldstate); +} + +inline uint8_t pcg_oneseq_8_rxs_m_xs_8_boundedrand_r(struct pcg_state_8* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_oneseq_8_rxs_m_xs_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t pcg_oneseq_16_rxs_m_xs_16_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_oneseq_16_step_r(rng); + return pcg_output_rxs_m_xs_16_16(oldstate); +} + +inline uint16_t +pcg_oneseq_16_rxs_m_xs_16_boundedrand_r(struct pcg_state_16* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_oneseq_16_rxs_m_xs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_oneseq_32_rxs_m_xs_32_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_oneseq_32_step_r(rng); + return pcg_output_rxs_m_xs_32_32(oldstate); +} + +inline uint32_t +pcg_oneseq_32_rxs_m_xs_32_boundedrand_r(struct pcg_state_32* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_oneseq_32_rxs_m_xs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint64_t pcg_oneseq_64_rxs_m_xs_64_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_oneseq_64_step_r(rng); + return pcg_output_rxs_m_xs_64_64(oldstate); +} + +inline uint64_t +pcg_oneseq_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_oneseq_64_rxs_m_xs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_oneseq_128_rxs_m_xs_128_random_r(struct pcg_state_128* rng) +{ + pcg_oneseq_128_step_r(rng); + return pcg_output_rxs_m_xs_128_128(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_oneseq_128_rxs_m_xs_128_boundedrand_r(struct pcg_state_128* rng, + pcg128_t bound) +{ + pcg128_t threshold = -bound % bound; + for (;;) { + pcg128_t r = pcg_oneseq_128_rxs_m_xs_128_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint16_t pcg_unique_16_rxs_m_xs_16_random_r(struct pcg_state_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_unique_16_step_r(rng); + return pcg_output_rxs_m_xs_16_16(oldstate); +} + +inline uint16_t +pcg_unique_16_rxs_m_xs_16_boundedrand_r(struct pcg_state_16* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_unique_16_rxs_m_xs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t pcg_unique_32_rxs_m_xs_32_random_r(struct pcg_state_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_unique_32_step_r(rng); + return pcg_output_rxs_m_xs_32_32(oldstate); +} + +inline uint32_t +pcg_unique_32_rxs_m_xs_32_boundedrand_r(struct pcg_state_32* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_unique_32_rxs_m_xs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint64_t pcg_unique_64_rxs_m_xs_64_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_unique_64_step_r(rng); + return pcg_output_rxs_m_xs_64_64(oldstate); +} + +inline uint64_t +pcg_unique_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_unique_64_rxs_m_xs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_unique_128_rxs_m_xs_128_random_r(struct pcg_state_128* rng) +{ + pcg_unique_128_step_r(rng); + return pcg_output_rxs_m_xs_128_128(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_unique_128_rxs_m_xs_128_boundedrand_r(struct pcg_state_128* rng, + pcg128_t bound) +{ + pcg128_t threshold = -bound % bound; + for (;;) { + pcg128_t r = pcg_unique_128_rxs_m_xs_128_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint8_t pcg_setseq_8_rxs_m_xs_8_random_r(struct pcg_state_setseq_8* rng) +{ + uint8_t oldstate = rng->state; + pcg_setseq_8_step_r(rng); + return pcg_output_rxs_m_xs_8_8(oldstate); +} + +inline uint8_t +pcg_setseq_8_rxs_m_xs_8_boundedrand_r(struct pcg_state_setseq_8* rng, + uint8_t bound) +{ + uint8_t threshold = ((uint8_t)(-bound)) % bound; + for (;;) { + uint8_t r = pcg_setseq_8_rxs_m_xs_8_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint16_t +pcg_setseq_16_rxs_m_xs_16_random_r(struct pcg_state_setseq_16* rng) +{ + uint16_t oldstate = rng->state; + pcg_setseq_16_step_r(rng); + return pcg_output_rxs_m_xs_16_16(oldstate); +} + +inline uint16_t +pcg_setseq_16_rxs_m_xs_16_boundedrand_r(struct pcg_state_setseq_16* rng, + uint16_t bound) +{ + uint16_t threshold = ((uint16_t)(-bound)) % bound; + for (;;) { + uint16_t r = pcg_setseq_16_rxs_m_xs_16_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint32_t +pcg_setseq_32_rxs_m_xs_32_random_r(struct pcg_state_setseq_32* rng) +{ + uint32_t oldstate = rng->state; + pcg_setseq_32_step_r(rng); + return pcg_output_rxs_m_xs_32_32(oldstate); +} + +inline uint32_t +pcg_setseq_32_rxs_m_xs_32_boundedrand_r(struct pcg_state_setseq_32* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_setseq_32_rxs_m_xs_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +inline uint64_t +pcg_setseq_64_rxs_m_xs_64_random_r(struct pcg_state_setseq_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_setseq_64_step_r(rng); + return pcg_output_rxs_m_xs_64_64(oldstate); +} + +inline uint64_t +pcg_setseq_64_rxs_m_xs_64_boundedrand_r(struct pcg_state_setseq_64* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_setseq_64_rxs_m_xs_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_setseq_128_rxs_m_xs_128_random_r(struct pcg_state_setseq_128* rng) +{ + pcg_setseq_128_step_r(rng); + return pcg_output_rxs_m_xs_128_128(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_setseq_128_rxs_m_xs_128_boundedrand_r(struct pcg_state_setseq_128* rng, + pcg128_t bound) +{ + pcg128_t threshold = -bound % bound; + for (;;) { + pcg128_t r = pcg_setseq_128_rxs_m_xs_128_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +/* Generation functions for XSL RR (only defined for "large" types) */ + +inline uint32_t pcg_oneseq_64_xsl_rr_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_oneseq_64_step_r(rng); + return pcg_output_xsl_rr_64_32(oldstate); +} + +inline uint32_t pcg_oneseq_64_xsl_rr_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_oneseq_64_xsl_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_oneseq_128_xsl_rr_64_random_r(struct pcg_state_128* rng) +{ + pcg_oneseq_128_step_r(rng); + return pcg_output_xsl_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_oneseq_128_xsl_rr_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_oneseq_128_xsl_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint32_t pcg_unique_64_xsl_rr_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_unique_64_step_r(rng); + return pcg_output_xsl_rr_64_32(oldstate); +} + +inline uint32_t pcg_unique_64_xsl_rr_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_unique_64_xsl_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_unique_128_xsl_rr_64_random_r(struct pcg_state_128* rng) +{ + pcg_unique_128_step_r(rng); + return pcg_output_xsl_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_unique_128_xsl_rr_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_unique_128_xsl_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint32_t +pcg_setseq_64_xsl_rr_32_random_r(struct pcg_state_setseq_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_setseq_64_step_r(rng); + return pcg_output_xsl_rr_64_32(oldstate); +} + +inline uint32_t +pcg_setseq_64_xsl_rr_32_boundedrand_r(struct pcg_state_setseq_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_setseq_64_xsl_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_setseq_128_xsl_rr_64_random_r(struct pcg_state_setseq_128* rng) +{ + pcg_setseq_128_step_r(rng); + return pcg_output_xsl_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t +pcg_setseq_128_xsl_rr_64_boundedrand_r(struct pcg_state_setseq_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_setseq_128_xsl_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint32_t pcg_mcg_64_xsl_rr_32_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_mcg_64_step_r(rng); + return pcg_output_xsl_rr_64_32(oldstate); +} + +inline uint32_t pcg_mcg_64_xsl_rr_32_boundedrand_r(struct pcg_state_64* rng, + uint32_t bound) +{ + uint32_t threshold = -bound % bound; + for (;;) { + uint32_t r = pcg_mcg_64_xsl_rr_32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_mcg_128_xsl_rr_64_random_r(struct pcg_state_128* rng) +{ + pcg_mcg_128_step_r(rng); + return pcg_output_xsl_rr_128_64(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline uint64_t pcg_mcg_128_xsl_rr_64_boundedrand_r(struct pcg_state_128* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_mcg_128_xsl_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +/* Generation functions for XSL RR RR (only defined for "large" types) */ + +inline uint64_t pcg_oneseq_64_xsl_rr_rr_64_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_oneseq_64_step_r(rng); + return pcg_output_xsl_rr_rr_64_64(oldstate); +} + +inline uint64_t +pcg_oneseq_64_xsl_rr_rr_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_oneseq_64_xsl_rr_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_oneseq_128_xsl_rr_rr_128_random_r(struct pcg_state_128* rng) +{ + pcg_oneseq_128_step_r(rng); + return pcg_output_xsl_rr_rr_128_128(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_oneseq_128_xsl_rr_rr_128_boundedrand_r(struct pcg_state_128* rng, + pcg128_t bound) +{ + pcg128_t threshold = -bound % bound; + for (;;) { + pcg128_t r = pcg_oneseq_128_xsl_rr_rr_128_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint64_t pcg_unique_64_xsl_rr_rr_64_random_r(struct pcg_state_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_unique_64_step_r(rng); + return pcg_output_xsl_rr_rr_64_64(oldstate); +} + +inline uint64_t +pcg_unique_64_xsl_rr_rr_64_boundedrand_r(struct pcg_state_64* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_unique_64_xsl_rr_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t pcg_unique_128_xsl_rr_rr_128_random_r(struct pcg_state_128* rng) +{ + pcg_unique_128_step_r(rng); + return pcg_output_xsl_rr_rr_128_128(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_unique_128_xsl_rr_rr_128_boundedrand_r(struct pcg_state_128* rng, + pcg128_t bound) +{ + pcg128_t threshold = -bound % bound; + for (;;) { + pcg128_t r = pcg_unique_128_xsl_rr_rr_128_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +inline uint64_t +pcg_setseq_64_xsl_rr_rr_64_random_r(struct pcg_state_setseq_64* rng) +{ + uint64_t oldstate = rng->state; + pcg_setseq_64_step_r(rng); + return pcg_output_xsl_rr_rr_64_64(oldstate); +} + +inline uint64_t +pcg_setseq_64_xsl_rr_rr_64_boundedrand_r(struct pcg_state_setseq_64* rng, + uint64_t bound) +{ + uint64_t threshold = -bound % bound; + for (;;) { + uint64_t r = pcg_setseq_64_xsl_rr_rr_64_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_setseq_128_xsl_rr_rr_128_random_r(struct pcg_state_setseq_128* rng) +{ + pcg_setseq_128_step_r(rng); + return pcg_output_xsl_rr_rr_128_128(rng->state); +} +#endif + +#if PCG_HAS_128BIT_OPS +inline pcg128_t +pcg_setseq_128_xsl_rr_rr_128_boundedrand_r(struct pcg_state_setseq_128* rng, + pcg128_t bound) +{ + pcg128_t threshold = -bound % bound; + for (;;) { + pcg128_t r = pcg_setseq_128_xsl_rr_rr_128_random_r(rng); + if (r >= threshold) + return r % bound; + } +} +#endif + +//// Typedefs +typedef struct pcg_state_setseq_64 pcg32_random_t; +typedef struct pcg_state_64 pcg32s_random_t; +typedef struct pcg_state_64 pcg32u_random_t; +typedef struct pcg_state_64 pcg32f_random_t; +//// random_r +#define pcg32_random_r pcg_setseq_64_xsh_rr_32_random_r +#define pcg32s_random_r pcg_oneseq_64_xsh_rr_32_random_r +#define pcg32u_random_r pcg_unique_64_xsh_rr_32_random_r +#define pcg32f_random_r pcg_mcg_64_xsh_rs_32_random_r +//// boundedrand_r +#define pcg32_boundedrand_r pcg_setseq_64_xsh_rr_32_boundedrand_r +#define pcg32s_boundedrand_r pcg_oneseq_64_xsh_rr_32_boundedrand_r +#define pcg32u_boundedrand_r pcg_unique_64_xsh_rr_32_boundedrand_r +#define pcg32f_boundedrand_r pcg_mcg_64_xsh_rs_32_boundedrand_r +//// srandom_r +#define pcg32_srandom_r pcg_setseq_64_srandom_r +#define pcg32s_srandom_r pcg_oneseq_64_srandom_r +#define pcg32u_srandom_r pcg_unique_64_srandom_r +#define pcg32f_srandom_r pcg_mcg_64_srandom_r +//// advance_r +#define pcg32_advance_r pcg_setseq_64_advance_r +#define pcg32s_advance_r pcg_oneseq_64_advance_r +#define pcg32u_advance_r pcg_unique_64_advance_r +#define pcg32f_advance_r pcg_mcg_64_advance_r + +#if PCG_HAS_128BIT_OPS +//// Typedefs +typedef struct pcg_state_setseq_128 pcg64_random_t; +typedef struct pcg_state_128 pcg64s_random_t; +typedef struct pcg_state_128 pcg64u_random_t; +typedef struct pcg_state_128 pcg64f_random_t; +//// random_r +#define pcg64_random_r pcg_setseq_128_xsl_rr_64_random_r +#define pcg64s_random_r pcg_oneseq_128_xsl_rr_64_random_r +#define pcg64u_random_r pcg_unique_128_xsl_rr_64_random_r +#define pcg64f_random_r pcg_mcg_128_xsl_rr_64_random_r +//// boundedrand_r +#define pcg64_boundedrand_r pcg_setseq_128_xsl_rr_64_boundedrand_r +#define pcg64s_boundedrand_r pcg_oneseq_128_xsl_rr_64_boundedrand_r +#define pcg64u_boundedrand_r pcg_unique_128_xsl_rr_64_boundedrand_r +#define pcg64f_boundedrand_r pcg_mcg_128_xsl_rr_64_boundedrand_r +//// srandom_r +#define pcg64_srandom_r pcg_setseq_128_srandom_r +#define pcg64s_srandom_r pcg_oneseq_128_srandom_r +#define pcg64u_srandom_r pcg_unique_128_srandom_r +#define pcg64f_srandom_r pcg_mcg_128_srandom_r +//// advance_r +#define pcg64_advance_r pcg_setseq_128_advance_r +#define pcg64s_advance_r pcg_oneseq_128_advance_r +#define pcg64u_advance_r pcg_unique_128_advance_r +#define pcg64f_advance_r pcg_mcg_128_advance_r +#endif + +//// Typedefs +typedef struct pcg_state_8 pcg8si_random_t; +typedef struct pcg_state_16 pcg16si_random_t; +typedef struct pcg_state_32 pcg32si_random_t; +typedef struct pcg_state_64 pcg64si_random_t; +//// random_r +#define pcg8si_random_r pcg_oneseq_8_rxs_m_xs_8_random_r +#define pcg16si_random_r pcg_oneseq_16_rxs_m_xs_16_random_r +#define pcg32si_random_r pcg_oneseq_32_rxs_m_xs_32_random_r +#define pcg64si_random_r pcg_oneseq_64_rxs_m_xs_64_random_r +//// boundedrand_r +#define pcg8si_boundedrand_r pcg_oneseq_8_rxs_m_xs_8_boundedrand_r +#define pcg16si_boundedrand_r pcg_oneseq_16_rxs_m_xs_16_boundedrand_r +#define pcg32si_boundedrand_r pcg_oneseq_32_rxs_m_xs_32_boundedrand_r +#define pcg64si_boundedrand_r pcg_oneseq_64_rxs_m_xs_64_boundedrand_r +//// srandom_r +#define pcg8si_srandom_r pcg_oneseq_8_srandom_r +#define pcg16si_srandom_r pcg_oneseq_16_srandom_r +#define pcg32si_srandom_r pcg_oneseq_32_srandom_r +#define pcg64si_srandom_r pcg_oneseq_64_srandom_r +//// advance_r +#define pcg8si_advance_r pcg_oneseq_8_advance_r +#define pcg16si_advance_r pcg_oneseq_16_advance_r +#define pcg32si_advance_r pcg_oneseq_32_advance_r +#define pcg64si_advance_r pcg_oneseq_64_advance_r + +#if PCG_HAS_128BIT_OPS +typedef struct pcg_state_128 pcg128si_random_t; +#define pcg128si_random_r pcg_oneseq_128_rxs_m_xs_128_random_r +#define pcg128si_boundedrand_r pcg_oneseq_128_rxs_m_xs_128_boundedrand_r +#define pcg128si_srandom_r pcg_oneseq_128_srandom_r +#define pcg128si_advance_r pcg_oneseq_128_advance_r +#endif + +//// Typedefs +typedef struct pcg_state_setseq_8 pcg8i_random_t; +typedef struct pcg_state_setseq_16 pcg16i_random_t; +typedef struct pcg_state_setseq_32 pcg32i_random_t; +typedef struct pcg_state_setseq_64 pcg64i_random_t; +//// random_r +#define pcg8i_random_r pcg_setseq_8_rxs_m_xs_8_random_r +#define pcg16i_random_r pcg_setseq_16_rxs_m_xs_16_random_r +#define pcg32i_random_r pcg_setseq_32_rxs_m_xs_32_random_r +#define pcg64i_random_r pcg_setseq_64_rxs_m_xs_64_random_r +//// boundedrand_r +#define pcg8i_boundedrand_r pcg_setseq_8_rxs_m_xs_8_boundedrand_r +#define pcg16i_boundedrand_r pcg_setseq_16_rxs_m_xs_16_boundedrand_r +#define pcg32i_boundedrand_r pcg_setseq_32_rxs_m_xs_32_boundedrand_r +#define pcg64i_boundedrand_r pcg_setseq_64_rxs_m_xs_64_boundedrand_r +//// srandom_r +#define pcg8i_srandom_r pcg_setseq_8_srandom_r +#define pcg16i_srandom_r pcg_setseq_16_srandom_r +#define pcg32i_srandom_r pcg_setseq_32_srandom_r +#define pcg64i_srandom_r pcg_setseq_64_srandom_r +//// advance_r +#define pcg8i_advance_r pcg_setseq_8_advance_r +#define pcg16i_advance_r pcg_setseq_16_advance_r +#define pcg32i_advance_r pcg_setseq_32_advance_r +#define pcg64i_advance_r pcg_setseq_64_advance_r + +#if PCG_HAS_128BIT_OPS +typedef struct pcg_state_setseq_128 pcg128i_random_t; +#define pcg128i_random_r pcg_setseq_128_rxs_m_xs_128_random_r +#define pcg128i_boundedrand_r pcg_setseq_128_rxs_m_xs_128_boundedrand_r +#define pcg128i_srandom_r pcg_setseq_128_srandom_r +#define pcg128i_advance_r pcg_setseq_128_advance_r +#endif + +extern uint32_t pcg32_random(); +extern uint32_t pcg32_boundedrand(uint32_t bound); +extern void pcg32_srandom(uint64_t seed, uint64_t seq); +extern void pcg32_advance(uint64_t delta); + +#if PCG_HAS_128BIT_OPS +extern uint64_t pcg64_random(); +extern uint64_t pcg64_boundedrand(uint64_t bound); +extern void pcg64_srandom(pcg128_t seed, pcg128_t seq); +extern void pcg64_advance(pcg128_t delta); +#endif + +/* + * Static initialization constants (if you can't call srandom for some + * bizarre reason). + */ + +#define PCG32_INITIALIZER PCG_STATE_SETSEQ_64_INITIALIZER +#define PCG32U_INITIALIZER PCG_STATE_UNIQUE_64_INITIALIZER +#define PCG32S_INITIALIZER PCG_STATE_ONESEQ_64_INITIALIZER +#define PCG32F_INITIALIZER PCG_STATE_MCG_64_INITIALIZER + +#if PCG_HAS_128BIT_OPS +#define PCG64_INITIALIZER PCG_STATE_SETSEQ_128_INITIALIZER +#define PCG64U_INITIALIZER PCG_STATE_UNIQUE_128_INITIALIZER +#define PCG64S_INITIALIZER PCG_STATE_ONESEQ_128_INITIALIZER +#define PCG64F_INITIALIZER PCG_STATE_MCG_128_INITIALIZER +#endif + +#if PCG_HAS_128BIT_OPS +#define PCG8SI_INITIALIZER PCG_STATE_ONESEQ_8_INITIALIZER +#define PCG16SI_INITIALIZER PCG_STATE_ONESEQ_16_INITIALIZER +#define PCG32SI_INITIALIZER PCG_STATE_ONESEQ_32_INITIALIZER +#define PCG64SI_INITIALIZER PCG_STATE_ONESEQ_64_INITIALIZER +#define PCG128SI_INITIALIZER PCG_STATE_ONESEQ_128_INITIALIZER +#endif + +#if PCG_HAS_128BIT_OPS +#define PCG8I_INITIALIZER PCG_STATE_SETSEQ_8_INITIALIZER +#define PCG16I_INITIALIZER PCG_STATE_SETSEQ_16_INITIALIZER +#define PCG32I_INITIALIZER PCG_STATE_SETSEQ_32_INITIALIZER +#define PCG64I_INITIALIZER PCG_STATE_SETSEQ_64_INITIALIZER +#define PCG128I_INITIALIZER PCG_STATE_SETSEQ_128_INITIALIZER +#endif + +#if __cplusplus +} +#endif + +#endif // PCG_VARIANTS_H_INCLUDED + diff --git a/pcg-random.cabal b/pcg-random.cabal new file mode 100644 --- /dev/null +++ b/pcg-random.cabal @@ -0,0 +1,73 @@ +name: pcg-random +version: 0.1.0.0 +synopsis: Haskell bindings to the PCG random number generator. +description: + PCG is a family of simple fast space-efficient statistically good + algorithms for random number generation with better-than-typical + cryptographic security. + . + This library implements bindings to the standard C implementation. + This includes the standard, unique, fast and single variants in the + pcg family. There is a pure implementation that can be used as a + generator with the random package as well as a faster primitive api + that includes functions for generating common types. + . + The generators in this module are suitable for use in parallel but + make sure threads don't share the same generator or things will go + horribly wrong. + +license: BSD3 +license-file: LICENSE +extra-source-files: README.md c/LICENSE.txt +author: Christopher Chalmers +maintainer: c.chalmers@me.com +Homepage: http://github.com/cchalmers/pcg-random +Bug-reports: http://github.com/cchalmers/pcg-random/issues +copyright: (c) 2014. Christopher Chalmers +category: System +build-type: Simple +cabal-version: >=1.10 + +source-repository head + type: git + location: git://github.com/cchalmers/pcg-random.git + +library + exposed-modules: + System.Random.PCG + System.Random.PCG.Class + System.Random.PCG.Fast + System.Random.PCG.Unique + System.Random.PCG.Single + hs-source-dirs: src + other-extensions: BangPatterns, CPP, ForeignFunctionInterface + build-depends: + base >=4.3 && < 5, + primitive >=0.4 && < 0.6, + random >=1.0 && < 2.0, + time + default-language: Haskell2010 + cc-options: + -DMEXP=19937 + -DNDEBUG + -finline-functions -fomit-frame-pointer + -fno-strict-aliasing + -std=c99 -fPIC + + c-sources: + c/pcg-advance-64.c + c/pcg-global-64.c + c/pcg-output-64.c + c/pcg-rngs-64.c + include-dirs: c + includes: pcg_variants.h + install-includes: pcg_variants.h + +test-suite doctests + default-language: Haskell2010 + type: exitcode-stdio-1.0 + main-is: doctest.hs + hs-source-dirs: test + build-depends: + base, + doctest diff --git a/src/System/Random/PCG.hs b/src/System/Random/PCG.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/PCG.hs @@ -0,0 +1,258 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE MultiParamTypeClasses #-} +#if __GLASGOW_HASKELL__ >= 707 +{-# LANGUAGE RoleAnnotations #-} +#endif +-------------------------------------------------------------------- +-- | +-- Module : System.Random.PCG +-- Copyright : Copyright (c) 2014, Christopher Chalmers +-- License : BSD3 +-- Maintainer : Christopher Chalmers +-- Stability : experimental +-- Portability: CPP, FFI +-- +-- Standard PCG Random Number Generator with chosen streams. See +-- for details. +-- +-- @ +-- import Control.Monad.ST +-- import System.Random.PCG +-- +-- three :: [Double] +-- three = runST $ do +-- g <- create +-- a <- uniform g +-- b <- uniform g +-- c <- uniform g +-- return [a,b,c] +-- @ + +module System.Random.PCG + ( -- * Generator + Gen, GenIO, GenST + , create, createSystemRandom, initialize, withSystemRandom + + -- * Getting random numbers + , Variate (..) + , advance, retract + + -- * Frozen generator + , FrozenGen + , save, restore, seed, initFrozen + + -- * Type restricted versions + , uniformW8, uniformW16, uniformW32, uniformW64 + , uniformI8, uniformI16, uniformI32, uniformI64 + , uniformF, uniformD, uniformBool + ) where + +import Control.Applicative +import Control.Monad.Primitive +import Foreign +import System.IO.Unsafe +import System.Random + +import System.Random.PCG.Class + +-- $setup +-- >>> import System.Random.PCG +-- >>> import System.Random.PCG.Class +-- >>> import Control.Monad + +------------------------------------------------------------------------ +-- State +------------------------------------------------------------------------ + +-- | Immutable snapshot of the state of a 'Gen'. +data FrozenGen = FrozenGen {-# UNPACK #-} !Word64 {-# UNPACK #-} !Word64 + deriving (Show, Eq, Ord) + +-- | Save the state of a 'Gen' in a 'FrozenGen'. +save :: PrimMonad m => Gen (PrimState m) -> m FrozenGen +save (Gen p) = unsafePrimToPrim (peek p) +{-# INLINE save #-} + +-- | Restore a 'Gen' from a 'FrozenGen'. +restore :: PrimMonad m => FrozenGen -> m (Gen (PrimState m)) +restore s = unsafePrimToPrim $ do + p <- malloc + poke p s + return (Gen p) +{-# INLINE restore #-} + +-- | Fixed seed. +seed :: FrozenGen +seed = FrozenGen 0x853c49e6748fea9b 0xda3e39cb94b95bdb + +-- | Generate a new seed using two 'Word64's. +-- +-- >>> initFrozen 0 0 +-- FrozenGen 6364136223846793006 1 +initFrozen :: Word64 -> Word64 -> FrozenGen +initFrozen w1 w2 = unsafeDupablePerformIO $ do + p <- malloc + pcg32_srandom_r p w1 w2 + peek p <* free p +{-# INLINE initFrozen #-} + +instance Storable FrozenGen where + sizeOf _ = 16 + {-# INLINE sizeOf #-} + alignment _ = 8 + {-# INLINE alignment #-} + poke ptr (FrozenGen x y) = poke ptr' x >> pokeElemOff ptr' 1 y + where ptr' = castPtr ptr + {-# INLINE poke #-} + peek ptr = FrozenGen <$> peek ptr' <*> peekElemOff ptr' 1 + where ptr' = castPtr ptr + {-# INLINE peek #-} + +------------------------------------------------------------------------ +-- PrimMonad interface +------------------------------------------------------------------------ + +-- | State of the random number generator +newtype Gen s = Gen (Ptr FrozenGen) + deriving (Eq, Ord) + +#if __GLASGOW_HASKELL__ >= 707 +type role Gen representational +#endif + +-- this should be type safe because the Gen cannot escape its PrimMonad + +-- | Type alias of 'Gen' specialized to 'IO'. +type GenIO = Gen RealWorld + +-- | Type alias of 'Gen' specialized to 'ST'. ( +type GenST s = Gen s +-- Note this doesn't force it to be in ST. You can write (STGen Realworld) +-- and it'll work in IO. Writing STGen s = Gen (PrimState (ST s)) doesn't +-- solve this. + +-- | Create a 'Gen' from a fixed initial 'seed'. +create :: PrimMonad m => m (Gen (PrimState m)) +create = restore seed + +-- | Initialize a generator with two words. +-- +-- >>> initialize 0 0 >>= save +-- FrozenGen 6364136223846793006 1 +initialize :: PrimMonad m => Word64 -> Word64 -> m (Gen (PrimState m)) +initialize a b = unsafePrimToPrim $ do + p <- malloc + pcg32_srandom_r p a b + return (Gen p) + +-- | Seed with system random number. (\"@\/dev\/urandom@\" on Unix-like +-- systems, time otherwise). +withSystemRandom :: PrimMonad m => (Gen (PrimState m) -> m a) -> IO a +withSystemRandom f = do + w1 <- sysRandom + w2 <- sysRandom + unsafePrimToIO $ initialize w1 w2 >>= f + +-- | Seed a PRNG with data from the system's fast source of pseudo-random +-- numbers. All the caveats of 'withSystemRandom' apply here as well. +createSystemRandom :: IO GenIO +createSystemRandom = withSystemRandom (return :: GenIO -> IO GenIO) + +-- -- | Generate a uniform 'Word32' bounded by the given bound. +-- uniformB :: PrimMonad m => Word32 -> Gen (PrimState m) -> m Word32 +-- uniformB u (Gen p) = unsafePrimToPrim $ pcg32_boundedrand_r p u +-- {-# INLINE uniformB #-} + +-- | Advance the given generator n steps in log(n) time. (Note that a +-- \"step\" is a single random 32-bit (or less) 'Variate'. Data types +-- such as 'Double' or 'Word64' require two \"steps\".) +-- +-- >>> create >>= \g -> replicateM_ 1000 (uniformW32 g) >> uniformW32 g +-- 3640764222 +-- >>> create >>= \g -> replicateM_ 500 (uniformD g) >> uniformW32 g +-- 3640764222 +-- >>> create >>= \g -> advance 1000 g >> uniformW32 g +-- 3640764222 +advance :: PrimMonad m => Word64 -> Gen (PrimState m) -> m () +advance u (Gen p) = unsafePrimToPrim $ pcg32_advance_r p u +{-# INLINE advance #-} + +-- | Retract the given generator n steps in log(2^64-n) time. This +-- is just @advance (-n)@. +-- +-- >>> create >>= \g -> replicateM 3 (uniformW32 g) +-- [355248013,41705475,3406281715] +-- >>> create >>= \g -> retract 1 g >> replicateM 3 (uniformW32 g) +-- [19683962,355248013,41705475] +retract :: PrimMonad m => Word64 -> Gen (PrimState m) -> m () +retract u g = advance (-u) g +{-# INLINE retract #-} + +------------------------------------------------------------------------ +-- Foreign calls +------------------------------------------------------------------------ + +-- It shouldn't be too hard to impliment the algorithm in pure haskell. +-- For now just use the c interface. + +-- For whatever reason, calling the #defined versions doesn't seem to work +-- so we need to call the low-level api directly + +foreign import ccall unsafe "pcg_setseq_64_srandom_r" + pcg32_srandom_r :: Ptr FrozenGen -> Word64 -> Word64 -> IO () + +foreign import ccall unsafe "pcg_setseq_64_xsh_rr_32_random_r" + pcg32_random_r :: Ptr FrozenGen -> IO Word32 + +-- foreign import ccall unsafe "pcg_setseq_64_xsh_rr_32_boundedrand_r" +-- pcg32_boundedrand_r :: Ptr FrozenGen -> Word32 -> IO Word32 + +foreign import ccall unsafe "pcg_setseq_64_advance_r" + pcg32_advance_r :: Ptr FrozenGen -> Word64 -> IO () + +------------------------------------------------------------------------ +-- Instances +------------------------------------------------------------------------ + +instance (PrimMonad m, s ~ PrimState m) => Generator (Gen s) m where + uniform1 f (Gen p) = unsafePrimToPrim $ f <$> pcg32_random_r p + {-# INLINE uniform1 #-} + + uniform2 f (Gen p) = unsafePrimToPrim $ do + w1 <- pcg32_random_r p + w2 <- pcg32_random_r p + return $ f w1 w2 + {-# INLINE uniform2 #-} + +instance RandomGen FrozenGen where + next s = unsafeDupablePerformIO $ do + p <- malloc + poke p s + w1 <- pcg32_random_r p + w2 <- pcg32_random_r p + s' <- peek p + free p + return (wordsTo64Bit w1 w2, s') + {-# INLINE next #-} + + split s = unsafeDupablePerformIO $ do + p <- malloc + poke p s + w1 <- pcg32_random_r p + w2 <- pcg32_random_r p + w3 <- pcg32_random_r p + w4 <- pcg32_random_r p + w5 <- pcg32_random_r p + w6 <- pcg32_random_r p + w7 <- pcg32_random_r p + w8 <- pcg32_random_r p + pcg32_srandom_r p (wordsTo64Bit w1 w2) (wordsTo64Bit w3 w4) + s1 <- peek p + pcg32_srandom_r p (wordsTo64Bit w5 w6) (wordsTo64Bit w7 w8) + s2 <- peek p + free p + return (s1,s2) + diff --git a/src/System/Random/PCG/Class.hs b/src/System/Random/PCG/Class.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/PCG/Class.hs @@ -0,0 +1,350 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UnboxedTuples #-} +-- | +-- Module : System.Random.PCG.Class +-- Copyright : Copyright (c) 2014, Christopher Chalmers +-- License : BSD3 +-- Maintainer : Christopher Chalmers +-- Stability : experimental +-- Portability: CPP +-- +-- Classes for working with random numbers along with utility functions. +-- In a future release this module may disappear and use another module +-- for this functionality. +module System.Random.PCG.Class + ( -- * Classes + Generator (..) + , Variate (..) + + -- * Type restricted versions + , uniformW8, uniformW16, uniformW32, uniformW64 + , uniformI8, uniformI16, uniformI32, uniformI64 + , uniformF, uniformD, uniformBool + + -- * Utilities + , Unsigned + , wordsTo64Bit + , wordToBool + , wordToFloat + , wordsToDouble + , sysRandom + ) where + +import Control.Monad +import Data.Bits +import Data.Int +import Data.Word +import Foreign (allocaBytes, peek) +import System.IO + +import qualified Control.Exception as E +import Data.IORef (atomicModifyIORef, newIORef) +import Data.Ratio (numerator, (%)) +import Data.Time.Clock.POSIX (getPOSIXTime) +import System.CPUTime (cpuTimePrecision, getCPUTime) +import System.IO.Unsafe (unsafePerformIO) + +class Monad m => Generator g m where + uniform1 :: (Word32 -> a) -> g -> m a + uniform2 :: (Word32 -> Word32 -> a) -> g -> m a + +class Variate a where + -- | Generate a uniformly distributed random vairate. + -- + -- * Use entire range for integral types. + -- + -- * Use (0,1] range for floating types. + uniform :: Generator g m => g -> m a + + -- | Generate a uniformly distributed random vairate in the given + -- range. + -- + -- * Use inclusive range for integral types. + -- + -- * Use (a,b] range for floating types. + uniformR :: Generator g m => (a,a) -> g -> m a + +------------------------------------------------------------------------ +-- Variate instances +------------------------------------------------------------------------ + +instance Variate Int8 where + uniform = uniform1 fromIntegral + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Int16 where + uniform = uniform1 fromIntegral + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Int32 where + uniform = uniform1 fromIntegral + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Int64 where + uniform = uniform2 wordsTo64Bit + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Word8 where + uniform = uniform1 fromIntegral + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Word16 where + uniform = uniform1 fromIntegral + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Word32 where + uniform = uniform1 fromIntegral + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Word64 where + uniform = uniform2 wordsTo64Bit + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Bool where + uniform = uniform1 wordToBool + {-# INLINE uniform #-} + uniformR (False,True) g = uniform g + uniformR (False,False) _ = return False + uniformR (True,True) _ = return True + uniformR (True,False) g = uniform g + {-# INLINE uniformR #-} + +instance Variate Float where + uniform = uniform1 wordToFloat + {-# INLINE uniform #-} + uniformR (x1,x2) = uniform1 (\w -> x1 + (x2-x1) * wordToFloat w) + {-# INLINE uniformR #-} + +instance Variate Double where + uniform = uniform2 wordsToDouble + {-# INLINE uniform #-} + uniformR (x1,x2) = uniform2 (\w1 w2 -> x1 + (x2-x1) * wordsToDouble w1 w2) + {-# INLINE uniformR #-} + +instance Variate Word where +#if WORD_SIZE_IN_BITS < 64 + uniform = uniform1 fromIntegral +#else + uniform = uniform2 wordsTo64Bit +#endif + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance Variate Int where +#if WORD_SIZE_IN_BITS < 64 + uniform = uniform1 fromIntegral +#else + uniform = uniform2 wordsTo64Bit +#endif + {-# INLINE uniform #-} + uniformR a g = uniformRange a g + {-# INLINE uniformR #-} + +instance (Variate a, Variate b) => Variate (a,b) where + uniform g = (,) `liftM` uniform g `ap` uniform g + {-# INLINE uniform #-} + uniformR ((x1,y1),(x2,y2)) g = (,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g + {-# INLINE uniformR #-} + +instance (Variate a, Variate b, Variate c) => Variate (a,b,c) where + uniform g = (,,) `liftM` uniform g `ap` uniform g `ap` uniform g + {-# INLINE uniform #-} + uniformR ((x1,y1,z1),(x2,y2,z2)) g = + (,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` uniformR (z1,z2) g + {-# INLINE uniformR #-} + +instance (Variate a, Variate b, Variate c, Variate d) => Variate (a,b,c,d) where + uniform g = (,,,) `liftM` uniform g `ap` uniform g `ap` uniform g `ap` uniform g + {-# INLINE uniform #-} + uniformR ((x1,y1,z1,t1),(x2,y2,z2,t2)) g = + (,,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` + uniformR (z1,z2) g `ap` uniformR (t1,t2) g + {-# INLINE uniformR #-} + +------------------------------------------------------------------------ +-- Type restricted versions +------------------------------------------------------------------------ + +uniformI8 :: Generator g m => g -> m Int8 +uniformI8 = uniform +{-# INLINE uniformI8 #-} + +uniformI16 :: Generator g m => g -> m Int16 +uniformI16 = uniform +{-# INLINE uniformI16 #-} + +uniformI32 :: Generator g m => g -> m Int32 +uniformI32 = uniform +{-# INLINE uniformI32 #-} + +uniformI64 :: Generator g m => g -> m Int64 +uniformI64 = uniform +{-# INLINE uniformI64 #-} + +uniformW8 :: Generator g m => g -> m Word8 +uniformW8 = uniform +{-# INLINE uniformW8 #-} + +uniformW16 :: Generator g m => g -> m Word16 +uniformW16 = uniform +{-# INLINE uniformW16 #-} + +uniformW32 :: Generator g m => g -> m Word32 +uniformW32 = uniform +{-# INLINE uniformW32 #-} + +uniformW64 :: Generator g m => g -> m Word64 +uniformW64 = uniform +{-# INLINE uniformW64 #-} + +uniformBool :: Generator g m => g -> m Bool +uniformBool = uniform +{-# INLINE uniformBool #-} + +uniformF :: Generator g m => g -> m Float +uniformF = uniform +{-# INLINE uniformF #-} + +uniformD :: Generator g m => g -> m Double +uniformD = uniform +{-# INLINE uniformD #-} + +------------------------------------------------------------------------ +-- Utilities +------------------------------------------------------------------------ + +sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a +sub x y = fromIntegral x - fromIntegral y +{-# INLINE sub #-} + +add :: (Integral a, Integral (Unsigned a)) => a -> Unsigned a -> a +add m x = m + fromIntegral x +{-# INLINE add #-} + +wordsTo64Bit :: Integral a => Word32 -> Word32 -> a +wordsTo64Bit x y = + fromIntegral ((fromIntegral x `shiftL` 32) .|. fromIntegral y :: Word64) +{-# INLINE wordsTo64Bit #-} + +wordToBool :: Word32 -> Bool +wordToBool i = (i .&. 1) /= 0 +{-# INLINE wordToBool #-} + +wordToFloat :: Word32 -> Float +wordToFloat x = (fromIntegral i * m_inv_32) + 0.5 + m_inv_33 + where m_inv_33 = 1.16415321826934814453125e-10 + m_inv_32 = 2.3283064365386962890625e-10 + i = fromIntegral x :: Int32 +{-# INLINE wordToFloat #-} + +wordsToDouble :: Word32 -> Word32 -> Double +wordsToDouble x y = (fromIntegral u * m_inv_32 + (0.5 + m_inv_53) + + fromIntegral (v .&. 0xFFFFF) * m_inv_52) + where m_inv_52 = 2.220446049250313080847263336181640625e-16 + m_inv_53 = 1.1102230246251565404236316680908203125e-16 + m_inv_32 = 2.3283064365386962890625e-10 + u = fromIntegral x :: Int32 + v = fromIntegral y :: Int32 +{-# INLINE wordsToDouble #-} + +-- IO randoms + +devRandom :: IO Word64 +devRandom = + allocaBytes 8 $ \buf -> do + nread <- withBinaryFile "/dev/urandom" ReadMode $ \h -> hGetBuf h buf 8 + when (nread /= 8) $ error "unable to read from /dev/urandom" + peek buf + +-- Acquire seed from current time. This is horrible fall-back for +-- Windows system. +acquireSeedTime :: IO Word64 +acquireSeedTime = do + c <- (numerator . (%cpuTimePrecision)) `liftM` getCPUTime + t <- toRational `liftM` getPOSIXTime + let n = fromIntegral (numerator t) :: Word64 + return $ wordsTo64Bit (fromIntegral c) (fromIntegral n) + +-- | Get a random number from system source. If \"@\/dev\/urandom@\" is +-- not found return inferior random number from time. +sysRandom :: IO Word64 +sysRandom = + devRandom `E.catch` \(_ :: E.IOException) -> do + seen <- atomicModifyIORef warned ((,) True) + unless seen $ E.handle (\(_::E.IOException) -> return ()) $ do + hPutStrLn stderr ("Warning: Couldn't open /dev/urandom") + hPutStrLn stderr ("Warning: using system clock for seed instead " ++ + "(quality will be lower)") + acquireSeedTime + where + warned = unsafePerformIO $ newIORef False + {-# NOINLINE warned #-} + +uniformRange :: ( Generator g m + , Integral a, Bounded a, Variate a + , Integral (Unsigned a), Bounded (Unsigned a), Variate (Unsigned a)) + => (a,a) -> g -> m a +uniformRange (x1,x2) g + | n == 0 = uniform g -- Abuse overflow in unsigned types + | otherwise = loop + where + -- Allow ranges where x2= 707 +{-# LANGUAGE RoleAnnotations #-} +#endif +-- | +-- Module : System.Random.PCG.Fast +-- Copyright : Copyright (c) 2014, Christopher Chalmers +-- License : BSD3 +-- Maintainer : Christopher Chalmers +-- Stability : experimental +-- Portability: CPP, FFI +-- +-- Fast variant of the PCG random number generator. This module performs +-- around 20% faster the multiple streams version but produces slightly +-- lower quality (still good) random numbers. +-- +-- See for details. +-- +-- @ +-- import Control.Monad.ST +-- import System.Random.PCG.Fast +-- +-- three :: [Double] +-- three = runST $ do +-- g <- create +-- a <- uniform g +-- b <- uniform g +-- c <- uniform g +-- return [a,b,c] +-- @ +module System.Random.PCG.Fast + ( -- * Gen + Gen, GenIO, GenST + , create, createSystemRandom, initialize, withSystemRandom + + -- * Getting random numbers + , Variate (..) + , advance, retract + + -- * Seeds + , FrozenGen, save, restore, seed, initFrozen + + -- * Type restricted versions + , uniformW8, uniformW16, uniformW32, uniformW64 + , uniformI8, uniformI16, uniformI32, uniformI64 + , uniformF, uniformD, uniformBool + ) where + +import Control.Applicative +import Control.Monad.Primitive +import Foreign +import System.IO.Unsafe +import System.Random + +import System.Random.PCG.Class + +-- $setup +-- >>> import System.Random.PCG.Fast +-- >>> import System.Random.PCG.Class +-- >>> import Control.Monad + +------------------------------------------------------------------------ +-- Seed +------------------------------------------------------------------------ + +newtype FrozenGen = FrozenGen Word64 + deriving (Show, Eq, Ord, Storable) + +-- | Save the state of a 'Gen' in a 'Seed'. +save :: PrimMonad m => Gen (PrimState m) -> m FrozenGen +save (Gen p) = unsafePrimToPrim (peek p) +{-# INLINE save #-} + +-- | Restore a 'Gen' from a 'Seed'. +restore :: PrimMonad m => FrozenGen -> m (Gen (PrimState m)) +restore s = unsafePrimToPrim $ do + p <- malloc + poke p s + return (Gen p) +{-# INLINE restore #-} + +-- | Generate a new seed using single 'Word64'. +-- +-- >>> initFrozen 0 +-- FrozenGen 1 +initFrozen :: Word64 -> FrozenGen +initFrozen w = unsafeDupablePerformIO $ do + p <- malloc + pcg32f_srandom_r p w + peek p <* free p +{-# INLINE initFrozen #-} + +-- | Standard initial seed. +seed :: FrozenGen +seed = FrozenGen 0xcafef00dd15ea5e5 + +-- | Create a 'Gen' from a fixed initial seed. +create :: PrimMonad m => m (Gen (PrimState m)) +create = restore seed + +------------------------------------------------------------------------ +-- Gen +------------------------------------------------------------------------ + +-- | State of the random number generator +newtype Gen s = Gen (Ptr FrozenGen) + deriving (Eq, Ord) +#if __GLASGOW_HASKELL__ >= 707 +type role Gen representational +#endif + +type GenIO = Gen RealWorld +type GenST = Gen + +-- | Initialize a generator a single word. +-- +-- >>> initialize 0 >>= save +-- FrozenGen 1 +initialize :: PrimMonad m => Word64 -> m (Gen (PrimState m)) +initialize a = unsafePrimToPrim $ do + p <- malloc + pcg32f_srandom_r p a + return (Gen p) + +-- | Seed with system random number. (\"@\/dev\/urandom@\" on Unix-like +-- systems, time otherwise). +withSystemRandom :: PrimMonad m => (Gen (PrimState m) -> m a) -> IO a +withSystemRandom f = do + w <- sysRandom + unsafePrimToIO $ initialize w >>= f + +-- | Seed a PRNG with data from the system's fast source of pseudo-random +-- numbers. All the caveats of 'withSystemRandom' apply here as well. +createSystemRandom :: IO GenIO +createSystemRandom = withSystemRandom (return :: GenIO -> IO GenIO) + +-- -- | Generate a uniform 'Word32' bounded by the given bound. +-- uniformB :: PrimMonad m => Word32 -> Gen (PrimState m) -> m Word32 +-- uniformB u (Gen p) = unsafePrimToPrim $ pcg32f_boundedrand_r p u +-- {-# INLINE uniformB #-} + +-- | Advance the given generator n steps in log(n) time. (Note that a +-- \"step\" is a single random 32-bit (or less) 'Variate'. Data types +-- such as 'Double' or 'Word64' require two \"steps\".) +-- +-- >>> create >>= \g -> replicateM_ 1000 (uniformW32 g) >> uniformW32 g +-- 3725702568 +-- >>> create >>= \g -> replicateM_ 500 (uniformD g) >> uniformW32 g +-- 3725702568 +-- >>> create >>= \g -> advance 1000 g >> uniformW32 g +-- 3725702568 +advance :: PrimMonad m => Word64 -> Gen (PrimState m) -> m () +advance u (Gen p) = unsafePrimToPrim $ pcg32f_advance_r p u +{-# INLINE advance #-} + +-- | Retract the given generator n steps in log(2^64-n) time. This +-- is just @advance (-n)@. +-- +-- >>> create >>= \g -> replicateM 3 (uniformW32 g) +-- [2951688802,2698927131,361549788] +-- >>> create >>= \g -> retract 1 g >> replicateM 3 (uniformW32 g) +-- [954135925,2951688802,2698927131] +retract :: PrimMonad m => Word64 -> Gen (PrimState m) -> m () +retract u g = advance (-u) g +{-# INLINE retract #-} + +------------------------------------------------------------------------ +-- Foreign calls +------------------------------------------------------------------------ + +foreign import ccall unsafe "pcg_mcg_64_srandom_r" + pcg32f_srandom_r :: Ptr FrozenGen -> Word64 -> IO () + +foreign import ccall unsafe "pcg_mcg_64_xsh_rs_32_random_r" + pcg32f_random_r :: Ptr FrozenGen -> IO Word32 + +-- foreign import ccall unsafe "pcg_mcg_64_xsh_rs_32_boundedrand_r" +-- pcg32f_boundedrand_r :: Ptr FrozenGen -> Word32 -> IO Word32 + +foreign import ccall unsafe "pcg_mcg_64_advance_r" + pcg32f_advance_r :: Ptr FrozenGen -> Word64 -> IO () + +------------------------------------------------------------------------ +-- Instances +------------------------------------------------------------------------ + +instance (PrimMonad m, s ~ PrimState m) => Generator (Gen s) m where + uniform1 f (Gen p) = unsafePrimToPrim $ f <$> pcg32f_random_r p + {-# INLINE uniform1 #-} + + uniform2 f (Gen p) = unsafePrimToPrim $ do + w1 <- pcg32f_random_r p + w2 <- pcg32f_random_r p + return $ f w1 w2 + {-# INLINE uniform2 #-} + +instance RandomGen FrozenGen where + next s = unsafeDupablePerformIO $ do + p <- malloc + poke p s + w1 <- pcg32f_random_r p + w2 <- pcg32f_random_r p + s' <- peek p + free p + return (wordsTo64Bit w1 w2, s') + {-# INLINE next #-} + + split s = unsafeDupablePerformIO $ do + p <- malloc + poke p s + w1 <- pcg32f_random_r p + w2 <- pcg32f_random_r p + w3 <- pcg32f_random_r p + w4 <- pcg32f_random_r p + pcg32f_srandom_r p (wordsTo64Bit w1 w2) + s1 <- peek p + pcg32f_srandom_r p (wordsTo64Bit w3 w4) + s2 <- peek p + free p + return (s1,s2) + diff --git a/src/System/Random/PCG/Single.hs b/src/System/Random/PCG/Single.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/PCG/Single.hs @@ -0,0 +1,219 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +#if __GLASGOW_HASKELL__ >= 707 +{-# LANGUAGE RoleAnnotations #-} +#endif +-------------------------------------------------------------------- +-- | +-- Module : System.Random.PCG.Single +-- Copyright : Copyright (c) 2014, Christopher Chalmers +-- License : BSD3 +-- Maintainer : Christopher Chalmers +-- Stability : experimental +-- Portability: CPP, FFI +-- +-- Single variant of the PCG random number generator. This module only +-- uses a single stream. See for details. +-- +-- @ +-- import Control.Monad.ST +-- import System.Random.PCG.Single +-- +-- three :: [Double] +-- three = runST $ do +-- g <- create +-- a <- uniform g +-- b <- uniform g +-- c <- uniform g +-- return [a,b,c] +-- @ +module System.Random.PCG.Single + ( -- * Gen + Gen, GenIO, GenST + , create, createSystemRandom, initialize, withSystemRandom + + -- * Getting random numbers + , Variate (..) + , advance, retract + + -- * Seeds + , FrozenGen + , save, restore, seed, initFrozen + + -- * Type restricted versions + , uniformW8, uniformW16, uniformW32, uniformW64 + , uniformI8, uniformI16, uniformI32, uniformI64 + , uniformF, uniformD, uniformBool + ) where + +import Control.Applicative +import Control.Monad.Primitive +import Foreign +import System.IO.Unsafe +import System.Random +import System.Random.PCG.Class + +-- $setup +-- >>> import System.Random.PCG.Single +-- >>> import System.Random.PCG.Class +-- >>> import Control.Monad + +------------------------------------------------------------------------ +-- Seed +------------------------------------------------------------------------ + +newtype FrozenGen = FrozenGen Word64 + deriving (Show, Eq, Ord, Storable) + +-- | Save the state of a 'Gen' in a 'Seed'. +save :: PrimMonad m => Gen (PrimState m) -> m FrozenGen +save (Gen p) = unsafePrimToPrim (peek p) +{-# INLINE save #-} + +-- | Restore a 'Gen' from a 'Seed'. +restore :: PrimMonad m => FrozenGen -> m (Gen (PrimState m)) +restore s = unsafePrimToPrim $ do + p <- malloc + poke p s + return (Gen p) +{-# INLINE restore #-} + +-- | Generate a new seed using single 'Word64'. +-- +-- >>> initFrozen 0 +-- FrozenGen 1876011003808476466 +initFrozen :: Word64 -> FrozenGen +initFrozen w = unsafeDupablePerformIO . alloca $ \p -> + pcg32s_srandom_r p w >> peek p +{-# INLINE initFrozen #-} + +-- | Standard initial seed. +seed :: FrozenGen +seed = FrozenGen 0x4d595df4d0f33173 + +-- | Create a 'Gen' from a fixed initial seed. +create :: PrimMonad m => m (Gen (PrimState m)) +create = restore seed + +------------------------------------------------------------------------ +-- Gen +------------------------------------------------------------------------ + +-- | State of the random number generator +newtype Gen s = Gen (Ptr FrozenGen) + deriving (Eq, Ord) +#if __GLASGOW_HASKELL__ >= 707 +type role Gen representational +#endif + +type GenIO = Gen RealWorld +type GenST = Gen + +-- | Initialize a generator a single word. +-- +-- >>> initialize 0 >>= save +-- FrozenGen 1876011003808476466 +initialize :: PrimMonad m => Word64 -> m (Gen (PrimState m)) +initialize a = unsafePrimToPrim $ do + p <- malloc + pcg32s_srandom_r p a + return (Gen p) + +-- | Seed with system random number. (\"@\/dev\/urandom@\" on Unix-like +-- systems, time otherwise). +withSystemRandom :: (GenIO -> IO a) -> IO a +withSystemRandom f = sysRandom >>= initialize >>= f + +-- | Seed a PRNG with data from the system's fast source of pseudo-random +-- numbers. All the caveats of 'withSystemRandom' apply here as well. +createSystemRandom :: IO GenIO +createSystemRandom = withSystemRandom return + +-- -- | Generate a uniform 'Word32' bounded by the given bound. +-- uniformB :: PrimMonad m => Word32 -> Gen (PrimState m) -> m Word32 +-- uniformB u (Gen p) = unsafePrimToPrim $ pcg32s_boundedrand_r p u +-- {-# INLINE uniformB #-} + +-- | Advance the given generator n steps in log(n) time. (Note that a +-- \"step\" is a single random 32-bit (or less) 'Variate'. Data types +-- such as 'Double' or 'Word64' require two \"steps\".) +-- +-- >>> create >>= \g -> replicateM_ 1000 (uniformW32 g) >> uniformW32 g +-- 787992099 +-- >>> create >>= \g -> replicateM_ 500 (uniformD g) >> uniformW32 g +-- 787992099 +-- >>> create >>= \g -> advance 1000 g >> uniformW32 g +-- 787992099 +advance :: PrimMonad m => Word64 -> Gen (PrimState m) -> m () +advance u (Gen p) = unsafePrimToPrim $ pcg32s_advance_r p u +{-# INLINE advance #-} + +-- | Retract the given generator n steps in log(2^64-n) time. This +-- is just @advance (-n)@. +-- +-- >>> create >>= \g -> replicateM 3 (uniformW32 g) +-- [1499251125,2252584727,2021901712] +-- >>> create >>= \g -> retract 1 g >> replicateM 3 (uniformW32 g) +-- [4031155981,1499251125,2252584727] +retract :: PrimMonad m => Word64 -> Gen (PrimState m) -> m () +retract u g = advance (-u) g +{-# INLINE retract #-} + +------------------------------------------------------------------------ +-- Foreign calls +------------------------------------------------------------------------ + +foreign import ccall unsafe "pcg_oneseq_64_srandom_r" + pcg32s_srandom_r :: Ptr FrozenGen -> Word64 -> IO () + +foreign import ccall unsafe "pcg_oneseq_64_xsh_rs_32_random_r" + pcg32s_random_r :: Ptr FrozenGen -> IO Word32 + +-- foreign import ccall unsafe "pcg_oneseq_64_xsh_rs_32_boundedrand_r" +-- pcg32s_boundedrand_r :: Ptr FrozenGen -> Word32 -> IO Word32 + +foreign import ccall unsafe "pcg_oneseq_64_advance_r" + pcg32s_advance_r :: Ptr FrozenGen -> Word64 -> IO () + +------------------------------------------------------------------------ +-- Instances +------------------------------------------------------------------------ + +instance (PrimMonad m, s ~ PrimState m) => Generator (Gen s) m where + uniform1 f (Gen p) = unsafePrimToPrim $ + f <$> pcg32s_random_r p + {-# INLINE uniform1 #-} + + uniform2 f (Gen p) = unsafePrimToPrim $ + f <$> pcg32s_random_r p <*> pcg32s_random_r p + {-# INLINE uniform2 #-} + +instance RandomGen FrozenGen where + next s = unsafeDupablePerformIO $ do + p <- malloc + poke p s + w1 <- pcg32s_random_r p + w2 <- pcg32s_random_r p + s' <- peek p + free p + return (wordsTo64Bit w1 w2, s') + {-# INLINE next #-} + + split s = unsafeDupablePerformIO $ do + p <- malloc + poke p s + w1 <- pcg32s_random_r p + w2 <- pcg32s_random_r p + w3 <- pcg32s_random_r p + w4 <- pcg32s_random_r p + pcg32s_srandom_r p (wordsTo64Bit w1 w2) + s1 <- peek p + pcg32s_srandom_r p (wordsTo64Bit w3 w4) + s2 <- peek p + free p + return (s1,s2) + diff --git a/src/System/Random/PCG/Unique.hs b/src/System/Random/PCG/Unique.hs new file mode 100644 --- /dev/null +++ b/src/System/Random/PCG/Unique.hs @@ -0,0 +1,142 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +{-# LANGUAGE MultiParamTypeClasses #-} +-------------------------------------------------------------------- +-- | +-- Module : System.Random.PCG.Unique +-- Copyright : Copyright (c) 2014, Christopher Chalmers +-- License : BSD3 +-- Maintainer : Christopher Chalmers +-- Stability : experimental +-- Portability: CPP, FFI +-- +-- Unique variant of the PCG random number generator. Guarantees the +-- sequence to be unique by using the pointer address to select the +-- output sequence. +-- +-- There is no way to freeze the state because then it wouldn't be +-- unique anymore. Also, generators can't be initialized in ST because +-- we don't know what pointer reference they'll get. +-- +-- See for details. +-- +-- @ +-- import System.Random.PCG.Unique +-- +-- three :: IO [Double] +-- three = do +-- g <- create +-- a <- uniform g +-- b <- uniform g +-- c <- uniform g +-- return [a,b,c] +-- @ + +module System.Random.PCG.Unique + ( -- * Gen + Gen + , create, createSystemRandom, initialize, withSystemRandom + + -- * Getting random numbers + , Variate (..) + , advance, retract + + -- * Type restricted versions + , uniformW8, uniformW16, uniformW32, uniformW64 + , uniformI8, uniformI16, uniformI32, uniformI64 + , uniformF, uniformD, uniformBool + ) where + +import Data.Functor +import Foreign + +import System.Random.PCG.Class + +-- | Standard initial seed. +seed :: Word64 +seed = 0x4d595df4d0f33173 + +-- | Create a 'Gen' from a fixed initial seed. +create :: IO Gen +create = initialize seed + -- do + -- p <- malloc + -- poke p seed + -- return (Gen p) + + -- Note that this does produce a unique sequence but if two generators + -- are created in sequence they'll have the similar pointer references + -- and the first couple of numbers are likely to be the same. Since + -- this is undesirable we run initialise to randomise it. + +------------------------------------------------------------------------ +-- Generator +------------------------------------------------------------------------ + +-- | State of the random number generator +newtype Gen = Gen (Ptr Word64) + deriving (Eq, Ord) + +-- | Create a generator from two words. Note: this is not the same as the +-- two words in a 'Seed'. +initialize :: Word64 -> IO Gen +initialize a = do + p <- malloc + pcg32u_srandom_r p a + return (Gen p) + +-- | Seed with system random number. (\"@\/dev\/urandom@\" on Unix-like +-- systems, time otherwise). +withSystemRandom :: (Gen -> IO a) -> IO a +withSystemRandom f = sysRandom >>= initialize >>= f + +-- | Seed a PRNG with data from the system's fast source of pseudo-random +-- numbers. All the caveats of 'withSystemRandom' apply here as well. +createSystemRandom :: IO Gen +createSystemRandom = withSystemRandom return + +-- -- | Generate a uniform 'Word32' bounded above by the given bound. +-- uniformB :: Word32 -> UGen -> IO Word32 +-- uniformB u (UGen p) = unsafePrimToPrim $ pcg32u_boundedrand_r p u +-- {-# INLINE uniformB #-} + +-- | Advance the given generator n steps in log(n) time. +advance :: Word64 -> Gen -> IO () +advance u (Gen p) = pcg32u_advance_r p u +{-# INLINE advance #-} + +-- | Retract the given generator n steps in log(2^64-n) time. This +-- is just @advance (-n)@. +retract :: Word64 -> Gen -> IO () +retract u g = advance (-u) g +{-# INLINE retract #-} + +------------------------------------------------------------------------ +-- Foreign calls +------------------------------------------------------------------------ + +foreign import ccall unsafe "pcg_unique_64_srandom_r" + pcg32u_srandom_r :: Ptr Word64 -> Word64 -> IO () + +foreign import ccall unsafe "pcg_unique_64_xsh_rs_32_random_r" + pcg32u_random_r :: Ptr Word64 -> IO Word32 + +-- foreign import ccall unsafe "pcg_unique_64_xsh_rs_32_boundedrand_r" +-- pcg32u_boundedrand_r :: Ptr Word64 -> Word32 -> IO Word32 + +foreign import ccall unsafe "pcg_unique_64_advance_r" + pcg32u_advance_r :: Ptr Word64 -> Word64 -> IO () + +------------------------------------------------------------------------ +-- Instances +------------------------------------------------------------------------ + +instance Generator Gen IO where + uniform1 f (Gen p) = f <$> pcg32u_random_r p + {-# INLINE uniform1 #-} + + uniform2 f (Gen p) = do + w1 <- pcg32u_random_r p + w2 <- pcg32u_random_r p + return $ f w1 w2 + {-# INLINE uniform2 #-} diff --git a/test/doctest.hs b/test/doctest.hs new file mode 100644 --- /dev/null +++ b/test/doctest.hs @@ -0,0 +1,11 @@ +import Test.DocTest + +main = doctest + [ "src/System/Random/PCG.hs" + , "src/System/Random/PCG/Class.hs" + , "src/System/Random/PCG/Fast.hs" + , "src/System/Random/PCG/Single.hs" + , "dist/build/c/pcg-advance-64.o" + , "dist/build/c/pcg-output-64.o" + , "dist/build/c/pcg-rngs-64.o" + ]