@@ -111,146 +111,49 @@ namespace bclibc
111111 BCLIBC_Curve BCLIBC_Curve_fromPylist (PyObject *data_points)
112112 {
113113 Py_ssize_t n = PyList_Size (data_points);
114- // Check for PyList_Size error or insufficient data
115114 if (n < 2 )
116115 {
117116 if (n < 0 )
118- {
119- // Error in Python List access
120117 return BCLIBC_Curve ();
121- }
122- // insufficient data; return empty curve
123118 PyErr_SetString (PyExc_ValueError, " BCLIBC_Curve requires at least 2 data points." );
124119 return BCLIBC_Curve ();
125120 }
126121
127- // --- Phase 1: Extract x (Mach) and y (CD) into vectors (RAII) ---
128- // Instead of malloc(x) and malloc(y), we use std::vector.
122+ // --- Phase 1: Prepare data ---
129123 std::vector<double > x ((size_t )n);
130124 std::vector<double > y ((size_t )n);
131125
132- // If memory allocation fails here, std::vector will throw std::bad_alloc,
133- // which Cython's 'except +' will catch.
134-
135126 for (Py_ssize_t i = 0 ; i < n; ++i)
136127 {
137- PyObject *item = PyList_GetItem (data_points, i); // borrowed
138- // Note: PyList_GetItem should not return nullptr here unless list size changed,
139- // but we keep the check for robustness.
128+ PyObject *item = PyList_GetItem (data_points, i);
140129 if (item == nullptr )
141130 {
142131 PyErr_SetString (PyExc_IndexError, " Curve generation: list item access failed." );
143132 return BCLIBC_Curve ();
144133 }
145134
146- // Get x (Mach) and y (CD) attributes
147- PyObject *xobj = PyObject_GetAttrString (item, " Mach" ); // new ref
148- PyObject *yobj = PyObject_GetAttrString (item, " CD" ); // new ref
135+ PyObject *xobj = PyObject_GetAttrString (item, " Mach" );
136+ PyObject *yobj = PyObject_GetAttrString (item, " CD" );
149137
150138 if (xobj == nullptr || yobj == nullptr )
151139 {
152- // Attribute access failed (PyErr_Occurred is true)
153140 Py_XDECREF (xobj);
154141 Py_XDECREF (yobj);
155142 return BCLIBC_Curve ();
156143 }
157144
158- // Convert and assign
159145 x[i] = PyFloat_AsDouble (xobj);
160146 y[i] = PyFloat_AsDouble (yobj);
161147
162148 Py_DECREF (xobj);
163149 Py_DECREF (yobj);
164150
165151 if (PyErr_Occurred ())
166- {
167- // Conversion failed (e.g., attribute was not a number)
168152 return BCLIBC_Curve ();
169- }
170- }
171-
172- // --- Phase 2: Calculate PCHIP Slopes and Coefficients ---
173-
174- Py_ssize_t nm1 = n - 1 ;
175-
176- // Temporary vectors for PCHIP calculation (RAII: no more free calls!)
177- std::vector<double > h ((size_t )nm1); // Steps h[i] = x[i+1] - x[i]
178- std::vector<double > d ((size_t )nm1); // Slopes d[i] = (y[i+1] - y[i]) / h[i]
179- std::vector<double > m ((size_t )n); // Final calculated slopes m[i]
180-
181- // Final result vector (n-1 segments)
182- BCLIBC_Curve curve_points ((size_t )nm1);
183-
184- for (Py_ssize_t i = 0 ; i < nm1; ++i)
185- {
186- h[i] = x[i + 1 ] - x[i];
187- d[i] = (y[i + 1 ] - y[i]) / h[i];
188- }
189-
190- if (n == 2 )
191- {
192- m[0 ] = d[0 ];
193- m[1 ] = d[0 ];
194- }
195- else
196- {
197- // Interior slopes (Fritsch–Carlson algorithm)
198- for (Py_ssize_t i = 1 ; i < n - 1 ; ++i)
199- {
200- if (d[i - 1 ] == 0.0 || d[i] == 0.0 || d[i - 1 ] * d[i] < 0.0 )
201- {
202- m[i] = 0.0 ;
203- }
204- else
205- {
206- double w1 = 2.0 * h[i] + h[i - 1 ];
207- double w2 = h[i] + 2.0 * h[i - 1 ];
208- m[i] = (w1 + w2) / (w1 / d[i - 1 ] + w2 / d[i]);
209- }
210- }
211-
212- // Endpoint m[0]
213- double m0 = ((2.0 * h[0 ] + h[1 ]) * d[0 ] - h[0 ] * d[1 ]) / (h[0 ] + h[1 ]);
214- if (m0 * d[0 ] <= 0.0 )
215- m0 = 0.0 ;
216- else if ((d[0 ] * d[1 ] < 0.0 ) && (std::fabs (m0) > 3.0 * std::fabs (d[0 ])))
217- m0 = 3.0 * d[0 ];
218- m[0 ] = m0;
219-
220- // Endpoint m[n-1]
221- double mn = ((2.0 * h[n - 2 ] + h[n - 3 ]) * d[n - 2 ] - h[n - 2 ] * d[n - 3 ]) / (h[n - 2 ] + h[n - 3 ]);
222- if (mn * d[n - 2 ] <= 0.0 )
223- mn = 0.0 ;
224- else if ((d[n - 2 ] * d[n - 3 ] < 0.0 ) && (std::fabs (mn) > 3.0 * std::fabs (d[n - 2 ])))
225- mn = 3.0 * d[n - 2 ];
226- m[n - 1 ] = mn;
227- }
228-
229- // Build per-segment cubic coefficients in dx=(x-x_i):
230- for (Py_ssize_t i = 0 ; i < nm1; ++i)
231- {
232- double H = h[i];
233- double yi = y[i];
234- double mi = m[i];
235- double mip1 = m[i + 1 ];
236-
237- // A and B helpers
238- double A = (y[i + 1 ] - yi - mi * H) / (H * H);
239- double B = (mip1 - mi) / H;
240-
241- double a = (B - 2.0 * A) / H;
242- double b = 3.0 * A - B;
243-
244- // Assign directly to the final curve vector element
245- curve_points[i].a = a;
246- curve_points[i].b = b;
247- curve_points[i].c = mi;
248- curve_points[i].d = yi;
249153 }
250154
251- // When we return curve_points, all temporary vectors (x, y, h, d, m)
252- // are automatically destroyed and memory freed.
253- return curve_points;
155+ // --- Phase 2: Call universal core ---
156+ return build_pchip_curve_from_arrays (x, y);
254157 }
255158}; // namespace bclibc
256159
0 commit comments