Skip to content

Commit 96ee075

Browse files
SwayamInSyncSwayamInSync
committed
adding ufuncs
1 parent 58fbc6f commit 96ee075

File tree

6 files changed

+384
-1
lines changed

6 files changed

+384
-1
lines changed

‎quaddtype/meson.build‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ py = py_mod.find_installation()
55

66
c = meson.get_compiler('c')
77

8-
sleef_dep = c.find_library('sleef')
8+
sleef_dep = c.find_library('sleef', dirs:['/usr/local/lib'])
99
sleefquad_dep = c.find_library('sleefquad')
1010

1111
incdir_numpy = run_command(py,

‎quaddtype/quaddtype/src/ops.hpp‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
#include <sleef.h>
2+
#include <sleefquad.h>
3+
4+
typedef int (*unary_op_def)(Sleef_quad *, Sleef_quad *);
5+
6+
static inline int
7+
quad_negative(Sleef_quad *op, Sleef_quad *out)
8+
{
9+
*out = Sleef_negq1(*op);
10+
return 0;
11+
}
12+
13+
static inline int
14+
quad_absolute(Sleef_quad *op, Sleef_quad *out)
15+
{
16+
*out = Sleef_fabsq1(*op);
17+
return 0;
18+
}
19+
20+
// binary ops
21+
typedef int (*binop_def)(Sleef_quad *, Sleef_quad *, Sleef_quad *);
22+
23+
static inline int quad_add(Sleef_quad *out, Sleef_quad *in1, Sleef_quad *in2)
24+
{
25+
*out = Sleef_addq1_u05(*in1, *in2);
26+
return 0;
27+
}
28+
29+
static inline int quad_sub(Sleef_quad *out, Sleef_quad *in1, Sleef_quad *in2)
30+
{
31+
*out = Sleef_subq1_u05(*in1, *in2);
32+
return 0;
33+
}

‎quaddtype/quaddtype/src/umath.cpp‎

Lines changed: 224 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,224 @@
1+
#define PY_ARRAY_UNIQUE_SYMBOL QuadPrecType_ARRAY_API
2+
#define PY_UFUNC_UNIQUE_SYMBOL QuadPrecType_UFUNC_API
3+
#define NPY_NO_DEPRECATED_API NPY_2_0_API_VERSION
4+
#define NPY_TARGET_VERSION NPY_2_0_API_VERSION
5+
#define NO_IMPORT_ARRAY
6+
#define NO_IMPORT_UFUNC
7+
8+
extern "C" {
9+
#include <Python.h>
10+
11+
#include "numpy/arrayobject.h"
12+
#include "numpy/ndarraytypes.h"
13+
#include "numpy/ufuncobject.h"
14+
15+
#include "numpy/dtype_api.h"
16+
}
17+
18+
#include "scalar.h"
19+
#include "dtype.h"
20+
#include "umath.h"
21+
#include "ops.hpp"
22+
23+
template <unary_op_def unary_op>
24+
int quad_generic_unary_op_strided_loop(PyArrayMethod_Context *context,
25+
char *const data[], npy_intp const dimensions[],
26+
npy_intp const strides[], NpyAuxData *auxdata)
27+
{
28+
npy_intp N = dimensions[0];
29+
char *in_ptr = data[0];
30+
char *out_ptr = data[1];
31+
npy_intp in_stride = strides[0];
32+
npy_intp out_stride = strides[1];
33+
34+
while (N--)
35+
{
36+
unary_op((Sleef_quad *)in_ptr, (Sleef_quad *)out_ptr);
37+
in_ptr += in_stride;
38+
out_ptr += out_stride;
39+
}
40+
return 0;
41+
}
42+
43+
static NPY_CASTING
44+
quad_unary_op_resolve_descriptors(PyObject *self,
45+
PyArray_DTypeMeta *dtypes[], QuadPrecDTypeObject *given_descrs[],
46+
QuadPrecDTypeObject *loop_descrs[], npy_intp *unused)
47+
{
48+
Py_INCREF(given_descrs[0]);
49+
loop_descrs[0] = given_descrs[0];
50+
51+
if (given_descrs[1] == NULL) {
52+
Py_INCREF(given_descrs[0]);
53+
loop_descrs[1] = given_descrs[0];
54+
return NPY_NO_CASTING;
55+
}
56+
Py_INCREF(given_descrs[1]);
57+
loop_descrs[1] = given_descrs[1];
58+
59+
return NPY_NO_CASTING; // Quad precision is always the same precision
60+
}
61+
62+
template <unary_op_def unary_op>
63+
int create_quad_unary_ufunc(PyObject *numpy, const char *ufunc_name)
64+
{
65+
PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name);
66+
if (ufunc == NULL) {
67+
return -1;
68+
}
69+
70+
PyArray_DTypeMeta *dtypes[2] = {
71+
&QuadPrecDType, &QuadPrecDType};
72+
73+
PyType_Slot slots[] = {
74+
{NPY_METH_resolve_descriptors, (void *)&quad_unary_op_resolve_descriptors},
75+
{NPY_METH_strided_loop, (void *)&quad_generic_unary_op_strided_loop<unary_op>},
76+
{0, NULL}
77+
};
78+
79+
PyArrayMethod_Spec Spec = {
80+
.name = "quad_unary_op",
81+
.nin = 1,
82+
.nout = 1,
83+
.casting = NPY_NO_CASTING,
84+
.flags = (NPY_ARRAYMETHOD_FLAGS)0,
85+
.dtypes = dtypes,
86+
.slots = slots,
87+
};
88+
89+
if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) {
90+
return -1;
91+
}
92+
93+
return 0;
94+
}
95+
96+
int init_quad_unary_ops(PyObject *numpy)
97+
{
98+
if (create_quad_unary_ufunc<quad_negative>(numpy, "negative") < 0) {
99+
return -1;
100+
}
101+
if (create_quad_unary_ufunc<quad_absolute>(numpy, "absolute") < 0) {
102+
return -1;
103+
}
104+
return 0;
105+
}
106+
107+
// Binary ufuncs
108+
109+
template <binop_def binop>
110+
int quad_generic_binop_strided_loop(PyArrayMethod_Context *context,
111+
char *const data[], npy_intp const dimensions[],
112+
npy_intp const strides[], NpyAuxData *auxdata)
113+
{
114+
npy_intp N = dimensions[0];
115+
char *in1_ptr = data[0], *in2_ptr = data[1];
116+
char *out_ptr = data[2];
117+
npy_intp in1_stride = strides[0];
118+
npy_intp in2_stride = strides[1];
119+
npy_intp out_stride = strides[2];
120+
121+
while (N--) {
122+
binop((Sleef_quad *)out_ptr, (Sleef_quad *)in1_ptr, (Sleef_quad *)in2_ptr);
123+
124+
in1_ptr += in1_stride;
125+
in2_ptr += in2_stride;
126+
out_ptr += out_stride;
127+
}
128+
return 0;
129+
}
130+
131+
static NPY_CASTING
132+
quad_binary_op_resolve_descriptors(PyObject *self,
133+
PyArray_DTypeMeta *dtypes[], QuadPrecDTypeObject *given_descrs[],
134+
QuadPrecDTypeObject *loop_descrs[], npy_intp *unused)
135+
{
136+
Py_INCREF(given_descrs[0]);
137+
loop_descrs[0] = given_descrs[0];
138+
Py_INCREF(given_descrs[1]);
139+
loop_descrs[1] = given_descrs[1];
140+
141+
if (given_descrs[2] == NULL) {
142+
Py_INCREF(given_descrs[0]);
143+
loop_descrs[2] = given_descrs[0];
144+
}
145+
else {
146+
Py_INCREF(given_descrs[2]);
147+
loop_descrs[2] = given_descrs[2];
148+
}
149+
150+
return NPY_NO_CASTING; // Quad precision is always the same precision
151+
}
152+
153+
// todo: skipping the promoter for now, since same type operation will be requried
154+
155+
template <binop_def binop>
156+
int create_quad_binary_ufunc(PyObject *numpy, const char *ufunc_name)
157+
{
158+
PyObject *ufunc = PyObject_GetAttrString(numpy, ufunc_name);
159+
if (ufunc == NULL) {
160+
return -1;
161+
}
162+
163+
PyArray_DTypeMeta *dtypes[3] = {
164+
&QuadPrecDType, &QuadPrecDType, &QuadPrecDType};
165+
166+
PyType_Slot slots[] = {
167+
{NPY_METH_resolve_descriptors,
168+
(void *)&quad_binary_op_resolve_descriptors},
169+
{NPY_METH_strided_loop,
170+
(void *)&quad_generic_binop_strided_loop<binop>},
171+
{0, NULL}
172+
};
173+
174+
PyArrayMethod_Spec Spec = {
175+
.name = "quad_binop",
176+
.nin = 2,
177+
.nout = 1,
178+
.casting = NPY_NO_CASTING,
179+
.flags = (NPY_ARRAYMETHOD_FLAGS)0,
180+
.dtypes = dtypes,
181+
.slots = slots,
182+
};
183+
184+
if (PyUFunc_AddLoopFromSpec(ufunc, &Spec) < 0) {
185+
return -1;
186+
}
187+
188+
return 0;
189+
}
190+
191+
int init_quad_binary_ops(PyObject *numpy)
192+
{
193+
if (create_quad_binary_ufunc<quad_add>(numpy, "add") < 0) {
194+
return -1;
195+
}
196+
if (create_quad_binary_ufunc<quad_sub>(numpy, "subtract") < 0) {
197+
return -1;
198+
}
199+
200+
return 0;
201+
}
202+
203+
int init_quad_umath(void)
204+
{
205+
PyObject * numpy = PyImport_ImportModule("numpy");
206+
if (!numpy)
207+
return -1;
208+
209+
if (init_quad_unary_ops(numpy) < 0) {
210+
goto err;
211+
}
212+
213+
if (init_quad_binary_ops(numpy) < 0) {
214+
goto err;
215+
}
216+
217+
Py_DECREF(numpy);
218+
return 0;
219+
220+
err:
221+
Py_DECREF(numpy);
222+
return -1;
223+
224+
}

‎quaddtype/quaddtype/src/umath.h‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
#ifndef _QUADDTYPE_UMATH_H
2+
#define _QUADDTYPE_UMATH_H
3+
4+
#ifdef __cplusplus
5+
extern "C" {
6+
#endif
7+
8+
int
9+
init_quad_umath(void);
10+
11+
#ifdef __cplusplus
12+
}
13+
#endif
14+
15+
#endif

‎quaddtype/tests/test_arr.py‎ renamed to ‎quaddtype/tests/test_quaddtype.py‎

File renamed without changes.

‎temp.py‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
import numpy as np
2+
from quaddtype import QuadPrecDType, QuadPrecision
3+
import matplotlib.pyplot as plt
4+
5+
def get_color(t, interior_t):
6+
epsilon = QuadPrecision("1e-10")
7+
8+
if abs(t - QuadPrecision(1.0)) < epsilon:
9+
value = int(255 * float(interior_t))
10+
return np.array([value, value, value], dtype=np.uint8)
11+
12+
t = np.power(t, 0.5)
13+
t = np.mod(t * 20, 1.0)
14+
15+
if t < 0.16:
16+
return np.array([0, int(255 * (t / 0.16)), int(128 + 127 * (t / 0.16))], dtype=np.uint8)
17+
elif t < 0.33:
18+
return np.array([0, 255, int(255 * (1 - (t - 0.16) / 0.17))], dtype=np.uint8)
19+
elif t < 0.5:
20+
return np.array([int(255 * ((t - 0.33) / 0.17)), 255, 0], dtype=np.uint8)
21+
elif t < 0.66:
22+
return np.array([255, int(255 * (1 - (t - 0.5) / 0.16)), 0], dtype=np.uint8)
23+
elif t < 0.83:
24+
return np.array([255, 0, int(255 * ((t - 0.66) / 0.17))], dtype=np.uint8)
25+
else:
26+
return np.array([int(255 * (1 - (t - 0.83) / 0.17)), 0, int(128 * ((t - 0.83) / 0.17))], dtype=np.uint8)
27+
28+
def iterate_and_compute_derivatives(c, max_iter):
29+
z = 0
30+
dz = 1
31+
dc = 0
32+
dzdz = 0
33+
34+
for _ in range(max_iter):
35+
dzdz = 2 * (z * dzdz + dz * dz)
36+
dz = 2 * z * dz + dc
37+
z = z * z + c
38+
dc = 1
39+
40+
return z, dz, dc, dzdz
41+
42+
def estimate_interior_distance(c, max_iter):
43+
z, dz, dc, dzdz = iterate_and_compute_derivatives(c, max_iter)
44+
45+
dz_abs_sq = np.abs(dz) ** 2
46+
numerator = 1 - dz_abs_sq
47+
48+
denominator = np.abs(dc * dz + dzdz * z * dc)
49+
50+
return numerator / denominator
51+
52+
def mandelbrot(c, max_iter, radius2):
53+
z = 0
54+
for i in range(max_iter):
55+
z = z * z + c
56+
if np.abs(z) ** 2 > radius2:
57+
log_zn = np.log(np.abs(z))
58+
nu = np.log(log_zn / np.log(2)) / np.log(2)
59+
return i + 1 - nu, z
60+
return max_iter, z
61+
62+
def mandelbrot_set(width, height, max_iter, center_r, center_i, zoom):
63+
radius = 2.0
64+
radius2 = radius * radius
65+
zoom_q = 1 / zoom
66+
67+
x = np.linspace(center_r - radius / zoom, center_r + radius / zoom, width)
68+
y = np.linspace(center_i - radius / zoom, center_i + radius / zoom, height)
69+
c = x[np.newaxis, :] + 1j * y[:, np.newaxis]
70+
71+
smooth_iter, final_z = np.frompyfunc(lambda c: mandelbrot(c, max_iter, radius2), 1, 2)(c)
72+
smooth_iter = smooth_iter.astype(np.float64)
73+
final_z = final_z.astype(np.complex128)
74+
75+
img = np.zeros((height, width, 3), dtype=np.uint8)
76+
77+
interior_mask = smooth_iter == max_iter
78+
interior_c = c[interior_mask]
79+
interior_distance = np.frompyfunc(lambda c: estimate_interior_distance(c, max_iter), 1, 1)(interior_c)
80+
interior_distance = interior_distance.astype(np.float64)
81+
interior_t = interior_distance - np.floor(interior_distance)
82+
83+
exterior_mask = ~interior_mask
84+
t = smooth_iter[exterior_mask] / max_iter
85+
86+
interior_colors = np.array(list(map(lambda t: get_color(1.0, t), interior_t)))
87+
exterior_colors = np.array(list(map(lambda t: get_color(t, 0.0), t)))
88+
89+
img[interior_mask] = interior_colors
90+
img[exterior_mask] = exterior_colors
91+
92+
return img
93+
94+
def plot_mandelbrot(width, height, max_iter, center_r, center_i, zoom):
95+
img_array = mandelbrot_set(width, height, max_iter, center_r, center_i, zoom)
96+
97+
plt.figure(figsize=(10, 10))
98+
plt.imshow(img_array)
99+
plt.axis('off')
100+
plt.title(f'Mandelbrot Set (zoom: {zoom}, center: {center_r} + {center_i}i, iterations: {max_iter}, dtype: numpy.float64)')
101+
plt.show()
102+
103+
if __name__ == "__main__":
104+
width = 800
105+
height = 800
106+
max_iter = 1000
107+
center_r = -0.75
108+
center_i = 0.0
109+
zoom = 1.0
110+
111+
plot_mandelbrot(width, height, max_iter, center_r, center_i, zoom)

0 commit comments

Comments
 (0)