source: git/libpolys/coeffs/flintcf_Zn.cc @ a507140

spielwiese
Last change on this file since a507140 was a507140, checked in by Hans Schoenemann <hannes@…>, 4 years ago
opt: nCoeffWrite/nCoeffString via nCoeffName
  • Property mode set to 100644
File size: 14.0 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: flint: nmod_poly_t
6*/
7#include <ctype.h> /* isdigit*/
8
9#include "misc/auxiliary.h"
10
11#ifdef HAVE_FLINT
12
13#include <flint/flint.h>
14#include <flint/nmod_poly.h>
15#include "factory/factory.h"
16
17#include "coeffs/coeffs.h"
18
19#include "coeffs/numbers.h"
20#include "coeffs/longrat.h"
21#include "coeffs/modulop.h"
22#include "coeffs/flintcf_Zn.h"
23
24typedef nmod_poly_struct *nmod_poly_ptr;
25
26/*2
27* extracts a long integer from s, returns the rest
28*/
29static const char* Eati(const char *s, int *i)
30{
31
32  if (((*s) >= '0') && ((*s) <= '9'))
33  {
34    unsigned long ii=0L;
35    do
36    {
37      ii *= 10;
38      ii += *s++ - '0';
39    }
40    while (((*s) >= '0') && ((*s) <= '9'));
41    *i=(int)ii;
42  }
43  else (*i) = 1;
44  return s;
45}
46
47static BOOLEAN CoeffIsEqual(const coeffs r, n_coeffType n, void * parameter)
48{
49  flintZn_struct *pp=(flintZn_struct*)parameter;
50  return (r->type==n) &&(r->ch==pp->ch)
51          &&(r->pParameterNames!=NULL)
52          &&(strcmp(r->pParameterNames[0],pp->name)==0);
53}
54static void KillChar(coeffs cf)
55{
56  omFree((ADDRESS)(cf->pParameterNames[0]));
57  omFreeSize(cf->pParameterNames,sizeof(char*));
58}
59static void SetChar(const coeffs r)
60{
61  // dummy
62}
63static number Mult(number a, number b, const coeffs c)
64{
65  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(*res));
66  nmod_poly_init(res,c->ch);
67  nmod_poly_mul(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
68  return (number)res;
69}
70static number Sub(number a, number b, const coeffs c)
71{
72  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
73  nmod_poly_init(res,c->ch);
74  nmod_poly_sub(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
75  return (number)res;
76}
77static number Add(number a, number b, const coeffs c)
78{
79  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
80  nmod_poly_init(res,c->ch);
81  nmod_poly_add(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
82  return (number)res;
83}
84static number Div(number a, number b, const coeffs c)
85{
86  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
87  nmod_poly_init(res,c->ch);
88  if(nmod_poly_is_zero((nmod_poly_ptr)b))
89  {
90     WerrorS(nDivBy0);
91  }
92  else
93  {
94    nmod_poly_div(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
95    nmod_poly_t mod;
96    nmod_poly_init(mod,c->ch);
97    nmod_poly_rem(mod,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
98    if (!nmod_poly_is_zero((nmod_poly_ptr)mod))
99    {
100      WerrorS("cannot divide");
101    }
102    nmod_poly_clear(mod);
103  }
104  return (number)res;
105}
106static number ExactDiv(number a, number b, const coeffs c)
107{
108  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
109  nmod_poly_init(res,c->ch);
110  if(nmod_poly_is_zero((nmod_poly_ptr)b))
111  {
112     WerrorS(nDivBy0);
113  }
114  else
115    nmod_poly_div(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
116  return (number)res;
117}
118static number IntMod(number a, number b, const coeffs c)
119{
120  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
121  nmod_poly_init(res,c->ch);
122  nmod_poly_rem(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
123  return (number)res;
124}
125static number Init (long i, const coeffs r)
126{
127  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
128  nmod_poly_init(res,r->ch);
129  i= i%r->ch;
130  if (i<0) i+=r->ch;
131  nmod_poly_set_coeff_ui(res,0,i);
132  return (number)res;
133}
134static number InitMPZ (mpz_t i, const coeffs r)
135{
136  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
137  nmod_poly_init(res,r->ch);
138  mpz_t tmp;
139  mpz_init(tmp);
140  slong ii=mpz_mod_ui(tmp,i,r->ch);
141  mpz_clear(tmp);
142  nmod_poly_set_coeff_ui(res,0,ii);
143  return (number)res;
144}
145static int Size (number n, const coeffs r)
146{
147  return nmod_poly_degree((nmod_poly_ptr)n);
148}
149static long Int (number &n, const coeffs r)
150{
151  if (nmod_poly_degree((nmod_poly_ptr)n)==0)
152  {
153    slong m;
154    m=nmod_poly_get_coeff_ui((nmod_poly_ptr)n,0);
155    return (long)m;
156  }
157  return 0;
158}
159static void MPZ(mpz_t result, number &n, const coeffs r)
160{
161  mpz_init(result);
162  if (nmod_poly_degree((nmod_poly_ptr)n)==0)
163  {
164    slong m;
165    m=nmod_poly_get_coeff_ui((nmod_poly_ptr)n,0);
166    mpz_set_ui(result,m);
167  }
168}
169static number Neg(number a, const coeffs r)
170{
171  nmod_poly_neg((nmod_poly_ptr)a,(nmod_poly_ptr)a);
172  return a;
173}
174static number Invers(number a, const coeffs r)
175{
176  if(nmod_poly_is_zero((nmod_poly_ptr)a))
177  {
178    WerrorS(nDivBy0);
179    return NULL;
180  }
181  if (nmod_poly_degree((nmod_poly_ptr)a)==0)
182  {
183    nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
184    nmod_poly_init(res,r->ch);
185    slong c=nmod_poly_get_coeff_ui((nmod_poly_ptr)a,0);
186    extern number  nvInvers (number c, const coeffs r);
187    c=(slong)nvInvers((number)c,r);
188    nmod_poly_set_coeff_ui((nmod_poly_ptr)a,0,c);
189    return (number)res;
190  }
191  else
192  {
193    WerrorS("not invertable");
194    return NULL;
195  }
196}
197static number Copy(number a, const coeffs r)
198{
199 nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
200 nmod_poly_init(res,r->ch);
201 nmod_poly_set(res,(nmod_poly_ptr)a);
202 return (number)res;
203}
204//static number RePart(number a, const coeffs r)
205//{
206//}
207//static number ImPart(number a, const coeffs r)
208//{
209//}
210static BOOLEAN IsOne (number a, const coeffs r);
211static BOOLEAN IsZero (number a, const coeffs r);
212//static void WriteLong(number &a, const coeffs r)
213//{
214//}
215static void WriteShort(number a, const coeffs r)
216{
217  //nmod_poly_print_pretty((nmod_poly_ptr)a,r->pParameterNames[0]);
218  if (IsOne(a,r)) StringAppendS("1");
219  else if (IsZero(a,r)) StringAppendS("0");
220  else
221  {
222    StringAppendS("(");
223    BOOLEAN need_plus=FALSE;
224    for(int i=nmod_poly_length((nmod_poly_ptr)a);i>=0;i--)
225    {
226      slong m=nmod_poly_get_coeff_ui((nmod_poly_ptr)a,i);
227      if (m!=0)
228      {
229        if (need_plus) StringAppendS("+");
230        need_plus=TRUE;
231        if (i>0)
232        {
233          if (m!=1) StringAppend("%d*",(int)m);
234          if (i>1)
235            StringAppend("%s^%d",r->pParameterNames[0],i);
236          else if (i==1)
237            StringAppend("%s",r->pParameterNames[0]);
238        }
239        else StringAppend("%d",(int)m);
240      }
241    }
242    StringAppendS(")");
243  }
244}
245static const char* Read(const char * st, number * a, const coeffs r)
246{
247// we only read "monomials" (i.e. [-][digits][parameter]),
248// everythings else (+,*,^,()) is left to the singular interpreter
249  const char *s=st;
250  *a=(number)omAlloc(sizeof(nmod_poly_t));
251  nmod_poly_init((nmod_poly_ptr)(*a),r->ch);
252  BOOLEAN neg=FALSE;
253  if (*s=='-') { neg=TRUE; s++;}
254  if (isdigit(*s))
255  {
256    int z;
257    s=Eati((char *)s, &z);
258    nmod_poly_set_coeff_ui((nmod_poly_ptr)(*a),0,z);
259  }
260  else if(strncmp(s,r->pParameterNames[0],strlen(r->pParameterNames[0]))==0)
261  {
262    nmod_poly_set_coeff_ui((nmod_poly_ptr)(*a),1,1);
263    s+=strlen(r->pParameterNames[0]);
264    if(isdigit(*s))
265    {
266      int i=1;
267      s=Eati(s,&i);
268      if (i!=1)
269      {
270        nmod_poly_set_coeff_ui((nmod_poly_ptr)(*a),1,0);
271        nmod_poly_set_coeff_ui((nmod_poly_ptr)(*a),i,1);
272      }
273    }
274  }
275  if (neg)
276    nmod_poly_neg((nmod_poly_ptr)(*a),(nmod_poly_ptr)(*a));
277  return s;
278}
279static void Normalize(number &a, const coeffs r)
280{
281}
282static BOOLEAN Greater (number a, number b, const coeffs r)
283{
284  if (nmod_poly_length((nmod_poly_ptr)a)>nmod_poly_length((nmod_poly_ptr)b))
285    return TRUE;
286   else if (nmod_poly_length((nmod_poly_ptr)a)<nmod_poly_length((nmod_poly_ptr)b))
287     return FALSE;
288   for(int i=nmod_poly_length((nmod_poly_ptr)a);i>=0;i--)
289   {
290     slong ac=nmod_poly_get_coeff_ui((nmod_poly_ptr)a,i);
291     slong bc=nmod_poly_get_coeff_ui((nmod_poly_ptr)b,i);
292     if (ac>bc) return TRUE;
293     else if (ac<bc) return FALSE;
294   }
295   return FALSE;
296}
297static BOOLEAN Equal (number a, number b, const coeffs r)
298{
299  return (nmod_poly_equal((nmod_poly_ptr)a,(nmod_poly_ptr)b));
300}
301static BOOLEAN IsZero (number a, const coeffs r)
302{
303  return nmod_poly_is_zero((nmod_poly_ptr)a);
304}
305static BOOLEAN IsOne (number a, const coeffs r)
306{
307  return nmod_poly_is_one((nmod_poly_ptr)a);
308}
309static BOOLEAN IsMOne (number k, const coeffs r)
310{
311  if (nmod_poly_length((nmod_poly_ptr)k)>0) return FALSE;
312  slong m=nmod_poly_get_coeff_ui((nmod_poly_ptr)k,0);
313  return (m+1==r->ch);
314}
315static BOOLEAN GreaterZero (number k, const coeffs r)
316{
317  // does it have a leading sign?
318  // no: 0 and 1 do not have, everything else is in (...)
319  return TRUE;
320}
321static void Power(number a, int i, number * result, const coeffs r)
322{
323  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
324  nmod_poly_init(res,r->ch);
325  *result=(number)res;
326  nmod_poly_pow((nmod_poly_ptr)(*result),(nmod_poly_ptr)a,i);
327}
328static number Gcd(number a, number b, const coeffs r)
329{
330  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
331  nmod_poly_init(res,r->ch);
332  nmod_poly_gcd(res,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
333  return (number)res;
334}
335static number ExtGcd(number a, number b, number *s, number *t,const coeffs r)
336{
337  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
338  nmod_poly_init(res,r->ch);
339  nmod_poly_init((nmod_poly_ptr)*s,r->ch);
340  nmod_poly_init((nmod_poly_ptr)*t,r->ch);
341  nmod_poly_xgcd(res,(nmod_poly_ptr)*s,(nmod_poly_ptr)*t,(nmod_poly_ptr)a,(nmod_poly_ptr)b);
342  return (number)res;
343}
344static number Lcm(number a, number b, const coeffs r)
345{
346  WerrorS("not yet: Lcm");
347}
348static void Delete(number * a, const coeffs r)
349{
350  if ((*a)!=NULL)
351  {
352    nmod_poly_clear((nmod_poly_ptr)*a);
353    omFree(*a);
354    *a=NULL;
355  }
356}
357static nMapFunc SetMap(const coeffs src, const coeffs dst)
358{
359  WerrorS("not yet: SetMap");
360  return NULL;
361}
362//static void InpMult(number &a, number b, const coeffs r)
363//{
364//}
365//static void InpAdd(number &a, number b, const coeffs r)
366//{
367//}
368static number Init_bigint(number i, const coeffs dummy, const coeffs dst)
369{
370  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
371  nmod_poly_init(res,dst->ch);
372  long ii;
373  if (SR_HDL(i) & SR_INT)
374  {
375    ii=SR_TO_INT(i) % dst->ch;
376  }
377  else
378  {
379    mpz_t tmp;
380    mpz_init(tmp);
381    ii=mpz_mod_ui(tmp,i->z,dst->ch);
382    mpz_clear(tmp);
383  }
384  if (ii<0) ii+=dst->ch;
385  nmod_poly_set_coeff_ui(res,0,ii);
386  return (number)res;
387}
388static number Farey(number p, number n, const coeffs)
389{
390  WerrorS("not yet: Farey");
391}
392static number ChineseRemainder(number *x, number *q,int rl, BOOLEAN sym,CFArray &inv_cache,const coeffs)
393{
394  WerrorS("not yet: ChineseRemainder");
395}
396static int ParDeg(number x,const coeffs r)
397{
398  return nmod_poly_degree((nmod_poly_ptr)x);
399}
400static number Parameter(const int i, const coeffs r)
401{
402  nmod_poly_ptr res=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
403  nmod_poly_init(res,r->ch);
404  nmod_poly_set_coeff_ui(res,1,1);
405  return (number)res;
406}
407// cfClearContent
408// cfClearDenominators
409static number ConvFactoryNSingN( const CanonicalForm n, const coeffs r)
410{
411}
412static CanonicalForm ConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs r )
413{
414  WerrorS("not yet: ConvSingNFactoryN");
415}
416static char * CoeffName(const coeffs r)
417{
418  STATIC_VAR char CoeffName_flint_Zn[20];
419  sprintf(CoeffName_flint_Zn,"flint:Z/%d[%s]",r->ch,r->pParameterNames[0]);
420  return (char*)CoeffName_flint_Zn;
421}
422coeffs flintZnInitCfByName(char *s,n_coeffType n)
423{
424  const char start[]="flint:Z/";
425  const int start_len=strlen(start);
426  if (strncmp(s,start,start_len)==0)
427  {
428    s+=start_len;
429    int p;
430    char st[10];
431    int l=sscanf(s,"%d[%s",&p,st);
432    if (l==2)
433    {
434      flintZn_struct info;
435      info.ch=p;
436      while (st[strlen(st)-1]==']') st[strlen(st)-1]='\0';
437      info.name=st;
438      return nInitChar(n,(void*)&info);
439    }
440  }
441  return NULL;
442}
443static void WriteFd(number a, const ssiInfo *d, const coeffs)
444{
445  // format: len a_len .. a_0
446  nmod_poly_ptr aa=(nmod_poly_ptr)a;
447  int l=nmod_poly_length(aa);
448  fprintf(d->f_write,"%d ",l);
449  for(int i=l; i>=0; i--)
450  {
451    ulong ul=nmod_poly_get_coeff_ui(aa,i);
452    fprintf(d->f_write,"%lu ", ul);
453  }
454}
455
456static number ReadFd(const ssiInfo *d, const coeffs r)
457{
458  // format: len a_len .. a_0
459  nmod_poly_ptr aa=(nmod_poly_ptr)omAlloc(sizeof(nmod_poly_t));
460  nmod_poly_init(aa,r->ch);
461  int l=s_readint(d->f_read);
462  unsigned long ul;
463  for (int i=l;i>=0;i--)
464  {
465    unsigned long ul=s_readlong(d->f_read);
466    nmod_poly_set_coeff_ui(aa,i,ul);
467  }
468  return (number)aa;
469}
470#ifdef LDEBUG
471static BOOLEAN DBTest(number a, const char *f, const int l, const coeffs r)
472{
473  return TRUE;
474}
475#endif
476BOOLEAN flintZn_InitChar(coeffs cf, void * infoStruct)
477{
478  flintZn_struct *pp=(flintZn_struct*)infoStruct;
479  cf->ch=pp->ch;
480
481  cf->cfCoeffName    = CoeffName;
482  cf->nCoeffIsEqual  = CoeffIsEqual;
483  cf->cfKillChar     = KillChar;
484  cf->cfSetChar      = SetChar;
485  cf->cfMult         = Mult;
486  cf->cfSub          = Sub;
487  cf->cfAdd          = Add;
488  cf->cfDiv          = Div;
489  cf->cfExactDiv     = ExactDiv; // ???
490  cf->cfInit         = Init;
491  cf->cfInitMPZ      = InitMPZ;
492  cf->cfSize         = Size;
493  cf->cfInt          = Int;
494  cf->cfMPZ          = MPZ;
495  cf->cfInpNeg       = Neg;
496  cf->cfInvers       = Invers;
497  cf->cfCopy         = Copy;
498  cf->cfRePart       = Copy;
499  // default: cf->cfImPart       = ndReturn0;
500  cf->cfWriteLong    = WriteShort; //WriteLong;
501  cf->cfWriteShort = WriteShort;
502  cf->cfRead         = Read;
503  cf->cfNormalize    = Normalize;
504
505  //cf->cfDivComp=
506  //cf->cfIsUnit=
507  //cf->cfGetUnit=
508  //cf->cfDivBy=
509
510  cf->cfGreater=Greater;
511  cf->cfEqual  =Equal;
512  cf->cfIsZero =IsZero;
513  cf->cfIsOne  =IsOne;
514  cf->cfIsMOne =IsMOne;
515  cf->cfGreaterZero=GreaterZero;
516
517  cf->cfPower        = Power;
518  //default: cf->cfGetDenom     = GetDenom;
519  //default: cf->cfGetNumerator = GetNumerator;
520  cf->cfGcd          = Gcd;
521  cf->cfExtGcd         = ExtGcd;
522  cf->cfLcm          = Lcm;
523  cf->cfDelete       = Delete;
524  cf->cfSetMap       = SetMap;
525  // default: cf->cfInpMult
526  // default: cf->cfInpAdd
527  cf->cfFarey        =Farey;
528  cf->cfChineseRemainder=ChineseRemainder;
529  cf->cfParDeg = ParDeg;
530  cf->cfParameter = Parameter;
531  //  cf->cfClearContent = ClearContent;
532  //  cf->cfClearDenominators = ClearDenominators;
533  cf->convFactoryNSingN=ConvFactoryNSingN;
534  cf->convSingNFactoryN=ConvSingNFactoryN;
535  cf->cfWriteFd = WriteFd;
536  cf->cfReadFd  = ReadFd;
537#ifdef LDEBUG
538  cf->cfDBTest       = DBTest;
539#endif
540
541  cf->iNumberOfParameters = 1;
542  char **pn=(char**)omAlloc0(sizeof(char*));
543  pn[0]=(char*)omStrDup(pp->name);
544  cf->pParameterNames = (const char **)pn;
545  cf->has_simple_Inverse= FALSE;
546  cf->has_simple_Alloc= FALSE;
547  cf->is_field=FALSE;
548
549  return FALSE;
550}
551#endif
Note: See TracBrowser for help on using the repository browser.