source: git/Singular/GMPrat.cc @ 584f84d

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