Context Navigation

source: git/kernel/fglm/fglmcomb.cc @ 066288

spielwiese

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

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

Download in other formats: