Context Navigation

source: git/kernel/fglmcomb.cc @ 16f511

spielwiese

Last change on this file since 16f511 was 16f511, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fixed the usage of "config.h" (if defined HAVE_CONFIG_H)
Property mode set to `100644`
File size: 15.6 KB

Line
1	// emacs edit mode for this file is -- C++ --
2
3	/****************************************
4	* Computer Algebra System SINGULAR *
5	****************************************/
6	/*
7	* ABSTRACT -
8	*/
9
10	#ifdef HAVE_CONFIG_H
11	#include "config.h"
12	#endif /* HAVE_CONFIG_H */
13	#include <kernel/mod2.h>
14
15	#ifdef HAVE_FACTORY
16	#define SI_DONT_HAVE_GLOBAL_VARS
17	#include <factory/factory.h>
18	#include <misc/options.h>
19	#include <kernel/polys.h>
20	#include <kernel/ideals.h>
21	#include <polys/monomials/ring.h>
22	#include <kernel/febase.h>
23	#include <polys/monomials/maps.h>
24	#include <omalloc/omalloc.h>
25	#include <kernel/fglmvec.h>
26	#include <kernel/fglmgauss.h>
27	#include <kernel/kstd1.h>
28
29	#include <kernel/fglm.h>
30
31	#ifndef NOSTREAMIO
32	# ifdef HAVE_IOSTREAM
33	# include <iostream>
34	# else
35	# include <iostream.h>
36	# endif
37	#endif
38
39	static void
40	fglmEliminateMonomials( poly * pptr, fglmVector & v, polyset monomials, int numMonoms )
41	{
42	poly temp = *pptr;
43	poly pretemp = NULL;
44	int point = 0;
45	int state;
46
47	while ( (temp != NULL) && (point < numMonoms) ) {
48	state= pCmp( temp, monomials[point] );
49	if ( state == 0 ) {
50	// Eliminate this monomial
51	poly todelete;
52	if ( pretemp == NULL ) {
53	todelete = temp;
54	pIter( *pptr );
55	temp= *pptr;
56	}
57	else {
58	todelete= temp;
59	pIter( temp );
60	pretemp->next= temp;
61	}
62	pGetCoeff( todelete )= nNeg( pGetCoeff( todelete ) );
63	number newelem = nAdd( pGetCoeff( todelete ), v.getconstelem( point+1 ) );
64	v.setelem( point+1, newelem );
65	nDelete( & pGetCoeff( todelete ) );
66	pLmFree( todelete );
67	point++;
68	}
69	else if ( state < 0 )
70	point++;
71	else {
72	pretemp= temp;
73	pIter( temp );
74	}
75	}
76	}
77
78	static BOOLEAN
79	fglmReductionStep( poly * pptr, ideal source, int * w )
80	{
81	// returns TRUE if the leading monomial was reduced
82	if ( *pptr == NULL ) return FALSE;
83	int k;
84	int best = 0;
85	for ( k= IDELEMS( source ) - 1; k >= 0; k-- ) {
86	if ( pDivisibleBy( (source->m)[k], *pptr ) ) {
87	if ( best == 0 ) {
88	best= k + 1;
89	}
90	else {
91	if ( w[k] < w[best-1] ) {
92	best= k + 1;
93	}
94	}
95	}
96	}
97	if ( best > 0 ) {
98	// OwnSPoly
99	poly p2 = (source->m)[best-1];
100	int i, diff;
101
102	poly m = pOne();
103	for ( i= (currRing->N); i > 0; i-- )
104	{
105	diff= pGetExp( *pptr, i ) - pGetExp( p2, i );
106	pSetExp( m, i, diff );
107	}
108	pSetm( m );
109	number n1 = nCopy( pGetCoeff( *pptr ) );
110	number n2 = pGetCoeff( p2 );
111
112	p2= pCopy( p2 );
113	pLmDelete(pptr);
114	pLmDelete( & p2 );
115	p2= pMult( m, p2 );
116
117	number temp = nDiv( n1, n2 );
118	n_Normalize( temp, currRing );
119	nDelete( & n1 );
120	n1= temp;
121	n1= nNeg( n1 );
122	pMult_nn( p2, n1 );
123	// pNormalize( p2 );
124	nDelete( & n1 );
125	pptr= pAdd( pptr, p2 );
126	}
127	return ( (best > 0) );
128	}
129
130	static void
131	fglmReduce( poly * pptr, fglmVector & v, polyset m, int numMonoms, ideal source, int * w )
132	{
133	BOOLEAN reduced = FALSE;
134	reduced= fglmReductionStep( pptr, source, w );
135	while ( reduced == TRUE ) {
136	fglmEliminateMonomials( pptr, v, m, numMonoms );
137	reduced= fglmReductionStep( pptr, source, w );
138	}
139	STICKYPROT( "<" );
140	poly temp = * pptr;
141	if ( temp != NULL ) {
142	while ( pNext( temp ) != NULL ) {
143	STICKYPROT( ">" );
144	reduced= fglmReductionStep( & pNext( temp ), source, w );
145	while ( reduced == TRUE ) {
146	fglmEliminateMonomials( & pNext( temp ), v, m, numMonoms );
147	reduced= fglmReductionStep( & pNext( temp ), source, w );
148	}
149	if ( pNext( temp ) != NULL ) {
150	pIter( temp );
151	}
152	}
153	}
154	}
155
156	static poly
157	fglmNewLinearCombination( ideal source, poly monset )
158	{
159	polyset m = NULL;
160	polyset nf = NULL;
161	fglmVector * mv = NULL;
162
163	fglmVector * v = NULL;
164	polyset basis = NULL;
165	int basisSize = 0;
166	int basisBS = 16;
167	int basisMax = basisBS;
168
169	int * weights = NULL;
170	int * lengthes = NULL;
171	int * order = NULL;
172
173	int numMonoms = 0;
174	int k;
175
176	// get number of monoms
177	numMonoms= pLength( monset );
178	STICKYPROT2( "%i monoms\n", numMonoms );
179
180	// Allcoate Memory and initialize sets
181	m= (polyset)omAlloc( numMonoms * sizeof( poly ) );
182	poly temp= monset;
183	for ( k= 0; k < numMonoms; k++ ) {
184	// m[k]= pOne();
185	// pSetExpV( m[k], temp->exp );
186	// pSetm( m[k] );
187	m[k]=pLmInit(temp);
188	pSetCoeff( m[k], nInit(1) );
189	pIter( temp );
190	}
191
192	nf= (polyset)omAlloc( numMonoms * sizeof( poly ) );
193
194	#ifndef HAVE_EXPLICIT_CONSTR
195	mv= new fglmVector[ numMonoms ];
196	#else
197	mv= (fglmVector )omAlloc( numMonoms sizeof( fglmVector ) );
198	#endif
199
200	#ifndef HAVE_EXPLICIT_CONSTR
201	v= new fglmVector[ numMonoms ];
202	#else
203	v= (fglmVector )omAlloc( numMonoms sizeof( fglmVector ) );
204	#endif
205
206	basisMax= basisBS;
207	basis= (polyset)omAlloc( basisMax * sizeof( poly ) );
208
209	weights= (int )omAlloc( IDELEMS( source ) sizeof( int ) );
210	STICKYPROT( "weights: " );
211	for ( k= 0; k < IDELEMS( source ); k++ ) {
212	poly temp= (source->m)[k];
213	int w = 0;
214	while ( temp != NULL ) {
215	w+= nSize( pGetCoeff( temp ) );
216	pIter( temp );
217	}
218	weights[k]= w;
219	STICKYPROT2( "%i ", w );
220	}
221	STICKYPROT( "\n" );
222	lengthes= (int )omAlloc( numMonoms sizeof( int ) );
223	order= (int )omAlloc( numMonoms sizeof( int ) );
224
225	// calculate the NormalForm in a special way
226	for ( k= 0; k < numMonoms; k++ )
227	{
228	STICKYPROT( "#" );
229	poly current = pCopy( m[k] );
230	fglmVector currV( numMonoms, k+1 );
231
232	fglmReduce( & current, currV, m, numMonoms, source, weights );
233	poly temp = current;
234	int b;
235	while ( temp != NULL )
236	{
237	BOOLEAN found = FALSE;
238	for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
239	{
240	if ( pLmEqual( temp, basis[b] ) )
241	{
242	found= TRUE;
243	}
244	}
245	if ( found == FALSE )
246	{
247	if ( basisSize == basisMax )
248	{
249	// Expand the basis
250	basis= (polyset)omReallocSize( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
251	basisMax+= basisBS;
252	}
253	// basis[basisSize]= pOne();
254	// pSetExpV( basis[basisSize], temp->exp );
255	// pSetm( basis[basisSize] );
256	basis[basisSize]=pLmInit(temp);
257	pSetCoeff( basis[basisSize], nInit(1) );
258	basisSize++;
259	}
260	pIter( temp );
261	}
262	nf[k]= current;
263	#ifndef HAVE_EXPLICIT_CONSTR
264	mv[k].mac_constr( currV );
265	#else
266	mv[k].fglmVector( currV );
267	#endif
268	STICKYPROT( "\n" );
269	}
270	// get the vector representation
271	for ( k= 0; k < numMonoms; k++ ) {
272	STICKYPROT( "." );
273
274	#ifndef HAVE_EXPLICIT_CONSTR
275	v[k].mac_constr_i( basisSize );
276	#else
277	v[k].fglmVector( basisSize );
278	#endif
279	poly mon= nf[k];
280	while ( mon != NULL ) {
281	BOOLEAN found = FALSE;
282	int b= 0;
283	while ( found == FALSE ) {
284	if ( pLmEqual( mon, basis[b] ) ) {
285	found= TRUE;
286	}
287	else {
288	b++;
289	}
290	}
291	number coeff = nCopy( pGetCoeff( mon ) );
292	v[k].setelem( b+1, coeff );
293	pIter( mon );
294	}
295	pDelete( nf + k );
296	}
297	// Free Memory not needed anymore
298	omFreeSize( (ADDRESS)nf, numMonoms * sizeof( poly ) );
299	omFreeSize( (ADDRESS)weights, IDELEMS( source ) * sizeof( int ) );
300
301	STICKYPROT2( "\nbasis size: %i\n", basisSize );
302	STICKYPROT( "(clear basis" );
303	for ( k= 0; k < basisSize; k++ )
304	pDelete( basis + k );
305	STICKYPROT( ")\n" );
306	// gauss reduce
307	gaussReducer gauss( basisSize );
308	BOOLEAN isZero = FALSE;
309	fglmVector p;
310
311	STICKYPROT( "sizes: " );
312	for ( k= 0; k < numMonoms; k++ ) {
313	lengthes[k]= v[k].numNonZeroElems();
314	STICKYPROT2( "%i ", lengthes[k] );
315	}
316	STICKYPROT( "\n" );
317
318	int act = 0;
319	while ( (isZero == FALSE) && (act < numMonoms) ) {
320	int best = 0;
321	for ( k= numMonoms - 1; k >= 0; k-- ) {
322	if ( lengthes[k] > 0 ) {
323	if ( best == 0 ) {
324	best= k+1;
325	}
326	else {
327	if ( lengthes[k] < lengthes[best-1] ) {
328	best= k+1;
329	}
330	}
331	}
332	}
333	lengthes[best-1]= 0;
334	order[act]= best-1;
335	STICKYPROT2( " (%i) ", best );
336	if ( ( isZero= gauss.reduce( v[best-1] )) == TRUE ) {
337	p= gauss.getDependence();
338	}
339	else {
340	STICKYPROT( "+" );
341	gauss.store();
342	act++;
343	}
344	#ifndef HAVE_EXPLICIT_CONSTR
345	v[best-1].clearelems();
346	#else
347	v[best-1].~fglmVector();
348	#endif
349	}
350	poly result = NULL;
351	if ( isZero == TRUE ) {
352	number gcd = p.gcd();
353	if ( ! nIsZero( gcd ) && ! ( nIsOne( gcd ) ) ) {
354	p/= gcd;
355	}
356	nDelete( & gcd );
357	fglmVector temp( numMonoms );
358	for ( k= 0; k < p.size(); k++ ) {
359	if ( ! p.elemIsZero( k+1 ) ) {
360	temp+= p.getconstelem( k+1 ) * mv[order[k]];
361	}
362	}
363	number denom = temp.clearDenom();
364	nDelete( & denom );
365	gcd = temp.gcd();
366	if ( ! nIsZero( gcd ) && ! nIsOne( gcd ) ) {
367	temp/= gcd;
368	}
369	nDelete( & gcd );
370
371	poly sum = NULL;
372	for ( k= 1; k <= numMonoms; k++ ) {
373	if ( ! temp.elemIsZero( k ) ) {
374	if ( result == NULL ) {
375	result= pCopy( m[k-1] );
376	sum= result;
377	}
378	else {
379	sum->next= pCopy( m[k-1] );
380	pIter( sum );
381	}
382	pSetCoeff( sum, nCopy( temp.getconstelem( k ) ) );
383	}
384	}
385	p_Content( result, currRing );
386	if ( ! nGreaterZero( pGetCoeff( result ) ) ) result= pNeg( result );
387	}
388	// Free Memory
389	omFreeSize( (ADDRESS)lengthes, numMonoms * sizeof( int ) );
390	omFreeSize( (ADDRESS)order, numMonoms * sizeof( int ) );
391	// for ( k= 0; k < numMonoms; k++ )
392	// v[k].~fglmVector();
393	#ifndef HAVE_EXPLICIT_CONSTR
394	delete [] v;
395	#else
396	omFreeSize( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
397	#endif
398
399	for ( k= 0; k < basisSize; k++ )
400	pDelete( basis + k );
401	omFreeSize( (ADDRESS)basis, basisMax * sizeof( poly ) );
402
403	#ifndef HAVE_EXPLICIT_CONSTR
404	delete [] mv;
405	#else
406	for ( k= 0; k < numMonoms; k++ )
407	mv[k].~fglmVector();
408	omFreeSize( (ADDRESS)mv, numMonoms * sizeof( fglmVector ) );
409	#endif
410
411	for ( k= 0; k < numMonoms; k++ )
412	pDelete( m + k );
413	omFreeSize( (ADDRESS)m, numMonoms * sizeof( poly ) );
414
415	STICKYPROT( "\n" );
416	return result;
417	}
418
419
420	static poly
421	fglmLinearCombination( ideal source, poly monset )
422	{
423	int k;
424	poly temp;
425
426	polyset m;
427	polyset nf;
428	fglmVector * v;
429
430	polyset basis;
431	int basisSize = 0;
432	int basisBS = 16;
433	int basisMax;
434	// get number of monomials in monset
435	int numMonoms = 0;
436	temp = monset;
437	while ( temp != NULL ) {
438	numMonoms++;
439	pIter( temp );
440	}
441	// Allocate Memory
442	m= (polyset)omAlloc( numMonoms * sizeof( poly ) );
443	nf= (polyset)omAlloc( numMonoms * sizeof( poly ) );
444	// Warning: The fglmVectors in v are yet not initialized
445	v= (fglmVector )omAlloc( numMonoms sizeof( fglmVector ) );
446	basisMax= basisBS;
447	basis= (polyset)omAlloc( basisMax * sizeof( poly ) );
448
449	// get the NormalForm and the basis monomials
450	temp= monset;
451	for ( k= 0; k < numMonoms; k++ ) {
452	poly mon= pHead( temp );
453	if ( ! nIsOne( pGetCoeff( mon ) ) ) {
454	nDelete( & pGetCoeff( mon ) );
455	pSetCoeff( mon, nInit( 1 ) );
456	}
457	STICKYPROT( "(" );
458	nf[k]= kNF( source, currRing->qideal, mon );
459	STICKYPROT( ")" );
460
461	// search through basis
462	STICKYPROT( "[" );
463	poly sm = nf[k];
464	while ( sm != NULL ) {
465	BOOLEAN found = FALSE;
466	int b;
467	for ( b= 0; (b < basisSize) && (found == FALSE); b++ )
468	if ( pLmEqual( sm, basis[b] ) ) found= TRUE;
469	if ( found == FALSE ) {
470	// Expand the basis
471	if ( basisSize == basisMax ) {
472	basis= (polyset)omReallocSize( basis, basisMax * sizeof( poly ), (basisMax + basisBS ) * sizeof( poly ) );
473	basisMax+= basisBS;
474	}
475	basis[basisSize]= pHead( sm );
476	nDelete( & pGetCoeff( basis[basisSize] ) );
477	pSetCoeff( basis[basisSize], nInit( 1 ) );
478	basisSize++;
479	}
480	pIter( sm );
481	}
482	STICKYPROT( "]" );
483	m[k]= mon;
484	pIter( temp );
485	}
486	// get the vector representation
487	STICKYPROT2( "(%i)", basisSize );
488	for ( k= 0; k < numMonoms; k++ ) {
489	#ifndef HAVE_EXPLICIT_CONSTR
490	v[k].mac_constr_i( basisSize );
491	#else
492	v[k].fglmVector( basisSize );
493	#endif
494	STICKYPROT( "(+" );
495	poly mon= nf[k];
496	while ( mon != NULL ) {
497	BOOLEAN found = FALSE;
498	int b= 0;
499	while ( found == FALSE ) {
500	if ( pLmEqual( mon, basis[b] ) )
501	found= TRUE;
502	else
503	b++;
504	}
505	number coeff = nCopy( pGetCoeff( mon ) );
506	v[k].setelem( b+1, coeff );
507	pIter( mon );
508	}
509	STICKYPROT( ")" );
510	}
511	// gauss reduce
512	gaussReducer gauss( basisSize );
513	BOOLEAN isZero = FALSE;
514	fglmVector p;
515	for ( k= 0; (k < numMonoms) && (isZero == FALSE); k++ ) {
516	STICKYPROT( "(-" );
517	if ( ( isZero= gauss.reduce( v[k] )) == TRUE )
518	p= gauss.getDependence();
519	else
520	gauss.store();
521	STICKYPROT( ")" );
522	}
523	poly comb = NULL;
524	if ( isZero == TRUE ) {
525	number gcd = p.gcd();
526	if ( ! nIsZero( gcd ) && ! ( nIsOne( gcd ) ) )
527	p/= gcd;
528	nDelete( & gcd );
529	for ( k= 1; k <= p.size(); k++ ) {
530	if ( ! p.elemIsZero( k ) ) {
531	poly temp = pCopy( m[k-1] );
532	pSetCoeff( temp, nCopy( p.getconstelem( k ) ) );
533	comb= pAdd( comb, temp );
534	}
535	}
536	if ( ! nGreaterZero( pGetCoeff( comb ) ) ) comb= pNeg( comb );
537	}
538
539	// Free Memory
540	for ( k= 0; k < numMonoms; k++ ) {
541	pDelete( m + k );
542	pDelete( nf + k );
543	}
544	omFreeSize( (ADDRESS)m, numMonoms * sizeof( poly ) );
545	omFreeSize( (ADDRESS)nf, numMonoms * sizeof( poly ) );
546	// Warning: At this point all Vectors in v have to be initialized
547	for ( k= numMonoms - 1; k >= 0; k-- ) v[k].~fglmVector();
548	omFreeSize( (ADDRESS)v, numMonoms * sizeof( fglmVector ) );
549	for ( k= 0; k < basisSize; k++ )
550	pDelete( basis + k );
551	omFreeSize( (ADDRESS)basis, basisMax * sizeof( poly ) );
552	STICKYPROT( "\n" );
553	return comb;
554	}
555
556	#endif
557	// Local Variables: ***
558	// compile-command: "make Singular" ***
559	// page-delimiter: "^\\(\\\|//!\\)" ***
560	// fold-internal-margins: nil ***
561	// End: ***

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

Download in other formats: