Skip to content

Commit 7087893

Browse files
committed
Add new cost() and sint() Cython functions that take turns
These return exact values for 90° points instead of very near values.
1 parent 828442d commit 7087893

File tree

6 files changed

+75
-46
lines changed

6 files changed

+75
-46
lines changed

src/flitter/language/functions.pyx

Lines changed: 13 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,14 @@ Flitter language functions
66

77
import cython
88

9-
from libc.math cimport floor, sin, cos, tan, asin, acos, sqrt, exp, atan2, log, log2, log10
9+
from libc.math cimport floor, sin, tan, asin, acos, sqrt, exp, atan2, log, log2, log10
1010
from libc.stdint cimport int64_t, uint64_t
1111
from cpython.object cimport PyObject
1212
from cpython.ref cimport Py_INCREF
1313
from cpython.tuple cimport PyTuple_New, PyTuple_GET_ITEM, PyTuple_SET_ITEM
1414

1515
from ..cache import SharedCache
16-
from ..model cimport Vector, Matrix33, Matrix44, Quaternion, Context, null_, true_, false_
16+
from ..model cimport Vector, Matrix33, Matrix44, Quaternion, Context, null_, true_, false_, cost, sint
1717

1818

1919
cdef double Pi = 3.141592653589793115997963468544185161590576171875
@@ -162,12 +162,12 @@ cdef class normal(uniform):
162162
if u1 == 0:
163163
u1, u2 = u2, u1
164164
self.R = sqrt(-2 * log(u1))
165-
self.th = Tau * u2
165+
self.th = u2
166166
self.i = i
167167
self.cached = True
168168
if odd:
169-
return self.R * sin(self.th)
170-
return self.R * cos(self.th)
169+
return self.R * sint(self.th)
170+
return self.R * cost(self.th)
171171

172172

173173
@context_func
@@ -241,7 +241,7 @@ def sinv(Vector theta not None):
241241
cdef Vector ys = Vector.__new__(Vector)
242242
cdef int64_t i, n = theta.length
243243
for i in range(ys.allocate_numbers(n)):
244-
ys.numbers[i] = sin(theta.numbers[i] * Tau)
244+
ys.numbers[i] = sint(theta.numbers[i])
245245
return ys
246246

247247

@@ -251,7 +251,7 @@ def cosv(Vector theta not None):
251251
cdef Vector ys = Vector.__new__(Vector)
252252
cdef int64_t i, n = theta.length
253253
for i in range(ys.allocate_numbers(n)):
254-
ys.numbers[i] = cos(theta.numbers[i] * Tau)
254+
ys.numbers[i] = cost(theta.numbers[i])
255255
return ys
256256

257257

@@ -292,8 +292,8 @@ def polar(Vector theta not None):
292292
cdef int64_t i, n = theta.length
293293
ys.allocate_numbers(n*2)
294294
for i in range(0, n):
295-
ys.numbers[i*2] = cos(theta.numbers[i] * Tau)
296-
ys.numbers[i*2+1] = sin(theta.numbers[i] * Tau)
295+
ys.numbers[i*2] = cost(theta.numbers[i])
296+
ys.numbers[i*2+1] = sint(theta.numbers[i])
297297
return ys
298298

299299

@@ -403,7 +403,7 @@ def sine(Vector xs not None):
403403
return null_
404404
cdef Vector ys = Vector.__new__(Vector)
405405
for i in range(ys.allocate_numbers(xs.length)):
406-
ys.numbers[i] = (1 - cos(Tau * xs.numbers[i])) / 2
406+
ys.numbers[i] = (1 - cost(xs.numbers[i])) / 2
407407
return ys
408408

409409

@@ -1071,11 +1071,11 @@ def oklch(Vector lch):
10711071
if lch.length != 3 or lch.objects is not None:
10721072
return null_
10731073
cdef Vector lab = Vector.__new__(Vector)
1074-
cdef double L=lch.numbers[0], C=lch.numbers[1], h=lch.numbers[2]*Tau
1074+
cdef double L=lch.numbers[0], C=lch.numbers[1], h=lch.numbers[2]
10751075
lab.allocate_numbers(3)
10761076
lab.numbers[0] = L
1077-
lab.numbers[1] = C * cos(h)
1078-
lab.numbers[2] = C * sin(h)
1077+
lab.numbers[1] = C * cost(h)
1078+
lab.numbers[2] = C * sint(h)
10791079
return XYZ_to_sRGB.vmul(OKLab_M1_inv.vmul(OKLab_M2_inv.vmul(lab).pow(Three)))
10801080

10811081

src/flitter/model.pxd

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,10 @@ from libc.stdint cimport int64_t, uint64_t
55
from cpython.unicode cimport PyUnicode_DATA, PyUnicode_GET_LENGTH, PyUnicode_KIND, PyUnicode_READ
66

77

8+
cdef double sint(double t) noexcept nogil
9+
cdef double cost(double t) noexcept nogil
10+
11+
812
# SplitMix64 algorithm [http://xoshiro.di.unimi.it/splitmix64.c]
913
#
1014
cdef uint64_t HASH_START

src/flitter/model.pyx

Lines changed: 33 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,32 @@ cdef dict SymbolTable = {}
3030
cdef dict ReverseSymbolTable = {}
3131

3232

33+
cdef double sint(double t) noexcept nogil:
34+
t = t - c_floor(t)
35+
if t == 0:
36+
return 0
37+
if t == 0.25:
38+
return 1
39+
if t == 0.5:
40+
return 0
41+
if t == 0.75:
42+
return -1
43+
return sin(Tau * t)
44+
45+
46+
cdef double cost(double t) noexcept nogil:
47+
t = t - c_floor(t)
48+
if t == 0:
49+
return 1
50+
if t == 0.25:
51+
return 0
52+
if t == 0.5:
53+
return -1
54+
if t == 0.75:
55+
return 0
56+
return cos(Tau * t)
57+
58+
3359
cdef inline int64_t vector_compare(Vector left, Vector right) noexcept:
3460
if left is right:
3561
return 0
@@ -1241,7 +1267,7 @@ cdef class Matrix33(Vector):
12411267
cdef Matrix33 _rotate(double turns):
12421268
if isnan(turns):
12431269
return None
1244-
cdef double theta = turns*Tau, cth = cos(theta), sth = sin(theta)
1270+
cdef double cth = cost(turns), sth = sint(turns)
12451271
cdef Matrix33 result = Matrix33.__new__(Matrix33)
12461272
cdef double* numbers = result._numbers
12471273
numbers[0] = cth
@@ -1573,7 +1599,7 @@ cdef class Matrix44(Vector):
15731599
cdef Matrix44 _rotate_x(double turns):
15741600
if isnan(turns):
15751601
return None
1576-
cdef double theta = turns*Tau, cth = cos(theta), sth = sin(theta)
1602+
cdef double cth = cost(turns), sth = sint(turns)
15771603
cdef Matrix44 result = Matrix44.__new__(Matrix44)
15781604
cdef double* numbers = result._numbers
15791605
numbers[0] = 1
@@ -1597,7 +1623,7 @@ cdef class Matrix44(Vector):
15971623
cdef Matrix44 _rotate_y(double turns):
15981624
if isnan(turns):
15991625
return None
1600-
cdef double theta = turns*Tau, cth = cos(theta), sth = sin(theta)
1626+
cdef double cth = cost(turns), sth = sint(turns)
16011627
cdef Matrix44 result = Matrix44.__new__(Matrix44)
16021628
cdef double* numbers = result._numbers
16031629
numbers[0] = cth
@@ -1623,7 +1649,7 @@ cdef class Matrix44(Vector):
16231649
cdef Matrix44 _rotate_z(double turns):
16241650
if isnan(turns):
16251651
return None
1626-
cdef double theta = turns*Tau, cth = cos(theta), sth = sin(theta)
1652+
cdef double cth = cost(turns), sth = sint(turns)
16271653
cdef Matrix44 result = Matrix44.__new__(Matrix44)
16281654
cdef double* numbers = result._numbers
16291655
numbers[0] = cth
@@ -2013,11 +2039,11 @@ cdef class Quaternion(Vector):
20132039
if axis is None or axis.numbers == NULL or axis.length != 3:
20142040
raise ValueError("Axis must be a numeric 3-vector")
20152041
cdef double x=axis.numbers[0], y=axis.numbers[1], z=axis.numbers[2]
2016-
cdef double half_theta = rotation * Pi
2017-
cdef double k = sin(half_theta) / sqrt(x*x + y*y + z*z)
2042+
cdef double half_rotation = rotation / 2
2043+
cdef double k = sint(half_rotation) / sqrt(x*x + y*y + z*z)
20182044
cdef Quaternion result = Quaternion.__new__(Quaternion)
20192045
cdef double* numbers = result._numbers
2020-
numbers[0] = cos(half_theta)
2046+
numbers[0] = cost(half_rotation)
20212047
numbers[1] = k * x
20222048
numbers[2] = k * y
20232049
numbers[3] = k * z

src/flitter/render/window/models.pyx

Lines changed: 10 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -6,14 +6,14 @@ import numpy as np
66
import trimesh
77
import trimesh.proximity
88

9-
from libc.math cimport cos, sin, sqrt, atan2, abs, floor as c_floor
9+
from libc.math cimport sqrt, atan2, abs, floor as c_floor
1010
from libc.stdint cimport int32_t, int64_t
1111
from cpython.object cimport PyObject
1212
from cpython.dict cimport PyDict_GetItem
1313

1414
from ... import name_patch
1515
from ...cache import SharedCache
16-
from ...model cimport true_, Context, Vector, StateDict, HASH_START, HASH_UPDATE, HASH_STRING, double_long, Matrix33
16+
from ...model cimport true_, Context, Vector, StateDict, HASH_START, HASH_UPDATE, HASH_STRING, double_long, Matrix33, sint, cost
1717
from ...timer cimport perf_counter
1818
from ...language.vm cimport VectorStack
1919

@@ -1294,8 +1294,8 @@ cdef class Sphere(PrimitiveModel):
12941294
r, z, v = 1, 0, 0.5
12951295
else:
12961296
th = hemisphere * (1 - <float>row / nrows) / 4
1297-
v, th = 2 * th + 0.5, th * Tau
1298-
r, z = cos(th), sin(th)
1297+
v = 2 * th + 0.5
1298+
r, z = cost(th), sint(th)
12991299
for col in range(row + 1):
13001300
if row == 0:
13011301
u = (side + 0.5) / 4
@@ -1310,8 +1310,7 @@ cdef class Sphere(PrimitiveModel):
13101310
y = r if side == 0 else -r if side == 2 else 0
13111311
else:
13121312
u = (side + (<float>col / row)) / 4
1313-
th = Tau * u
1314-
x, y = r * cos(th), r * sin(th)
1313+
x, y = r * cost(u), r * sint(u)
13151314
vertices[i, 0], vertices[i, 1], vertices[i, 2] = x, y, z
13161315
vertices[i, 3], vertices[i, 4], vertices[i, 5] = x, y, z
13171316
vertices[i, 6], vertices[i, 7] = u, v
@@ -1369,16 +1368,15 @@ cdef class Cylinder(PrimitiveModel):
13691368
cdef float[:, :] vertices = vertices_array
13701369
cdef object faces_array = np.empty((n*4, 3), dtype='i4')
13711370
cdef int32_t[:, :] faces = faces_array
1372-
cdef float x, y, th, u, u_
1371+
cdef float x, y, u, u_
13731372
for i in range(n+1):
13741373
j = k = i * 6
13751374
u = <float>i / n
13761375
u_ = (i+0.5) / n
13771376
if i == 0 or i == n:
13781377
x, y = 1, 0
13791378
else:
1380-
th = Tau * u
1381-
x, y = cos(th), sin(th)
1379+
x, y = cost(u), sint(u)
13821380
# bottom centre (k):
13831381
vertices[j, 0], vertices[j, 1], vertices[j, 2] = 0, 0, -0.5
13841382
vertices[j, 3], vertices[j, 4], vertices[j, 5] = 0, 0, -1
@@ -1459,17 +1457,15 @@ cdef class Cone(PrimitiveModel):
14591457
cdef float[:, :] vertices = vertices_array
14601458
cdef object faces_array = np.empty((n*2, 3), dtype='i4')
14611459
cdef int32_t[:, :] faces = faces_array
1462-
cdef float x, y, th, u, u_
1460+
cdef float x, y, u, u_
14631461
for i in range(n+1):
14641462
j = k = i * 4
14651463
u = <double>i / n
14661464
u_ = (i+0.5) / n
1467-
th_ = Tau * u_
14681465
if i == 0 or i == n:
14691466
x, y = 1, 0
14701467
else:
1471-
th = Tau * u
1472-
x, y = cos(th), sin(th)
1468+
x, y = cost(u), sint(u)
14731469
# bottom centre (k):
14741470
vertices[j, 0], vertices[j, 1], vertices[j, 2] = 0, 0, -0.5
14751471
vertices[j, 3], vertices[j, 4], vertices[j, 5] = 0, 0, -1
@@ -1487,7 +1483,7 @@ cdef class Cone(PrimitiveModel):
14871483
j += 1
14881484
# side top (k+3):
14891485
vertices[j, 0], vertices[j, 1], vertices[j, 2] = 0, 0, 0.5
1490-
vertices[j, 3], vertices[j, 4], vertices[j, 5] = cos(th_)*RootHalf, sin(th_)*RootHalf, RootHalf
1486+
vertices[j, 3], vertices[j, 4], vertices[j, 5] = cost(u_)*RootHalf, sint(u_)*RootHalf, RootHalf
14911487
vertices[j, 6], vertices[j, 7] = u_, 1
14921488
if i < n:
14931489
j = i * 2

tests/test_functions.py

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -495,7 +495,8 @@ def setUp(self):
495495
def test_cos(self):
496496
self.assertEqual(cosv(null), null)
497497
self.assertEqual(cosv(Vector('hello')), null)
498-
theta = Vector.range(0, 1, 0.01)
498+
self.assertEqual(cosv(Vector([0, 0.25, 0.5, 0.75, 1])), Vector([1, 0, -1, 0, 1]))
499+
theta = Vector.range(0.001, 1, 0.01)
499500
values = [math.cos(th) for th in theta*Tau]
500501
for i in range(len(values)):
501502
self.assertEqual(cosv(theta.item(i)), values[i])
@@ -513,7 +514,8 @@ def test_acos(self):
513514
def test_sin(self):
514515
self.assertEqual(sinv(null), null)
515516
self.assertEqual(sinv(Vector('hello')), null)
516-
theta = Vector.range(0, 1, 0.01)
517+
self.assertEqual(sinv(Vector([0, 0.25, 0.5, 0.75, 1])), Vector([0, 1, 0, -1, 0]))
518+
theta = Vector.range(0.001, 1, 0.01)
517519
values = [math.sin(th) for th in theta*Tau]
518520
for i in range(len(values)):
519521
self.assertEqual(sinv(theta.item(i)), values[i])
@@ -571,7 +573,8 @@ def test_dot(self):
571573
def test_polar(self):
572574
self.assertEqual(polar(null), null)
573575
self.assertEqual(polar(Vector('hello')), null)
574-
theta = Vector.range(0, 1, 0.01)
576+
self.assertEqual(polar(Vector([0, 0.25, 0.5, 0.75, 1])), Vector([1, 0, 0, 1, -1, 0, 0, -1, 1, 0]))
577+
theta = Vector.range(0.001, 1, 0.01)
575578
values = [(math.cos(th), math.sin(th)) for th in theta*Tau]
576579
for i in range(len(values)):
577580
self.assertEqual(polar(theta.item(i)), values[i])

tests/test_quaternion.py

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -128,18 +128,18 @@ def test_multiply(self):
128128
self.assertAllAlmostEqual(q @ q @ q @ q @ q @ q, Quaternion(1))
129129

130130
def test_inverse(self):
131-
self.assertAlmostEqual(Quaternion.euler([1, 0, 0], 0.125).inverse(), Quaternion.euler([1, 0, 0], -0.125))
132-
self.assertAlmostEqual(Quaternion.euler([0, 1, 0], 0.125).inverse(), Quaternion.euler([0, 1, 0], -0.125))
133-
self.assertAlmostEqual(Quaternion.euler([0, 0, 1], 0.125).inverse(), Quaternion.euler([0, 0, 1], -0.125))
131+
self.assertAllAlmostEqual(Quaternion.euler([1, 0, 0], 0.125).inverse(), Quaternion.euler([1, 0, 0], -0.125))
132+
self.assertAllAlmostEqual(Quaternion.euler([0, 1, 0], 0.125).inverse(), Quaternion.euler([0, 1, 0], -0.125))
133+
self.assertAllAlmostEqual(Quaternion.euler([0, 0, 1], 0.125).inverse(), Quaternion.euler([0, 0, 1], -0.125))
134134
q = Quaternion.euler([1, 1, 1], 1/3)
135-
self.assertAlmostEqual(q @ q.inverse(), Quaternion())
136-
self.assertAlmostEqual(q.inverse() @ q, Quaternion())
137-
self.assertAlmostEqual(q @ q.inverse() @ q, q)
135+
self.assertAllAlmostEqual(q @ q.inverse(), Quaternion())
136+
self.assertAllAlmostEqual(q.inverse() @ q, Quaternion())
137+
self.assertAllAlmostEqual(q @ q.inverse() @ q, q)
138138

139139
def test_normalize(self):
140-
self.assertAlmostEqual(Quaternion([0.5, 0.5, 0.5, 0.5]).normalize(), [0.5, 0.5, 0.5, 0.5])
141-
self.assertAlmostEqual(Quaternion([1, 1, 1, 1]).normalize(), [0.5, 0.5, 0.5, 0.5])
142-
self.assertAlmostEqual(Quaternion(10).normalize(), Quaternion(1))
140+
self.assertAllAlmostEqual(Quaternion([0.5, 0.5, 0.5, 0.5]).normalize(), [0.5, 0.5, 0.5, 0.5])
141+
self.assertAllAlmostEqual(Quaternion([1, 1, 1, 1]).normalize(), [0.5, 0.5, 0.5, 0.5])
142+
self.assertAllAlmostEqual(Quaternion(10).normalize(), Quaternion(1))
143143

144144
def test_conjugate(self):
145145
qx = Quaternion.euler([1, 0, 0], 0.25)

0 commit comments

Comments
 (0)