source: git/Singular/GMPrat.cc @ f9aada

spielwiese
Last change on this file since f9aada was f9aada, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: simplified GMPrat git-svn-id: file:///usr/local/Singular/svn/trunk@5566 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.6 KB
Line 
1// ----------------------------------------------------------------------------
2//  GMPrat.cc
3//  begin of file
4//  originally written by Gerd Sussner, sussner@mi.uni-erlangen.de
5//  copied by Stephan Endrass, endrass@mathematik.uni-mainz.de
6//  23.7.99
7// ----------------------------------------------------------------------------
8
9#define  GMPRAT_CC
10
11#include "mod2.h"
12
13#ifdef HAVE_SPECTRUM
14
15#ifdef   GMPRAT_PRINT
16#include <iostream.h>
17#ifndef  GMPRAT_IOSTREAM
18#include <stdio.h>
19#endif
20#endif
21
22#include <stdlib.h>
23#include <float.h>
24#include <math.h>
25#include <ctype.h>
26#include <string.h>
27
28#include "GMPrat.h"
29
30// ----------------------------------------------------------------------------
31//  Miscellaneous
32// ----------------------------------------------------------------------------
33
34Rational Rational::save;    // dummy variable
35
36// ----------------------------------------------------------------------------
37//  disconnect a rational from its reference
38// ----------------------------------------------------------------------------
39
40void    Rational::disconnect( )
41{
42    if( p->n>1)
43    {
44        p->n--;
45        p = new rep;
46    }
47    else
48    {
49        mpq_clear(p->rat);
50    }
51    mpq_init(p->rat);
52}
53
54// ----------------------------------------------------------------------------
55//  Constructors
56// ----------------------------------------------------------------------------
57
58Rational::Rational( )
59{
60    p = new rep;
61    mpq_init( p->rat );
62}
63
64Rational::Rational( long int a )
65{
66    p = new rep;
67    mpq_init( p->rat );
68    mpq_set_si( p->rat,a,1 );
69}
70
71Rational::Rational( unsigned long int a )
72{
73    p = new rep;
74    mpq_init( p->rat );
75    mpq_set_ui( p->rat,a,1 );
76}
77
78Rational::Rational( int a )
79{
80    p = new rep;
81    mpq_init( p->rat );
82    mpq_set_si( p->rat,(long int)a,1 );
83}
84
85Rational::Rational( unsigned int a )
86{
87    p = new rep;
88    mpq_init( p->rat );
89    mpq_set_ui( p->rat,(unsigned long)a,1 );
90}
91
92Rational::Rational( const Rational& a )
93{
94    a.p->n++;
95    p=a.p;
96}
97
98Rational::Rational( double a )
99{
100    mpz_t   h1,h2;
101    int     i=0;
102    double  aa=a;
103
104    p = new rep;
105    mpq_init( p->rat );
106    mpz_init_set_ui( h1,1 );
107
108    while( fabs( 10.0*aa ) < DBL_MAX && i<DBL_DIG )
109    {
110        aa *= 10.;
111        mpz_mul_ui( h1,h1,10 );
112        i++;
113    }
114    mpz_init_set_d( h2,aa );
115    mpq_set_num( p->rat,h2 );
116    mpq_set_den( p->rat,h1 );
117    mpq_canonicalize( p->rat );
118    mpz_clear( h1 );
119    mpz_clear( h2 );
120}
121
122Rational::Rational(float a)
123{
124    mpz_t h1,h2;
125    int i=0;
126    double aa=a;
127
128    p=new rep;
129    mpq_init(p->rat);
130    mpz_init_set_ui(h1,1);
131    while(fabs(10.*aa) < FLT_MAX && i<FLT_DIG){
132        aa*=10.;
133        mpz_mul_ui(h1,h1,10);
134        i++;
135    }
136    mpz_init_set_d(h2,aa);
137    mpq_set_num(p->rat,h2);
138    mpq_set_den(p->rat,h1);
139    mpq_canonicalize(p->rat);
140    mpz_clear(h1);
141    mpz_clear(h2);
142}
143
144// ----------------------------------------------------------------------------
145//  Create a Rational from a string like "+1234.5678e-1234"
146// ----------------------------------------------------------------------------
147
148Rational::Rational( char *s )
149{
150    mpz_t   h1;
151    int     exp=0,i=0;
152    char    *s1,*s2,*ss;
153
154    ss = new char[strlen(s)+1];
155    strcpy( ss,s );
156    s1 = ss;
157    p = new rep;
158    mpq_init( p->rat );
159    if( isdigit(s1[0]) || s1[0]=='-' || s1[0]=='+' )
160    {
161        if (s1[0]=='+') ++s1;
162        if (s1[0]=='-') ++i;
163
164        while( isdigit( s1[i] ) )
165        {
166            ++i;
167        }
168        if (s1[i]=='.')
169        {
170            ++i;
171            while( isdigit( s1[i] ) )
172            {
173                s1[i-1]=s1[i];
174                ++i;
175                --exp;
176            }
177            s1[i-1]='\0';
178        }
179        if (s1[i]=='e' || s1[i]=='E')
180        {
181            s2=s1+i+1;
182        }
183        else
184            s2="";
185
186        s1[i]='\0';
187        i=exp+atoi(s2);
188        mpz_init_set_str(h1,s1,10);
189        delete[] ss;
190        mpq_set_z(p->rat,h1);
191        mpq_set_ui(save.p->rat,10,1);
192        if (i>0)
193        {
194            for(int j=0;j<i;j++)
195                mpq_mul(p->rat,p->rat,save.p->rat);
196        }
197        else if (i<0)
198        {
199            for(int j=0;j>i;j--)
200                mpq_div(p->rat,p->rat,save.p->rat);
201        }
202        mpq_canonicalize(p->rat);
203        mpz_clear(h1);
204    }
205}
206
207// ----------------------------------------------------------------------------
208//  Constructors with two arguments: numerator and denominator
209// ----------------------------------------------------------------------------
210
211Rational::Rational(long int a, unsigned long int b)
212{
213    p = new rep;
214    mpq_init( p->rat );
215    mpq_set_si( p->rat,a,b );
216    mpq_canonicalize( p->rat );
217}
218
219Rational::Rational(unsigned long int a, unsigned long int b)
220{
221    p=new rep;
222    mpq_init(p->rat);
223    mpq_set_ui(p->rat,a,b);
224    mpq_canonicalize(p->rat);
225}
226
227Rational::Rational(int a, unsigned int b)
228{
229    p=new rep;
230    mpq_init(p->rat);
231    mpq_set_si(p->rat,(long int) a,(unsigned long int) b);
232    mpq_canonicalize(p->rat);
233}
234
235Rational::Rational(unsigned int a, unsigned int b)
236{
237    p=new rep;
238    mpq_init(p->rat);
239    mpq_set_ui(p->rat,(unsigned long) a,(unsigned long int) b);
240    mpq_canonicalize(p->rat);
241}
242
243Rational::Rational(const Rational& a, const Rational& b)
244{
245    p=new rep;
246    mpq_init(p->rat);
247    mpq_div(p->rat, a.p->rat, b.p->rat);
248}
249
250Rational::Rational(long int a, long int b)
251{
252    if (b<0) a=-a;
253    p=new rep;
254    mpq_init(p->rat);
255    mpq_set_si(p->rat,a,abs(b));
256    mpq_canonicalize(p->rat);
257}
258
259Rational::Rational(int a, int b)
260{
261    if (b<0) a=-a;
262    p=new rep;
263    mpq_init(p->rat);
264    mpq_set_si(p->rat,(long int) a,(unsigned long int) abs(b));
265    mpq_canonicalize(p->rat);
266}
267
268Rational::Rational(char *sn, char *sd)
269{
270  Rational
271    h=sd;
272
273  p=new rep;
274  mpq_init(p->rat);
275  *this=sn;
276  mpq_div(p->rat,p->rat,h.p->rat);
277}
278
279// ----------------------------------------------------------------------------
280//  Destructor
281// ----------------------------------------------------------------------------
282
283Rational::~Rational()
284{
285  if (--p->n==0){
286    mpq_clear(p->rat);
287    delete p;
288  }
289}
290
291// ----------------------------------------------------------------------------
292//  Assignment operators
293// ----------------------------------------------------------------------------
294Rational& Rational::operator=(int a)
295{
296  disconnect();
297  mpq_set_si(p->rat,(long int) a,1);
298  return *this;
299}
300
301Rational& Rational::operator=(double a)
302{
303  mpz_t
304    h1,
305    h2;
306  int
307    i=0;
308  double
309    aa=a;
310
311  disconnect();
312  mpz_init_set_ui(h1,1);
313  while(fabs(10.*aa) < DBL_MAX && i<DBL_DIG){
314    aa*=10.;
315    mpz_mul_ui(h1,h1,10);
316    i++;
317  }
318  mpz_init_set_d(h2,aa);
319  mpq_set_num(p->rat,h2);
320  mpq_set_den(p->rat,h1);
321  mpq_canonicalize(p->rat);
322  mpz_clear(h1);
323  mpz_clear(h2);
324  return *this;
325}
326
327Rational& Rational::operator=(float a)
328{
329  mpz_t
330    h1,
331    h2;
332  int
333    i=0;
334  double
335    aa=a;
336
337  disconnect();
338  mpz_init_set_ui(h1,1);
339  while(fabs(10.*aa) < FLT_MAX && i<FLT_DIG){
340    aa*=10.;
341    mpz_mul_ui(h1,h1,10);
342    i++;
343  }
344  mpz_init_set_d(h2,aa);
345  mpq_set_num(p->rat,h2);
346  mpq_set_den(p->rat,h1);
347  mpq_canonicalize(p->rat);
348  mpz_clear(h1);
349  mpz_clear(h2);
350  return *this;
351}
352
353Rational& Rational::operator=(char *s)
354{
355  mpz_t
356    h1;
357  int
358    exp=0,
359    i=0;
360  char
361    *s1=s,
362    *s2,
363    *ss;
364
365  ss=new char[strlen(s)+1];
366  strcpy(ss,s);
367  s1=ss;
368  disconnect();
369  if (isdigit(s1[0]) || s1[0]=='-' || s1[0]=='+'){
370    if (s1[0]=='+') ++s1;
371    if (s1[0]=='-') ++i;
372    while(isdigit(s1[i]))
373      ++i;
374    if (s1[i]=='.'){
375      ++i;
376      while(isdigit(s1[i])){
377        s1[i-1]=s1[i];
378        ++i;
379        --exp;
380      }
381      s1[i-1]='\0';
382    }
383    if (s1[i]=='e' || s1[i]=='E'){
384      s2=s1+i+1;
385    }
386    else
387      s2="";
388    s1[i]='\0';
389    i=exp+atoi(s2);
390    mpz_init_set_str(h1,s1,10);
391    delete[] ss;
392    mpq_set_z(p->rat,h1);
393    mpq_set_ui(save.p->rat,10,1);
394    if (i>0){
395      for(int j=0;j<i;j++)
396        mpq_mul(p->rat,p->rat,save.p->rat);
397    }
398    else if (i<0){
399      for(int j=0;j>i;j--)
400        mpq_div(p->rat,p->rat,save.p->rat);
401    }
402    mpq_canonicalize(p->rat);
403    mpz_clear(h1);
404  }
405  else
406    mpq_set_ui(p->rat,0,1);
407  return *this;
408}
409
410Rational& Rational::operator=(const Rational& a)
411{
412  a.p->n++;
413  if (--p->n==0){
414    mpq_clear(p->rat);
415    delete p;
416  }
417  p=a.p;
418  return *this;
419}
420
421// ----------------------------------------------------------------------------
422//  Numerator and denominator
423// ----------------------------------------------------------------------------
424
425Rational Rational::get_num( )
426{
427    Rational erg;
428
429    mpq_set_num( erg.p->rat,mpq_numref( p->rat ) );
430
431    return  erg;
432}
433
434int Rational::get_num_si( )
435{
436    return  mpz_get_si( mpq_numref( p->rat ) );
437}
438
439Rational Rational::get_den( )
440{
441    Rational erg;
442
443    mpq_set_num( erg.p->rat,mpq_denref( p->rat ) );
444
445    return  erg;
446}
447
448int Rational::get_den_si( )
449{
450    return  mpz_get_si( mpq_denref( p->rat ) );
451}
452
453// ----------------------------------------------------------------------------
454//  Casting
455// ----------------------------------------------------------------------------
456
457Rational::operator bool()
458{
459    if (mpq_sgn(p->rat)) return true;
460    return false;
461}
462
463Rational::operator long int()
464{
465  mpz_t
466    h;
467  long int
468    ret_val;
469
470  mpz_init(h);
471  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
472  ret_val=mpz_get_si(h);
473  mpz_clear(h);
474
475  return ret_val;
476}
477
478Rational::operator unsigned long int()
479{
480  mpz_t
481    h;
482  unsigned long int
483    ret_val;
484
485  mpz_init(h);
486  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
487  ret_val=mpz_get_ui(h);
488  mpz_clear(h);
489  return ret_val;
490}
491
492Rational::operator int()
493{
494  mpz_t
495    h;
496  long int
497    ret_val;
498
499  mpz_init(h);
500  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
501  ret_val=mpz_get_si(h);
502  mpz_clear(h);
503
504  return ret_val;
505}
506
507Rational::operator unsigned int()
508{
509  mpz_t
510    h;
511  unsigned long int
512    ret_val;
513
514  mpz_init(h);
515  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
516  ret_val=mpz_get_ui(h);
517  mpz_clear(h);
518  return ret_val;
519}
520
521Rational::operator short int()
522{
523  mpz_t
524    h;
525  short int
526    ret_val;
527
528  mpz_init(h);
529  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
530  ret_val=mpz_get_si(h);
531  mpz_clear(h);
532  return ret_val;
533}
534
535Rational::operator unsigned short int()
536{
537  mpz_t
538    h;
539  unsigned short int
540    ret_val;
541
542  mpz_init(h);
543  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
544  ret_val=mpz_get_ui(h);
545  mpz_clear(h);
546  return ret_val;
547}
548
549Rational::operator char()
550{
551  mpz_t
552    h;
553  char
554    ret_val;
555
556  mpz_init(h);
557  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
558  ret_val=mpz_get_si(h);
559  mpz_clear(h);
560  return ret_val;
561}
562
563Rational::operator unsigned char()
564{
565  mpz_t
566    h;
567  unsigned char
568    ret_val;
569
570  mpz_init(h);
571  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
572  ret_val=mpz_get_ui(h);
573  mpz_clear(h);
574  return ret_val;
575}
576
577Rational::operator double()
578{
579  return mpq_get_d(p->rat);
580}
581
582Rational::operator float()
583{
584  return mpq_get_d(p->rat);
585}
586
587// ----------------------------------------------------------------------------
588//  Unary minus
589// ----------------------------------------------------------------------------
590
591Rational
592Rational::operator-()
593{
594  Rational
595    erg;
596
597  mpq_neg(erg.p->rat,p->rat);
598  return erg;
599}
600
601Rational operator - ( const Rational &r )
602{
603  Rational
604    erg;
605
606  mpq_neg(erg.p->rat,r.p->rat);
607  return erg;
608}
609
610// ----------------------------------------------------------------------------
611//  Inverse
612// ----------------------------------------------------------------------------
613
614Rational
615Rational::operator~()
616{
617  Rational
618    erg;
619
620  mpq_inv(erg.p->rat,p->rat);
621  return erg;
622}
623
624// ----------------------------------------------------------------------------
625//  +=, -= ...
626// ----------------------------------------------------------------------------
627
628Rational&
629Rational::operator+=(const Rational &a)
630{
631  mpq_set(save.p->rat,p->rat);
632  disconnect();
633  mpq_add(p->rat,save.p->rat,a.p->rat);
634  return *this;
635}
636
637Rational&
638Rational::operator-=(const Rational &a)
639{
640  mpq_set(save.p->rat,p->rat);
641  disconnect();
642  mpq_sub(p->rat,save.p->rat,a.p->rat);
643  return *this;
644}
645
646Rational&
647Rational::operator*=(const Rational &a)
648{
649  mpq_set(save.p->rat,p->rat);
650  disconnect();
651  mpq_mul(p->rat,save.p->rat,a.p->rat);
652  return *this;
653}
654
655Rational&
656Rational::operator/=(const Rational &a)
657{
658  mpq_set(save.p->rat,p->rat);
659  disconnect();
660  mpq_div(p->rat,save.p->rat,a.p->rat);
661  return *this;
662}
663
664// ----------------------------------------------------------------------------
665//  Increment and decrement
666// ----------------------------------------------------------------------------
667
668Rational&
669Rational::operator++()
670{
671  mpq_set(save.p->rat,p->rat);
672  *this=1;
673  mpq_add(p->rat,p->rat,save.p->rat);
674  return *this;
675}
676
677Rational
678Rational::operator++(int)
679{
680  Rational
681    erg(*this);
682
683  mpq_set(save.p->rat,p->rat);
684  *this=1;
685  mpq_add(p->rat,p->rat,save.p->rat);
686  return erg;
687}
688
689Rational&
690Rational::operator--()
691{
692  mpq_set(save.p->rat,p->rat);
693  *this=1;
694  mpq_sub(p->rat,save.p->rat,p->rat);
695  return *this;
696}
697
698Rational
699Rational::operator--(int)
700{
701  Rational
702    erg(*this);
703
704  mpq_set(save.p->rat,p->rat);
705  *this=1;
706  mpq_sub(p->rat,save.p->rat,p->rat);
707  return erg;
708}
709
710// ----------------------------------------------------------------------------
711//  Relational operators
712// ----------------------------------------------------------------------------
713
714bool operator<(const Rational& a,const Rational& b)
715{
716  if (mpq_cmp(a.p->rat,b.p->rat)<0) return true;
717  return false;
718}
719
720bool operator<=(const Rational& a,const Rational& b)
721{
722  if (mpq_cmp(a.p->rat,b.p->rat)>0) return false;
723  return true;
724}
725
726bool operator>(const Rational& a,const Rational& b)
727{
728  if (mpq_cmp(a.p->rat,b.p->rat)>0) return true;
729  return false;
730}
731
732bool operator>=(const Rational& a,const Rational& b)
733{
734  if (mpq_cmp(a.p->rat,b.p->rat)<0) return false;
735  return true;
736}
737
738bool operator==(const Rational& a,const Rational& b)
739{
740  if (mpq_equal(a.p->rat,b.p->rat)) return true;
741  return false;
742}
743
744bool operator!=(const Rational& a,const Rational& b)
745{
746  if (mpq_equal(a.p->rat,b.p->rat)) return false;
747  return true;
748}
749
750// ----------------------------------------------------------------------------
751//  Ostream
752// ----------------------------------------------------------------------------
753
754#ifdef GMPRAT_PRINT
755ostream& operator<< (ostream& s,const Rational& a)
756{
757    char *snum,*sdenom;
758
759    snum   = mpz_get_str( NULL,10,mpq_numref(a.p->rat) );
760    sdenom = mpz_get_str( NULL,10,mpq_denref(a.p->rat) );
761
762    if( sdenom[0] == '1' && sdenom[1] == '\0' )
763    {
764        #ifdef GMPRAT_IOSTREAM
765            s << snum;
766        #else
767            fprintf( stdout,snum );
768        #endif
769    }
770    else
771    {
772        #ifdef GMPRAT_IOSTREAM
773            s << snum << "/" << sdenom;
774        #else
775            fprintf( stdout,snum );
776            fprintf( stdout,"/" );
777            fprintf( stdout,sdenom );
778        #endif
779    }
780
781    //free( snum );
782    //free( sdenom );
783
784    return s;
785}
786#endif
787
788unsigned int Rational::length( ) const
789{
790    char *snum = (char*)NULL;
791    char *sden = (char*)NULL;
792
793    snum = mpz_get_str( snum,10,mpq_numref( p->rat ) );
794    sden = mpz_get_str( sden,10,mpq_denref( p->rat ) );
795
796    int length = strlen( snum );
797
798    if( sden[0] != '1' || sden[1] != '\0' ) length += strlen( sden ) + 1;
799
800    free( snum );
801    free( sden );
802
803    return  length;
804}
805
806// ----------------------------------------------------------------------------
807//  Operators
808// ----------------------------------------------------------------------------
809
810Rational
811operator+(const Rational& a,const Rational &b)
812{
813  Rational
814    erg(a);
815
816  return erg+=b;
817}
818
819Rational
820operator-(const Rational& a,const Rational &b)
821{
822  Rational
823    erg(a);
824
825  return erg-=b;
826}
827
828Rational
829operator*(const Rational& a,const Rational &b)
830{
831  Rational
832    erg(a);
833
834  return erg*=b;
835}
836
837Rational pow( const Rational& a,int e )
838{
839    Rational erg(1);
840
841    for( int i=0; i<e; i++ )
842    {
843        erg *= a;
844    }
845    return erg;
846}
847
848Rational operator/(const Rational& a,const Rational &b)
849{
850  Rational
851    erg(a);
852
853  return erg/=b;
854}
855
856int sgn(const Rational& a)
857{
858  return mpq_sgn(a.p->rat);
859}
860
861Rational
862abs(const Rational& a)
863{
864  Rational
865    erg;
866
867  if (mpq_sgn(a.p->rat)<0)
868    mpq_neg(erg.p->rat,a.p->rat);
869  else
870    mpq_set(erg.p->rat,a.p->rat);
871  return erg;
872}
873
874Rational gcd( const Rational &a,const Rational &b )
875{
876    if( a == 0 )
877    {
878        if( b == 0 )
879        {
880            return  (Rational)1;
881        }
882        else
883        {
884            return  abs( b );
885        }
886    }
887    else if( b == 0 )
888    {
889        return  abs( a );
890    }
891
892    Rational erg;
893
894    mpz_gcd( mpq_numref( erg.p->rat ),
895            mpq_numref( a.p->rat ),mpq_numref( b.p->rat ) );
896    mpz_gcd( mpq_denref( erg.p->rat ),
897            mpq_denref( a.p->rat ),mpq_denref( b.p->rat ) );
898
899    //mpq_canonicalize( erg.p->rat );
900
901    return  abs( erg );
902}
903
904Rational gcd( Rational *a,int n )
905{
906    if( n == 1 )
907    {
908        return  a[0];
909    }
910
911    Rational g = gcd( a[0],a[1] );
912
913    for( int i=2; i<n; i++ )
914    {
915        g = gcd( g,a[i] );
916    }
917
918    return  g;
919}
920
921Rational lcm( const Rational &a,const Rational &b )
922{
923    if( a == 0 )
924    {
925        return b;
926    }
927    else if( b == 0 )
928    {
929        return a;
930    }
931
932    return a*b/gcd(a,b);
933}
934
935Rational lcm( Rational *a,int n )
936{
937    if( n == 1 )
938    {
939        return  a[0];
940    }
941
942    Rational g = lcm( a[0],a[1] );
943
944    for( int i=2; i<n; i++ )
945    {
946        g = lcm( g,a[i] );
947    }
948
949    return  g;
950}
951
952double  Rational::complexity( ) const
953{
954    double num = mpz_get_d( mpq_numref( p->rat ) );
955    double den = mpz_get_d( mpq_denref( p->rat ) );
956
957    if( num < 0 ) num = -num;
958    if( den < 0 ) den = -den;
959
960    return  ( num > den ? num : den );
961}
962
963#endif /* HAVE_SPECTRUM */
964// ----------------------------------------------------------------------------
965//  GMPrat.cc
966//  end of file
967// ----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.