Skip to content

Commit 825ec2b

Browse files
v923zv923z
authored
re-name integration constants (#714)
1 parent 66daa89 commit 825ec2b

File tree

3 files changed

+58
-52
lines changed

3 files changed

+58
-52
lines changed

‎code/scipy/integrate/integrate.c

Lines changed: 51 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -41,21 +41,21 @@
4141

4242
#if MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_DOUBLE
4343
ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-14), 0x283424dcUL, 0x3e901b2b29a4692bULL);
44-
#define MACHEPS MICROPY_FLOAT_CONST(1e-17)
44+
#define ULAB_MACHEPS MICROPY_FLOAT_CONST(1e-17)
4545
#else
4646
ULAB_DEFINE_FLOAT_CONST(etolerance, MICROPY_FLOAT_CONST(1e-8), 0x358637cfUL, 0x3e7010c6f7d42d18ULL);
47-
#define MACHEPS MICROPY_FLOAT_CONST(1e-8)
47+
#define ULAB_MACHEPS MICROPY_FLOAT_CONST(1e-8)
4848
#endif
4949

50-
#define ZERO MICROPY_FLOAT_CONST(0.0)
51-
#define POINT_TWO_FIVE MICROPY_FLOAT_CONST(0.25)
52-
#define ONE MICROPY_FLOAT_CONST(1.0)
53-
#define TWO MICROPY_FLOAT_CONST(2.0)
54-
#define FOUR MICROPY_FLOAT_CONST(4.0)
55-
#define SIX MICROPY_FLOAT_CONST(6.0)
56-
#define TEN MICROPY_FLOAT_CONST(10.0)
57-
#define FIFTEEN MICROPY_FLOAT_CONST(15.0)
58-
#define EPS_5 MICROPY_FLOAT_CONST(1e-5)
50+
#define ULAB_ZERO MICROPY_FLOAT_CONST(0.0)
51+
#define ULAB_POINT_TWO_FIVE MICROPY_FLOAT_CONST(0.25)
52+
#define ULAB_ONE MICROPY_FLOAT_CONST(1.0)
53+
#define ULAB_TWO MICROPY_FLOAT_CONST(2.0)
54+
#define ULAB_FOUR MICROPY_FLOAT_CONST(4.0)
55+
#define ULAB_SIX MICROPY_FLOAT_CONST(6.0)
56+
#define ULAB_TEN MICROPY_FLOAT_CONST(10.0)
57+
#define ULAB_FIFTEEN MICROPY_FLOAT_CONST(15.0)
58+
#define ULAB_EPSILON_5 MICROPY_FLOAT_CONST(1e-5)
5959

6060

6161
static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun, mp_float_t x, mp_obj_t *fargs, uint8_t nparams) {
@@ -68,7 +68,7 @@ static mp_float_t integrate_python_call(const mp_obj_type_t *type, mp_obj_t fun,
6868

6969
// sign helper function
7070
int sign(mp_float_t x) {
71-
if (x >= ZERO)
71+
if (x >= ULAB_ZERO)
7272
return 1;
7373
else
7474
return -1;
@@ -85,7 +85,7 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_
8585
mp_obj_t fargs[1];
8686
mp_float_t h2 = integrate_python_call(type, fun, a + d/2, fargs, 0) - integrate_python_call(type, fun, (a + d*2)*4, fargs, 0);
8787
int i = 1, j = 32; // j=32 is optimal to find r
88-
if (isfinite(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > EPS_5) { // if |h2| > 2^-16
88+
if (isfinite(h2) && MICROPY_FLOAT_C_FUN(fabs)(h2) > ULAB_EPSILON_5) { // if |h2| > 2^-16
8989
mp_float_t r, fl, fr, h, s = 0, lfl, lfr, lr = 2;
9090
do { // find max j such that fl and fr are finite
9191
j /= 2;
@@ -118,8 +118,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_
118118
if (s > eps) { // if sum of |h| > eps
119119
h = lfl - lfr; // use last fl and fr before the sign change
120120
r = lr; // use last r before the sign change
121-
if (h != ZERO) // if last diff != 0, back up r by one step
122-
r /= TWO;
121+
if (h != ULAB_ZERO) // if last diff != 0, back up r by one step
122+
r /= ULAB_TWO;
123123
if (MICROPY_FLOAT_C_FUN(fabs)(lfl) < MICROPY_FLOAT_C_FUN(fabs)(lfr))
124124
d /= r; // move d closer to the finite endpoint
125125
else
@@ -135,8 +135,8 @@ mp_float_t exp_sinh_opt_d(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_
135135
mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint16_t n, mp_float_t eps, mp_float_t *e) {
136136
const mp_obj_type_t *type = mp_obj_get_type(fun);
137137
mp_obj_t fargs[1];
138-
const mp_float_t tol = TEN*eps;
139-
mp_float_t c = ZERO, d = ONE, s, sign = ONE, v, h = TWO;
138+
const mp_float_t tol = ULAB_TEN * eps;
139+
mp_float_t c = ULAB_ZERO, d = ULAB_ONE, s, sign = ULAB_ONE, v, h = ULAB_TWO;
140140
int k = 0, mode = 0; // Tanh-Sinh = 0, Exp-Sinh = 1, Sinh-Sinh = 2
141141
if (b < a) { // swap bounds
142142
v = b;
@@ -145,8 +145,8 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
145145
sign = -1;
146146
}
147147
if (isfinite(a) && isfinite(b)) {
148-
c = (a+b)/TWO;
149-
d = (b-a)/TWO;
148+
c = (a+b) / ULAB_TWO;
149+
d = (b-a) / ULAB_TWO;
150150
v = c;
151151
}
152152
else if (isfinite(a)) {
@@ -165,20 +165,20 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
165165
}
166166
else {
167167
mode = 2; // Sinh-Sinh
168-
v = ZERO;
168+
v = ULAB_ZERO;
169169
}
170170
s = integrate_python_call(type, fun, v, fargs, 0);
171171
do {
172-
mp_float_t p = ZERO, q, fp = ZERO, fm = ZERO, t, eh;
173-
h /= TWO;
172+
mp_float_t p = ULAB_ZERO, q, fp = ULAB_ZERO, fm = ULAB_ZERO, t, eh;
173+
h /= ULAB_TWO;
174174
t = eh = MICROPY_FLOAT_C_FUN(exp)(h);
175-
if (k > ZERO)
175+
if (k > ULAB_ZERO)
176176
eh *= eh;
177177
if (mode == 0) { // Tanh-Sinh
178178
do {
179-
mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(ONE / t - t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
180-
mp_float_t r = TWO * u / (ONE + u); // = 1 - tanh(sinh(j*h))
181-
mp_float_t w = (t + ONE / t) * r / (ONE + u); // = cosh(j*h)/cosh(sinh(j*h))^2
179+
mp_float_t u = MICROPY_FLOAT_C_FUN(exp)(ULAB_ONE / t - t); // = exp(-2*sinh(j*h)) = 1/exp(sinh(j*h))^2
180+
mp_float_t r = ULAB_TWO * u / (ULAB_ONE + u); // = 1 - tanh(sinh(j*h))
181+
mp_float_t w = (t + ULAB_ONE / t) * r / (ULAB_ONE + u); // = cosh(j*h)/cosh(sinh(j*h))^2
182182
mp_float_t x = d*r;
183183
if (a+x > a) { // if too close to a then reuse previous fp
184184
mp_float_t y = integrate_python_call(type, fun, a+x, fargs, 0);
@@ -196,11 +196,11 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
196196
} while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p));
197197
}
198198
else {
199-
t /= TWO;
199+
t /= ULAB_TWO;
200200
do {
201-
mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t - POINT_TWO_FIVE / t); // = exp(sinh(j*h))
201+
mp_float_t r = MICROPY_FLOAT_C_FUN(exp)(t - ULAB_POINT_TWO_FIVE / t); // = exp(sinh(j*h))
202202
mp_float_t x, y, w = r;
203-
q = ZERO;
203+
q = ULAB_ZERO;
204204
if (mode == 1) { // Exp-Sinh
205205
x = c + d/r;
206206
if (x == c) // if x hit the finite endpoint then break
@@ -210,8 +210,8 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
210210
q += y/w;
211211
}
212212
else { // Sinh-Sinh
213-
r = (r - ONE / r) / TWO; // = sinh(sinh(j*h))
214-
w = (w + ONE / w) / TWO; // = cosh(sinh(j*h))
213+
r = (r - ULAB_ONE / r) / ULAB_TWO; // = sinh(sinh(j*h))
214+
w = (w + ULAB_ONE / w) / ULAB_TWO; // = cosh(sinh(j*h))
215215
x = c - d*r;
216216
y = integrate_python_call(type, fun, x, fargs, 0);
217217
if (isfinite(y)) // if f(x) is finite, add to local sum
@@ -221,7 +221,7 @@ mp_float_t tanhsinh(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, u
221221
y = integrate_python_call(type, fun, x, fargs, 0);
222222
if (isfinite(y)) // if f(x) is finite, add to local sum
223223
q += y*w;
224-
q *= t + POINT_TWO_FIVE / t; // q *= cosh(j*h)
224+
q *= t + ULAB_POINT_TWO_FIVE / t; // q *= cosh(j*h)
225225
p += q;
226226
t *= eh;
227227
} while (MICROPY_FLOAT_C_FUN(fabs)(q) > eps*MICROPY_FLOAT_C_FUN(fabs)(p));
@@ -319,12 +319,12 @@ mp_float_t qromb(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, uint
319319
for (i = 1; i < n; ++i) {
320320
unsigned long long k = 1UL << i;
321321
unsigned long long s = 1;
322-
mp_float_t sum = ZERO;
322+
mp_float_t sum = ULAB_ZERO;
323323
mp_float_t *Rt;
324-
h /= TWO;
324+
h /= ULAB_TWO;
325325
for (j = 1; j < k; j += 2)
326326
sum += integrate_python_call(type, fun, a+j*h, fargs, 0);
327-
Ru[0] = h*sum + Ro[0] / TWO;
327+
Ru[0] = h*sum + Ro[0] / ULAB_TWO;
328328
for (j = 1; j <= i; ++j) {
329329
s <<= 2;
330330
Ru[j] = (s*Ru[j-1] - Ro[j-1])/(s-1);
@@ -408,17 +408,17 @@ mp_float_t as(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, mp_floa
408408
mp_float_t fb, mp_float_t v, mp_float_t eps, int n, mp_float_t t) {
409409
const mp_obj_type_t *type = mp_obj_get_type(fun);
410410
mp_obj_t fargs[1];
411-
mp_float_t h = (b-a) / TWO;
412-
mp_float_t f1 = integrate_python_call(type, fun, a + h / TWO, fargs, 0);
413-
mp_float_t f2 = integrate_python_call(type, fun, b - h / TWO, fargs, 0);
414-
mp_float_t sl = h*(fa + FOUR * f1 + fm) / SIX;
415-
mp_float_t sr = h*(fm + FOUR * f2 + fb) / SIX;
411+
mp_float_t h = (b-a) / ULAB_TWO;
412+
mp_float_t f1 = integrate_python_call(type, fun, a + h / ULAB_TWO, fargs, 0);
413+
mp_float_t f2 = integrate_python_call(type, fun, b - h / ULAB_TWO, fargs, 0);
414+
mp_float_t sl = h*(fa + ULAB_FOUR * f1 + fm) / ULAB_SIX;
415+
mp_float_t sr = h*(fm + ULAB_FOUR * f2 + fb) / ULAB_SIX;
416416
mp_float_t s = sl+sr;
417-
mp_float_t d = (s-v) / FIFTEEN;
417+
mp_float_t d = (s-v) / ULAB_FIFTEEN;
418418
mp_float_t m = a+h;
419419
if (n <= 0 || MICROPY_FLOAT_C_FUN(fabs)(d) < eps)
420420
return t + s + d; // note: fabs(d) can be used as error estimate
421-
eps /= TWO;
421+
eps /= ULAB_TWO;
422422
--n;
423423
t = as(fun, a, m, fa, f1, fm, sl, eps, n, t);
424424
return as(fun, m, b, fm, f2, fb, sr, eps, n, t);
@@ -430,7 +430,7 @@ mp_float_t qasi(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n
430430
mp_float_t fa = integrate_python_call(type, fun, a, fargs, 0);
431431
mp_float_t fm = integrate_python_call(type, fun, (a+b)/2, fargs, 0);
432432
mp_float_t fb = integrate_python_call(type, fun, b, fargs, 0);
433-
mp_float_t v = (fa + FOUR * fm + fb) * (b-a) / SIX;
433+
mp_float_t v = (fa + ULAB_FOUR * fm + fb) * (b-a) / ULAB_SIX;
434434
return as(fun, a, b, fa, fm, fb, v, eps, n, 0);
435435
}
436436

@@ -562,8 +562,8 @@ mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_floa
562562
};
563563
const mp_obj_type_t *type = mp_obj_get_type(fun);
564564
mp_obj_t fargs[1];
565-
mp_float_t p = ZERO; // kronrod quadrature sum
566-
mp_float_t q = ZERO; // gauss quadrature sum
565+
mp_float_t p = ULAB_ZERO; // kronrod quadrature sum
566+
mp_float_t q = ULAB_ZERO; // gauss quadrature sum
567567
mp_float_t fp, fm;
568568
mp_float_t e;
569569
int i;
@@ -581,24 +581,24 @@ mp_float_t gk(mp_float_t (*fun)(mp_float_t), mp_float_t c, mp_float_t d, mp_floa
581581
p += (fp + fm) * weights[i];
582582
}
583583
*err = MICROPY_FLOAT_C_FUN(fabs)(p - q);
584-
e = MICROPY_FLOAT_C_FUN(fabs)(2 * p * MACHEPS); // optional, to take 1e-17 MachEps prec. into account
584+
e = MICROPY_FLOAT_C_FUN(fabs)(2 * p * ULAB_MACHEPS); // optional, to take 1e-17 MachEps prec. into account
585585
if (*err < e)
586586
*err = e;
587587
return p;
588588
}
589589

590590
mp_float_t qakro(mp_float_t (*fun)(mp_float_t), mp_float_t a, mp_float_t b, int n, mp_float_t tol, mp_float_t eps, mp_float_t *err) {
591-
mp_float_t c = (a+b) / TWO;
592-
mp_float_t d = (b-a) / TWO;
591+
mp_float_t c = (a+b) / ULAB_TWO;
592+
mp_float_t d = (b-a) / ULAB_TWO;
593593
mp_float_t e;
594594
mp_float_t r = gk(fun, c, d, &e);
595595
mp_float_t s = d*r;
596596
mp_float_t t = MICROPY_FLOAT_C_FUN(fabs)(s*tol);
597-
if (tol == ZERO)
597+
if (tol == ULAB_ZERO)
598598
tol = t;
599599
if (n > 0 && t < e && tol < e) {
600-
s = qakro(fun, a, c, n-1, t / TWO, eps, err);
601-
s += qakro(fun, c, b, n-1, t / TWO, eps, &e);
600+
s = qakro(fun, a, c, n-1, t / ULAB_TWO, eps, err);
601+
s += qakro(fun, c, b, n-1, t / ULAB_TWO, eps, &e);
602602
*err += e;
603603
return s;
604604
}

‎code/ulab.c

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@
3333
#include "user/user.h"
3434
#include "utils/utils.h"
3535

36-
#define ULAB_VERSION 6.7.3
36+
#define ULAB_VERSION 6.7.4
3737
#define xstr(s) str(s)
3838
#define str(s) #s
3939

‎docs/ulab-change-log.md

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,9 @@
1+
Sun, 16 Mar 2025
2+
3+
version 6.7.4
4+
5+
re-name integration constants to avoid name clash with EPS ports
6+
17
Sun, 26 Jan 2025
28

39
version 6.7.3

0 commit comments

Comments
 (0)