2Y_cpoly.c @ 6c2e3aa

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 6c2e3aa was 79492222, checked in by krzywon, 10 years ago
Changed the file and folder names to remove all SANS references.
Property mode set to `100644`
File size: 17.1 KB

Line
1	/*
2	*******************************************************************************
3	*
4	*
5	* Copyright (c) 2002
6	* Henrik Vestermark
7	* Denmark
8	*
9	* All Rights Reserved
10	*
11	* This source file is subject to the terms and conditions of the
12	* Henrik Vestermark Software License Agreement which restricts the manner
13	* in which it may be used.
14	*
15	*******************************************************************************
16	*
17	*
18	* Module name : cpoly.cpp
19	* Module ID Nbr :
20	* Description : cpoly.cpp -- Jenkins-Traub real polynomial root finder.
21	* Translation of TOMS493 from FORTRAN to C. This
22	* implementation of Jenkins-Traub partially adapts
23	* the original code to a C environment by restruction
24	* many of the 'goto' controls to better fit a block
25	* structured form. It also eliminates the global memory
26	* allocation in favor of local, dynamic memory management.
27	*
28	* The calling conventions are slightly modified to return
29	* the number of roots found as the function value.
30	*
31	* INPUT:
32	* opr - double precision vector of real coefficients in order of
33	* decreasing powers.
34	* opi - double precision vector of imaginary coefficients in order of
35	* decreasing powers.
36	* degree - integer degree of polynomial
37	*
38	* OUTPUT:
39	* zeror,zeroi - output double precision vectors of the
40	* real and imaginary parts of the zeros.
41	* to be consistent with rpoly.cpp the zeros is inthe index
42	* [0..max_degree-1]
43	*
44	* RETURN:
45	* returnval: -1 if leading coefficient is zero, otherwise
46	* number of roots found.
47	* --------------------------------------------------------------------------
48	* Change Record :
49	*
50	* Version Author/Date Description of changes
51	* ------- ----------- ----------------------
52	* 01.01 HVE/021101 Initial release
53	*
54	* -- SRK changed name of pi to piw (and all occurances in this file) to avoid the conflict with pi in Apple's fp.h
55	*
56	* End of Change Record
57	* --------------------------------------------------------------------------
58	*/
59
60
61	/* define version string */
62
63	#include <stdlib.h>
64	#include <math.h>
65	#include <float.h>
66
67	static double sr, si, tr, ti, pvr, pvi, are, mre, eta, infin;
68	static int nn;
69	static double pr, piw, hr, hi, qpr, qpi, qhr, qhi, shr, shi;
70
71	static void noshft( const int l1 );
72	static void fxshft( const int l2, double zr, double zi, int *conv );
73	static void vrshft( const int l3, double zr, double zi, int *conv );
74	static void calct( int *bol );
75	static void nexth( const int bol );
76	static void polyev( const int nn, const double sr, const double si, const double pr[], const double piw[], double qr[], double qi[], double pvr, double pvi );
77	static double errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, const double are, const double mre );
78	static void cauchy( const int nn, double pt[], double q[], double *fn_val );
79	static double scale( const int nn, const double pt[], const double eta, const double infin, const double smalno, const double base );
80	static void cdivid( const double ar, const double ai, const double br, const double bi, double cr, double ci );
81	static double cmod( const double r, const double i );
82	static void mcon( double eta, double infiny, double smalno, double base );
83
84	int cpoly( const double opr, const double opi, int degree, double zeror, double zeroi )
85	{
86	int cnt1, cnt2, idnn2, i, conv;
87	double xx, yy, cosr, sinr, smalno, base, xxx, zr, zi, bnd;
88
89	mcon( &eta, &infin, &smalno, &base );
90	are = eta;
91	mre = 2.0 * sqrt( 2.0 ) * eta;
92	xx = 0.70710678;
93	yy = -xx;
94	cosr = -0.060756474;
95	sinr = -0.99756405;
96	nn = degree;
97
98	// Algorithm fails if the leading coefficient is zero
99	if( opr[ 0 ] == 0 && opi[ 0 ] == 0 )
100	return -1;
101
102	// Allocate arrays
103	pr = malloc( sizeof( double ) * ( degree + 1 ) );
104	piw = malloc( sizeof( double ) * ( degree + 1 ) );
105	hr =malloc( sizeof( double ) * ( degree + 1 ) );
106	hi = malloc( sizeof( double ) * ( degree + 1 ) );
107	qpr= malloc( sizeof( double ) * ( degree + 1 ) );
108	qpi= malloc( sizeof( double ) * ( degree + 1 ) );
109	qhr= malloc( sizeof( double ) * ( degree + 1 ) );
110	qhi= malloc( sizeof( double ) * ( degree + 1 ) );
111	shr= malloc( sizeof( double ) * ( degree + 1 ) );
112	shi= malloc( sizeof( double ) * ( degree + 1 ) );
113
114	// Remove the zeros at the origin if any
115	while( opr[ nn ] == 0 && opi[ nn ] == 0 )
116	{
117	idnn2 = degree - nn;
118	zeror[ idnn2 ] = 0;
119	zeroi[ idnn2 ] = 0;
120	nn--;
121	}
122
123	// Make a copy of the coefficients
124	for( i = 0; i <= nn; i++ )
125	{
126	pr[ i ] = opr[ i ];
127	piw[ i ] = opi[ i ];
128	shr[ i ] = cmod( pr[ i ], piw[ i ] );
129	}
130
131	// Scale the polynomial
132	bnd = scale( nn, shr, eta, infin, smalno, base );
133	if( bnd != 1 )
134	for( i = 0; i <= nn; i++ )
135	{
136	pr[ i ] *= bnd;
137	piw[ i ] *= bnd;
138	}
139
140	search:
141	if( nn <= 1 )
142	{
143	cdivid( -pr[ 1 ], -piw[ 1 ], pr[ 0 ], piw[ 0 ], &zeror[ degree-1 ], &zeroi[ degree-1 ] );
144	goto finish;
145	}
146
147	// Calculate bnd, alower bound on the modulus of the zeros
148	for( i = 0; i<= nn; i++ )
149	shr[ i ] = cmod( pr[ i ], piw[ i ] );
150
151	cauchy( nn, shr, shi, &bnd );
152
153	// Outer loop to control 2 Major passes with different sequences of shifts
154	for( cnt1 = 1; cnt1 <= 2; cnt1++ )
155	{
156	// First stage calculation , no shift
157	noshft( 5 );
158
159	// Inner loop to select a shift
160	for( cnt2 = 1; cnt2 <= 9; cnt2++ )
161	{
162	// Shift is chosen with modulus bnd and amplitude rotated by 94 degree from the previous shif
163	xxx = cosr * xx - sinr * yy;
164	yy = sinr * xx + cosr * yy;
165	xx = xxx;
166	sr = bnd * xx;
167	si = bnd * yy;
168
169	// Second stage calculation, fixed shift
170	fxshft( 10 * cnt2, &zr, &zi, &conv );
171	if( conv )
172	{
173	// The second stage jumps directly to the third stage ieration
174	// If successful the zero is stored and the polynomial deflated
175	idnn2 = degree - nn;
176	zeror[ idnn2 ] = zr;
177	zeroi[ idnn2 ] = zi;
178	nn--;
179	for( i = 0; i <= nn; i++ )
180	{
181	pr[ i ] = qpr[ i ];
182	piw[ i ] = qpi[ i ];
183	}
184	goto search;
185	}
186	// If the iteration is unsuccessful another shift is chosen
187	}
188	// if 9 shifts fail, the outer loop is repeated with another sequence of shifts
189	}
190
191	// The zerofinder has failed on two major passes
192	// return empty handed with the number of roots found (less than the original degree)
193	degree -= nn;
194
195	finish:
196	// Deallocate arrays
197	free( pr );
198	free( piw );
199	free( hr );
200	free( hi );
201	free( qpr );
202	free( qpi );
203	free( qhr );
204	free( qhi );
205	free( shr );
206	free( shi );
207
208	return degree;
209	}
210
211
212	// COMPUTES THE DERIVATIVE POLYNOMIAL AS THE INITIAL H
213	// POLYNOMIAL AND COMPUTES L1 NO-SHIFT H POLYNOMIALS.
214	//
215	static void noshft( const int l1 )
216	{
217	int i, j, jj, n, nm1;
218	double xni, t1, t2;
219
220	n = nn;
221	nm1 = n - 1;
222	for( i = 0; i < n; i++ )
223	{
224	xni = nn - i;
225	hr[ i ] = xni * pr[ i ] / n;
226	hi[ i ] = xni * piw[ i ] / n;
227	}
228	for( jj = 1; jj <= l1; jj++ )
229	{
230	if( cmod( hr[ n - 1 ], hi[ n - 1 ] ) > eta * 10 * cmod( pr[ n - 1 ], piw[ n - 1 ] ) )
231	{
232	cdivid( -pr[ nn ], -piw[ nn ], hr[ n - 1 ], hi[ n - 1 ], &tr, &ti );
233	for( i = 0; i < nm1; i++ )
234	{
235	j = nn - i - 1;
236	t1 = hr[ j - 1 ];
237	t2 = hi[ j - 1 ];
238	hr[ j ] = tr * t1 - ti * t2 + pr[ j ];
239	hi[ j ] = tr * t2 + ti * t1 + piw[ j ];
240	}
241	hr[ 0 ] = pr[ 0 ];
242	hi[ 0 ] = piw[ 0 ];
243	}
244	else
245	{
246	// If the constant term is essentially zero, shift H coefficients
247	for( i = 0; i < nm1; i++ )
248	{
249	j = nn - i - 1;
250	hr[ j ] = hr[ j - 1 ];
251	hi[ j ] = hi[ j - 1 ];
252	}
253	hr[ 0 ] = 0;
254	hi[ 0 ] = 0;
255	}
256	}
257	}
258
259	// COMPUTES L2 FIXED-SHIFT H POLYNOMIALS AND TESTS FOR CONVERGENCE.
260	// INITIATES A VARIABLE-SHIFT ITERATION AND RETURNS WITH THE
261	// APPROXIMATE ZERO IF SUCCESSFUL.
262	// L2 - LIMIT OF FIXED SHIFT STEPS
263	// ZR,ZI - APPROXIMATE ZERO IF CONV IS .TRUE.
264	// CONV - LOGICAL INDICATING CONVERGENCE OF STAGE 3 ITERATION
265	//
266	static void fxshft( const int l2, double zr, double zi, int *conv )
267	{
268	int i, j, n;
269	int test, pasd, bol;
270	double otr, oti, svsr, svsi;
271
272	n = nn;
273	polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
274	test = 1;
275	pasd = 0;
276
277	// Calculate first T = -P(S)/H(S)
278	calct( &bol );
279
280	// Main loop for second stage
281	for( j = 1; j <= l2; j++ )
282	{
283	otr = tr;
284	oti = ti;
285
286	// Compute the next H Polynomial and new t
287	nexth( bol );
288	calct( &bol );
289	*zr = sr + tr;
290	*zi = si + ti;
291
292	// Test for convergence unless stage 3 has failed once or this
293	// is the last H Polynomial
294	if( !( bol \|\| !test \|\| j == 12 ) )
295	if( cmod( tr - otr, ti - oti ) < 0.5 * cmod( zr, zi ) )
296	{
297	if( pasd )
298	{
299	// The weak convergence test has been passwed twice, start the third stage
300	// Iteration, after saving the current H polynomial and shift
301	for( i = 0; i < n; i++ )
302	{
303	shr[ i ] = hr[ i ];
304	shi[ i ] = hi[ i ];
305	}
306	svsr = sr;
307	svsi = si;
308	vrshft( 10, zr, zi, conv );
309	if( *conv ) return;
310
311	//The iteration failed to converge. Turn off testing and restore h,s,pv and T
312	test = 0;
313	for( i = 0; i < n; i++ )
314	{
315	hr[ i ] = shr[ i ];
316	hi[ i ] = shi[ i ];
317	}
318	sr = svsr;
319	si = svsi;
320	polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
321	calct( &bol );
322	continue;
323	}
324	pasd = 1;
325	}
326	else
327	pasd = 0;
328	}
329
330	// Attempt an iteration with final H polynomial from second stage
331	vrshft( 10, zr, zi, conv );
332	}
333
334	// CARRIES OUT THE THIRD STAGE ITERATION.
335	// L3 - LIMIT OF STEPS IN STAGE 3.
336	// ZR,ZI - ON ENTRY CONTAINS THE INITIAL ITERATE, IF THE
337	// ITERATION CONVERGES IT CONTAINS THE FINAL ITERATE ON EXIT.
338	// CONV - .TRUE. IF ITERATION CONVERGES
339	//
340	static void vrshft( const int l3, double zr, double zi, int *conv )
341	{
342	int b, bol;
343	int i, j;
344	double mp, ms, omp, relstp, r1, r2, tp;
345
346	*conv = 0;
347	b = 0;
348	sr = *zr;
349	si = *zi;
350
351	// Main loop for stage three
352	for( i = 1; i <= l3; i++ )
353	{
354	// Evaluate P at S and test for convergence
355	polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
356	mp = cmod( pvr, pvi );
357	ms = cmod( sr, si );
358	if( mp <= 20 * errev( nn, qpr, qpi, ms, mp, are, mre ) )
359	{
360	// Polynomial value is smaller in value than a bound onthe error
361	// in evaluationg P, terminate the ietartion
362	*conv = 1;
363	*zr = sr;
364	*zi = si;
365	return;
366	}
367	if( i != 1 )
368	{
369	if( !( b \|\| mp < omp \|\| relstp >= 0.05 ) )
370	{
371	// Iteration has stalled. Probably a cluster of zeros. Do 5 fixed
372	// shift steps into the cluster to force one zero to dominate
373	tp = relstp;
374	b = 1;
375	if( relstp < eta ) tp = eta;
376	r1 = sqrt( tp );
377	r2 = sr * ( 1 + r1 ) - si * r1;
378	si = sr * r1 + si * ( 1 + r1 );
379	sr = r2;
380	polyev( nn, sr, si, pr, piw, qpr, qpi, &pvr, &pvi );
381	for( j = 1; j <= 5; j++ )
382	{
383	calct( &bol );
384	nexth( bol );
385	}
386	omp = infin;
387	goto _20;
388	}
389
390	// Exit if polynomial value increase significantly
391	if( mp *0.1 > omp ) return;
392	}
393
394	omp = mp;
395
396	// Calculate next iterate
397	_20: calct( &bol );
398	nexth( bol );
399	calct( &bol );
400	if( !bol )
401	{
402	relstp = cmod( tr, ti ) / cmod( sr, si );
403	sr += tr;
404	si += ti;
405	}
406	}
407	}
408
409	// COMPUTES T = -P(S)/H(S).
410	// BOOL - LOGICAL, SET TRUE IF H(S) IS ESSENTIALLY ZERO.
411	static void calct( int *bol )
412	{
413	int n;
414	double hvr, hvi;
415
416	n = nn;
417
418	// evaluate h(s)
419	polyev( n - 1, sr, si, hr, hi, qhr, qhi, &hvr, &hvi );
420	bol = cmod( hvr, hvi ) <= are 10 * cmod( hr[ n - 1 ], hi[ n - 1 ] ) ? 1 : 0;
421	if( !*bol )
422	{
423	cdivid( -pvr, -pvi, hvr, hvi, &tr, &ti );
424	return;
425	}
426
427	tr = 0;
428	ti = 0;
429	}
430
431	// CALCULATES THE NEXT SHIFTED H POLYNOMIAL.
432	// BOOL - LOGICAL, IF .TRUE. H(S) IS ESSENTIALLY ZERO
433	//
434	static void nexth( const int bol )
435	{
436	int j, n;
437	double t1, t2;
438
439	n = nn;
440	if( !bol )
441	{
442	for( j = 1; j < n; j++ )
443	{
444	t1 = qhr[ j - 1 ];
445	t2 = qhi[ j - 1 ];
446	hr[ j ] = tr * t1 - ti * t2 + qpr[ j ];
447	hi[ j ] = tr * t2 + ti * t1 + qpi[ j ];
448	}
449	hr[ 0 ] = qpr[ 0 ];
450	hi[ 0 ] = qpi[ 0 ];
451	return;
452	}
453
454	// If h[s] is zero replace H with qh
455	for( j = 1; j < n; j++ )
456	{
457	hr[ j ] = qhr[ j - 1 ];
458	hi[ j ] = qhi[ j - 1 ];
459	}
460	hr[ 0 ] = 0;
461	hi[ 0 ] = 0;
462	}
463
464	// EVALUATES A POLYNOMIAL P AT S BY THE HORNER RECURRENCE
465	// PLACING THE PARTIAL SUMS IN Q AND THE COMPUTED VALUE IN PV.
466	//
467	static void polyev( const int nn, const double sr, const double si, const double pr[], const double piw[], double qr[], double qi[], double pvr, double pvi )
468	{
469	int i;
470	double t;
471
472	qr[ 0 ] = pr[ 0 ];
473	qi[ 0 ] = piw[ 0 ];
474	*pvr = qr[ 0 ];
475	*pvi = qi[ 0 ];
476
477	for( i = 1; i <= nn; i++ )
478	{
479	t = ( pvr ) sr - ( pvi ) si + pr[ i ];
480	pvi = ( pvr ) * si + ( pvi ) sr + piw[ i ];
481	*pvr = t;
482	qr[ i ] = *pvr;
483	qi[ i ] = *pvi;
484	}
485	}
486
487	// BOUNDS THE ERROR IN EVALUATING THE POLYNOMIAL BY THE HORNER RECURRENCE.
488	// QR,QI - THE PARTIAL SUMS
489	// MS -MODULUS OF THE POINT
490	// MP -MODULUS OF POLYNOMIAL VALUE
491	// ARE, MRE -ERROR BOUNDS ON COMPLEX ADDITION AND MULTIPLICATION
492	//
493	static double errev( const int nn, const double qr[], const double qi[], const double ms, const double mp, const double are, const double mre )
494	{
495	int i;
496	double e;
497
498	e = cmod( qr[ 0 ], qi[ 0 ] ) * mre / ( are + mre );
499	for( i = 0; i <= nn; i++ )
500	e = e * ms + cmod( qr[ i ], qi[ i ] );
501
502	return e * ( are + mre ) - mp * mre;
503	}
504
505	// CAUCHY COMPUTES A LOWER BOUND ON THE MODULI OF THE ZEROS OF A
506	// POLYNOMIAL - PT IS THE MODULUS OF THE COEFFICIENTS.
507	//
508	static void cauchy( const int nn, double pt[], double q[], double *fn_val )
509	{
510	int i, n;
511	double x, xm, f, dx, df;
512
513	pt[ nn ] = -pt[ nn ];
514
515	// Compute upper estimate bound
516	n = nn;
517	x = exp( log( -pt[ nn ] ) - log( pt[ 0 ] ) ) / n;
518	if( pt[ n - 1 ] != 0 )
519	{
520	// Newton step at the origin is better, use it
521	xm = -pt[ nn ] / pt[ n - 1 ];
522	if( xm < x ) x = xm;
523	}
524
525	// Chop the interval (0,x) until f < 0
526	while(1)
527	{
528	xm = x * 0.1;
529	f = pt[ 0 ];
530	for( i = 1; i <= nn; i++ )
531	f = f * xm + pt[ i ];
532	if( f <= 0 )
533	break;
534	x = xm;
535	}
536	dx = x;
537
538	// Do Newton iteration until x converges to two decimal places
539	while( fabs( dx / x ) > 0.005 )
540	{
541	q[ 0 ] = pt[ 0 ];
542	for( i = 1; i <= nn; i++ )
543	q[ i ] = q[ i - 1 ] * x + pt[ i ];
544	f = q[ nn ];
545	df = q[ 0 ];
546	for( i = 1; i < n; i++ )
547	df = df * x + q[ i ];
548	dx = f / df;
549	x -= dx;
550	}
551
552	*fn_val = x;
553	}
554
555	// RETURNS A SCALE FACTOR TO MULTIPLY THE COEFFICIENTS OF THE POLYNOMIAL.
556	// THE SCALING IS DONE TO AVOID OVERFLOW AND TO AVOID UNDETECTED UNDERFLOW
557	// INTERFERING WITH THE CONVERGENCE CRITERION. THE FACTOR IS A POWER OF THE
558	// BASE.
559	// PT - MODULUS OF COEFFICIENTS OF P
560	// ETA, INFIN, SMALNO, BASE - CONSTANTS DESCRIBING THE FLOATING POINT ARITHMETIC.
561	//
562	static double scale( const int nn, const double pt[], const double eta, const double infin, const double smalno, const double base )
563	{
564	int i, l;
565	double hi, lo, max, min, x, sc;
566	double fn_val;
567
568	// Find largest and smallest moduli of coefficients
569	hi = sqrt( infin );
570	lo = smalno / eta;
571	max = 0;
572	min = infin;
573
574	for( i = 0; i <= nn; i++ )
575	{
576	x = pt[ i ];
577	if( x > max ) max = x;
578	if( x != 0 && x < min ) min = x;
579	}
580
581	// Scale only if there are very large or very small components
582	fn_val = 1;
583	if( min >= lo && max <= hi ) return fn_val;
584	x = lo / min;
585	if( x <= 1 )
586	sc = 1 / ( sqrt( max )* sqrt( min ) );
587	else
588	{
589	sc = x;
590	if( infin / sc > max ) sc = 1;
591	}
592	l = (int)( log( sc ) / log(base ) + 0.5 );
593	fn_val = pow( base , l );
594	return fn_val;
595	}
596
597	// COMPLEX DIVISION C = A/B, AVOIDING OVERFLOW.
598	//
599	static void cdivid( const double ar, const double ai, const double br, const double bi, double cr, double ci )
600	{
601	double r, d, t, infin;
602
603	if( br == 0 && bi == 0 )
604	{
605	// Division by zero, c = infinity
606	mcon( &t, &infin, &t, &t );
607	*cr = infin;
608	*ci = infin;
609	return;
610	}
611
612	if( fabs( br ) < fabs( bi ) )
613	{
614	r = br/ bi;
615	d = bi + r * br;
616	cr = ( ar r + ai ) / d;
617	ci = ( ai r - ar ) / d;
618	return;
619	}
620
621	r = bi / br;
622	d = br + r * bi;
623	cr = ( ar + ai r ) / d;
624	ci = ( ai - ar r ) / d;
625	}
626
627	// MODULUS OF A COMPLEX NUMBER AVOIDING OVERFLOW.
628	//
629	static double cmod( const double r, const double i )
630	{
631	double ar, ai;
632
633	ar = fabs( r );
634	ai = fabs( i );
635	if( ar < ai )
636	return ai * sqrt( 1.0 + pow( ( ar / ai ), 2.0 ) );
637
638	if( ar > ai )
639	return ar * sqrt( 1.0 + pow( ( ai / ar ), 2.0 ) );
640
641	return ar * sqrt( 2.0 );
642	}
643
644	// MCON PROVIDES MACHINE CONSTANTS USED IN VARIOUS PARTS OF THE PROGRAM.
645	// THE USER MAY EITHER SET THEM DIRECTLY OR USE THE STATEMENTS BELOW TO
646	// COMPUTE THEM. THE MEANING OF THE FOUR CONSTANTS ARE -
647	// ETA THE MAXIMUM RELATIVE REPRESENTATION ERROR WHICH CAN BE DESCRIBED
648	// AS THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
649	// 1.0_dp + ETA > 1.0.
650	// INFINY THE LARGEST FLOATING-POINT NUMBER
651	// SMALNO THE SMALLEST POSITIVE FLOATING-POINT NUMBER
652	// BASE THE BASE OF THE FLOATING-POINT NUMBER SYSTEM USED
653	//
654	static void mcon( double eta, double infiny, double smalno, double base )
655	{
656	*base = 10;//DBL_RADIX;
657	*eta = DBL_EPSILON;
658	*infiny = DBL_MAX;
659	*smalno = DBL_MIN;
660	}

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/src/sas/models/c_extension/libigor/2Y_cpoly.c @ 6c2e3aa

Download in other formats: