1 #include "ml_math.h"
2 #include "ml_macros.h"
3 #include <math.h>
4 #include <float.h>
5 #ifdef ML_COMPLEX
6 #include <complex.h>
7 #endif
8
9 #define MATH_REAL(NAME, CNAME) \
10 ML_METHOD_DECL(NAME ## Method, NULL); \
11 \
12 ML_METHOD(NAME ## Method, MLRealT) { \
13 return ml_real(CNAME(ml_real_value(Args[0]))); \
14 }
15
16 #ifdef ML_COMPLEX
17
18 #define MATH_NUMBER(NAME, CNAME) \
19 ML_METHOD_DECL(NAME ## Method, NULL); \
20 \
21 ML_METHOD(NAME ## Method, MLRealT) { \
22 return ml_real(CNAME(ml_real_value(Args[0]))); \
23 } \
24 \
25 ML_METHOD(NAME ## Method, MLComplexT) { \
26 complex double Result = c ## CNAME(ml_complex_value(Args[0])); \
27 if (fabs(cimag(Result)) <= DBL_EPSILON) { \
28 return ml_real(creal(Result)); \
29 } else { \
30 return ml_complex(Result); \
31 } \
32 }
33
34 #else
35
36 #define MATH_NUMBER(NAME, CNAME) \
37 ML_METHOD_DECL(NAME ## Method, NULL); \
38 \
39 ML_METHOD(NAME ## Method, MLRealT) { \
40 return ml_real(CNAME(ml_real_value(Args[0]))); \
41 }
42
43 #endif
44
45 #define MATH_REAL_REAL(NAME, CNAME) \
46 ML_METHOD_DECL(NAME ## Method, NULL); \
47 \
48 ML_METHOD(NAME ## Method, MLRealT, MLRealT) { \
49 return ml_real(CNAME(ml_real_value(Args[0]), ml_real_value(Args[1]))); \
50 }
51
52 ML_METHOD("%", MLRealT, MLRealT) {
53 return ml_real(fmod(ml_real_value(Args[0]), ml_real_value(Args[1])));
54 }
55
56 ML_METHOD("^", MLIntegerT, MLIntegerT) {
57 int64_t Base = ml_integer_value_fast(Args[0]);
58 int64_t Exponent = ml_integer_value_fast(Args[1]);
59 if (Exponent >= 0) {
60 int64_t N = 1;
61 while (Exponent) {
62 if (Exponent & 1) N *= Base;
63 Base *= Base;
64 Exponent >>= 1;
65 }
66 return ml_integer(N);
67 } else {
68 return ml_real(pow(Base, Exponent));
69 }
70 }
71
72 ML_METHOD("^", MLRealT, MLIntegerT) {
73 return ml_real(pow(ml_double_value_fast(Args[0]), ml_integer_value_fast(Args[1])));
74 }
75
76 ML_METHOD("^", MLRealT, MLRealT) {
77 return ml_real(pow(ml_double_value_fast(Args[0]), ml_double_value_fast(Args[1])));
78 }
79
80 #ifdef ML_COMPLEX
81
82 ML_METHOD("^", MLComplexT, MLIntegerT) {
83 complex double Base = ml_complex_value_fast(Args[0]);
84 int64_t Power = ml_integer_value_fast(Args[1]);
85 if (Power == 0) return ml_real(0);
86 complex double Result;
87 if (Power > 0 && Power < 10) {
88 Result = Base;
89 while (--Power > 0) Result *= Base;
90 } else {
91 Result = cpow(Base, Power);
92 }
93 if (fabs(cimag(Result)) <= DBL_EPSILON) {
94 return ml_real(creal(Result));
95 } else {
96 return ml_complex(Result);
97 }
98 }
99
100 ML_METHOD("^", MLComplexT, MLNumberT) {
101 complex double V = cpow(ml_complex_value_fast(Args[0]), ml_complex_value(Args[1]));
102 if (fabs(cimag(V)) <= DBL_EPSILON) {
103 return ml_real(creal(V));
104 } else {
105 return ml_complex(V);
106 }
107 }
108
109 ML_METHOD("^", MLNumberT, MLComplexT) {
110 complex double V = cpow(ml_complex_value(Args[0]), ml_complex_value_fast(Args[1]));
111 if (fabs(cimag(V)) < DBL_EPSILON) {
112 return ml_real(creal(V));
113 } else {
114 return ml_complex(V);
115 }
116 }
117
118 #endif
119
120 ML_METHOD("!", MLIntegerT) {
121 int N = ml_integer_value(Args[0]);
122 if (N > 20) return ml_error("RangeError", "Factorials over 20 are not supported yet");
123 int64_t F = N;
124 while (--N > 1) F *= N;
125 return ml_integer(F);
126 }
127
128 #define MIN(X, Y) (X < Y) ? X : Y
129
130 ML_METHOD("!", MLIntegerT, MLIntegerT) {
131 int N = ml_integer_value(Args[0]);
132 int K = ml_integer_value(Args[1]);
133 if (K > 20) return ml_error("RangeError", "Factorials over 20 are not supported yet");
134 int64_t C = 1;
135 if (K > N - K) K = N - K;
136 for (int I = 0; I < K; ++I) {
137 C *= (N - I);
138 C /= (I + 1);
139 }
140 return ml_integer(C);
141 }
142
143 MATH_NUMBER(Acos, acos);
144 MATH_NUMBER(Asin, asin);
145 MATH_NUMBER(Atan, atan);
ML_METHOD(AtanMethod,MLRealT,MLRealT)146 ML_METHOD(AtanMethod, MLRealT, MLRealT) {
147 //@atan
148 //>number
149 return ml_real(atan2(ml_real_value(Args[0]), ml_real_value(Args[1])));
150 }
151 MATH_REAL(Ceil, ceil);
152 MATH_NUMBER(Cos, cos);
153 MATH_NUMBER(Cosh, cosh);
154 MATH_NUMBER(Exp, exp);
155 MATH_REAL(Abs, fabs);
156 MATH_REAL(Floor, floor);
157 MATH_NUMBER(Log, log);
158 MATH_NUMBER(Log10, log10);
159 MATH_NUMBER(Sin, sin);
160 MATH_NUMBER(Sinh, sinh);
161 MATH_NUMBER(Sqrt, sqrt);
ML_METHOD(SqrtMethod,MLIntegerT)162 ML_METHOD(SqrtMethod, MLIntegerT) {
163 //@sqrt
164 //>number
165 int64_t N = ml_integer_value(Args[0]);
166 if (N < 0) return ml_real(-NAN);
167 if (N <= 1) return Args[0];
168 int64_t X = N >> 1;
169 for (;;) {
170 int64_t X1 = (X + N / X) >> 1;
171 if (X1 >= X) break;
172 X = X1;
173 }
174 if (X * X == N) return ml_integer(X);
175 return ml_real(sqrt(N));
176 }
177 MATH_NUMBER(Tan, tan);
178 MATH_NUMBER(Tanh, tanh);
179 MATH_REAL(Erf, erf);
180 MATH_REAL(Erfc, erfc);
181 MATH_REAL_REAL(Hypot, hypot);
182 MATH_REAL(Gamma, lgamma);
183 MATH_NUMBER(Acosh, acosh);
184 MATH_NUMBER(Asinh, asinh);
185 MATH_NUMBER(Atanh, atanh);
186 MATH_REAL(Cbrt, cbrt);
187 MATH_REAL(Expm1, expm1);
188 MATH_REAL(Log1p, log1p);
189 MATH_REAL_REAL(Rem, remainder);
190 MATH_REAL(Round, round);
191
ML_FUNCTION(IntegerRandom)192 ML_FUNCTION(IntegerRandom) {
193 //@integer::random
194 //<Min?:number
195 //<Max?:number
196 //>integer
197 // Returns a random integer between :mini:`Min` and :mini:`Max` (where :mini:`Max <= 2³² - 1`.
198 // If omitted, :mini:`Min` defaults to :mini:`0` and :mini:`Max` defaults to :mini:`2³² - 1`.
199 if (Count == 2) {
200 ML_CHECK_ARG_TYPE(0, MLRealT);
201 ML_CHECK_ARG_TYPE(1, MLRealT);
202 int Base = ml_integer_value(Args[0]);
203 int Limit = ml_integer_value(Args[1]) + 1 - Base;
204 if (Limit <= 0) return Args[0];
205 int Divisor = RAND_MAX / Limit;
206 int Random;
207 do Random = random() / Divisor; while (Random > Limit);
208 return ml_integer(Base + Random);
209 } else if (Count == 1) {
210 ML_CHECK_ARG_TYPE(0, MLRealT);
211 int Limit = ml_integer_value(Args[0]);
212 if (Limit <= 0) return Args[0];
213 int Divisor = RAND_MAX / Limit;
214 int Random;
215 do Random = random() / Divisor; while (Random > Limit);
216 return ml_integer(Random + 1);
217 } else {
218 return ml_integer(random());
219 }
220 }
221
ML_FUNCTION(IntegerRandomPermutation)222 ML_FUNCTION(IntegerRandomPermutation) {
223 //@integer::random_permutation
224 //<Max:number
225 ML_CHECK_ARG_TYPE(0, MLIntegerT);
226 int Limit = ml_integer_value(Args[0]);
227 if (Limit <= 0) return ml_error("ValueError", "Permutation requires positive size");
228 ml_value_t *Permutation = ml_list();
229 ml_list_put(Permutation, ml_integer(1));
230 for (int I = 2; I <= Limit; ++I) {
231 int Divisor = RAND_MAX / I, J;
232 do J = random() / Divisor; while (J > I);
233 ++J;
234 if (J == I) {
235 ml_list_put(Permutation, ml_integer(I));
236 } else {
237 ml_value_t *Old = ml_list_get(Permutation, J);
238 ml_list_set(Permutation, J, ml_integer(I));
239 ml_list_put(Permutation, Old);
240 }
241 }
242 return Permutation;
243 }
244
ML_FUNCTION(IntegerRandomCycle)245 ML_FUNCTION(IntegerRandomCycle) {
246 //@integer::random_cycle
247 //<Max:number
248 ML_CHECK_ARG_TYPE(0, MLIntegerT);
249 int Limit = ml_integer_value(Args[0]);
250 if (Limit <= 0) return ml_error("ValueError", "Permutation requires positive size");
251 ml_value_t *Permutation = ml_list();
252 ml_list_put(Permutation, ml_integer(1));
253 if (Limit == 1) return Permutation;
254 ml_list_push(Permutation, ml_integer(2));
255 for (int I = 2; I < Limit; ++I) {
256 int Divisor = RAND_MAX / I, J;
257 do J = random() / Divisor; while (J > I);
258 ++J;
259 ml_value_t *Old = ml_list_get(Permutation, J);
260 ml_list_set(Permutation, J, ml_integer(I + 1));
261 ml_list_put(Permutation, Old);
262 }
263 return Permutation;
264 }
265
ML_FUNCTION(RealRandom)266 ML_FUNCTION(RealRandom) {
267 //@real::random
268 //<Min?:number
269 //<Max?:number
270 //>real
271 // Returns a random real between :mini:`Min` and :mini:`Max`.
272 // If omitted, :mini:`Min` defaults to :mini:`0` and :mini:`Max` defaults to :mini:`1`.
273 if (Count == 2) {
274 ML_CHECK_ARG_TYPE(0, MLRealT);
275 ML_CHECK_ARG_TYPE(1, MLRealT);
276 double Base = ml_real_value(Args[0]);
277 double Limit = ml_real_value(Args[1]) - Base;
278 if (Limit <= 0) return Args[0];
279 double Scale = Limit / (double)RAND_MAX;
280 return ml_real(Base + random() * Scale);
281 } else if (Count == 1) {
282 double Limit = ml_real_value(Args[0]);
283 if (Limit <= 0) return Args[0];
284 double Scale = Limit / (double)RAND_MAX;
285 return ml_real(random() * Scale);
286 } else {
287 return ml_real(random() / (double)RAND_MAX);
288 }
289 }
290
ml_math_init(stringmap_t * Globals)291 void ml_math_init(stringmap_t *Globals) {
292 srandom(time(NULL));
293 #include "ml_math_init.c"
294 stringmap_insert(MLIntegerT->Exports, "random", IntegerRandom);
295 stringmap_insert(MLIntegerT->Exports, "permutation", IntegerRandomPermutation);
296 stringmap_insert(MLIntegerT->Exports, "cycle", IntegerRandomCycle);
297 stringmap_insert(MLRealT->Exports, "random", RealRandom);
298 if (Globals) {
299 stringmap_insert(Globals, "math", ml_module("math",
300 "acos", AcosMethod,
301 "asin", AsinMethod,
302 "atan", AtanMethod,
303 "ceil", CeilMethod,
304 "cos", CosMethod,
305 "cosh", CoshMethod,
306 "exp", ExpMethod,
307 "abs", AbsMethod,
308 "floor", FloorMethod,
309 "log", LogMethod,
310 "log10", Log10Method,
311 "sin", SinMethod,
312 "sinh", SinhMethod,
313 "sqrt", SqrtMethod,
314 "√", SqrtMethod,
315 "tan", TanMethod,
316 "tanh", TanhMethod,
317 "erf", ErfMethod,
318 "erfc", ErfcMethod,
319 "hypot", HypotMethod,
320 "gamma", GammaMethod,
321 "acosh", AcoshMethod,
322 "asinh", AsinhMethod,
323 "atanh", AtanhMethod,
324 "cbrt", CbrtMethod,
325 "∛", CbrtMethod,
326 "expm1", Expm1Method,
327 "log1p", Log1pMethod,
328 "rem", RemMethod,
329 "round", RoundMethod,
330 "pi", ml_real(M_PI),
331 "π", ml_real(M_PI),
332 "e", ml_real(M_E),
333 "ℯ", ml_real(M_E),
334 NULL));
335 }
336 }
337