2Y_cpoly.c @ 79492222

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 79492222 was 79492222, checked in by krzywon, 9 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 @ 79492222

Download in other formats: