1 /*
2 * Copyright (c) 2021 Calvin Rose
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to
6 * deal in the Software without restriction, including without limitation the
7 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
8 * sell copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
19 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
20 * IN THE SOFTWARE.
21 */
22
23 #ifndef JANET_AMALG
24 #include "features.h"
25 #include <janet.h>
26 #include "state.h"
27 #include "util.h"
28 #endif
29
30 #include <math.h>
31
32 static int janet_rng_get(void *p, Janet key, Janet *out);
33 static Janet janet_rng_next(void *p, Janet key);
34
janet_rng_marshal(void * p,JanetMarshalContext * ctx)35 static void janet_rng_marshal(void *p, JanetMarshalContext *ctx) {
36 JanetRNG *rng = (JanetRNG *)p;
37 janet_marshal_abstract(ctx, p);
38 janet_marshal_int(ctx, (int32_t) rng->a);
39 janet_marshal_int(ctx, (int32_t) rng->b);
40 janet_marshal_int(ctx, (int32_t) rng->c);
41 janet_marshal_int(ctx, (int32_t) rng->d);
42 janet_marshal_int(ctx, (int32_t) rng->counter);
43 }
44
janet_rng_unmarshal(JanetMarshalContext * ctx)45 static void *janet_rng_unmarshal(JanetMarshalContext *ctx) {
46 JanetRNG *rng = janet_unmarshal_abstract(ctx, sizeof(JanetRNG));
47 rng->a = (uint32_t) janet_unmarshal_int(ctx);
48 rng->b = (uint32_t) janet_unmarshal_int(ctx);
49 rng->c = (uint32_t) janet_unmarshal_int(ctx);
50 rng->d = (uint32_t) janet_unmarshal_int(ctx);
51 rng->counter = (uint32_t) janet_unmarshal_int(ctx);
52 return rng;
53 }
54
55 const JanetAbstractType janet_rng_type = {
56 "core/rng",
57 NULL,
58 NULL,
59 janet_rng_get,
60 NULL,
61 janet_rng_marshal,
62 janet_rng_unmarshal,
63 NULL, /* tostring */
64 NULL, /* compare */
65 NULL, /* hash */
66 janet_rng_next,
67 JANET_ATEND_NEXT
68 };
69
janet_default_rng(void)70 JanetRNG *janet_default_rng(void) {
71 return &janet_vm.rng;
72 }
73
janet_rng_seed(JanetRNG * rng,uint32_t seed)74 void janet_rng_seed(JanetRNG *rng, uint32_t seed) {
75 rng->a = seed;
76 rng->b = 0x97654321u;
77 rng->c = 123871873u;
78 rng->d = 0xf23f56c8u;
79 rng->counter = 0u;
80 /* First several numbers aren't that random. */
81 for (int i = 0; i < 16; i++) janet_rng_u32(rng);
82 }
83
janet_rng_longseed(JanetRNG * rng,const uint8_t * bytes,int32_t len)84 void janet_rng_longseed(JanetRNG *rng, const uint8_t *bytes, int32_t len) {
85 uint8_t state[16] = {0};
86 for (int32_t i = 0; i < len; i++)
87 state[i & 0xF] ^= bytes[i];
88 rng->a = state[0] + (state[1] << 8) + (state[2] << 16) + (state[3] << 24);
89 rng->b = state[4] + (state[5] << 8) + (state[6] << 16) + (state[7] << 24);
90 rng->c = state[8] + (state[9] << 8) + (state[10] << 16) + (state[11] << 24);
91 rng->d = state[12] + (state[13] << 8) + (state[14] << 16) + (state[15] << 24);
92 rng->counter = 0u;
93 /* a, b, c, d can't all be 0 */
94 if (rng->a == 0) rng->a = 1u;
95 for (int i = 0; i < 16; i++) janet_rng_u32(rng);
96 }
97
janet_rng_u32(JanetRNG * rng)98 uint32_t janet_rng_u32(JanetRNG *rng) {
99 /* Algorithm "xorwow" from p. 5 of Marsaglia, "Xorshift RNGs" */
100 uint32_t t = rng->d;
101 uint32_t const s = rng->a;
102 rng->d = rng->c;
103 rng->c = rng->b;
104 rng->b = s;
105 t ^= t >> 2;
106 t ^= t << 1;
107 t ^= s ^ (s << 4);
108 rng->a = t;
109 rng->counter += 362437;
110 return t + rng->counter;
111 }
112
janet_rng_double(JanetRNG * rng)113 double janet_rng_double(JanetRNG *rng) {
114 uint32_t hi = janet_rng_u32(rng);
115 uint32_t lo = janet_rng_u32(rng);
116 uint64_t big = (uint64_t)(lo) | (((uint64_t) hi) << 32);
117 return ldexp((double)(big >> (64 - 52)), -52);
118 }
119
120 JANET_CORE_FN(cfun_rng_make,
121 "(math/rng &opt seed)",
122 "Creates a Psuedo-Random number generator, with an optional seed. "
123 "The seed should be an unsigned 32 bit integer or a buffer. "
124 "Do not use this for cryptography. Returns a core/rng abstract type."
125 ) {
126 janet_arity(argc, 0, 1);
127 JanetRNG *rng = janet_abstract(&janet_rng_type, sizeof(JanetRNG));
128 if (argc == 1) {
129 if (janet_checkint(argv[0])) {
130 uint32_t seed = (uint32_t)(janet_getinteger(argv, 0));
131 janet_rng_seed(rng, seed);
132 } else {
133 JanetByteView bytes = janet_getbytes(argv, 0);
134 janet_rng_longseed(rng, bytes.bytes, bytes.len);
135 }
136 } else {
137 janet_rng_seed(rng, 0);
138 }
139 return janet_wrap_abstract(rng);
140 }
141
142 JANET_CORE_FN(cfun_rng_uniform,
143 "(math/rng-uniform rng)",
144 "Extract a random number in the range [0, 1) from the RNG."
145 ) {
146 janet_fixarity(argc, 1);
147 JanetRNG *rng = janet_getabstract(argv, 0, &janet_rng_type);
148 return janet_wrap_number(janet_rng_double(rng));
149 }
150
151 JANET_CORE_FN(cfun_rng_int,
152 "(math/rng-int rng &opt max)",
153 "Extract a random random integer in the range [0, max] from the RNG. If "
154 "no max is given, the default is 2^31 - 1."
155 ) {
156 janet_arity(argc, 1, 2);
157 JanetRNG *rng = janet_getabstract(argv, 0, &janet_rng_type);
158 if (argc == 1) {
159 uint32_t word = janet_rng_u32(rng) >> 1;
160 return janet_wrap_integer(word);
161 } else {
162 int32_t max = janet_optnat(argv, argc, 1, INT32_MAX);
163 if (max == 0) return janet_wrap_number(0.0);
164 uint32_t modulo = (uint32_t) max;
165 uint32_t maxgen = INT32_MAX;
166 uint32_t maxword = maxgen - (maxgen % modulo);
167 uint32_t word;
168 do {
169 word = janet_rng_u32(rng) >> 1;
170 } while (word > maxword);
171 return janet_wrap_integer(word % modulo);
172 }
173 }
174
rng_get_4bytes(JanetRNG * rng,uint8_t * buf)175 static void rng_get_4bytes(JanetRNG *rng, uint8_t *buf) {
176 uint32_t word = janet_rng_u32(rng);
177 buf[0] = word & 0xFF;
178 buf[1] = (word >> 8) & 0xFF;
179 buf[2] = (word >> 16) & 0xFF;
180 buf[3] = (word >> 24) & 0xFF;
181 }
182
183 JANET_CORE_FN(cfun_rng_buffer,
184 "(math/rng-buffer rng n &opt buf)",
185 "Get n random bytes and put them in a buffer. Creates a new buffer if no buffer is "
186 "provided, otherwise appends to the given buffer. Returns the buffer."
187 ) {
188 janet_arity(argc, 2, 3);
189 JanetRNG *rng = janet_getabstract(argv, 0, &janet_rng_type);
190 int32_t n = janet_getnat(argv, 1);
191 JanetBuffer *buffer = janet_optbuffer(argv, argc, 2, n);
192
193 /* Split into first part (that is divisible by 4), and rest */
194 int32_t first_part = n & ~3;
195 int32_t second_part = n - first_part;
196
197 /* Get first part in chunks of 4 bytes */
198 janet_buffer_extra(buffer, n);
199 uint8_t *buf = buffer->data + buffer->count;
200 for (int32_t i = 0; i < first_part; i += 4) rng_get_4bytes(rng, buf + i);
201 buffer->count += first_part;
202
203 /* Get remaining 0 - 3 bytes */
204 if (second_part) {
205 uint8_t wordbuf[4] = {0};
206 rng_get_4bytes(rng, wordbuf);
207 janet_buffer_push_bytes(buffer, wordbuf, second_part);
208 }
209
210 return janet_wrap_buffer(buffer);
211 }
212
213 static const JanetMethod rng_methods[] = {
214 {"uniform", cfun_rng_uniform},
215 {"int", cfun_rng_int},
216 {"buffer", cfun_rng_buffer},
217 {NULL, NULL}
218 };
219
janet_rng_get(void * p,Janet key,Janet * out)220 static int janet_rng_get(void *p, Janet key, Janet *out) {
221 (void) p;
222 if (!janet_checktype(key, JANET_KEYWORD)) return 0;
223 return janet_getmethod(janet_unwrap_keyword(key), rng_methods, out);
224 }
225
janet_rng_next(void * p,Janet key)226 static Janet janet_rng_next(void *p, Janet key) {
227 (void) p;
228 return janet_nextmethod(rng_methods, key);
229 }
230
231 /* Get a random number */
232 JANET_CORE_FN(janet_rand,
233 "(math/random)",
234 "Returns a uniformly distributed random number between 0 and 1") {
235 (void) argv;
236 janet_fixarity(argc, 0);
237 return janet_wrap_number(janet_rng_double(&janet_vm.rng));
238 }
239
240 /* Seed the random number generator */
241 JANET_CORE_FN(janet_srand,
242 "(math/seedrandom seed)",
243 "Set the seed for the random number generator. seed should be "
244 "an integer or a buffer."
245 ) {
246 janet_fixarity(argc, 1);
247 if (janet_checkint(argv[0])) {
248 uint32_t seed = (uint32_t)(janet_getinteger(argv, 0));
249 janet_rng_seed(&janet_vm.rng, seed);
250 } else {
251 JanetByteView bytes = janet_getbytes(argv, 0);
252 janet_rng_longseed(&janet_vm.rng, bytes.bytes, bytes.len);
253 }
254 return janet_wrap_nil();
255 }
256
257 #define JANET_DEFINE_MATHOP(name, fop, doc)\
258 JANET_CORE_FN(janet_##name, "(math/" #name " x)", doc) {\
259 janet_fixarity(argc, 1); \
260 double x = janet_getnumber(argv, 0); \
261 return janet_wrap_number(fop(x)); \
262 }
263
264 JANET_DEFINE_MATHOP(acos, acos, "Returns the arccosize of x.")
265 JANET_DEFINE_MATHOP(asin, asin, "Returns the arcsin of x.")
266 JANET_DEFINE_MATHOP(atan, atan, "Returns the arctangent of x.")
267 JANET_DEFINE_MATHOP(cos, cos, "Returns the cosine of x.")
268 JANET_DEFINE_MATHOP(cosh, cosh, "Returns the hyperbolic cosine of x.")
269 JANET_DEFINE_MATHOP(acosh, acosh, "Returns the hyperbolic arccosine of x.")
270 JANET_DEFINE_MATHOP(sin, sin, "Returns the sine of x.")
271 JANET_DEFINE_MATHOP(sinh, sinh, "Returns the hyperbolic sine of x.")
272 JANET_DEFINE_MATHOP(asinh, asinh, "Returns the hypberbolic arcsine of x.")
273 JANET_DEFINE_MATHOP(tan, tan, "Returns the tangent of x.")
274 JANET_DEFINE_MATHOP(tanh, tanh, "Returns the hyperbolic tangent of x.")
275 JANET_DEFINE_MATHOP(atanh, atanh, "Returns the hyperbolic arctangent of x.")
276 JANET_DEFINE_MATHOP(exp, exp, "Returns e to the power of x.")
277 JANET_DEFINE_MATHOP(exp2, exp2, "Returns 2 to the power of x.")
278 JANET_DEFINE_MATHOP(expm1, expm1, "Returns e to the power of x minus 1.")
279 JANET_DEFINE_MATHOP(log, log, "Returns the natural logarithm of x.")
280 JANET_DEFINE_MATHOP(log10, log10, "Returns the log base 10 of x.")
281 JANET_DEFINE_MATHOP(log2, log2, "Returns the log base 2 of x.")
282 JANET_DEFINE_MATHOP(sqrt, sqrt, "Returns the square root of x.")
283 JANET_DEFINE_MATHOP(cbrt, cbrt, "Returns the cube root of x.")
284 JANET_DEFINE_MATHOP(ceil, ceil, "Returns the smallest integer value number that is not less than x.")
285 JANET_DEFINE_MATHOP(fabs, fabs, "Return the absolute value of x.")
286 JANET_DEFINE_MATHOP(floor, floor, "Returns the largest integer value number that is not greater than x.")
287 JANET_DEFINE_MATHOP(trunc, trunc, "Returns the integer between x and 0 nearest to x.")
288 JANET_DEFINE_MATHOP(round, round, "Returns the integer nearest to x.")
289 JANET_DEFINE_MATHOP(gamma, tgamma, "Returns gamma(x).")
290 JANET_DEFINE_MATHOP(lgamma, lgamma, "Returns log-gamma(x).")
291 JANET_DEFINE_MATHOP(log1p, log1p, "Returns (log base e of x) + 1 more accurately than (+ (math/log x) 1)")
292 JANET_DEFINE_MATHOP(erf, erf, "Returns the error function of x.")
293 JANET_DEFINE_MATHOP(erfc, erfc, "Returns the complementary error function of x.")
294
295 #define JANET_DEFINE_MATH2OP(name, fop, signature, doc)\
296 JANET_CORE_FN(janet_##name, signature, doc) {\
297 janet_fixarity(argc, 2); \
298 double lhs = janet_getnumber(argv, 0); \
299 double rhs = janet_getnumber(argv, 1); \
300 return janet_wrap_number(fop(lhs, rhs)); \
301 }
302
303 JANET_DEFINE_MATH2OP(atan2, atan2, "(math/atan2 y x)", "Returns the arctangent of y/x. Works even when x is 0.")
304 JANET_DEFINE_MATH2OP(pow, pow, "(math/pow a x)", "Returns a to the power of x.")
305 JANET_DEFINE_MATH2OP(hypot, hypot, "(math/hypot a b)", "Returns c from the equation c^2 = a^2 + b^2.")
306 JANET_DEFINE_MATH2OP(nextafter, nextafter, "(math/next x y)", "Returns the next representable floating point vaue after x in the direction of y.")
307
308 JANET_CORE_FN(janet_not, "(not x)", "Returns the boolean inverse of x.") {
309 janet_fixarity(argc, 1);
310 return janet_wrap_boolean(!janet_truthy(argv[0]));
311 }
312
janet_gcd(double x,double y)313 static double janet_gcd(double x, double y) {
314 if (isnan(x) || isnan(y)) {
315 #ifdef NAN
316 return NAN;
317 #else
318 return 0.0 \ 0.0;
319 #endif
320 }
321 if (isinf(x) || isinf(y)) return INFINITY;
322 while (y != 0) {
323 double temp = y;
324 y = fmod(x, y);
325 x = temp;
326 }
327 return x;
328 }
329
janet_lcm(double x,double y)330 static double janet_lcm(double x, double y) {
331 return (x / janet_gcd(x, y)) * y;
332 }
333
334 JANET_CORE_FN(janet_cfun_gcd, "(math/gcd x y)",
335 "Returns the greatest common divisor between x and y.") {
336 janet_fixarity(argc, 2);
337 double x = janet_getnumber(argv, 0);
338 double y = janet_getnumber(argv, 1);
339 return janet_wrap_number(janet_gcd(x, y));
340 }
341
342 JANET_CORE_FN(janet_cfun_lcm, "(math/lcm x y)",
343 "Returns the least common multiple of x and y.") {
344 janet_fixarity(argc, 2);
345 double x = janet_getnumber(argv, 0);
346 double y = janet_getnumber(argv, 1);
347 return janet_wrap_number(janet_lcm(x, y));
348 }
349
350 /* Module entry point */
janet_lib_math(JanetTable * env)351 void janet_lib_math(JanetTable *env) {
352 JanetRegExt math_cfuns[] = {
353 JANET_CORE_REG("not", janet_not),
354 JANET_CORE_REG("math/random", janet_rand),
355 JANET_CORE_REG("math/seedrandom", janet_srand),
356 JANET_CORE_REG("math/cos", janet_cos),
357 JANET_CORE_REG("math/sin", janet_sin),
358 JANET_CORE_REG("math/tan", janet_tan),
359 JANET_CORE_REG("math/acos", janet_acos),
360 JANET_CORE_REG("math/asin", janet_asin),
361 JANET_CORE_REG("math/atan", janet_atan),
362 JANET_CORE_REG("math/exp", janet_exp),
363 JANET_CORE_REG("math/log", janet_log),
364 JANET_CORE_REG("math/log10", janet_log10),
365 JANET_CORE_REG("math/log2", janet_log2),
366 JANET_CORE_REG("math/sqrt", janet_sqrt),
367 JANET_CORE_REG("math/cbrt", janet_cbrt),
368 JANET_CORE_REG("math/floor", janet_floor),
369 JANET_CORE_REG("math/ceil", janet_ceil),
370 JANET_CORE_REG("math/pow", janet_pow),
371 JANET_CORE_REG("math/abs", janet_fabs),
372 JANET_CORE_REG("math/sinh", janet_sinh),
373 JANET_CORE_REG("math/cosh", janet_cosh),
374 JANET_CORE_REG("math/tanh", janet_tanh),
375 JANET_CORE_REG("math/atanh", janet_atanh),
376 JANET_CORE_REG("math/asinh", janet_asinh),
377 JANET_CORE_REG("math/acosh", janet_acosh),
378 JANET_CORE_REG("math/atan2", janet_atan2),
379 JANET_CORE_REG("math/rng", cfun_rng_make),
380 JANET_CORE_REG("math/rng-uniform", cfun_rng_uniform),
381 JANET_CORE_REG("math/rng-int", cfun_rng_int),
382 JANET_CORE_REG("math/rng-buffer", cfun_rng_buffer),
383 JANET_CORE_REG("math/hypot", janet_hypot),
384 JANET_CORE_REG("math/exp2", janet_exp2),
385 JANET_CORE_REG("math/log1p", janet_log1p),
386 JANET_CORE_REG("math/gamma", janet_gamma),
387 JANET_CORE_REG("math/log-gamma", janet_lgamma),
388 JANET_CORE_REG("math/erfc", janet_erfc),
389 JANET_CORE_REG("math/erf", janet_erf),
390 JANET_CORE_REG("math/expm1", janet_expm1),
391 JANET_CORE_REG("math/trunc", janet_trunc),
392 JANET_CORE_REG("math/round", janet_round),
393 JANET_CORE_REG("math/next", janet_nextafter),
394 JANET_CORE_REG("math/gcd", janet_cfun_gcd),
395 JANET_CORE_REG("math/lcm", janet_cfun_lcm),
396 JANET_REG_END
397 };
398 janet_core_cfuns_ext(env, NULL, math_cfuns);
399 janet_register_abstract_type(&janet_rng_type);
400 #ifdef JANET_BOOTSTRAP
401 JANET_CORE_DEF(env, "math/pi", janet_wrap_number(3.1415926535897931),
402 "The value pi.");
403 JANET_CORE_DEF(env, "math/e", janet_wrap_number(2.7182818284590451),
404 "The base of the natural log.");
405 JANET_CORE_DEF(env, "math/inf", janet_wrap_number(INFINITY),
406 "The number representing positive infinity");
407 JANET_CORE_DEF(env, "math/-inf", janet_wrap_number(-INFINITY),
408 "The number representing negative infinity");
409 JANET_CORE_DEF(env, "math/int32-min", janet_wrap_number(INT32_MIN),
410 "The minimum contiguous integer representable by a 32 bit signed integer");
411 JANET_CORE_DEF(env, "math/int32-max", janet_wrap_number(INT32_MAX),
412 "The maximum contiguous integer represtenable by a 32 bit signed integer");
413 JANET_CORE_DEF(env, "math/int-min", janet_wrap_number(JANET_INTMIN_DOUBLE),
414 "The minimum contiguous integer representable by a double (2^53)");
415 JANET_CORE_DEF(env, "math/int-max", janet_wrap_number(JANET_INTMAX_DOUBLE),
416 "The maximum contiguous integer represtenable by a double (-(2^53))");
417 #ifdef NAN
418 JANET_CORE_DEF(env, "math/nan", janet_wrap_number(NAN), "Not a number (IEEE-754 NaN)");
419 #else
420 JANET_CORE_DEF(env, "math/nan", janet_wrap_number(0.0 / 0.0), "Not a number (IEEE-754 NaN)");
421 #endif
422 #endif
423 }
424