source: git/Singular/GMPrat.cc @ 545bf8

spielwiese
Last change on this file since 545bf8 was 545bf8, checked in by Hans Schönemann <hannes@…>, 23 years ago
*hannes: GPMrat- global destruktor git-svn-id: file:///usr/local/Singular/svn/trunk@5479 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 22.5 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// ----------------------------------------------------------------------------
33extern int mmInit();
34static int GMPIsInitialized=mmInit();
35
36Rational Rational::save;    // dummy variable
37Rational RAT_EPS("1e-15");  // epsilon used in transcendental functions
38
39// ----------------------------------------------------------------------------
40//  the number Pi up to 25, 50, 75 and 100 digits
41// ----------------------------------------------------------------------------
42
43//Rational RAT_PI_25("3.141592653589793238462643");
44//Rational RAT_PI_50("3.1415926535897932384626433832795028841971693993751");
45//Rational RAT_PI_75("3.141592653589793238462643383279502884197169399375105820974944592307816406286");
46//Rational RAT_PI_100("3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068");
47//Rational RAT_PI(RAT_PI_25);
48
49
50// ----------------------------------------------------------------------------
51//  disconnect a rational from its reference
52// ----------------------------------------------------------------------------
53
54void    Rational::disconnect( )
55{
56    if( p->n>1)
57    {
58        p->n--;
59        p = new rep;
60    }
61    else
62    {
63        mpq_clear(p->rat);
64    }
65    mpq_init(p->rat);
66}
67
68// ----------------------------------------------------------------------------
69//  Constructors
70// ----------------------------------------------------------------------------
71
72Rational::Rational( )
73{
74    p = new rep;
75    mpq_init( p->rat );
76}
77
78Rational::Rational( long int a )
79{
80    p = new rep;
81    mpq_init( p->rat );
82    mpq_set_si( p->rat,a,1 );
83}
84
85Rational::Rational( unsigned long int a )
86{
87    p = new rep;
88    mpq_init( p->rat );
89    mpq_set_ui( p->rat,a,1 );
90}
91
92Rational::Rational( int a )
93{
94    p = new rep;
95    mpq_init( p->rat );
96    mpq_set_si( p->rat,(long int)a,1 );
97}
98
99Rational::Rational( unsigned int a )
100{
101    p = new rep;
102    mpq_init( p->rat );
103    mpq_set_ui( p->rat,(unsigned long)a,1 );
104}
105
106Rational::Rational( short int a )
107{
108    p = new rep;
109    mpq_init( p->rat );
110    mpq_set_si( p->rat,(long int)a,1 );
111}
112
113Rational::Rational( unsigned short int a )
114{
115    p = new rep;
116    mpq_init( p->rat );
117    mpq_set_ui( p->rat,(unsigned long int)a,1 );
118}
119
120Rational::Rational( char a )
121{
122    p = new rep;
123    mpq_init( p->rat );
124    mpq_set_si( p->rat,(long int)a,1 );
125}
126
127Rational::Rational( unsigned char a )
128{
129    p = new rep;
130    mpq_init( p->rat );
131    mpq_set_ui( p->rat,(unsigned long int)a,1 );
132}
133
134Rational::Rational( const Rational& a )
135{
136    a.p->n++;
137    p=a.p;
138}
139
140Rational::Rational( double a )
141{
142    mpz_t   h1,h2;
143    int     i=0;
144    double  aa=a;
145
146    p = new rep;
147    mpq_init( p->rat );
148    mpz_init_set_ui( h1,1 );
149
150    while( fabs( 10.0*aa ) < DBL_MAX && i<DBL_DIG )
151    {
152        aa *= 10.;
153        mpz_mul_ui( h1,h1,10 );
154        i++;
155    }
156    mpz_init_set_d( h2,aa );
157    mpq_set_num( p->rat,h2 );
158    mpq_set_den( p->rat,h1 );
159    mpq_canonicalize( p->rat );
160    mpz_clear( h1 );
161    mpz_clear( h2 );
162}
163
164Rational::Rational(float a)
165{
166    mpz_t h1,h2;
167    int i=0;
168    double aa=a;
169
170    p=new rep;
171    mpq_init(p->rat);
172    mpz_init_set_ui(h1,1);
173    while(fabs(10.*aa) < FLT_MAX && i<FLT_DIG){
174        aa*=10.;
175        mpz_mul_ui(h1,h1,10);
176        i++;
177    }
178    mpz_init_set_d(h2,aa);
179    mpq_set_num(p->rat,h2);
180    mpq_set_den(p->rat,h1);
181    mpq_canonicalize(p->rat);
182    mpz_clear(h1);
183    mpz_clear(h2);
184}
185
186// ----------------------------------------------------------------------------
187//  Create a Rational from a string like "+1234.5678e-1234"
188// ----------------------------------------------------------------------------
189
190Rational::Rational( char *s )
191{
192    mpz_t   h1;
193    int     exp=0,i=0;
194    char    *s1,*s2,*ss;
195
196    ss = new char[strlen(s)+1];
197    strcpy( ss,s );
198    s1 = ss;
199    p = new rep;
200    mpq_init( p->rat );
201    if( isdigit(s1[0]) || s1[0]=='-' || s1[0]=='+' )
202    {
203        if (s1[0]=='+') ++s1;
204        if (s1[0]=='-') ++i;
205
206        while( isdigit( s1[i] ) )
207        {
208            ++i;
209        }
210        if (s1[i]=='.')
211        {
212            ++i;
213            while( isdigit( s1[i] ) )
214            {
215                s1[i-1]=s1[i];
216                ++i;
217                --exp;
218            }
219            s1[i-1]='\0';
220        }
221        if (s1[i]=='e' || s1[i]=='E')
222        {
223            s2=s1+i+1;
224        }
225        else
226            s2="";
227
228        s1[i]='\0';
229        i=exp+atoi(s2);
230        mpz_init_set_str(h1,s1,10);
231        delete[] ss;
232        mpq_set_z(p->rat,h1);
233        mpq_set_ui(save.p->rat,10,1);
234        if (i>0)
235        {
236            for(int j=0;j<i;j++)
237                mpq_mul(p->rat,p->rat,save.p->rat);
238        }
239        else if (i<0)
240        {
241            for(int j=0;j>i;j--)
242                mpq_div(p->rat,p->rat,save.p->rat);
243        }
244        mpq_canonicalize(p->rat);
245        mpz_clear(h1);
246    }
247}
248
249// ----------------------------------------------------------------------------
250//  Constructors with two arguments: numerator and denominator
251// ----------------------------------------------------------------------------
252
253Rational::Rational(long int a, unsigned long int b)
254{
255    p = new rep;
256    mpq_init( p->rat );
257    mpq_set_si( p->rat,a,b );
258    mpq_canonicalize( p->rat );
259}
260
261Rational::Rational(unsigned long int a, unsigned long int b)
262{
263    p=new rep;
264    mpq_init(p->rat);
265    mpq_set_ui(p->rat,a,b);
266    mpq_canonicalize(p->rat);
267}
268
269Rational::Rational(int a, unsigned int b)
270{
271    p=new rep;
272    mpq_init(p->rat);
273    mpq_set_si(p->rat,(long int) a,(unsigned long int) b);
274    mpq_canonicalize(p->rat);
275}
276
277Rational::Rational(unsigned int a, unsigned int b)
278{
279    p=new rep;
280    mpq_init(p->rat);
281    mpq_set_ui(p->rat,(unsigned long) a,(unsigned long int) b);
282    mpq_canonicalize(p->rat);
283}
284
285Rational::Rational(short int a, unsigned short int b)
286{
287    p=new rep;
288    mpq_init(p->rat);
289    mpq_set_si(p->rat,(long int) a,(unsigned long int) b);
290    mpq_canonicalize(p->rat);
291}
292
293Rational::Rational(unsigned short int a, unsigned short int b)
294{
295    p=new rep;
296    mpq_init(p->rat);
297    mpq_set_ui(p->rat,(unsigned long int) a,(unsigned long int) b);
298    mpq_canonicalize(p->rat);
299}
300
301Rational::Rational(char a, unsigned char b)
302{
303    p=new rep;
304    mpq_init(p->rat);
305    mpq_set_si(p->rat,(long int) a,(unsigned long int) b);
306    mpq_canonicalize(p->rat);
307}
308
309Rational::Rational(unsigned char a, unsigned char b)
310{
311    p=new rep;
312    mpq_init(p->rat);
313    mpq_set_ui(p->rat,(unsigned long int) a,(unsigned long int) b);
314    mpq_canonicalize(p->rat);
315}
316
317Rational::Rational(const Rational& a, const Rational& b)
318{
319    p=new rep;
320    mpq_init(p->rat);
321    mpq_div(p->rat, a.p->rat, b.p->rat);
322}
323
324Rational::Rational(long int a, long int b)
325{
326    if (b<0) a=-a;
327    p=new rep;
328    mpq_init(p->rat);
329    mpq_set_si(p->rat,a,abs(b));
330    mpq_canonicalize(p->rat);
331}
332
333Rational::Rational(int a, int b)
334{
335    if (b<0) a=-a;
336    p=new rep;
337    mpq_init(p->rat);
338    mpq_set_si(p->rat,(long int) a,(unsigned long int) abs(b));
339    mpq_canonicalize(p->rat);
340}
341
342Rational::Rational(short int a, short int b)
343{
344    if (b<0) a=-a;
345    p=new rep;
346    mpq_init(p->rat);
347    mpq_set_si(p->rat,(long int) a,(unsigned long int) abs(b));
348    mpq_canonicalize(p->rat);
349}
350
351Rational::Rational(char a, char b)
352{
353    if (b<0) a=-a;
354    p=new rep;
355    mpq_init(p->rat);
356    mpq_set_si(p->rat,(long int) a,(unsigned long int) abs(b));
357    mpq_canonicalize(p->rat);
358}
359
360Rational::Rational(char *sn, char *sd)
361{
362  Rational
363    h=sd;
364
365  p=new rep;
366  mpq_init(p->rat);
367  *this=sn;
368  mpq_div(p->rat,p->rat,h.p->rat);
369}
370
371// ----------------------------------------------------------------------------
372//  Destructor
373// ----------------------------------------------------------------------------
374
375Rational::~Rational()
376{
377  if (--p->n==0){
378    mpq_clear(p->rat);
379    delete p;
380  }
381}
382
383// ----------------------------------------------------------------------------
384//  Assignment operators
385// ----------------------------------------------------------------------------
386
387Rational& Rational::operator=(long int a)
388{
389  disconnect();
390  mpq_set_si(p->rat,a,1);
391  return *this;
392}
393
394Rational& Rational::operator=(unsigned long int a)
395{
396  disconnect();
397  mpq_set_ui(p->rat,a,1);
398  return *this;
399}
400
401Rational& Rational::operator=(int a)
402{
403  disconnect();
404  mpq_set_si(p->rat,(long int) a,1);
405  return *this;
406}
407
408Rational& Rational::operator=(unsigned int a)
409{
410  disconnect();
411  mpq_set_ui(p->rat,(unsigned long int) a,1);
412  return *this;
413}
414
415Rational& Rational::operator=(short int a)
416{
417  disconnect();
418  mpq_set_si(p->rat,(long int) a,1);
419  return *this;
420}
421
422Rational& Rational::operator=(unsigned short int a)
423{
424  disconnect();
425  mpq_set_ui(p->rat,(unsigned long int) a,1);
426  return *this;
427}
428
429Rational& Rational::operator=(char a)
430{
431  disconnect();
432  mpq_set_si(p->rat,(long int) a,1);
433  return *this;
434}
435
436Rational& Rational::operator=(unsigned char a)
437{
438  disconnect();
439  mpq_set_ui(p->rat,(unsigned long int) a,1);
440  return *this;
441}
442
443Rational& Rational::operator=(double a)
444{
445  mpz_t
446    h1,
447    h2;
448  int
449    i=0;
450  double
451    aa=a;
452
453  disconnect();
454  mpz_init_set_ui(h1,1);
455  while(fabs(10.*aa) < DBL_MAX && i<DBL_DIG){
456    aa*=10.;
457    mpz_mul_ui(h1,h1,10);
458    i++;
459  }
460  mpz_init_set_d(h2,aa);
461  mpq_set_num(p->rat,h2);
462  mpq_set_den(p->rat,h1);
463  mpq_canonicalize(p->rat);
464  mpz_clear(h1);
465  mpz_clear(h2);
466  return *this;
467}
468
469Rational& Rational::operator=(float a)
470{
471  mpz_t
472    h1,
473    h2;
474  int
475    i=0;
476  double
477    aa=a;
478
479  disconnect();
480  mpz_init_set_ui(h1,1);
481  while(fabs(10.*aa) < FLT_MAX && i<FLT_DIG){
482    aa*=10.;
483    mpz_mul_ui(h1,h1,10);
484    i++;
485  }
486  mpz_init_set_d(h2,aa);
487  mpq_set_num(p->rat,h2);
488  mpq_set_den(p->rat,h1);
489  mpq_canonicalize(p->rat);
490  mpz_clear(h1);
491  mpz_clear(h2);
492  return *this;
493}
494
495Rational& Rational::operator=(char *s)
496{
497  mpz_t
498    h1;
499  int
500    exp=0,
501    i=0;
502  char
503    *s1=s,
504    *s2,
505    *ss;
506
507  ss=new char[strlen(s)+1];
508  strcpy(ss,s);
509  s1=ss;
510  disconnect();
511  if (isdigit(s1[0]) || s1[0]=='-' || s1[0]=='+'){
512    if (s1[0]=='+') ++s1;
513    if (s1[0]=='-') ++i;
514    while(isdigit(s1[i]))
515      ++i;
516    if (s1[i]=='.'){
517      ++i;
518      while(isdigit(s1[i])){
519        s1[i-1]=s1[i];
520        ++i;
521        --exp;
522      }
523      s1[i-1]='\0';
524    }
525    if (s1[i]=='e' || s1[i]=='E'){
526      s2=s1+i+1;
527    }
528    else
529      s2="";
530    s1[i]='\0';
531    i=exp+atoi(s2);
532    mpz_init_set_str(h1,s1,10);
533    delete[] ss;
534    mpq_set_z(p->rat,h1);
535    mpq_set_ui(save.p->rat,10,1);
536    if (i>0){
537      for(int j=0;j<i;j++)
538        mpq_mul(p->rat,p->rat,save.p->rat);
539    }
540    else if (i<0){
541      for(int j=0;j>i;j--)
542        mpq_div(p->rat,p->rat,save.p->rat);
543    }
544    mpq_canonicalize(p->rat);
545    mpz_clear(h1);
546  }
547  else
548    mpq_set_ui(p->rat,0,1);
549  return *this;
550}
551
552Rational& Rational::operator=(const Rational& a)
553{
554  a.p->n++;
555  if (--p->n==0){
556    mpq_clear(p->rat);
557    delete p;
558  }
559  p=a.p;
560  return *this;
561}
562
563// ----------------------------------------------------------------------------
564//  Numerator and denominator
565// ----------------------------------------------------------------------------
566
567Rational Rational::get_num( )
568{
569    Rational erg;
570
571    mpq_set_num( erg.p->rat,mpq_numref( p->rat ) );
572
573    return  erg;
574}
575
576int Rational::get_num_si( )
577{
578    return  mpz_get_si( mpq_numref( p->rat ) );
579}
580
581Rational Rational::get_den( )
582{
583    Rational erg;
584
585    mpq_set_num( erg.p->rat,mpq_denref( p->rat ) );
586
587    return  erg;
588}
589
590int Rational::get_den_si( )
591{
592    return  mpz_get_si( mpq_denref( p->rat ) );
593}
594
595// ----------------------------------------------------------------------------
596//  Casting
597// ----------------------------------------------------------------------------
598
599Rational::operator bool()
600{
601    if (mpq_sgn(p->rat)) return true;
602    return false;
603}
604
605Rational::operator long int()
606{
607  mpz_t
608    h;
609  long int
610    ret_val;
611
612  mpz_init(h);
613  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
614  ret_val=mpz_get_si(h);
615  mpz_clear(h);
616
617  return ret_val;
618}
619
620Rational::operator unsigned long int()
621{
622  mpz_t
623    h;
624  unsigned long int
625    ret_val;
626
627  mpz_init(h);
628  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
629  ret_val=mpz_get_ui(h);
630  mpz_clear(h);
631  return ret_val;
632}
633
634Rational::operator int()
635{
636  mpz_t
637    h;
638  long int
639    ret_val;
640
641  mpz_init(h);
642  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
643  ret_val=mpz_get_si(h);
644  mpz_clear(h);
645
646  return ret_val;
647}
648
649Rational::operator unsigned int()
650{
651  mpz_t
652    h;
653  unsigned long int
654    ret_val;
655
656  mpz_init(h);
657  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
658  ret_val=mpz_get_ui(h);
659  mpz_clear(h);
660  return ret_val;
661}
662
663Rational::operator short int()
664{
665  mpz_t
666    h;
667  short int
668    ret_val;
669
670  mpz_init(h);
671  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
672  ret_val=mpz_get_si(h);
673  mpz_clear(h);
674  return ret_val;
675}
676
677Rational::operator unsigned short int()
678{
679  mpz_t
680    h;
681  unsigned short int
682    ret_val;
683
684  mpz_init(h);
685  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
686  ret_val=mpz_get_ui(h);
687  mpz_clear(h);
688  return ret_val;
689}
690
691Rational::operator char()
692{
693  mpz_t
694    h;
695  char
696    ret_val;
697
698  mpz_init(h);
699  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
700  ret_val=mpz_get_si(h);
701  mpz_clear(h);
702  return ret_val;
703}
704
705Rational::operator unsigned char()
706{
707  mpz_t
708    h;
709  unsigned char
710    ret_val;
711
712  mpz_init(h);
713  mpz_tdiv_q(h,mpq_numref(p->rat),mpq_denref(p->rat));
714  ret_val=mpz_get_ui(h);
715  mpz_clear(h);
716  return ret_val;
717}
718
719Rational::operator double()
720{
721  return mpq_get_d(p->rat);
722}
723
724Rational::operator float()
725{
726  return mpq_get_d(p->rat);
727}
728
729// ----------------------------------------------------------------------------
730//  Unary minus
731// ----------------------------------------------------------------------------
732
733Rational
734Rational::operator-()
735{
736  Rational
737    erg;
738
739  mpq_neg(erg.p->rat,p->rat);
740  return erg;
741}
742
743Rational operator - ( const Rational &r )
744{
745  Rational
746    erg;
747
748  mpq_neg(erg.p->rat,r.p->rat);
749  return erg;
750}
751
752// ----------------------------------------------------------------------------
753//  Inverse
754// ----------------------------------------------------------------------------
755
756Rational
757Rational::operator~()
758{
759  Rational
760    erg;
761
762  mpq_inv(erg.p->rat,p->rat);
763  return erg;
764}
765
766// ----------------------------------------------------------------------------
767//  +=, -= ...
768// ----------------------------------------------------------------------------
769
770Rational&
771Rational::operator+=(const Rational &a)
772{
773  mpq_set(save.p->rat,p->rat);
774  disconnect();
775  mpq_add(p->rat,save.p->rat,a.p->rat);
776  return *this;
777}
778
779Rational&
780Rational::operator-=(const Rational &a)
781{
782  mpq_set(save.p->rat,p->rat);
783  disconnect();
784  mpq_sub(p->rat,save.p->rat,a.p->rat);
785  return *this;
786}
787
788Rational&
789Rational::operator*=(const Rational &a)
790{
791  mpq_set(save.p->rat,p->rat);
792  disconnect();
793  mpq_mul(p->rat,save.p->rat,a.p->rat);
794  return *this;
795}
796
797Rational&
798Rational::operator/=(const Rational &a)
799{
800  mpq_set(save.p->rat,p->rat);
801  disconnect();
802  mpq_div(p->rat,save.p->rat,a.p->rat);
803  return *this;
804}
805
806// ----------------------------------------------------------------------------
807//  Increment and decrement
808// ----------------------------------------------------------------------------
809
810Rational&
811Rational::operator++()
812{
813  mpq_set(save.p->rat,p->rat);
814  *this=1;
815  mpq_add(p->rat,p->rat,save.p->rat);
816  return *this;
817}
818
819Rational
820Rational::operator++(int)
821{
822  Rational
823    erg(*this);
824
825  mpq_set(save.p->rat,p->rat);
826  *this=1;
827  mpq_add(p->rat,p->rat,save.p->rat);
828  return erg;
829}
830
831Rational&
832Rational::operator--()
833{
834  mpq_set(save.p->rat,p->rat);
835  *this=1;
836  mpq_sub(p->rat,save.p->rat,p->rat);
837  return *this;
838}
839
840Rational
841Rational::operator--(int)
842{
843  Rational
844    erg(*this);
845
846  mpq_set(save.p->rat,p->rat);
847  *this=1;
848  mpq_sub(p->rat,save.p->rat,p->rat);
849  return erg;
850}
851
852// ----------------------------------------------------------------------------
853//  Relational operators
854// ----------------------------------------------------------------------------
855
856bool operator<(const Rational& a,const Rational& b)
857{
858  if (mpq_cmp(a.p->rat,b.p->rat)<0) return true;
859  return false;
860}
861
862bool operator<=(const Rational& a,const Rational& b)
863{
864  if (mpq_cmp(a.p->rat,b.p->rat)>0) return false;
865  return true;
866}
867
868bool operator>(const Rational& a,const Rational& b)
869{
870  if (mpq_cmp(a.p->rat,b.p->rat)>0) return true;
871  return false;
872}
873
874bool operator>=(const Rational& a,const Rational& b)
875{
876  if (mpq_cmp(a.p->rat,b.p->rat)<0) return false;
877  return true;
878}
879
880bool operator==(const Rational& a,const Rational& b)
881{
882  if (mpq_equal(a.p->rat,b.p->rat)) return true;
883  return false;
884}
885
886bool operator!=(const Rational& a,const Rational& b)
887{
888  if (mpq_equal(a.p->rat,b.p->rat)) return false;
889  return true;
890}
891
892// ----------------------------------------------------------------------------
893//  Ostream
894// ----------------------------------------------------------------------------
895
896#ifdef GMPRAT_PRINT
897ostream& operator<< (ostream& s,const Rational& a)
898{
899    char *snum,*sdenom;
900
901    snum   = mpz_get_str( NULL,10,mpq_numref(a.p->rat) );
902    sdenom = mpz_get_str( NULL,10,mpq_denref(a.p->rat) );
903
904    if( sdenom[0] == '1' && sdenom[1] == '\0' )
905    {
906        #ifdef GMPRAT_IOSTREAM
907            s << snum;
908        #else
909            fprintf( stdout,snum );
910        #endif
911    }
912    else
913    {
914        #ifdef GMPRAT_IOSTREAM
915            s << snum << "/" << sdenom;
916        #else
917            fprintf( stdout,snum );
918            fprintf( stdout,"/" );
919            fprintf( stdout,sdenom );
920        #endif
921    }
922
923    //free( snum );
924    //free( sdenom );
925
926    return s;
927}
928#endif
929
930unsigned int Rational::length( ) const
931{
932    char *snum = (char*)NULL;
933    char *sden = (char*)NULL;
934
935    snum = mpz_get_str( snum,10,mpq_numref( p->rat ) );
936    sden = mpz_get_str( sden,10,mpq_denref( p->rat ) );
937
938    int length = strlen( snum );
939
940    if( sden[0] != '1' || sden[1] != '\0' ) length += strlen( sden ) + 1;
941
942    free( snum );
943    free( sden );
944
945    return  length;
946}
947
948// ----------------------------------------------------------------------------
949//  Operators
950// ----------------------------------------------------------------------------
951
952Rational
953operator+(const Rational& a,const Rational &b)
954{
955  Rational
956    erg(a);
957
958  return erg+=b;
959}
960
961Rational
962operator-(const Rational& a,const Rational &b)
963{
964  Rational
965    erg(a);
966
967  return erg-=b;
968}
969
970Rational
971operator*(const Rational& a,const Rational &b)
972{
973  Rational
974    erg(a);
975
976  return erg*=b;
977}
978
979Rational pow( const Rational& a,int e )
980{
981    Rational erg(1);
982
983    for( int i=0; i<e; i++ )
984    {
985        erg *= a;
986    }
987    return erg;
988}
989
990Rational operator/(const Rational& a,const Rational &b)
991{
992  Rational
993    erg(a);
994
995  return erg/=b;
996}
997
998int sgn(const Rational& a)
999{
1000  return mpq_sgn(a.p->rat);
1001}
1002
1003Rational
1004abs(const Rational& a)
1005{
1006  Rational
1007    erg;
1008
1009  if (mpq_sgn(a.p->rat)<0)
1010    mpq_neg(erg.p->rat,a.p->rat);
1011  else
1012    mpq_set(erg.p->rat,a.p->rat);
1013  return erg;
1014}
1015
1016Rational
1017sqrt(const Rational& a)
1018{
1019  Rational
1020    erg;
1021  mpz_t
1022    h1,
1023    h2;
1024
1025  mpz_init_set(h1,mpq_numref(a.p->rat));
1026  mpz_init(h2);
1027  mpz_set_ui(h2,1);
1028  mpz_mul(h1,h1,mpq_denref(a.p->rat));
1029  mpz_mul(h1,h1,mpq_denref(RAT_EPS.p->rat));
1030  mpz_mul(h1,h1,mpq_denref(RAT_EPS.p->rat));
1031  mpz_div(h1,h1,mpq_numref(RAT_EPS.p->rat));
1032  mpz_div(h1,h1,mpq_numref(RAT_EPS.p->rat));
1033  mpz_sqrt(h1,h1);
1034  mpz_mul(h2,mpq_denref(a.p->rat),mpq_denref(RAT_EPS.p->rat));
1035  mpz_div(h2,h2,mpq_numref(RAT_EPS.p->rat));
1036  mpq_set_num(erg.p->rat,h1);
1037  mpq_set_den(erg.p->rat,h2);
1038  mpq_canonicalize(erg.p->rat);
1039  mpz_clear(h1);
1040  mpz_clear(h2);
1041  return erg;
1042}
1043
1044Rational
1045exp(const Rational& a)
1046{
1047  Rational
1048    erg(1),
1049    pow(a),
1050    fak(1),
1051    rem,
1052    i(2);
1053
1054  do{
1055    erg+=pow/fak;
1056    fak*=i++;
1057    pow*=a;
1058    rem=(((Rational)2)*abs(pow))/fak;
1059    if (rem<RAT_EPS) break;
1060  } while (true);
1061  return erg;
1062}
1063
1064Rational
1065sin(const Rational& a)
1066{
1067  Rational
1068    erg(a),
1069    pow(a),
1070    fak(6),
1071    rem,
1072    i(4),
1073    aq(a);
1074
1075  aq*=a;
1076  do {
1077    pow*=-aq;
1078    erg+=pow/fak;
1079    fak*=i++;
1080    fak*=i++;
1081    rem=(abs(pow*aq))/fak;
1082    if (rem<RAT_EPS) break;
1083  } while(true);
1084  return erg;
1085}
1086
1087Rational
1088cos(const Rational& a)
1089{
1090  Rational
1091    erg(1),
1092    pow(1),
1093    fak(2),
1094    rem,
1095    i(3),
1096    aq(a);
1097
1098  aq*=a;
1099  do {
1100    pow*=-aq;
1101    erg+=pow/fak;
1102    fak*=i++;
1103    fak*=i++;
1104    rem=(abs(pow*aq))/fak;
1105    if (rem<RAT_EPS) break;
1106  } while(true);
1107  return erg;
1108}
1109
1110Rational
1111tan(const Rational& a)
1112{
1113  Rational
1114    ergsin(a),
1115    ergcos(1),
1116    pow(a),
1117    fak(2),
1118    i(3),
1119    rem,
1120    eps2(((Rational)2)*RAT_EPS);
1121
1122  do{
1123    pow*=a*(-1);
1124    ergcos+=pow/fak;
1125    pow*=a;
1126    fak*=i++;
1127    ergsin+=pow/fak;
1128    fak*=i++;
1129    rem=(abs(pow*a*a))/fak;
1130    if (rem<eps2) break;
1131  } while(true);
1132  return ergsin/ergcos;
1133}
1134
1135Rational gcd( const Rational &a,const Rational &b )
1136{
1137    if( a == 0 )
1138    {
1139        if( b == 0 )
1140        {
1141            return  (Rational)1;
1142        }
1143        else
1144        {
1145            return  abs( b );
1146        }
1147    }
1148    else if( b == 0 )
1149    {
1150        return  abs( a );
1151    }
1152
1153    Rational erg;
1154
1155    mpz_gcd( mpq_numref( erg.p->rat ),
1156            mpq_numref( a.p->rat ),mpq_numref( b.p->rat ) );
1157    mpz_gcd( mpq_denref( erg.p->rat ),
1158            mpq_denref( a.p->rat ),mpq_denref( b.p->rat ) );
1159
1160    //mpq_canonicalize( erg.p->rat );
1161
1162    return  abs( erg );
1163}
1164
1165Rational gcd( Rational *a,int n )
1166{
1167    if( n == 1 )
1168    {
1169        return  a[0];
1170    }
1171
1172    Rational g = gcd( a[0],a[1] );
1173
1174    for( int i=2; i<n; i++ )
1175    {
1176        g = gcd( g,a[i] );
1177    }
1178
1179    return  g;
1180}
1181
1182Rational lcm( const Rational &a,const Rational &b )
1183{
1184    if( a == 0 )
1185    {
1186        return b;
1187    }
1188    else if( b == 0 )
1189    {
1190        return a;
1191    }
1192
1193    return a*b/gcd(a,b);
1194}
1195
1196Rational lcm( Rational *a,int n )
1197{
1198    if( n == 1 )
1199    {
1200        return  a[0];
1201    }
1202
1203    Rational g = lcm( a[0],a[1] );
1204
1205    for( int i=2; i<n; i++ )
1206    {
1207        g = lcm( g,a[i] );
1208    }
1209
1210    return  g;
1211}
1212
1213double  Rational::complexity( ) const
1214{
1215    double num = mpz_get_d( mpq_numref( p->rat ) );
1216    double den = mpz_get_d( mpq_denref( p->rat ) );
1217
1218    if( num < 0 ) num = -num;
1219    if( den < 0 ) den = -den;
1220
1221    return  ( num > den ? num : den );
1222}
1223
1224#endif /* HAVE_SPECTRUM */
1225// ----------------------------------------------------------------------------
1226//  GMPrat.cc
1227//  end of file
1228// ----------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.