1 /*
2     Multipack project.
3     This file is generated by setmodules.py. Do not modify it.
4  */
5 
6 #include <Python.h>
7 #include "numpy/arrayobject.h"
8 
9 #define PyInt_AsLong PyLong_AsLong
10 
11 static PyObject *fitpack_error;
12 #include "__fitpack.h"
13 
14 #ifdef HAVE_ILP64
15 
16 #define F_INT npy_int64
17 #define F_INT_NPY NPY_INT64
18 
19 #if NPY_BITSOF_SHORT == 64
20 #define F_INT_PYFMT   "h"
21 #elif NPY_BITSOF_INT == 64
22 #define F_INT_PYFMT   "i"
23 #elif NPY_BITSOF_LONG == 64
24 #define F_INT_PYFMT   "l"
25 #elif NPY_BITSOF_LONGLONG == 64
26 #define F_INT_PYFMT   "L"
27 #else
28 #error No compatible 64-bit integer size. \
29        Please contact NumPy maintainers and give detailed information about your \
30        compiler and platform, or set NPY_USE_BLAS64_=0
31 #endif
32 
33 #else
34 
35 #define F_INT int
36 #define F_INT_NPY NPY_INT
37 #define F_INT_PYFMT   "i"
38 
39 #endif
40 
41 
42 /*
43  * Functions moved verbatim from __fitpack.h
44  */
45 
46 
47 /*
48  * Python-C wrapper of FITPACK (by P. Dierckx) (in netlib known as dierckx)
49  * Author: Pearu Peterson <pearu@ioc.ee>
50  * June 1.-4., 1999
51  * June 7. 1999
52  * $Revision$
53  * $Date$
54  */
55 
56 /*  module_methods:
57  * {"_curfit", fitpack_curfit, METH_VARARGS, doc_curfit},
58  * {"_spl_", fitpack_spl_, METH_VARARGS, doc_spl_},
59  * {"_splint", fitpack_splint, METH_VARARGS, doc_splint},
60  * {"_sproot", fitpack_sproot, METH_VARARGS, doc_sproot},
61  * {"_spalde", fitpack_spalde, METH_VARARGS, doc_spalde},
62  * {"_parcur", fitpack_parcur, METH_VARARGS, doc_parcur},
63  * {"_surfit", fitpack_surfit, METH_VARARGS, doc_surfit},
64  * {"_bispev", fitpack_bispev, METH_VARARGS, doc_bispev},
65  * {"_insert", fitpack_insert, METH_VARARGS, doc_insert},
66  */
67 
68 /* link libraries: (one item per line)
69    ddierckx
70  */
71 /* python files: (to be imported to Multipack.py)
72    fitpack.py
73  */
74 
75 #if defined(UPPERCASE_FORTRAN)
76 	#if defined(NO_APPEND_FORTRAN)
77 	/* nothing to do */
78 	#else
79 		#define CURFIT CURFIT_
80 		#define PERCUR PERCUR_
81 		#define SPALDE SPALDE_
82 		#define SPLDER SPLDER_
83 		#define SPLEV  SPLEV_
84 		#define SPLINT SPLINT_
85 		#define SPROOT SPROOT_
86 		#define PARCUR PARCUR_
87 		#define CLOCUR CLOCUR_
88 		#define SURFIT SURFIT_
89 		#define BISPEV BISPEV_
90 		#define PARDER PARDER_
91 		#define INSERT INSERT_
92 	#endif
93 #else
94 	#if defined(NO_APPEND_FORTRAN)
95 		#define CURFIT curfit
96 		#define PERCUR percur
97 		#define SPALDE spalde
98 		#define SPLDER splder
99 		#define SPLEV splev
100 		#define SPLINT splint
101 		#define SPROOT sproot
102 		#define PARCUR parcur
103 		#define CLOCUR clocur
104 		#define SURFIT surfit
105 		#define BISPEV bispev
106 		#define PARDER parder
107 		#define INSERT insert
108 	#else
109 		#define CURFIT curfit_
110 		#define PERCUR percur_
111 		#define SPALDE spalde_
112 		#define SPLDER splder_
113 		#define SPLEV splev_
114 		#define SPLINT splint_
115 		#define SPROOT sproot_
116 		#define PARCUR parcur_
117 		#define CLOCUR clocur_
118 		#define SURFIT surfit_
119 		#define BISPEV bispev_
120 		#define PARDER parder_
121 		#define INSERT insert_
122 	#endif
123 #endif
124 
125 void CURFIT(F_INT*,F_INT*,double*,double*,double*,double*,
126         double*,F_INT*,double*,F_INT*,F_INT*,double*,double*,
127         double*,double*,F_INT*,F_INT*,F_INT*);
128 void PERCUR(F_INT*,F_INT*,double*,double*,double*,F_INT*,
129         double*,F_INT*,F_INT*,double*,double*,double*,
130         double*,F_INT*,F_INT*,F_INT*);
131 void SPALDE(double*,F_INT*,double*,F_INT*,double*,double*,F_INT*);
132 void SPLDER(double*,F_INT*,double*,F_INT*,F_INT*,double*,
133         double*,F_INT*,F_INT*,double*,F_INT*);
134 void SPLEV(double*,F_INT*,double*,F_INT*,double*,double*,F_INT*,F_INT*,F_INT*);
135 double SPLINT(double*,F_INT*,double*,F_INT*,double*,double*,double*);
136 void SPROOT(double*,F_INT*,double*,double*,F_INT*,F_INT*,F_INT*);
137 void PARCUR(F_INT*,F_INT*,F_INT*,F_INT*,double*,F_INT*,double*,
138         double*,double*,double*,F_INT*,double*,F_INT*,F_INT*,
139         double*,F_INT*,double*,double*,double*,F_INT*,F_INT*,F_INT*);
140 void CLOCUR(F_INT*,F_INT*,F_INT*,F_INT*,double*,F_INT*,double*,
141         double*,F_INT*,double*,F_INT*,F_INT*,double*,F_INT*,
142         double*,double*,double*,F_INT*,F_INT*,F_INT*);
143 void SURFIT(F_INT*,F_INT*,double*,double*,double*,double*,
144         double*,double*,double*,double*,F_INT*,F_INT*,double*,
145         F_INT*,F_INT*,F_INT*,double*,F_INT*,double*,F_INT*,double*,
146         double*,double*,double*,F_INT*,double*,F_INT*,F_INT*,F_INT*,F_INT*);
147 void BISPEV(double*,F_INT*,double*,F_INT*,double*,F_INT*,F_INT*,
148         double*,F_INT*,double*,F_INT*,double*,double*,F_INT*,
149         F_INT*,F_INT*,F_INT*);
150 void PARDER(double*,F_INT*,double*,F_INT*,double*,F_INT*,F_INT*,
151         F_INT*,F_INT*,double*,F_INT*,double*,F_INT*,double*,
152         double*,F_INT*,F_INT*,F_INT*,F_INT*);
153 void INSERT(F_INT*,double*,F_INT*,double*,F_INT*,double*,double*,
154         F_INT*,double*,F_INT*,F_INT*);
155 
156 /* Note that curev, cualde need no interface. */
157 
158 static char doc_bispev[] = " [z,ier] = _bispev(tx,ty,c,kx,ky,x,y,nux,nuy)";
159 static PyObject *
fitpack_bispev(PyObject * dummy,PyObject * args)160 fitpack_bispev(PyObject *dummy, PyObject *args)
161 {
162     F_INT nx, ny, kx, ky, mx, my, lwrk, *iwrk, kwrk, ier, lwa, nux, nuy;
163     npy_intp mxy;
164     double *tx, *ty, *c, *x, *y, *z, *wrk, *wa = NULL;
165     PyArrayObject *ap_x = NULL, *ap_y = NULL, *ap_z = NULL, *ap_tx = NULL;
166     PyArrayObject *ap_ty = NULL, *ap_c = NULL;
167     PyObject *x_py = NULL, *y_py = NULL, *c_py = NULL, *tx_py = NULL, *ty_py = NULL;
168 
169     if (!PyArg_ParseTuple(args, ("OOO" F_INT_PYFMT F_INT_PYFMT "OO" F_INT_PYFMT F_INT_PYFMT),
170                           &tx_py,&ty_py,&c_py,&kx,&ky,&x_py,&y_py,&nux,&nuy)) {
171         return NULL;
172     }
173     ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1);
174     ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, NPY_DOUBLE, 0, 1);
175     ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1);
176     ap_tx = (PyArrayObject *)PyArray_ContiguousFromObject(tx_py, NPY_DOUBLE, 0, 1);
177     ap_ty = (PyArrayObject *)PyArray_ContiguousFromObject(ty_py, NPY_DOUBLE, 0, 1);
178     if (ap_x == NULL
179             || ap_y == NULL
180             || ap_c == NULL
181             || ap_tx == NULL
182             || ap_ty == NULL) {
183         goto fail;
184     }
185     x = (double *) ap_x->data;
186     y = (double *) ap_y->data;
187     c = (double *) ap_c->data;
188     tx = (double *) ap_tx->data;
189     ty = (double *) ap_ty->data;
190     nx = ap_tx->dimensions[0];
191     ny = ap_ty->dimensions[0];
192     mx = ap_x->dimensions[0];
193     my = ap_y->dimensions[0];
194     mxy = mx*my;
195     if (my != 0 && mxy/my != mx) {
196         /* Integer overflow */
197         PyErr_Format(PyExc_RuntimeError,
198                      "Cannot produce output of size %dx%d (size too large)",
199                      mx, my);
200         goto fail;
201     }
202     ap_z = (PyArrayObject *)PyArray_SimpleNew(1,&mxy,NPY_DOUBLE);
203     if (ap_z == NULL) {
204         goto fail;
205     }
206     z = (double *) ap_z->data;
207     if (nux || nuy) {
208         lwrk = mx*(kx + 1 - nux) + my*(ky + 1 - nuy) + (nx - kx - 1)*(ny - ky - 1);
209     }
210     else {
211         lwrk = mx*(kx + 1) + my*(ky + 1);
212     }
213     kwrk = mx + my;
214     lwa = lwrk + kwrk;
215     if ((wa = malloc(lwa*sizeof(double))) == NULL) {
216         PyErr_NoMemory();
217         goto fail;
218     }
219     wrk = wa;
220     iwrk = (int *)(wrk + lwrk);
221     if (nux || nuy) {
222         PARDER(tx, &nx, ty, &ny, c, &kx, &ky, &nux, &nuy, x, &mx, y, &my, z,
223                 wrk, &lwrk, iwrk, &kwrk, &ier);
224     }
225     else {
226         BISPEV(tx, &nx, ty, &ny, c, &kx, &ky, x, &mx, y, &my, z, wrk, &lwrk,
227                 iwrk, &kwrk, &ier);
228     }
229 
230     free(wa);
231     Py_DECREF(ap_x);
232     Py_DECREF(ap_y);
233     Py_DECREF(ap_c);
234     Py_DECREF(ap_tx);
235     Py_DECREF(ap_ty);
236     return Py_BuildValue("Ni",PyArray_Return(ap_z),ier);
237 
238 fail:
239     free(wa);
240     Py_XDECREF(ap_x);
241     Py_XDECREF(ap_y);
242     Py_XDECREF(ap_z);
243     Py_XDECREF(ap_c);
244     Py_XDECREF(ap_tx);
245     Py_XDECREF(ap_ty);
246     return NULL;
247 }
248 
249 static char doc_surfit[] = " [tx,ty,c,o] = _surfit(x, y, z, w, xb, xe, yb, ye,"\
250       " kx,ky,iopt,s,eps,tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)";
251 static PyObject *
fitpack_surfit(PyObject * dummy,PyObject * args)252 fitpack_surfit(PyObject *dummy, PyObject *args)
253 {
254     F_INT iopt, m, kx, ky, nxest, nyest, lwrk1, lwrk2, *iwrk, kwrk, ier;
255     F_INT lwa, nxo, nyo, i, lcest, nmax, nx, ny, lc;
256     npy_intp dims[1];
257     double *x, *y, *z, *w, xb, xe, yb, ye, s, *tx, *ty, *c, fp;
258     double *wrk1, *wrk2, *wa = NULL, eps;
259     PyArrayObject *ap_x = NULL, *ap_y = NULL, *ap_z, *ap_w = NULL;
260     PyArrayObject *ap_tx = NULL,*ap_ty = NULL,*ap_c = NULL, *ap_wrk = NULL;
261     PyObject *x_py = NULL, *y_py = NULL, *z_py = NULL, *w_py = NULL;
262     PyObject *tx_py = NULL, *ty_py = NULL, *wrk_py = NULL;
263 
264     nx = ny = ier = nxo = nyo = 0;
265     if (!PyArg_ParseTuple(args, ("OOOOdddd" F_INT_PYFMT F_INT_PYFMT F_INT_PYFMT
266                                  "ddOO" F_INT_PYFMT F_INT_PYFMT "O"
267                                  F_INT_PYFMT F_INT_PYFMT),
268                 &x_py, &y_py, &z_py, &w_py, &xb, &xe, &yb, &ye,
269                 &kx, &ky, &iopt, &s, &eps, &tx_py, &ty_py, &nxest,
270                 &nyest, &wrk_py, &lwrk1, &lwrk2)) {
271         return NULL;
272     }
273     ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1);
274     ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, NPY_DOUBLE, 0, 1);
275     ap_z = (PyArrayObject *)PyArray_ContiguousFromObject(z_py, NPY_DOUBLE, 0, 1);
276     ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, NPY_DOUBLE, 0, 1);
277     ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, NPY_DOUBLE, 0, 1);
278     /*ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, F_INT_NPY, 0, 1);*/
279     if (ap_x == NULL
280             || ap_y == NULL
281             || ap_z == NULL
282             || ap_w == NULL
283             || ap_wrk == NULL) {
284         goto fail;
285     }
286     x = (double *) ap_x->data;
287     y = (double *) ap_y->data;
288     z = (double *) ap_z->data;
289     w = (double *) ap_w->data;
290     m = ap_x->dimensions[0];
291     nmax = nxest;
292     if (nmax < nyest) {
293         nmax = nyest;
294     }
295     lcest = (nxest - kx - 1)*(nyest - ky - 1);
296     kwrk = m + (nxest - 2*kx - 1)*(nyest - 2*ky - 1);
297     lwa = 2*nmax + lcest + lwrk1 + lwrk2 + kwrk;
298     if ((wa = malloc(lwa*sizeof(double))) == NULL) {
299         PyErr_NoMemory();
300         goto fail;
301     }
302     /*
303      * NOTE: the work arrays MUST be aligned on double boundaries, as Fortran
304      *       compilers (e.g. gfortran) may assume that.  Malloc gives suitable
305      *       alignment, so we are here careful not to mess it up.
306      */
307     tx = wa;
308     ty = tx + nmax;
309     c = ty + nmax;
310     wrk1 = c + lcest;
311     iwrk = (F_INT *)(wrk1 + lwrk1);
312     wrk2 = ((double *)iwrk) + kwrk;
313     if (iopt) {
314         ap_tx = (PyArrayObject *)PyArray_ContiguousFromObject(tx_py, NPY_DOUBLE, 0, 1);
315         ap_ty = (PyArrayObject *)PyArray_ContiguousFromObject(ty_py, NPY_DOUBLE, 0, 1);
316         if (ap_tx == NULL || ap_ty == NULL) {
317             goto fail;
318         }
319         nx = nxo = ap_tx->dimensions[0];
320         ny = nyo = ap_ty->dimensions[0];
321         memcpy(tx,ap_tx->data,nx*sizeof(double));
322         memcpy(ty,ap_ty->data,ny*sizeof(double));
323     }
324     if (iopt==1) {
325         lc = (nx - kx - 1)*(ny - ky - 1);
326         memcpy(wrk1, ap_wrk->data, lc*sizeof(double));
327         /*memcpy(iwrk,ap_iwrk->data,n*sizeof(int));*/
328     }
329     SURFIT(&iopt, &m, x, y, z, w, &xb, &xe, &yb, &ye, &kx, &ky,
330             &s, &nxest, &nyest, &nmax, &eps, &nx, tx, &ny, ty,
331             c, &fp, wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, &ier);
332     i = 0;
333     while ((ier > 10) && (i++ < 5)) {
334         lwrk2 = ier;
335         if ((wrk2 = malloc(lwrk2*sizeof(double))) == NULL) {
336             PyErr_NoMemory();
337             goto fail;
338         }
339         SURFIT(&iopt, &m, x, y, z, w, &xb, &xe, &yb, &ye, &kx, &ky,
340                 &s, &nxest, &nyest, &nmax, &eps, &nx, tx, &ny, ty,
341                 c, &fp, wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, &ier);
342         free(wrk2);
343     }
344     if (ier == 10) {
345         PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
346         goto fail;
347     }
348     lc = (nx - kx - 1)*(ny - ky - 1);
349     Py_XDECREF(ap_tx);
350     Py_XDECREF(ap_ty);
351     dims[0] = nx;
352     ap_tx = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
353     dims[0] = ny;
354     ap_ty = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
355     dims[0] = lc;
356     ap_c = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
357     if (ap_tx == NULL
358             || ap_ty == NULL
359             || ap_c == NULL) {
360         goto fail;
361     }
362     if ((iopt == 0)||(nx > nxo)||(ny > nyo)) {
363         Py_XDECREF(ap_wrk);
364         dims[0] = lc;
365         ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
366         if (ap_wrk == NULL) {
367             goto fail;
368         }
369         /*ap_iwrk = (PyArrayObject *)PyArray_SimpleNew(1,&n,F_INT_NPY);*/
370     }
371     if (ap_wrk->dimensions[0] < lc) {
372         Py_XDECREF(ap_wrk);
373         dims[0] = lc;
374         ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
375         if (ap_wrk == NULL) {
376             goto fail;
377         }
378     }
379     memcpy(ap_tx->data, tx, nx*sizeof(double));
380     memcpy(ap_ty->data, ty, ny*sizeof(double));
381     memcpy(ap_c->data, c, lc*sizeof(double));
382     memcpy(ap_wrk->data, wrk1, lc*sizeof(double));
383     /*memcpy(ap_iwrk->data,iwrk,n*sizeof(int));*/
384     free(wa);
385     Py_DECREF(ap_x);
386     Py_DECREF(ap_y);
387     Py_DECREF(ap_z);
388     Py_DECREF(ap_w);
389     return Py_BuildValue("NNN{s:N,s:i,s:d}",PyArray_Return(ap_tx),
390             PyArray_Return(ap_ty),PyArray_Return(ap_c),
391             "wrk",PyArray_Return(ap_wrk),
392             "ier",ier,"fp",fp);
393 
394 fail:
395     free(wa);
396     Py_XDECREF(ap_x);
397     Py_XDECREF(ap_y);
398     Py_XDECREF(ap_z);
399     Py_XDECREF(ap_w);
400     Py_XDECREF(ap_tx);
401     Py_XDECREF(ap_ty);
402     Py_XDECREF(ap_wrk);
403     /*Py_XDECREF(ap_iwrk);*/
404     if (!PyErr_Occurred()) {
405         PyErr_SetString(PyExc_ValueError,
406                 "An error occurred.");
407     }
408     return NULL;
409 }
410 
411 
412 static char doc_parcur[] = " [t,c,o] = _parcur(x,w,u,ub,ue,k,iopt,ipar,s,t,nest,wrk,iwrk,per)";
413 static PyObject *
fitpack_parcur(PyObject * dummy,PyObject * args)414 fitpack_parcur(PyObject *dummy, PyObject *args)
415 {
416     F_INT k, iopt, ipar, nest, *iwrk, idim, m, mx, no=0, nc, ier, lwa, lwrk, i, per;
417     F_INT n=0,  lc;
418     npy_intp dims[1];
419     double *x, *w, *u, *c, *t, *wrk, *wa=NULL, ub, ue, fp, s;
420     PyObject *x_py = NULL, *u_py = NULL, *w_py = NULL, *t_py = NULL;
421     PyObject *wrk_py=NULL, *iwrk_py=NULL;
422     PyArrayObject *ap_x = NULL, *ap_u = NULL, *ap_w = NULL, *ap_t = NULL, *ap_c = NULL;
423     PyArrayObject *ap_wrk = NULL, *ap_iwrk = NULL;
424 
425     if (!PyArg_ParseTuple(args, ("OOOdd" F_INT_PYFMT F_INT_PYFMT F_INT_PYFMT
426                                  "dO" F_INT_PYFMT "OO" F_INT_PYFMT),
427                           &x_py, &w_py, &u_py, &ub, &ue, &k, &iopt, &ipar,
428                           &s, &t_py, &nest, &wrk_py, &iwrk_py, &per)) {
429         return NULL;
430     }
431     ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1);
432     ap_u = (PyArrayObject *)PyArray_ContiguousFromObject(u_py, NPY_DOUBLE, 0, 1);
433     ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, NPY_DOUBLE, 0, 1);
434     ap_wrk=(PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, NPY_DOUBLE, 0, 1);
435     ap_iwrk=(PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, F_INT_NPY, 0, 1);
436     if (ap_x == NULL
437             || ap_u == NULL
438             || ap_w == NULL
439             || ap_wrk == NULL
440             || ap_iwrk == NULL) {
441         goto fail;
442     }
443     x = (double *) ap_x->data;
444     u = (double *) ap_u->data;
445     w = (double *) ap_w->data;
446     m = ap_w->dimensions[0];
447     mx = ap_x->dimensions[0];
448     idim = mx/m;
449     if (per) {
450         lwrk = m*(k + 1) + nest*(7 + idim + 5*k);
451     }
452     else {
453         lwrk = m*(k + 1) + nest*(6 + idim + 3*k);
454     }
455     nc = idim*nest;
456     lwa = nc + 2*nest + lwrk;
457     if ((wa = malloc(lwa*sizeof(double))) == NULL) {
458         PyErr_NoMemory();
459         goto fail;
460     }
461     t = wa;
462     c = t + nest;
463     wrk = c + nc;
464     iwrk = (F_INT *)(wrk + lwrk);
465     if (iopt) {
466         ap_t=(PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
467         if (ap_t == NULL) {
468             goto fail;
469         }
470         n = no = ap_t->dimensions[0];
471         memcpy(t, ap_t->data, n*sizeof(double));
472         Py_DECREF(ap_t);
473         ap_t = NULL;
474     }
475     if (iopt == 1) {
476         memcpy(wrk, ap_wrk->data, n*sizeof(double));
477         memcpy(iwrk, ap_iwrk->data, n*sizeof(F_INT));
478     }
479     if (per) {
480         CLOCUR(&iopt, &ipar, &idim, &m, u, &mx, x, w, &k, &s, &nest,
481                 &n, t, &nc, c, &fp, wrk, &lwrk, iwrk, &ier);
482     }
483     else {
484         PARCUR(&iopt, &ipar, &idim, &m, u, &mx, x, w, &ub, &ue, &k,
485                 &s, &nest, &n, t, &nc, c, &fp, wrk, &lwrk, iwrk, &ier);
486     }
487     if (ier == 10) {
488         PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
489         goto fail;
490     }
491     if (ier > 0 && n == 0) {
492         n = 1;
493     }
494     lc = (n - k - 1)*idim;
495     dims[0] = n;
496     ap_t = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
497     dims[0] = lc;
498     ap_c = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
499     if (ap_t == NULL || ap_c == NULL) {
500         goto fail;
501     }
502     if (iopt != 1|| n > no) {
503         Py_XDECREF(ap_wrk);
504         ap_wrk = NULL;
505         Py_XDECREF(ap_iwrk);
506         ap_iwrk = NULL;
507 
508         dims[0] = n;
509         ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
510         if (ap_wrk == NULL) {
511             goto fail;
512         }
513         ap_iwrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, F_INT_NPY);
514         if (ap_iwrk == NULL) {
515             goto fail;
516         }
517     }
518     memcpy(ap_t->data, t, n*sizeof(double));
519     for (i = 0; i < idim; i++)
520         memcpy((double *)ap_c->data + i*(n - k - 1), c + i*n, (n - k - 1)*sizeof(double));
521     memcpy(ap_wrk->data, wrk, n*sizeof(double));
522     memcpy(ap_iwrk->data, iwrk, n*sizeof(F_INT));
523     free(wa);
524     Py_DECREF(ap_x);
525     Py_DECREF(ap_w);
526     return Py_BuildValue(("NN{s:N,s:d,s:d,s:N,s:N,s:" F_INT_PYFMT ",s:d}"), PyArray_Return(ap_t),
527             PyArray_Return(ap_c), "u", PyArray_Return(ap_u), "ub", ub, "ue", ue,
528             "wrk", PyArray_Return(ap_wrk), "iwrk", PyArray_Return(ap_iwrk),
529             "ier", ier, "fp",fp);
530 fail:
531     free(wa);
532     Py_XDECREF(ap_x);
533     Py_XDECREF(ap_u);
534     Py_XDECREF(ap_w);
535     Py_XDECREF(ap_t);
536     Py_XDECREF(ap_wrk);
537     Py_XDECREF(ap_iwrk);
538     return NULL;
539 }
540 
541 static char doc_curfit[] = " [t,c,o] = _curfit(x,y,w,xb,xe,k,iopt,s,t,nest,wrk,iwrk,per)";
542 static PyObject *
fitpack_curfit(PyObject * dummy,PyObject * args)543 fitpack_curfit(PyObject *dummy, PyObject *args)
544 {
545     F_INT iopt, m, k, nest, lwrk, *iwrk, ier, lwa, no=0, per;
546     F_INT n, lc;
547     npy_intp dims[1];
548     double *x, *y, *w, xb, xe, s, *t, *c, fp, *wrk, *wa = NULL;
549     PyArrayObject *ap_x = NULL, *ap_y = NULL, *ap_w = NULL, *ap_t = NULL, *ap_c = NULL;
550     PyArrayObject *ap_wrk = NULL, *ap_iwrk = NULL;
551     PyObject *x_py = NULL, *y_py = NULL, *w_py = NULL, *t_py = NULL;
552     PyObject *wrk_py=NULL, *iwrk_py=NULL;
553 
554     if (!PyArg_ParseTuple(args, ("OOOdd" F_INT_PYFMT F_INT_PYFMT
555                                  "dO" F_INT_PYFMT "OO" F_INT_PYFMT),
556                           &x_py, &y_py, &w_py, &xb, &xe, &k, &iopt,
557                           &s, &t_py, &nest, &wrk_py, &iwrk_py, &per)) {
558         return NULL;
559     }
560     ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1);
561     ap_y = (PyArrayObject *)PyArray_ContiguousFromObject(y_py, NPY_DOUBLE, 0, 1);
562     ap_w = (PyArrayObject *)PyArray_ContiguousFromObject(w_py, NPY_DOUBLE, 0, 1);
563     ap_wrk = (PyArrayObject *)PyArray_ContiguousFromObject(wrk_py, NPY_DOUBLE, 0, 1);
564     ap_iwrk = (PyArrayObject *)PyArray_ContiguousFromObject(iwrk_py, F_INT_NPY, 0, 1);
565     if (ap_x == NULL
566             || ap_y == NULL
567             || ap_w == NULL
568             || ap_wrk == NULL
569             || ap_iwrk == NULL) {
570         goto fail;
571     }
572     x = (double *) ap_x->data;
573     y = (double *) ap_y->data;
574     w = (double *) ap_w->data;
575     m = ap_x->dimensions[0];
576     if (per) {
577         lwrk = m*(k + 1) + nest*(8 + 5*k);
578     }
579     else {
580         lwrk = m*(k + 1) + nest*(7 + 3*k);
581     }
582     lwa = 3*nest + lwrk;
583     if ((wa = malloc(lwa*sizeof(double))) == NULL) {
584         PyErr_NoMemory();
585         goto fail;
586     }
587     t = wa;
588     c = t + nest;
589     wrk = c + nest;
590     iwrk = (F_INT *)(wrk + lwrk);
591     if (iopt) {
592         ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
593         if (ap_t == NULL) {
594             goto fail;
595         }
596         n = no = ap_t->dimensions[0];
597         memcpy(t, ap_t->data, n*sizeof(double));
598     }
599     if (iopt == 1) {
600         memcpy(wrk, ap_wrk->data, n*sizeof(double));
601         memcpy(iwrk, ap_iwrk->data, n*sizeof(F_INT));
602     }
603     if (per)
604         PERCUR(&iopt, &m, x, y, w, &k, &s, &nest, &n, t, c, &fp, wrk,
605                 &lwrk, iwrk, &ier);
606     else {
607         CURFIT(&iopt, &m, x, y, w, &xb, &xe, &k, &s, &nest, &n, t, c,
608                 &fp, wrk, &lwrk, iwrk, &ier);
609     }
610     if (ier == 10) {
611         PyErr_SetString(PyExc_ValueError, "Invalid inputs.");
612         goto fail;
613     }
614     lc = n - k - 1;
615     if (!iopt) {
616         dims[0] = n;
617         ap_t = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
618         if (ap_t == NULL) {
619             goto fail;
620         }
621     }
622     dims[0] = lc;
623     ap_c = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
624     if (ap_c == NULL) {
625         goto fail;
626     }
627     if (iopt == 0 || n > no) {
628         Py_XDECREF(ap_wrk);
629         Py_XDECREF(ap_iwrk);
630         dims[0] = n;
631         ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
632         ap_iwrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, F_INT_NPY);
633         if (ap_wrk == NULL || ap_iwrk == NULL) {
634             goto fail;
635         }
636     }
637     memcpy(ap_t->data, t, n*sizeof(double));
638     memcpy(ap_c->data, c, lc*sizeof(double));
639     memcpy(ap_wrk->data, wrk, n*sizeof(double));
640     memcpy(ap_iwrk->data, iwrk, n*sizeof(F_INT));
641     free(wa);
642     Py_DECREF(ap_x);
643     Py_DECREF(ap_y);
644     Py_DECREF(ap_w);
645     return Py_BuildValue(("NN{s:N,s:N,s:" F_INT_PYFMT ",s:d}"), PyArray_Return(ap_t),
646             PyArray_Return(ap_c), "wrk", PyArray_Return(ap_wrk),
647             "iwrk", PyArray_Return(ap_iwrk), "ier", ier, "fp", fp);
648 
649 fail:
650     free(wa);
651     Py_XDECREF(ap_x);
652     Py_XDECREF(ap_y);
653     Py_XDECREF(ap_w);
654     Py_XDECREF(ap_t);
655     Py_XDECREF(ap_wrk);
656     Py_XDECREF(ap_iwrk);
657     return NULL;
658 }
659 
660 static char doc_spl_[] = " [y,ier] = _spl_(x,nu,t,c,k,e)";
661 static PyObject *
fitpack_spl_(PyObject * dummy,PyObject * args)662 fitpack_spl_(PyObject *dummy, PyObject *args)
663 {
664     F_INT n, nu, ier, k, m, e=0;
665     npy_intp dims[1];
666     double *x, *y, *t, *c, *wrk = NULL;
667     PyArrayObject *ap_x = NULL, *ap_y = NULL, *ap_t = NULL, *ap_c = NULL;
668     PyObject *x_py = NULL, *t_py = NULL, *c_py = NULL;
669 
670     if (!PyArg_ParseTuple(args, ("O" F_INT_PYFMT "OO" F_INT_PYFMT F_INT_PYFMT),
671                           &x_py, &nu, &t_py, &c_py, &k, &e)) {
672         return NULL;
673     }
674     ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, NPY_DOUBLE, 0, 1);
675     ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
676     ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1);
677     if ((ap_x == NULL || ap_t == NULL || ap_c == NULL)) {
678         goto fail;
679     }
680     x = (double *)ap_x->data;
681     m = ap_x->dimensions[0];
682     t = (double *)ap_t->data;
683     c = (double *)ap_c->data;
684     n = ap_t->dimensions[0];
685     dims[0] = m;
686     ap_y = (PyArrayObject *)PyArray_SimpleNew(1,dims,NPY_DOUBLE);
687     if (ap_y == NULL) {
688         goto fail;
689     }
690     y = (double *)ap_y->data;
691     if ((wrk = malloc(n*sizeof(double))) == NULL) {
692         PyErr_NoMemory();
693         goto fail;
694     }
695     if (nu) {
696         SPLDER(t, &n, c, &k, &nu, x, y, &m, &e, wrk, &ier);
697     }
698     else {
699         SPLEV(t, &n, c, &k, x, y, &m, &e, &ier);
700     }
701     free(wrk);
702     Py_DECREF(ap_x);
703     Py_DECREF(ap_c);
704     Py_DECREF(ap_t);
705     return Py_BuildValue(("N" F_INT_PYFMT), PyArray_Return(ap_y), ier);
706 
707 fail:
708     free(wrk);
709     Py_XDECREF(ap_x);
710     Py_XDECREF(ap_c);
711     Py_XDECREF(ap_t);
712     return NULL;
713 }
714 
715 static char doc_splint[] = " [aint,wrk] = _splint(t,c,k,a,b)";
716 static PyObject *
fitpack_splint(PyObject * dummy,PyObject * args)717 fitpack_splint(PyObject *dummy, PyObject *args)
718 {
719     F_INT k, n;
720     npy_intp dims[1];
721     double *t, *c, *wrk = NULL, a, b, aint;
722     PyArrayObject *ap_t = NULL, *ap_c = NULL;
723     PyArrayObject *ap_wrk = NULL;
724     PyObject *t_py = NULL, *c_py = NULL;
725 
726     if (!PyArg_ParseTuple(args, ("OO" F_INT_PYFMT "dd"),&t_py,&c_py,&k,&a,&b)) {
727         return NULL;
728     }
729     ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
730     ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1);
731     if ((ap_t == NULL || ap_c == NULL)) {
732         goto fail;
733     }
734     t = (double *)ap_t->data;
735     c = (double *)ap_c->data;
736     n = ap_t->dimensions[0];
737     dims[0] = n;
738     ap_wrk = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
739     if (ap_wrk == NULL) {
740         goto fail;
741     }
742     wrk = (double *)ap_wrk->data;
743     aint = SPLINT(t,&n,c,&k,&a,&b,wrk);
744     Py_DECREF(ap_c);
745     Py_DECREF(ap_t);
746     return Py_BuildValue("dN", aint, PyArray_Return(ap_wrk));
747 
748 fail:
749     Py_XDECREF(ap_c);
750     Py_XDECREF(ap_t);
751     return NULL;
752 }
753 
754 static char doc_sproot[] = " [z,ier] = _sproot(t,c,k,mest)";
755 static PyObject *
fitpack_sproot(PyObject * dummy,PyObject * args)756 fitpack_sproot(PyObject *dummy, PyObject *args)
757 {
758     F_INT n, k, m, mest, ier;
759     npy_intp dims[1];
760     double *t, *c, *z = NULL;
761     PyArrayObject *ap_t = NULL, *ap_c = NULL;
762     PyArrayObject *ap_z = NULL;
763     PyObject *t_py = NULL, *c_py = NULL;
764 
765     if (!PyArg_ParseTuple(args, ("OO" F_INT_PYFMT F_INT_PYFMT),
766                           &t_py,&c_py,&k,&mest)) {
767         return NULL;
768     }
769     ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
770     ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1);
771     if ((ap_t == NULL || ap_c == NULL)) {
772         goto fail;
773     }
774     t = (double *)ap_t->data;
775     c = (double *)ap_c->data;
776     n = ap_t->dimensions[0];
777     if ((z = malloc(mest*sizeof(double))) == NULL) {
778         PyErr_NoMemory();
779         goto fail;
780     }
781     m = 0;
782     SPROOT(t,&n,c,z,&mest,&m,&ier);
783     if (ier==10) {
784         m = 0;
785     }
786     dims[0] = m;
787     ap_z = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
788     if (ap_z == NULL) {
789         goto fail;
790     }
791     memcpy(ap_z->data,z,m*sizeof(double));
792     free(z);
793     Py_DECREF(ap_c);
794     Py_DECREF(ap_t);
795     return Py_BuildValue(("N" F_INT_PYFMT), PyArray_Return(ap_z), ier);
796 
797 fail:
798     free(z);
799     Py_XDECREF(ap_c);
800     Py_XDECREF(ap_t);
801     return NULL;
802 }
803 
804 static char doc_spalde[] = " [d,ier] = _spalde(t,c,k,x)";
805 static PyObject *
fitpack_spalde(PyObject * dummy,PyObject * args)806 fitpack_spalde(PyObject *dummy, PyObject *args)
807 {
808     F_INT n, k, ier, k1;
809     npy_intp dims[1];
810     double *t, *c, *d = NULL, x;
811     PyArrayObject *ap_t = NULL, *ap_c = NULL, *ap_d = NULL;
812     PyObject *t_py = NULL, *c_py = NULL;
813 
814     if (!PyArg_ParseTuple(args, ("OO" F_INT_PYFMT "d"),
815                           &t_py,&c_py,&k,&x)) {
816         return NULL;
817     }
818     ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
819     ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1);
820     if ((ap_t == NULL || ap_c == NULL)) {
821         goto fail;
822     }
823     t = (double *)ap_t->data;
824     c = (double *)ap_c->data;
825     n = ap_t->dimensions[0];
826     k1 = k + 1;
827     dims[0] = k1;
828     ap_d = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
829     if (ap_d == NULL) {
830         goto fail;
831     }
832     d = (double *)ap_d->data;
833     SPALDE(t, &n, c, &k1, &x, d, &ier);
834     Py_DECREF(ap_c);
835     Py_DECREF(ap_t);
836     return Py_BuildValue(("N" F_INT_PYFMT), PyArray_Return(ap_d), ier);
837 
838 fail:
839     Py_XDECREF(ap_c);
840     Py_XDECREF(ap_t);
841     return NULL;
842 }
843 
844 static char doc_insert[] = " [tt,cc,ier] = _insert(iopt,t,c,k,x,m)";
845 static PyObject *
fitpack_insert(PyObject * dummy,PyObject * args)846 fitpack_insert(PyObject *dummy, PyObject*args)
847 {
848     F_INT iopt, n, nn, k, ier, m, nest;
849     npy_intp dims[1];
850     double x;
851     double *t_in, *c_in, *t_out, *c_out, *t_buf = NULL, *c_buf = NULL, *p;
852     double *t1, *t2, *c1, *c2;
853     PyArrayObject *ap_t_in = NULL, *ap_c_in = NULL, *ap_t_out = NULL, *ap_c_out = NULL;
854     PyObject *t_py = NULL, *c_py = NULL;
855     PyObject *ret = NULL;
856 
857     if (!PyArg_ParseTuple(args,(F_INT_PYFMT "OO" F_INT_PYFMT "d" F_INT_PYFMT),
858                           &iopt,&t_py,&c_py,&k, &x, &m)) {
859         return NULL;
860     }
861     ap_t_in = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, NPY_DOUBLE, 0, 1);
862     ap_c_in = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, NPY_DOUBLE, 0, 1);
863     if (ap_t_in == NULL || ap_c_in == NULL) {
864         goto fail;
865     }
866     t_in = (double *)ap_t_in->data;
867     c_in = (double *)ap_c_in->data;
868     n = ap_t_in->dimensions[0];
869     nest = n + m;
870     dims[0] = nest;
871     ap_t_out = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
872     ap_c_out = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
873     if (ap_t_out == NULL || ap_c_out == NULL) {
874         goto fail;
875     }
876     t_out = (double *)ap_t_out->data;
877     c_out = (double *)ap_c_out->data;
878 
879     /*
880      * Call the INSERT routine m times to insert m-multiplicity knot, ie.:
881      *
882      *     for _ in range(n, nest):
883      *         t, c = INSERT(t, c)
884      *     return t, c
885      *
886      * We need to ensure that input and output buffers given to INSERT routine
887      * do not point to same memory, which is not allowed by Fortran. For this,
888      * we use temporary storage, and cycle between it and the output buffers.
889      */
890     t2 = t_in;
891     c2 = c_in;
892     t1 = t_out;
893     c1 = c_out;
894 
895     for ( ; n < nest; n++) {
896         /* Swap buffers */
897         p = t2; t2 = t1; t1 = p;
898         p = c2; c2 = c1; c1 = p;
899 
900         /* Allocate temporary buffer (needed for m > 1) */
901         if (t2 == t_in) {
902             if (t_buf == NULL) {
903                 t_buf = calloc(nest, sizeof(double));
904                 c_buf = calloc(nest, sizeof(double));
905                 if (t_buf == NULL || c_buf == NULL) {
906                     PyErr_NoMemory();
907                     goto fail;
908                 }
909             }
910             t2 = t_buf;
911             c2 = c_buf;
912         }
913 
914         INSERT(&iopt, t1, &n, c1, &k, &x, t2, &nn, c2, &nest, &ier);
915 
916         if (ier) {
917             break;
918         }
919     }
920 
921     /* Ensure output ends up in output buffers */
922     if (t2 != t_out) {
923         memcpy(t_out, t2, nest * sizeof(double));
924         memcpy(c_out, c2, nest * sizeof(double));
925     }
926 
927     Py_DECREF(ap_c_in);
928     Py_DECREF(ap_t_in);
929     free(t_buf);
930     free(c_buf);
931     ret = Py_BuildValue(("NN" F_INT_PYFMT), PyArray_Return(ap_t_out), PyArray_Return(ap_c_out), ier);
932     return ret;
933 
934 fail:
935     Py_XDECREF(ap_c_out);
936     Py_XDECREF(ap_t_out);
937     Py_XDECREF(ap_c_in);
938     Py_XDECREF(ap_t_in);
939     free(t_buf);
940     free(c_buf);
941     return NULL;
942 }
943 
944 
945 /*
946  * Given a set of (N+1) samples:  A default set of knots is constructed
947  * using the samples xk plus 2*(K-1) additional knots where
948  * K = max(order,1) and the knots are chosen so that distances
949  * are symmetric around the first and last samples: x_0 and x_N.
950  *
951  * There should be a vector of N+K coefficients for the spline
952  * curve in coef.  These coefficients form the curve as
953  *
954  * s(x) = sum(c_j B_{j,K}(x), j=-K..N-1)
955  *
956  * The spline function is evaluated at all points xx.
957  * The approximation interval is from xk[0] to xk[-1]
958  * Any xx outside that interval is set automatically to 0.0
959  */
960 static char doc_bspleval[] = "y = _bspleval(xx,xk,coef,k,{deriv (0)})\n"
961 "\n"
962 "The spline is defined by the approximation interval xk[0] to xk[-1],\n"
963 "the length of xk (N+1), the order of the spline, k, and \n"
964 "the number of coeficients N+k.  The coefficients range from xk_{-K}\n"
965 "to xk_{N-1} inclusive and are all the coefficients needed to define\n"
966 "an arbitrary spline of order k, on the given approximation interval\n"
967 "\n"
968 "Extra knot points are internally added using knot-point symmetry \n"
969 "around xk[0] and xk[-1]";
970 
971 static PyObject *
_bspleval(PyObject * dummy,PyObject * args)972 _bspleval(PyObject *dummy, PyObject *args)
973 {
974     int k, kk, N, i, ell, dk, deriv = 0;
975     PyObject *xx_py = NULL, *coef_py = NULL, *x_i_py = NULL;
976     PyArrayObject *xx = NULL, *coef = NULL, *x_i = NULL, *yy = NULL;
977     PyArrayIterObject *xx_iter;
978     double *t = NULL, *h = NULL, *ptr;
979     double x0, xN, xN1, arg, sp, cval;
980 
981     if (!PyArg_ParseTuple(args, "OOOi|i", &xx_py, &x_i_py, &coef_py, &k, &deriv)) {
982         return NULL;
983     }
984     if (k < 0) {
985         PyErr_Format(PyExc_ValueError,
986                 "order (%d) must be >=0", k);
987         return NULL;
988     }
989     if (deriv > k) {
990         PyErr_Format(PyExc_ValueError,
991                 "derivative (%d) must be <= order (%d)", deriv, k);
992         return NULL;
993     }
994     kk = k;
995     if (k == 0) {
996         kk = 1;
997     }
998     dk = (k == 0 ? 0 : 1);
999     x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
1000     coef = (PyArrayObject *)PyArray_FROMANY(coef_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
1001     xx = (PyArrayObject *)PyArray_FROMANY(xx_py, NPY_DOUBLE, 0, 0, NPY_ALIGNED);
1002     if (x_i == NULL || coef == NULL || xx == NULL) {
1003         goto fail;
1004     }
1005 
1006     N = PyArray_DIM(x_i, 0) - 1;
1007     if (PyArray_DIM(coef, 0) < N + k) {
1008         PyErr_Format(PyExc_ValueError,
1009                 "too few coefficients (have %d need at least %d)",
1010                 (int)PyArray_DIM(coef, 0), N + k);
1011         goto fail;
1012     }
1013 
1014     /* create output values */
1015     yy = (PyArrayObject *)PyArray_EMPTY(xx->nd, xx->dimensions, NPY_DOUBLE, 0);
1016     if (yy == NULL) {
1017         goto fail;
1018     }
1019     /*
1020      * create dummy knot array with new knots inserted at the end
1021      * selected as mirror symmetric versions of the old knots
1022      */
1023     t = malloc(sizeof(double)*(N + 2*kk - 1));
1024     if (t == NULL) {
1025         PyErr_NoMemory();
1026         goto fail;
1027     }
1028     x0 = *((double *)PyArray_DATA(x_i));
1029     xN = *((double *)PyArray_DATA(x_i) + N);
1030     for (i = 0; i < kk - 1; i++) { /* fill in ends if kk > 1*/
1031         t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i, kk - 1 - i)));
1032         t[kk+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i, N - 1 - i)));
1033     }
1034     ptr = t + (kk - 1);
1035     for (i = 0; i <= N; i++) {
1036         *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)));
1037     }
1038 
1039     /*
1040      * Create work array to hold computed non-zero values for
1041      * the spline for a value of x.
1042      */
1043     h = malloc(sizeof(double)*(2*kk+1));
1044     if (h == NULL) {
1045         PyErr_NoMemory();
1046         goto fail;
1047     }
1048 
1049     /* Determine the spline for each value of x */
1050     xx_iter = (PyArrayIterObject *)PyArray_IterNew((PyObject *)xx);
1051     if (xx_iter == NULL) {
1052         goto fail;
1053     }
1054     ptr = PyArray_DATA(yy);
1055 
1056     while(PyArray_ITER_NOTDONE(xx_iter)) {
1057         arg = *((double *)PyArray_ITER_DATA(xx_iter));
1058         if ((arg < x0) || (arg > xN)) {
1059             /*
1060              * If we are outside the interpolation region,
1061              * fill with zeros
1062              */
1063             *ptr++ = 0.0;
1064         }
1065         else {
1066             /*
1067              * Find the interval that arg lies between in the set of knots
1068              * t[ell] <= arg < t[ell+1] (last-knot use the previous interval)
1069              */
1070             xN1 = *((double *)PyArray_DATA(x_i) + N-1);
1071             if (arg >= xN1) {
1072                 ell = N + kk - 2;
1073             }
1074             else {
1075                 ell = kk - 1;
1076                 while ((arg > t[ell])) {
1077                     ell++;
1078                 }
1079                 if (arg != t[ell]) {
1080                     ell--;
1081                 }
1082             }
1083 
1084             _deBoor_D(t, arg, k, ell, deriv, h);
1085 
1086             sp = 0.0;
1087             for (i = 0; i <= k; i++) {
1088                 cval = *((double *)(PyArray_GETPTR1(coef, ell - i + dk)));
1089                 sp += cval * h[k - i];
1090             }
1091             *ptr++ = sp;
1092         }
1093         PyArray_ITER_NEXT(xx_iter);
1094     }
1095     Py_DECREF(xx_iter);
1096     Py_DECREF(x_i);
1097     Py_DECREF(coef);
1098     Py_DECREF(xx);
1099     free(t);
1100     free(h);
1101     return PyArray_Return(yy);
1102 
1103 fail:
1104     Py_XDECREF(xx);
1105     Py_XDECREF(coef);
1106     Py_XDECREF(x_i);
1107     Py_XDECREF(yy);
1108     free(t);
1109     free(h);
1110     return NULL;
1111 }
1112 
1113 
1114 /*
1115  * Given a set of (N+1) sample positions:
1116  * Construct the diagonals of the (N+1) x (N+K) matrix that is needed to find
1117  * the coefficients of a spline fit of order K.
1118  * Note that K>=2 because for K=0,1, the coefficients are just the
1119  * sample values themselves.
1120  *
1121  * The equation that expresses the constraints is
1122  *
1123  * s(x_i) = sum(c_j B_{j,K}(x_i), j=-K..N-1) = w_i   for i=0..N
1124  *
1125  * This is equivalent to
1126  *
1127  * w = B*c   where c.T = [c_{-K}, c{-K+1}, ..., c_{N-1}] and
1128  *
1129  * Therefore B is an (N+1) times (N+K) matrix with entries
1130  *
1131  * B_{j,K}(x_i)  for column j=-K..N-1
1132  * and row i=0..N
1133  *
1134  * This routine takes the N+1 sample positions and the order k and
1135  * constructs the banded constraint matrix B (with k non-zero diagonals)
1136  *
1137  * The returned array is (N+1) times (N+K) ready to be either used
1138  * to compute a minimally Kth-order derivative discontinuous spline
1139  * or to be expanded with an additional K-1 constraints to be used in
1140  * an exact spline specification.
1141  */
1142 static char doc_bsplmat[] = "B = _bsplmat(order,xk)\n"
1143 "Construct the constraint matrix for spline fitting of order k\n"
1144 "given sample positions in xk.\n"
1145 "\n"
1146 "If xk is an integer (N+1), then the result is equivalent to\n"
1147 "xk=arange(N+1)+x0 for any value of x0.   This produces the\n"
1148 "integer-spaced, or cardinal spline matrix a bit faster.";
1149 static PyObject *
_bsplmat(PyObject * dummy,PyObject * args)1150 _bsplmat(PyObject *dummy, PyObject *args) {
1151     int k, N, i, numbytes, j, equal;
1152     npy_intp dims[2];
1153     PyObject *x_i_py = NULL;
1154     PyArrayObject *x_i = NULL, *BB = NULL;
1155     double *t = NULL, *h = NULL, *ptr;
1156     double x0, xN, arg;
1157 
1158     if (!PyArg_ParseTuple(args, "iO", &k, &x_i_py)) {
1159         return NULL;
1160     }
1161     if (k < 2) {
1162         PyErr_Format(PyExc_ValueError, "order (%d) must be >=2", k);
1163         return NULL;
1164     }
1165 
1166     equal = 0;
1167     N = PySequence_Length(x_i_py);
1168     if (N == -1 && PyErr_Occurred()) {
1169         PyErr_Clear();
1170         N = PyInt_AsLong(x_i_py);
1171         if (N == -1 && PyErr_Occurred()) {
1172             goto fail;
1173         }
1174         equal = 1;
1175     }
1176     N -= 1;
1177 
1178     /* create output matrix */
1179     dims[0] = N + 1;
1180     dims[1] = N + k;
1181     BB = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
1182     if (BB == NULL) {
1183         goto fail;
1184     }
1185 
1186     t = malloc(sizeof(double)*(N + 2*k - 1));
1187     if (t == NULL) {
1188         PyErr_NoMemory();
1189         goto fail;
1190     }
1191 
1192     /*
1193      * Create work array to hold computed non-zero values for
1194      * the spline for a value of x.
1195      */
1196     h = malloc(sizeof(double)*(2*k + 1));
1197     if (h == NULL) {
1198         PyErr_NoMemory();
1199         goto fail;
1200     }
1201 
1202     numbytes = k*sizeof(double);
1203 
1204     if (equal) {
1205         /*
1206          * points equally spaced by 1
1207          * we run deBoor's algorithm one time with artificially created knots
1208          * Then, we keep copying the result to every row
1209          */
1210 
1211         /* Create knots at equally-spaced locations from -(K-1) to N+K-1 */
1212         ptr = t;
1213         for (i = -k + 1; i < N + k; i++) {
1214             *ptr++ = i;
1215         }
1216         j = k - 1;
1217         _deBoor_D(t, 0, k, k-1, 0, h);
1218         ptr = PyArray_DATA(BB);
1219         N = N+1;
1220         for (i = 0; i < N; i++) {
1221             memcpy(ptr, h, numbytes);
1222             ptr += N + k;
1223         }
1224         goto finish;
1225     }
1226 
1227     /* Not-equally spaced */
1228     x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
1229     if (x_i == NULL) {
1230         goto fail;
1231     }
1232     /*
1233      * create dummy knot array with new knots inserted at the end
1234      * selected as mirror symmetric versions of the old knots
1235      */
1236     x0 = *((double *)PyArray_DATA(x_i));
1237     xN = *((double *)PyArray_DATA(x_i) + N);
1238     for (i = 0; i < k - 1; i++) {
1239         /* fill in ends if k > 1*/
1240         t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i, k - 1 - i)));
1241         t[k+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i, N - 1 - i)));
1242     }
1243     ptr = t + (k - 1);
1244     for (i = 0; i <= N; i++) {
1245         *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)));
1246     }
1247 
1248 
1249     /*
1250      * Determine the K+1 non-zero values of the spline and place them in the
1251      * correct location in the matrix for each row (along the diagonals).
1252      * In fact, the last member is always zero so only K non-zero values
1253      * are present.
1254      */
1255     ptr = PyArray_DATA(BB);
1256     for (i = 0, j = k - 1; i < N; i++, j++) {
1257         arg = *((double *)PyArray_DATA(x_i) + i);
1258         _deBoor_D(t, arg, k, j, 0, h);
1259         memcpy(ptr, h, numbytes);
1260         /* advance to next row shifted over one */
1261         ptr += (N + k + 1);
1262     }
1263     /* Last row is different the first coefficient is zero.*/
1264     _deBoor_D(t, xN, k, j - 1, 0, h);
1265     memcpy(ptr, h + 1, numbytes);
1266 
1267 finish:
1268     Py_XDECREF(x_i);
1269     free(t);
1270     free(h);
1271     return (PyObject *)BB;
1272 
1273 fail:
1274     Py_XDECREF(x_i);
1275     Py_XDECREF(BB);
1276     free(t);
1277     free(h);
1278     return NULL;
1279 }
1280 
1281 
1282 
1283 /*
1284  * Given a set of (N+1) sample positions:
1285  * Construct the (N-1) x (N+K) error matrix J_{ij} such that
1286  *
1287  * for i=1..N-1,
1288  *
1289  * e_i = sum(J_{ij}c_{j},j=-K..N-1)
1290  *
1291  * is the discontinuity of the Kth derivative at the point i in the spline.
1292  *
1293  * This routine takes the N+1 sample positions and the order k and
1294  * constructs the banded matrix J
1295  *
1296  * The returned array is (N+1) times (N+K) ready to be either used
1297  * to compute a minimally Kth-order derivative discontinuous spline
1298  * or to be expanded with an additional K-1 constraints to be used in
1299  * an exact reconstruction approach.
1300  */
1301 static char doc_bspldismat[] = "B = _bspldismat(order,xk)\n"
1302 "Construct the kth derivative discontinuity jump constraint matrix \n"
1303 "for spline fitting of order k given sample positions in xk.\n"
1304 "\n"
1305 "If xk is an integer (N+1), then the result is equivalent to\n"
1306 "xk=arange(N+1)+x0 for any value of x0.   This produces the\n"
1307 "integer-spaced matrix a bit faster.  If xk is a 2-tuple (N+1,dx)\n"
1308 "then it produces the result as if the sample distance were dx";
1309 static PyObject *
_bspldismat(PyObject * dummy,PyObject * args)1310 _bspldismat(PyObject *dummy, PyObject *args)
1311 {
1312     int k, N, i, j, equal, m;
1313     npy_intp dims[2];
1314     PyObject *x_i_py = NULL;
1315     PyArrayObject *x_i = NULL, *BB = NULL;
1316     double *t = NULL, *h = NULL, *ptr, *dptr;
1317     double x0, xN, dx;
1318 
1319     if (!PyArg_ParseTuple(args, "iO", &k, &x_i_py)) {
1320         return NULL;
1321     }
1322     if (k < 2) {
1323         PyErr_Format(PyExc_ValueError, "order (%d) must be >=2", k);
1324         return NULL;
1325     }
1326 
1327     equal = 0;
1328     N = PySequence_Length(x_i_py);
1329     if (N == 2 || (N == -1 && PyErr_Occurred())) {
1330         PyErr_Clear();
1331         if (PyTuple_Check(x_i_py)) {
1332             /* x_i_py = (N+1, dx) */
1333             N = PyInt_AsLong(PyTuple_GET_ITEM(x_i_py, 0));
1334             dx = PyFloat_AsDouble(PyTuple_GET_ITEM(x_i_py, 1));
1335         }
1336         else {
1337             N = PyInt_AsLong(x_i_py);
1338             if (N == -1 && PyErr_Occurred()) {
1339                 goto fail;
1340             }
1341             dx = 1.0;
1342         }
1343         equal = 1;
1344     }
1345     N -= 1;
1346 
1347     if (N < 2) {
1348         PyErr_Format(PyExc_ValueError, "too few samples (%d)", N);
1349         return NULL;
1350     }
1351     /* create output matrix */
1352     dims[0] = N - 1;
1353     dims[1] = N + k;
1354     BB = (PyArrayObject *)PyArray_ZEROS(2, dims, NPY_DOUBLE, 0);
1355     if (BB == NULL) {
1356         goto fail;
1357     }
1358     t = malloc(sizeof(double)*(N+2*k-1));
1359     if (t == NULL) {
1360         PyErr_NoMemory();
1361         goto fail;
1362     }
1363 
1364     /*
1365      * Create work array to hold computed non-zero values for
1366      * the spline for a value of x.
1367      */
1368     h = malloc(sizeof(double)*(2*k + 1));
1369     if (h == NULL) {
1370         PyErr_NoMemory();
1371         goto fail;
1372     }
1373 
1374     if (equal) {
1375         /*
1376          * points equally spaced by 1
1377          * we run deBoor's full derivative algorithm twice, subtract the results
1378          * offset by one and then copy the result one time with artificially created knots
1379          * Then, we keep copying the result to every row
1380          */
1381 
1382         /* Create knots at equally-spaced locations from -(K-1) to N+K-1 */
1383         double *tmp, factor;
1384         int numbytes;
1385         numbytes = (k + 2)*sizeof(double);
1386         tmp = malloc(numbytes);
1387         if (tmp == NULL) {
1388             PyErr_NoMemory();
1389             goto fail;
1390         }
1391         ptr = t;
1392         for (i = -k + 1; i < N + k; i++) {
1393             *ptr++ = i;
1394         }
1395         j = k - 1;
1396         _deBoor_D(t, 0, k, k-1, k, h);
1397         ptr = tmp;
1398         for (m = 0; m <= k; m++) {
1399             *ptr++ = -h[m];
1400         }
1401         _deBoor_D(t, 0, k, k, k, h);
1402         ptr = tmp + 1;
1403         for (m = 0; m <= k; m++) {
1404             *ptr++ += h[m];
1405         }
1406         if (dx != 1.0) {
1407             factor = pow(dx, (double)k);
1408             for (m = 0; m < k + 2; m++) {
1409                 tmp[m] /= factor;
1410             }
1411         }
1412         ptr = PyArray_DATA(BB);
1413         for (i = 0; i < N - 1; i++) {
1414             memcpy(ptr, tmp, numbytes);
1415             ptr += N + k + 1;
1416         }
1417         free(tmp);
1418         goto finish;
1419     }
1420 
1421     /* Not-equally spaced */
1422     x_i = (PyArrayObject *)PyArray_FROMANY(x_i_py, NPY_DOUBLE, 1, 1, NPY_ALIGNED);
1423     if (x_i == NULL) {
1424         goto fail;
1425     }
1426     /*
1427      * create dummy knot array with new knots inserted at the end
1428      * selected as mirror symmetric versions of the old knots
1429      */
1430     x0 = *((double *)PyArray_DATA(x_i));
1431     xN = *((double *)PyArray_DATA(x_i) + N);
1432     for (i = 0; i < k - 1; i++) {
1433         /* fill in ends if k > 1*/
1434         t[i] = 2*x0 - *((double *)(PyArray_GETPTR1(x_i, k - 1 - i)));
1435         t[k+N+i] = 2*xN - *((double *)(PyArray_GETPTR1(x_i, N - 1 - i)));
1436     }
1437     ptr = t + (k - 1);
1438     for (i = 0; i <= N; i++) {
1439         *ptr++ = *((double *)(PyArray_GETPTR1(x_i, i)));
1440     }
1441 
1442 
1443     /*
1444      * Determine the K+1 non-zero values of the discontinuity jump matrix
1445      * and place them in the correct location in the matrix for each row
1446      * (along the diagonals).
1447      *
1448      * The matrix is
1449      *
1450      * J_{ij} = b^{(k)}_{j,k}(x^{+}_i) - b^{(k)}_{j,k}(x^{-}_i)
1451      */
1452     ptr = PyArray_DATA(BB);
1453     dptr = ptr;
1454     for (i = 0, j = k - 1; i < N - 1; i++, j++) {
1455         _deBoor_D(t, 0, k, j, k, h);
1456         /* We need to copy over but negate the terms */
1457         for (m = 0; m <= k; m++) {
1458             *ptr++ = -h[m];
1459         }
1460         /*
1461          * If we are past the first row, then we need to also add the current
1462          * values result to the previous row
1463          */
1464         if (i > 0) {
1465             for (m = 0; m <= k; m++) {
1466                 *dptr++ += h[m];
1467             }
1468         }
1469         /* store location of last start position plus one.*/
1470         dptr = ptr - k;
1471         /* advance to next row shifted over one */
1472         ptr += N;
1473     }
1474     /* We need to finish the result for the last row. */
1475     _deBoor_D(t, 0, k, j, k, h);
1476     for (m = 0; m <= k; m++) {
1477         *dptr++ += h[m];
1478     }
1479 
1480 finish:
1481     Py_XDECREF(x_i);
1482     free(t);
1483     free(h);
1484     return (PyObject *)BB;
1485 
1486 fail:
1487     Py_XDECREF(x_i);
1488     Py_XDECREF(BB);
1489     free(t);
1490     free(h);
1491     return NULL;
1492 }
1493 
1494 /* End of functions moved verbatim from __fitpack.h */
1495 
1496 
1497 
1498 static struct PyMethodDef fitpack_module_methods[] = {
1499 {"_curfit",
1500     fitpack_curfit,
1501     METH_VARARGS, doc_curfit},
1502 {"_spl_",
1503     fitpack_spl_,
1504     METH_VARARGS, doc_spl_},
1505 {"_splint",
1506     fitpack_splint,
1507     METH_VARARGS, doc_splint},
1508 {"_sproot",
1509     fitpack_sproot,
1510     METH_VARARGS, doc_sproot},
1511 {"_spalde",
1512     fitpack_spalde,
1513     METH_VARARGS, doc_spalde},
1514 {"_parcur",
1515     fitpack_parcur,
1516     METH_VARARGS, doc_parcur},
1517 {"_surfit",
1518     fitpack_surfit,
1519     METH_VARARGS, doc_surfit},
1520 {"_bispev",
1521     fitpack_bispev,
1522     METH_VARARGS, doc_bispev},
1523 {"_insert",
1524     fitpack_insert,
1525     METH_VARARGS, doc_insert},
1526 {"_bspleval",
1527     _bspleval,
1528     METH_VARARGS, doc_bspleval},
1529 {"_bsplmat",
1530     _bsplmat,
1531     METH_VARARGS, doc_bsplmat},
1532 {"_bspldismat",
1533     _bspldismat,
1534     METH_VARARGS, doc_bspldismat},
1535 {NULL, NULL, 0, NULL}
1536 };
1537 
1538 static struct PyModuleDef moduledef = {
1539     PyModuleDef_HEAD_INIT,
1540     "_fitpack",
1541     NULL,
1542     -1,
1543     fitpack_module_methods,
1544     NULL,
1545     NULL,
1546     NULL,
1547     NULL
1548 };
1549 
PyInit__fitpack(void)1550 PyObject *PyInit__fitpack(void)
1551 {
1552     PyObject *m, *d, *s;
1553 
1554     m = PyModule_Create(&moduledef);
1555     import_array();
1556 
1557     d = PyModule_GetDict(m);
1558 
1559     s = PyUnicode_FromString(" 1.7 ");
1560     PyDict_SetItemString(d, "__version__", s);
1561     fitpack_error = PyErr_NewException ("fitpack.error", NULL, NULL);
1562     Py_DECREF(s);
1563     if (PyErr_Occurred()) {
1564         Py_FatalError("can't initialize module fitpack");
1565     }
1566 
1567     return m;
1568 }
1569