source: git/libpolys/coeffs/flintcf_Q.cc @ 14a7ef

spielwiese
Last change on this file since 14a7ef was 14a7ef, checked in by Hans Schoenemann <hannes@…>, 4 years ago
flintQ: rational functions
  • Property mode set to 100644
File size: 15.7 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/*
5* ABSTRACT: flint: fmpq_poly
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/fmpq_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
22typedef fmpq_poly_struct *fmpq_poly_ptr;
23typedef fmpz *fmpz_ptr;
24/*2
25* extracts a long integer from s, returns the rest
26*/
27static char * nlEatLong(char *s, mpz_ptr i)
28{
29  const char * start=s;
30
31  while (*s >= '0' && *s <= '9') s++;
32  if (*s=='\0')
33  {
34    mpz_set_str(i,start,10);
35  }
36  else
37  {
38    char c=*s;
39    *s='\0';
40    mpz_set_str(i,start,10);
41    *s=c;
42  }
43  return s;
44}
45
46static BOOLEAN CoeffIsEqual(const coeffs r, n_coeffType n, void * parameter)
47{
48  return (r->type==n);
49}
50static void SetChar(const coeffs r)
51{
52  // dummy
53}
54static number Mult(number a, number b, const coeffs c)
55{
56  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
57  fmpq_poly_init(res);
58  fmpq_poly_mul(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
59  return (number)res;
60}
61static number Sub(number a, number b, const coeffs c)
62{
63  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
64  fmpq_poly_init(res);
65  fmpq_poly_sub(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
66  return (number)res;
67}
68static number Add(number a, number b, const coeffs c)
69{
70  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
71  fmpq_poly_init(res);
72  fmpq_poly_add(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
73  return (number)res;
74}
75static number Div(number a, number b, const coeffs c)
76{
77  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
78  fmpq_poly_init(res);
79  if(fmpq_poly_is_zero((fmpq_poly_ptr)b))
80  {
81     WerrorS(nDivBy0);
82  }
83  else
84  {
85    fmpq_poly_div(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
86    fmpq_poly_t mod;
87    fmpq_poly_init(mod);
88    fmpq_poly_rem(mod,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
89    if (!fmpq_poly_is_zero((fmpq_poly_ptr)mod))
90    {
91      WerrorS("cannot divide");
92    }
93    fmpq_poly_clear(mod);
94  }
95  return (number)res;
96}
97static number ExactDiv(number a, number b, const coeffs c)
98{
99  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
100  fmpq_poly_init(res);
101  if(fmpq_poly_is_zero((fmpq_poly_ptr)b))
102  {
103     WerrorS(nDivBy0);
104  }
105  else
106    fmpq_poly_div(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
107  return (number)res;
108}
109static number IntMod(number a, number b, const coeffs c)
110{
111  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
112  fmpq_poly_init(res);
113  fmpq_poly_rem(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
114  return (number)res;
115}
116static number Init (long i, const coeffs r)
117{
118  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
119  fmpq_poly_init(res);
120  fmpq_poly_set_si(res,i);
121  return (number)res;
122}
123static number InitMPZ (mpz_t i, const coeffs r)
124{
125  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
126  fmpq_poly_init(res);
127  fmpq_poly_set_mpz(res,i);
128  return (number)res;
129}
130static int Size (number n, const coeffs r)
131{
132  return fmpq_poly_degree((fmpq_poly_ptr)n);
133}
134static long Int (number &n, const coeffs r)
135{
136  if (fmpq_poly_degree((fmpq_poly_ptr)n)==0)
137  {
138    mpq_t m;
139    mpq_init(m);
140    fmpq_poly_get_coeff_mpq(m,(fmpq_poly_ptr)n,0);
141    mpz_t num,den;
142    mpz_init(num);
143    mpz_init(den);
144    mpq_get_num(num,m);
145    mpq_get_den(den,m);
146    long nl=mpz_get_si(num);
147    if (mpz_cmp_si(num,nl)!=0) nl=0;
148    long dl=mpz_get_si(den);
149    if ((dl!=1)||(mpz_cmp_si(den,dl)!=0)) nl=0;
150    mpz_clear(num);
151    mpz_clear(den);
152    mpq_clear(m);
153    return nl;
154  }
155  return 0;
156}
157static void MPZ(mpz_t result, number &n, const coeffs r)
158{
159  mpz_init(result);
160  if (fmpq_poly_degree((fmpq_poly_ptr)n)==0)
161  {
162    mpq_t m;
163    mpq_init(m);
164    fmpq_poly_get_coeff_mpq(m,(fmpq_poly_ptr)n,0);
165    mpz_t den;
166    mpz_init(den);
167    mpq_get_num(result,m);
168    mpq_get_den(den,m);
169    int dl=(int)mpz_get_si(den);
170    if ((dl!=1)||(mpz_cmp_si(den,(long)dl)!=0)) mpz_set_ui(result,0);
171    mpz_clear(den);
172    mpq_clear(m);
173  }
174}
175static number Neg(number a, const coeffs r)
176{
177  fmpq_poly_neg((fmpq_poly_ptr)a,(fmpq_poly_ptr)a);
178  return a;
179}
180static number Invers(number a, const coeffs r)
181{
182  if(fmpq_poly_is_zero((fmpq_poly_ptr)a))
183  {
184    WerrorS(nDivBy0);
185    return NULL;
186  }
187  if (fmpq_poly_degree((fmpq_poly_ptr)a)==0)
188  {
189    fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
190    fmpq_poly_init(res);
191    fmpq_poly_inv(res,(fmpq_poly_ptr)a);
192    return (number)res;
193  }
194  else
195  {
196    WerrorS("not invertable");
197    return NULL;
198  }
199}
200static number Copy(number a, const coeffs r)
201{
202 fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
203 fmpq_poly_init(res);
204 fmpq_poly_set(res,(fmpq_poly_ptr)a);
205 return (number)res;
206}
207//static number RePart(number a, const coeffs r)
208//{
209//}
210//static number ImPart(number a, const coeffs r)
211//{
212//}
213static BOOLEAN IsOne (number a, const coeffs r);
214static BOOLEAN IsZero (number a, const coeffs r);
215//static void WriteLong(number &a, const coeffs r)
216//{
217//}
218static void WriteShort(number a, const coeffs r)
219{
220  //fmpq_poly_print_pretty((fmpq_poly_ptr)a,r->pParameterNames[0]);
221  if (IsOne(a,r)) StringAppendS("1");
222  else if (IsZero(a,r)) StringAppendS("0");
223  else
224  {
225  StringAppendS("(");
226  mpq_t m;
227  mpq_init(m);
228  mpz_t num,den;
229  mpz_init(num);
230  mpz_init(den);
231  BOOLEAN need_plus=FALSE;
232  for(int i=fmpq_poly_length((fmpq_poly_ptr)a);i>=0;i--)
233  {
234    fmpq_poly_get_coeff_mpq(m,(fmpq_poly_ptr)a,i);
235    mpq_get_num(num,m);
236    mpq_get_den(den,m);
237    if (mpz_sgn1(num)!=0)
238    {
239      if (need_plus && (mpz_sgn1(num)>0))
240        StringAppendS("+");
241      need_plus=TRUE;
242      int l=mpz_sizeinbase(num,10);
243      l=si_max(l,(int)mpz_sizeinbase(den,10));
244      l+=2;
245      char *s=(char*)omAlloc(l);
246      char *z=mpz_get_str(s,10,num);
247      if ((i==0)
248      ||(mpz_cmp_si(num,1)!=0)
249      ||(mpz_cmp_si(den,1)!=0))
250      {
251        StringAppendS(z);
252        if (mpz_cmp_si(den,1)!=0)
253        {
254          StringAppendS("/");
255          z=mpz_get_str(s,10,den);
256          StringAppendS(z);
257        }
258        if (i!=0) StringAppendS("*");
259      }
260      if (i>1)
261        StringAppend("%s^%d",r->pParameterNames[0],i);
262      else if (i==1)
263        StringAppend("%s",r->pParameterNames[0]);
264    }
265  }
266  mpz_clear(den);
267  mpz_clear(num);
268  mpq_clear(m);
269  StringAppendS(")");
270  }
271}
272static const char* Read(const char * st, number * a, const coeffs r)
273{
274// we only read "monomials" (i.e. [-][digits][parameter]),
275// everythings else (+,*,^,()) is left to the singular interpreter
276  char *s=(char *)st;
277  *a=(number)omAlloc(sizeof(fmpq_poly_t));
278  fmpq_poly_init((fmpq_poly_ptr)(*a));
279  BOOLEAN neg=FALSE;
280  if (*s=='-') { neg=TRUE; s++;}
281  if (isdigit(*s))
282  {
283    mpz_t z;
284    mpz_init(z);
285    s=nlEatLong((char *)s, z);
286    fmpq_poly_set_mpz((fmpq_poly_ptr)(*a),z);
287    if (*s == '/')
288    {
289      s++;
290      s=nlEatLong((char *)s, z);
291      fmpq_poly_scalar_div_mpz((fmpq_poly_ptr)(*a),(fmpq_poly_ptr)(*a),z);
292    }
293    mpz_clear(z);
294  }
295  else if(strncmp(s,r->pParameterNames[0],strlen(r->pParameterNames[0]))==0)
296  {
297    fmpq_poly_set_coeff_si((fmpq_poly_ptr)(*a),1,1);
298    s+=strlen(r->pParameterNames[0]);
299    if(isdigit(*s))
300    {
301      int i=1;
302      s=nEati(s,&i,0);
303      if (i!=1)
304      {
305        fmpq_poly_set_coeff_si((fmpq_poly_ptr)(*a),1,0);
306        fmpq_poly_set_coeff_si((fmpq_poly_ptr)(*a),i,1);
307      }
308    }
309  }
310  if (neg)
311    fmpq_poly_neg((fmpq_poly_ptr)(*a),(fmpq_poly_ptr)(*a));
312  return s;
313}
314static void Normalize(number &a, const coeffs r)
315{
316  fmpq_poly_canonicalise((fmpq_poly_ptr)a);
317}
318static BOOLEAN Greater (number a, number b, const coeffs r)
319{
320  return (fmpq_poly_cmp((fmpq_poly_ptr)a,(fmpq_poly_ptr)b)>0);
321}
322static BOOLEAN Equal (number a, number b, const coeffs r)
323{
324  return (fmpq_poly_equal((fmpq_poly_ptr)a,(fmpq_poly_ptr)b));
325}
326static BOOLEAN IsZero (number a, const coeffs r)
327{
328  return fmpq_poly_is_zero((fmpq_poly_ptr)a);
329}
330static BOOLEAN IsOne (number a, const coeffs r)
331{
332  return fmpq_poly_is_one((fmpq_poly_ptr)a);
333}
334static BOOLEAN IsMOne (number k, const coeffs r)
335{
336  if (fmpq_poly_length((fmpq_poly_ptr)k)>0) return FALSE;
337  fmpq_poly_canonicalise((fmpq_poly_ptr)k);
338  mpq_t m;
339  mpq_init(m);
340  fmpq_poly_get_coeff_mpq(m,(fmpq_poly_ptr)k,0);
341  mpz_t num,den;
342  mpz_init(num);
343  mpq_get_num(num,m);
344  BOOLEAN result=TRUE;
345  if (mpz_cmp_si(num,(long)-1)!=0) result=FALSE;
346  else
347  {
348    mpz_init(den);
349    mpq_get_den(den,m);
350    int dl=(int)mpz_get_si(den);
351    if ((dl!=1)||(mpz_cmp_si(den,(long)dl)!=0)) result=FALSE;
352    mpz_clear(den);
353  }
354  mpz_clear(num);
355  mpq_clear(m);
356  return (result);
357}
358static BOOLEAN GreaterZero (number k, const coeffs r)
359{
360  // does it have a leading sign?
361  // no: 0 and 1 do not have, everything else is in (...)
362  return TRUE;
363}
364static void Power(number a, int i, number * result, const coeffs r)
365{
366  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
367  fmpq_poly_init(res);
368  *result=(number)res;
369  fmpq_poly_pow((fmpq_poly_ptr)(*result),(fmpq_poly_ptr)a,i);
370}
371static number GetDenom(number &n, const coeffs r)
372{
373 fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
374 fmpq_poly_init(res);
375 fmpz_ptr den=fmpq_poly_denref(res);
376 fmpq_poly_set_fmpz(res,den);
377 return (number)res;
378}
379static number GetNumerator(number &n, const coeffs r)
380{
381 fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
382 fmpq_poly_init(res);
383 fmpq_poly_set(res,(fmpq_poly_ptr)n);
384 fmpz_ptr den=fmpq_poly_denref(res);
385 fmpq_poly_scalar_mul_fmpz(res,res,den);
386 return (number)res;
387}
388static number Gcd(number a, number b, const coeffs r)
389{
390  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
391  fmpq_poly_init(res);
392  fmpq_poly_gcd(res,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
393  return (number)res;
394}
395static number ExtGcd(number a, number b, number *s, number *t,const coeffs r)
396{
397  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
398  fmpq_poly_init(res);
399  fmpq_poly_init((fmpq_poly_ptr)*s);
400  fmpq_poly_init((fmpq_poly_ptr)*t);
401  fmpq_poly_xgcd(res,(fmpq_poly_ptr)*s,(fmpq_poly_ptr)*t,(fmpq_poly_ptr)a,(fmpq_poly_ptr)b);
402  return (number)res;
403}
404static number Lcm(number a, number b, const coeffs r)
405{
406  WerrorS("not yet: Lcm");
407}
408static void Delete(number * a, const coeffs r)
409{
410  if ((*a)!=NULL)
411  {
412    fmpq_poly_clear((fmpq_poly_ptr)*a);
413    omFree(*a);
414    *a=NULL;
415  }
416}
417static nMapFunc SetMap(const coeffs src, const coeffs dst)
418{
419  WerrorS("not yet: SetMap");
420  return NULL;
421}
422//static void InpMult(number &a, number b, const coeffs r)
423//{
424//}
425//static void InpAdd(number &a, number b, const coeffs r)
426//{
427//}
428static number Init_bigint(number i, const coeffs dummy, const coeffs dst)
429{
430  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
431  fmpq_poly_init(res);
432  if (SR_HDL(i) & SR_INT)
433  {
434    fmpq_poly_set_si(res,SR_TO_INT(i));
435  }
436  else
437    fmpq_poly_set_mpz(res,i->z);
438  return (number)res;
439}
440static number Farey(number p, number n, const coeffs)
441{
442  WerrorS("not yet: Farey");
443}
444static number ChineseRemainder(number *x, number *q,int rl, BOOLEAN sym,CFArray &inv_cache,const coeffs)
445{
446  WerrorS("not yet: ChineseRemainder");
447}
448static int ParDeg(number x,const coeffs r)
449{
450  return fmpq_poly_degree((fmpq_poly_ptr)x);
451}
452static number Parameter(const int i, const coeffs r)
453{
454  fmpq_poly_ptr res=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
455  fmpq_poly_init(res);
456  fmpq_poly_set_coeff_si(res,1,1);
457  return (number)res;
458}
459static void WriteFd(number a, const ssiInfo *d, const coeffs)
460{
461  // format: len a_len(num den) .. a_0
462  fmpq_poly_ptr aa=(fmpq_poly_ptr)a;
463  int l=fmpq_poly_length(aa);
464  fprintf(d->f_write,"%d ",l);
465  mpq_t m;
466  mpq_init(m);
467  mpz_t num,den;
468  mpz_init(num);
469  mpz_init(den);
470  for(int i=l; i>=0; i--)
471  {
472    fmpq_poly_get_coeff_mpq(m,(fmpq_poly_ptr)a,i);
473    mpq_get_num(num,m);
474    mpq_get_den(den,m);
475    mpz_out_str (d->f_write,SSI_BASE, num);
476    fputc(' ',d->f_write);
477    mpz_out_str (d->f_write,SSI_BASE, den);
478    fputc(' ',d->f_write);
479  }
480  mpz_clear(den);
481  mpz_clear(num);
482  mpq_clear(m);
483}
484static number ReadFd(const ssiInfo *d, const coeffs)
485{
486  // format: len a_len .. a_0
487  fmpq_poly_ptr aa=(fmpq_poly_ptr)omAlloc(sizeof(fmpq_poly_t));
488  fmpq_poly_init(aa);
489  int l=s_readint(d->f_read);
490  mpz_t nm;
491  mpz_init(nm);
492  mpq_t m;
493  mpq_init(m);
494  for (int i=l;i>=0;i--)
495  {
496
497    s_readmpz_base (d->f_read,nm, SSI_BASE);
498    mpq_set_num(m,nm);
499    s_readmpz_base (d->f_read,nm, SSI_BASE);
500    mpq_set_den(m,nm);
501    fmpq_poly_set_coeff_mpq(aa,i,m);
502  }
503  mpz_clear(nm);
504  mpq_clear(m);
505  return (number)aa;
506}
507// cfClearContent
508// cfClearDenominators
509static number ConvFactoryNSingN( const CanonicalForm n, const coeffs r)
510{
511  WerrorS("not yet: ConvFactoryNSingN");
512}
513static CanonicalForm ConvSingNFactoryN( number n, BOOLEAN setChar, const coeffs r )
514{
515  WerrorS("not yet: ConvSingNFactoryN");
516}
517char * CoeffName(const coeffs r)
518{
519  STATIC_VAR char CoeffName_flint_Q[20];
520  sprintf(CoeffName_flint_Q,"flintQp[%s]",r->pParameterNames[0]);
521  return (char*)CoeffName_flint_Q;
522
523}
524static char* CoeffString(const coeffs r)
525{
526  char *buf=(char*)omAlloc(12+strlen(r->pParameterNames[0]));
527  sprintf(buf,"flintQp(\"%s\")",r->pParameterNames[0]);
528  return buf;
529}
530static void CoeffWrite(const coeffs r, BOOLEAN details)
531{
532  PrintS(CoeffName(r));
533}
534coeffs flintQInitCfByName(char *s,n_coeffType n)
535{
536  const char start[]="flintQp[";
537  const int start_len=strlen(start);
538  if (strncmp(s,start,start_len)==0)
539  {
540    s+=start_len;
541    char st[10];
542    int l=sscanf(s,"%s",st);
543    if (l==1)
544    {
545      while (st[strlen(st)-1]==']') st[strlen(st)-1]='\0';
546      return nInitChar(n,(void*)st);
547    }
548  }
549  return NULL;
550}
551#ifdef LDEBUG
552static BOOLEAN DBTest(number a, const char *f, const int l, const coeffs r)
553{
554  return TRUE;
555}
556#endif
557static void KillChar(coeffs cf)
558{
559  omFree((ADDRESS)(cf->pParameterNames[0]));
560  omFreeSize(cf->pParameterNames,sizeof(char*));
561}
562BOOLEAN flintQ_InitChar(coeffs cf, void * infoStruct)
563{
564  char *pp=(char*)infoStruct;
565  cf->cfCoeffString  = CoeffString;
566  cf->cfCoeffName    = CoeffName;
567  cf->cfCoeffWrite   = CoeffWrite;
568  cf->nCoeffIsEqual  = CoeffIsEqual;
569  cf->cfKillChar     = KillChar;
570  cf->cfSetChar      = SetChar;
571  cf->ch=0; //char 0
572  cf->cfMult         = Mult;
573  cf->cfSub          = Sub;
574  cf->cfAdd          = Add;
575  cf->cfDiv          = Div;
576  cf->cfExactDiv     = ExactDiv; // ???
577  cf->cfInit         =Init;
578  cf->cfInitMPZ      =InitMPZ;
579  cf->cfSize         = Size;
580  cf->cfInt          = Int;
581  cf->cfMPZ          = MPZ;
582  cf->cfInpNeg       = Neg;
583  cf->cfInvers       = Invers;
584  cf->cfCopy         = Copy;
585  cf->cfRePart       = Copy;
586  // default: cf->cfImPart       = ndReturn0;
587  cf->cfWriteLong    = WriteShort; //WriteLong;
588  cf->cfWriteShort = WriteShort;
589  cf->cfRead         = Read;
590  cf->cfNormalize    = Normalize;
591
592  //cf->cfDivComp=
593  //cf->cfIsUnit=
594  //cf->cfGetUnit=
595  //cf->cfDivBy=
596
597  cf->cfGreater=Greater;
598  cf->cfEqual  =Equal;
599  cf->cfIsZero =IsZero;
600  cf->cfIsOne  =IsOne;
601  cf->cfIsMOne =IsMOne;
602  cf->cfGreaterZero=GreaterZero;
603
604  cf->cfPower        = Power;
605  cf->cfGetDenom     = GetDenom;
606  cf->cfGetNumerator = GetNumerator;
607  cf->cfGcd          = Gcd;
608  cf->cfExtGcd         = ExtGcd;
609  cf->cfLcm          = Lcm;
610  cf->cfDelete       = Delete;
611  cf->cfSetMap       = SetMap;
612  // default: cf->cfInpMult
613  // default: cf->cfInpAdd
614  cf->cfFarey        =Farey;
615  cf->cfChineseRemainder=ChineseRemainder;
616  cf->cfParDeg = ParDeg;
617  cf->cfParameter = Parameter;
618  //  cf->cfClearContent = ClearContent;
619  //  cf->cfClearDenominators = ClearDenominators;
620  cf->convFactoryNSingN=ConvFactoryNSingN;
621  cf->convSingNFactoryN=ConvSingNFactoryN;
622  cf->cfWriteFd  = WriteFd;
623  cf->cfReadFd = ReadFd;
624#ifdef LDEBUG
625  cf->cfDBTest       = DBTest;
626#endif
627
628  cf->iNumberOfParameters = 1;
629  char **pn=(char**)omAlloc0(sizeof(char*));
630  pn[0]=omStrDup(pp);
631  cf->pParameterNames = (const char **)pn;
632  cf->has_simple_Inverse= FALSE;
633  cf->has_simple_Alloc= FALSE;
634  cf->is_field=FALSE;
635
636  return FALSE;
637}
638#endif
Note: See TracBrowser for help on using the repository browser.