Skip to content

Commit 04480ce

Browse files
SwayamInSyncSwayamInSync
committed
np.dot WIP
1 parent 7b43f3b commit 04480ce

File tree

3 files changed

+308
-2
lines changed

3 files changed

+308
-2
lines changed

‎quaddtype/numpy_quaddtype/src/scalar.c‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -226,12 +226,12 @@ PyTypeObject QuadPrecision_Type = {
226226
.tp_repr = (reprfunc)QuadPrecision_repr_dragon4,
227227
.tp_str = (reprfunc)QuadPrecision_str_dragon4,
228228
.tp_as_number = &quad_as_scalar,
229-
.tp_richcompare = (richcmpfunc)quad_richcompare
230-
229+
.tp_richcompare = (richcmpfunc)quad_richcompare,
231230
};
232231

233232
int
234233
init_quadprecision_scalar(void)
235234
{
235+
QuadPrecision_Type.tp_base = &PyFloat_Type; // this is not working (subclassing to np.floating)
236236
return PyType_Ready(&QuadPrecision_Type);
237237
}

‎quaddtype/numpy_quaddtype/src/scalar_ops.cpp‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -207,6 +207,27 @@ quad_richcompare(QuadPrecisionObject *self, PyObject *other, int cmp_op)
207207
return PyBool_FromLong(cmp);
208208
}
209209

210+
static PyObject *
211+
QuadPrecision_float(QuadPrecisionObject *self)
212+
{
213+
if (self->backend == BACKEND_SLEEF) {
214+
return PyFloat_FromDouble(Sleef_cast_to_doubleq1(self->value.sleef_value));
215+
} else {
216+
return PyFloat_FromDouble((double)self->value.longdouble_value);
217+
}
218+
}
219+
220+
static PyObject *
221+
QuadPrecision_int(QuadPrecisionObject *self)
222+
{
223+
if (self->backend == BACKEND_SLEEF) {
224+
return PyLong_FromLongLong(Sleef_cast_to_int64q1(self->value.sleef_value));
225+
} else {
226+
return PyLong_FromLongLong((long long)self->value.longdouble_value);
227+
}
228+
}
229+
230+
210231
PyNumberMethods quad_as_scalar = {
211232
.nb_add = (binaryfunc)quad_binary_func<quad_add, ld_add>,
212233
.nb_subtract = (binaryfunc)quad_binary_func<quad_sub, ld_sub>,
@@ -216,5 +237,7 @@ PyNumberMethods quad_as_scalar = {
216237
.nb_positive = (unaryfunc)quad_unary_func<quad_positive, ld_positive>,
217238
.nb_absolute = (unaryfunc)quad_unary_func<quad_absolute, ld_absolute>,
218239
.nb_bool = (inquiry)quad_nonzero,
240+
.nb_int = (unaryfunc)QuadPrecision_int,
241+
.nb_float = (unaryfunc)QuadPrecision_float,
219242
.nb_true_divide = (binaryfunc)quad_binary_func<quad_div, ld_div>,
220243
};

‎temp.py‎

Lines changed: 283 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,283 @@
1+
import numpy_quaddtype as npq
2+
import numpy as np
3+
4+
5+
def test_scalar_ops(backend):
6+
print(f"\nTesting scalar operations for {backend} backend:")
7+
8+
# Create QuadPrecision instances
9+
q1 = npq.QuadPrecision(
10+
"3.14159265358979323846264338327950288", backend=backend)
11+
q2 = npq.QuadPrecision(
12+
"-2.71828182845904523536028747135266250", backend=backend)
13+
14+
# Test unary operations
15+
print("\nUnary operations:")
16+
print(f" Negation of q1: {-q1}")
17+
print(f" Absolute value of q2: {abs(q2)}")
18+
19+
# Test binary operations
20+
print("\nBinary operations:")
21+
print(f" Addition: {q1 + q2}")
22+
print(f" Subtraction: {q1 - q2}")
23+
print(f" Multiplication: {q1 * q2}")
24+
print(f" Division: {q1 / q2}")
25+
26+
# Test comparison operations
27+
print("\nComparison operations:")
28+
print(f" q1 == q2: {q1 == q2}")
29+
print(f" q1 != q2: {q1 != q2}")
30+
print(f" q1 < q2: {q1 < q2}")
31+
print(f" q1 <= q2: {q1 <= q2}")
32+
print(f" q1 > q2: {q1 > q2}")
33+
print(f" q1 >= q2: {q1 >= q2}")
34+
35+
# Test operations with Python numbers
36+
print("\nOperations with Python numbers:")
37+
print(f" q1 + 1: {q1 + 1}")
38+
print(f" q1 - 2.5: {q1 - 2.5}")
39+
print(f" q1 * 3: {q1 * 3}")
40+
print(f" q1 / 2: {q1 / 2}")
41+
42+
# Test boolean conversion
43+
print("\nBoolean conversion:")
44+
print(f" bool(q1): {np.bool(q1)}")
45+
print(
46+
f" bool(npq.QuadPrecision('0', backend=backend)): {np.bool(npq.QuadPrecision('0', backend=backend))}")
47+
48+
49+
def test_casting(backend):
50+
print(f"\nTesting {backend} backend:")
51+
52+
# Create QuadPrecision instances
53+
q1 = npq.QuadPrecision(
54+
"3.14159265358979323846264338327950288", backend=backend)
55+
q2 = npq.QuadPrecision(
56+
"-2.71828182845904523536028747135266250", backend=backend)
57+
58+
# Test casting from QuadPrecision to numpy dtypes
59+
print("Casting from QuadPrecision to numpy dtypes:")
60+
print(f" float32: {np.float32(q1)}")
61+
print(f" float64: {np.float64(q1)}")
62+
print(f" int64: {np.int64(q1)}")
63+
print(f" uint64: {np.uint64(q1)}")
64+
65+
# Test casting from numpy dtypes to QuadPrecision
66+
print("\nCasting from numpy dtypes to QuadPrecision:")
67+
print(
68+
f" float32: {np.float32(3.14159).astype(npq.QuadPrecDType(backend=backend))}")
69+
print(
70+
f" float64: {np.float64(2.71828182845904).astype(npq.QuadPrecDType(backend=backend))}")
71+
print(
72+
f" int64: {np.int64(-1234567890).astype(npq.QuadPrecDType(backend=backend))}")
73+
print(
74+
f" uint64: {np.uint64(9876543210).astype(npq.QuadPrecDType(backend=backend))}")
75+
76+
# Test array operations
77+
print("\nArray operations:")
78+
q_array = np.array([q1, q2], dtype=npq.QuadPrecDType(backend=backend))
79+
print(f" QuadPrecision array: {q_array}")
80+
81+
np_array = np.array([3.14, -2.71, 1.41, -1.73], dtype=np.float64)
82+
q_from_np = np_array.astype(npq.QuadPrecDType(backend=backend))
83+
print(f" Numpy to QuadPrecision: {q_from_np}")
84+
85+
back_to_np = np.array(q_from_np, dtype=np.float64)
86+
print(f" QuadPrecision to Numpy: {back_to_np}")
87+
88+
# Test precision maintenance
89+
large_int = 12345678901234567890
90+
q_large = np.array([large_int], dtype=np.uint64).astype(
91+
npq.QuadPrecDType(backend=backend))[0]
92+
print(f"\nPrecision test:")
93+
print(f" Original large int: {large_int}")
94+
print(f" QuadPrecision: {q_large}")
95+
print(f" Back to int: {np.uint64(q_large)}")
96+
97+
# Test edge cases
98+
99+
100+
def test_edge_cases(backend):
101+
print(f"\nTesting negative numbers for {backend} backend:")
102+
103+
# Test various negative numbers
104+
test_values = [
105+
-1.0,
106+
-1e10,
107+
-1e100,
108+
-1e300,
109+
np.nextafter(np.finfo(np.float64).min, 0),
110+
np.finfo(np.float64).min
111+
]
112+
113+
for value in test_values:
114+
q_value = npq.QuadPrecision(str(value), backend=backend)
115+
print(f" Original: {value}")
116+
print(f" QuadPrecision: {q_value}")
117+
print(f" Back to float64: {np.float64(q_value)}")
118+
print()
119+
120+
# Test value beyond float64 precision
121+
beyond_float64_precision = "1.7976931348623157081452742373170435e+308"
122+
q_beyond = npq.QuadPrecision(beyond_float64_precision, backend=backend)
123+
print(f" Beyond float64 precision: {q_beyond}")
124+
q_float64_max = npq.QuadPrecision(
125+
str(np.finfo(np.float64).max), backend=backend)
126+
diff = q_beyond - q_float64_max
127+
print(f" Difference from float64 max: {diff}")
128+
print(
129+
f" Difference is positive: {diff > npq.QuadPrecision('0', backend=backend)}")
130+
131+
# Test epsilon (smallest representable difference between two numbers)
132+
q_epsilon = npq.QuadPrecision(
133+
str(np.finfo(np.float64).eps), backend=backend)
134+
print(f" Float64 epsilon in QuadPrecision: {q_epsilon}")
135+
q_one = npq.QuadPrecision("1", backend=backend)
136+
q_one_plus_epsilon = q_one + q_epsilon
137+
print(f" 1 + epsilon != 1: {q_one_plus_epsilon != q_one}")
138+
print(f" (1 + epsilon) - 1: {q_one_plus_epsilon - q_one}")
139+
140+
141+
def test_ufuncs(backend):
142+
print(f"\nTesting ufuncs for {backend} backend:")
143+
144+
# Create QuadPrecision arrays
145+
q_array1 = np.array([1, 2, 3], dtype=npq.QuadPrecDType(backend=backend))
146+
q_array2 = np.array([1, 2, 3], dtype=npq.QuadPrecDType(backend=backend))
147+
148+
# Test unary ufuncs
149+
print("\nUnary unfuncs:")
150+
print(f" negative: {np.negative(q_array1)}")
151+
print(f" absolute: {np.absolute(q_array1)}")
152+
print(f" rint: {np.rint(q_array1)}")
153+
print(f" floor: {np.floor(q_array1)}")
154+
print(f" ceil: {np.ceil(q_array1)}")
155+
print(f" trunc: {np.trunc(q_array1)}")
156+
print(f" sqrt: {np.sqrt(q_array1)}")
157+
print(f" square: {np.square(q_array1)}")
158+
print(f" log: {np.log(q_array1)}")
159+
print(f" log2: {np.log2(q_array1)}")
160+
print(f" log10: {np.log10(q_array1)}")
161+
print(f" exp: {np.exp(q_array1)}")
162+
print(f" exp2: {np.exp2(q_array1)}")
163+
164+
# Test binary ufuncs
165+
print("\nBinary ufuncs:")
166+
print(f" add: {np.add(q_array1, q_array2)}")
167+
print(f" subtract: {np.subtract(q_array1, q_array2)}")
168+
print(f" multiply: {np.multiply(q_array1, q_array2)}")
169+
print(f" divide: {np.divide(q_array1, q_array2)}")
170+
print(f" power: {np.power(q_array1, q_array2)}")
171+
print(f" mod: {np.mod(q_array1, q_array2)}")
172+
print(f" minimum: {np.minimum(q_array1, q_array2)}")
173+
print(f" maximum: {np.maximum(q_array1, q_array2)}")
174+
175+
# Test comparison ufuncs
176+
print("\nComparison ufuncs:")
177+
print(f" equal: {np.equal(q_array1, q_array2)}")
178+
print(f" not_equal: {np.not_equal(q_array1, q_array2)}")
179+
print(f" less: {np.less(q_array1, q_array2)}")
180+
print(f" less_equal: {np.less_equal(q_array1, q_array2)}")
181+
print(f" greater: {np.greater(q_array1, q_array2)}")
182+
print(f" greater_equal: {np.greater_equal(q_array1, q_array2)}")
183+
184+
# Test mixed operations with numpy arrays
185+
print(f"Testing backend: {backend}")
186+
print("\nMixed operations with numpy arrays:")
187+
np_array = np.array([1.0, 2.0, 3.0], dtype=np.float64)
188+
print(f" add: {np.add(q_array1, np_array)}")
189+
print(f" multiply: {np.multiply(q_array1, np_array)}")
190+
print(f" divide: {np.divide(q_array1, np_array)}")
191+
192+
# Test reduction operations
193+
print("\nReduction operations:")
194+
print(f" sum: {np.sum(q_array1)}")
195+
print(f" prod: {np.prod(q_array1)}")
196+
print(f" min: {np.min(q_array1)}")
197+
print(f" max: {np.max(q_array1)}")
198+
199+
from numpy_quaddtype import QuadPrecision, QuadPrecDType
200+
201+
def test_quad_precision():
202+
print("Testing QuadPrecision scalar:")
203+
204+
# Test different initializations
205+
values = [
206+
0,
207+
1,
208+
-1,
209+
3.14159265358979323846,
210+
1e100,
211+
1e-100,
212+
float('inf'),
213+
float('-inf'),
214+
float('nan')
215+
]
216+
217+
for val in values:
218+
q = QuadPrecision(val)
219+
print(f"Value: {val}")
220+
print(f" str: {str(q)}")
221+
print(f" repr: {repr(q)}")
222+
223+
# Test different backends
224+
print("\nTesting backends:")
225+
q_sleef = QuadPrecision(3.14159265358979323846, backend='sleef')
226+
q_longdouble = QuadPrecision(3.14159265358979323846, backend='longdouble')
227+
print(f"Sleef: {q_sleef}")
228+
print(f"Long double: {q_longdouble}")
229+
230+
def test_quad_dtype():
231+
print("\nTesting QuadPrecDType:")
232+
233+
# Create an array with QuadPrecDType
234+
arr = np.array([0, 1, -1, 3.14159265358979323846, 1e100, 1e-100],
235+
dtype=QuadPrecDType())
236+
237+
print("Array elements:")
238+
for elem in arr:
239+
print(f" {elem}")
240+
241+
print("\nFull array:")
242+
print(arr)
243+
244+
# Test different backends
245+
print("\nTesting backends in arrays:")
246+
arr_sleef = np.array([3.14159265358979323846], dtype=QuadPrecDType(backend='sleef'))
247+
arr_longdouble = np.array([3.14159265358979323846], dtype=QuadPrecDType(backend='longdouble'))
248+
print(f"Sleef array: {arr_sleef}")
249+
print(f"Long double array: {arr_longdouble}")
250+
251+
def test_operations():
252+
print("\nTesting basic operations:")
253+
a = QuadPrecision(3.14159265358979323846)
254+
b = QuadPrecision(2.71828182845904523536)
255+
256+
print(f"a = {a}")
257+
print(f"b = {b}")
258+
print(f"a + b = {a + b}")
259+
print(f"a - b = {a - b}")
260+
print(f"a * b = {a * b}")
261+
print(f"a / b = {a / b}")
262+
263+
def test():
264+
# Run tests for both backends
265+
for backend in ['sleef', 'longdouble']:
266+
test_scalar_ops(backend)
267+
test_casting(backend)
268+
test_edge_cases(backend)
269+
test_ufuncs(backend)
270+
test_quad_precision()
271+
test_quad_dtype()
272+
test_operations()
273+
print("*"*50)
274+
275+
print("All tests completed successfully")
276+
277+
def dot(a, b):
278+
r = np.dot(a, b)
279+
return r
280+
281+
if __name__ == "__main__":
282+
a = np.array([1, 2, 3], dtype=QuadPrecDType())
283+
print(dot(a, a))

0 commit comments

Comments
 (0)