source: sage/rings/rational.pyx @ 6224:0ac6643cec6a

Revision 6224:0ac6643cec6a, 45.3 KB checked in by Robert Bradshaw <robertwb@…>, 6 years ago (diff)

integer, rational powering
Line 
1r"""
2Rational Numbers
3
4AUTHORS:
5    -- William Stein (2005): first version
6    -- William Stein (2006-02-22): floor and ceil (pure fast GMP versions).
7    -- Gonzalo Tornaria and William Stein (2006-03-02): greatly improved
8                    python/GMP conversion; hashing
9    -- William Stein and Naqi Jaffery (2006-03-06): height, sqrt examples,
10          and improve behavior of sqrt.
11    -- David Harvey (2006-09-15): added nth_root
12    -- Pablo De Napoli (2007-04-01): corrected the implementations of
13       multiplicative_order, is_one; optimized __nonzero__ ; documented: lcm,gcd
14
15TESTS:
16    sage: a = -2/3
17    sage: a == loads(dumps(a))
18    True
19
20
21"""
22
23
24###########################################################################
25#       Copyright (C) 2004, 2006 William Stein <wstein@gmail.com>
26#
27#  Distributed under the terms of the GNU General Public License (GPL)
28#                  http://www.gnu.org/licenses/
29###########################################################################
30
31include "../ext/interrupt.pxi"  # ctrl-c interrupt block support
32include "../ext/gmp.pxi"
33include "../ext/stdsage.pxi"
34
35
36import operator
37
38from sage.misc.mathml import mathml
39
40import sage.misc.misc as misc
41import sage.rings.rational_field
42import sage.libs.pari.all
43
44cimport integer
45import integer
46
47from integer_ring import ZZ
48
49from sage.structure.element cimport Element, RingElement, ModuleElement
50from sage.structure.element import bin_op
51from sage.categories.morphism cimport Morphism
52
53import sage.rings.real_mpfr
54
55cimport sage.ext.arith
56import  sage.ext.arith
57
58cdef sage.ext.arith.arith_int ai
59ai = sage.ext.arith.arith_int()
60
61import sys
62MAX_UNSIGNED_LONG = 2 * sys.maxint  #TODO: this was copied from integer.pyx -- should this be defined somewhere more general?
63
64cdef extern from "mpz_pylong.h":
65    cdef mpz_get_pylong(mpz_t src)
66    cdef int mpz_set_pylong(mpz_t dst, src) except -1
67    cdef long mpz_pythonhash(mpz_t src)
68
69cdef class Rational(sage.structure.element.FieldElement)
70
71cdef public void set_from_mpq(Rational self, mpq_t value):
72    mpq_set(self.value, value)
73   
74cdef public void set_from_Rational(Rational self, Rational other):
75    mpq_set(self.value, other.value)
76
77cdef public void set_from_Integer(Rational self, integer.Integer other):
78    mpq_set_z(self.value, other.value)
79
80cdef object Rational_mul_(Rational a, Rational b):
81    cdef Rational x
82    x = <Rational> PY_NEW(Rational)
83
84    _sig_on
85    mpq_mul(x.value, a.value, b.value)
86    _sig_off   
87   
88    return x
89
90cdef object Rational_div_(Rational a, Rational b):
91    cdef Rational x
92    x = <Rational> PY_NEW(Rational)
93
94    _sig_on
95    mpq_div(x.value, a.value, b.value)
96    _sig_off
97   
98    return x
99
100cdef Rational_add_(Rational self, Rational other):
101    cdef Rational x
102    x = <Rational> PY_NEW(Rational)   
103    _sig_on
104    mpq_add(x.value, self.value, other.value)
105    _sig_off
106    return x
107
108cdef Rational_sub_(Rational self, Rational other):
109    cdef Rational x
110    x = <Rational> PY_NEW(Rational)   
111
112    _sig_on
113    mpq_sub(x.value, self.value, other.value)
114    _sig_off
115   
116    return x
117
118cdef object the_rational_ring
119the_rational_ring = sage.rings.rational_field.Q
120
121# make sure zero/one elements are set
122cdef set_zero_one_elements():
123    global the_rational_ring
124    the_rational_ring._zero_element = Rational(0)
125    the_rational_ring._one_element = Rational(1)
126
127set_zero_one_elements()
128   
129cdef class Rational(sage.structure.element.FieldElement):
130    """
131    A Rational number.
132
133    Rational numbers are implemented using the GMP C library.
134
135    EXAMPLES:
136        sage: a = -2/3
137        sage: type(a)
138        <type 'sage.rings.rational.Rational'>
139        sage: parent(a)
140        Rational Field
141        sage: Rational('1/0')
142        Traceback (most recent call last):
143        ...
144        TypeError: unable to convert 1/0 to a rational       
145        sage: Rational(1.5)
146        3/2
147        sage: Rational('9/6')
148        3/2
149        sage: Rational((2^99,2^100))
150        1/2
151        sage: Rational(("2", "10"), 16)
152        1/8
153
154    """
155    def __new__(self, x=None, int base=0):
156        global the_rational_ring
157        mpq_init(self.value)
158        self._parent = the_rational_ring
159
160    def __init__(self, x=None, unsigned int base=0):
161        if not (x is None):
162            self.__set_value(x, base)
163
164    def __reduce__(self):
165        return sage.rings.rational.make_rational, (self.str(32),)
166
167    def __index__(self):
168        """
169        Needed so integers can be used as list indices.
170
171        EXAMPLES:
172            sage: v = [1,2,3,4,5]
173            sage: v[3/1]
174            4
175            sage: v[3/2]
176            Traceback (most recent call last):
177            ...
178            TypeError: rational is not an integer
179        """
180        if self.denominator() == 1:
181            return int(self)
182        raise TypeError, "rational is not an integer"
183
184    def _reduce_set(self, s):
185        mpq_set_str(self.value, s, 32)
186
187    def __set_value(self, x, unsigned int base):
188        cdef int n
189        cdef Rational temp_rational
190       
191        if isinstance(x, Rational):
192            set_from_Rational(self, x)
193           
194        elif isinstance(x, int):
195            i = x
196            mpq_set_si(self.value, i, 1)
197
198        elif isinstance(x, long):
199            mpz_set_pylong(mpq_numref(self.value), x)
200           
201        elif isinstance(x, integer.Integer):
202            set_from_Integer(self, x)
203
204        elif isinstance(x, sage.rings.real_mpfr.RealNumber):
205           
206            if x == 0:
207                mpq_set_si(self.value, 0, 1)
208                return
209            if not base:
210                set_from_Rational(self, x.simplest_rational())
211            else:
212                xstr = x.str(base)
213                if '.' in xstr:
214                    exp = (len(xstr) - (xstr.index('.') +1))
215                    p = base**exp
216                    pstr = '1'+'0'*exp
217                    s = xstr.replace('.','') +'/'+pstr
218                    n = mpq_set_str( self.value, s, base)
219                    if n or mpz_cmp_si(mpq_denref(self.value), 0) == 0:
220                        raise TypeError, "unable to convert %s to a rational"%x
221                    mpq_canonicalize(self.value)                   
222                else:
223                    n = mpq_set_str(self.value, xstr, base)
224                    if n or mpz_cmp_si(mpq_denref(self.value), 0) == 0:
225                        raise TypeError, "unable to convert %s to a rational"%x
226                    mpq_canonicalize(self.value)
227           
228        elif isinstance(x, str):
229            n = mpq_set_str(self.value, x, base)
230            if n or mpz_cmp_si(mpq_denref(self.value), 0) == 0:
231                raise TypeError, "unable to convert %s to a rational"%x
232            mpq_canonicalize(self.value)
233           
234        elif hasattr(x, "_rational_"):
235            set_from_Rational(self, x._rational_())
236           
237        elif isinstance(x, tuple) and len(x) == 2:
238            num = x[0]
239            denom = x[1]
240            if isinstance(num, int) and isinstance(denom, int):
241                mpq_set_si(self.value, num, denom)
242            else:
243                if not isinstance(num, integer.Integer):
244                    num = integer.Integer(num, base)
245                if not isinstance(denom, integer.Integer):
246                    denom = integer.Integer(denom, base)
247                mpz_set(mpq_numref(self.value), (<integer.Integer>num).value)
248                mpz_set(mpq_denref(self.value), (<integer.Integer>denom).value)
249            if mpz_sgn(mpq_denref(self.value)) == 0:
250                raise ValueError, "denominator must not be 0"
251            mpq_canonicalize(self.value)
252
253        elif isinstance(x, list) and len(x) == 1:
254            self.__set_value(x[0], base)
255           
256        elif isinstance(x, sage.libs.pari.all.pari_gen):
257            # TODO: figure out how to convert to pari integer in base 16
258            s = str(x)
259            n = mpq_set_str(self.value, s, 0)
260            if n or mpz_cmp_si(mpq_denref(self.value), 0) == 0:
261                raise TypeError, "Unable to coerce %s (%s) to Rational"%(x,type(x))
262           
263        elif hasattr(x, 'rational_reconstruction'):
264            temp_rational = x.rational_reconstruction()
265            mpq_set(self.value, temp_rational.value)
266       
267        else:
268           
269            raise TypeError, "Unable to coerce %s (%s) to Rational"%(x,type(x))
270
271    cdef void set_from_mpq(Rational self, mpq_t value):
272        mpq_set(self.value, value)
273
274
275    def list(self):
276        """
277        Return a list with the rational element in it, to be
278        compatible with the method for number fields.
279
280        EXAMPLES:
281        sage: m = 5/3
282        sage: m.list()
283        [5/3]
284        """
285        return [ self ]
286
287
288    def __richcmp__(left, right, int op):
289        """
290        EXAMPLES:
291            sage: 1/3 < 2/3
292            True
293            sage: 2/3 < 1/3
294            False
295            sage: 4/5 < 2.0
296            True
297            sage: 4/5 < 0.8
298            False
299        """
300        return (<sage.structure.element.Element>left)._richcmp(right, op)
301
302    cdef int _cmp_c_impl(left, sage.structure.element.Element right) except -2:
303        cdef int i
304        i = mpq_cmp((<Rational>left).value, (<Rational>right).value)
305        if i < 0: return -1
306        elif i == 0: return 0
307        else: return 1
308
309    def copy(self):
310        cdef Rational z
311        z = <Rational> PY_NEW(Rational)           
312        mpq_set(z.value, self.value)
313        return z
314
315    def  __dealloc__(self):
316        mpq_clear(self.value)
317       
318    def __repr__(self):
319        return self.str()
320
321    def _latex_(self):
322        if self.denom() == 1:
323            return str(self.numer())
324        else:
325            if self < 0:
326                return "-\\frac{%s}{%s}"%(-self.numer(), self.denom())
327            else:
328               return "\\frac{%s}{%s}"%(self.numer(), self.denom())
329
330    def _mathml_(self):
331        if self.denom() == 1:
332            return '<mn>%s</mn>'%(self.numer())
333        else:
334            t = ''
335            if self < 0:
336                t = t + '<mo>-</mo>'
337            t = t + '<mfrac><mrow>%s</mrow><mrow>%s</mrow></mfrac>'%(
338                mathml(abs(self.numer())), mathml(self.denom()))
339            return t
340
341    def _im_gens_(self, codomain, im_gens):
342        return codomain._coerce_(self)
343
344    def lcm(self, Rational other):
345        r"""
346        Return the least common multiple of self and other.
347       
348        One way to define this notion is the following:
349       
350        Note that each rational positive rational number can be written
351        as a product of primes with integer (positive or negative)
352        powers in a unique way.
353       
354        Then, the LCM of two rational numbers x,y can be defined by
355        specifying that the  exponent of every prime p in lcm(x,y)
356        is the supremum of the exponents of p in x, 
357        and the exponent of p  in y
358        (The primes that does not appear in the decomposition of x
359        or y are considered to have exponent zero).
360       
361        This definition  is consistent with the definition of the LCM
362        in the rational integers. Our hopefully interesting notion of LCM
363        for rational numbers is illustrated in the examples below.
364
365        EXAMPLES:
366       
367            sage: lcm(2/3,1/5)
368            2
369           
370            This is consistent with the definition above, since:
371            $$
372               2/3 = 2^1 * 3^{-1}*5^0
373            $$
374            $$
375               1/5 = 2^0 * 3^0   *5^{-1}
376            $$
377            and hence,
378            $$
379               lcm(2/3,1/5)= 2^1*3^0*5^0 = 2.
380            $$
381           
382            sage: lcm(2/3,7/5)
383            14
384           
385            In this example:
386            $$
387               2/3 = 2^1*3^{-1}*5^0    * 7^0
388            $$
389            $$
390               7/5 = 2^0*3^0   *5^{-1} * 7^1
391            $$
392            $$
393               lcm(2/3,7/5) = 2^1*3^0*5^0*7^1 = 14
394            $$
395           
396            sage: lcm(1/3,1/5)
397            1
398           
399            In this example:
400            $$
401               1/3 = 3^{-1}*5^0
402            $$
403            $$
404               1/5 = 3^0 * 5^{-1}
405            $$
406            $$
407               lcm(1/3,1/5)=3^0*5^0=1
408            $$
409           
410            sage: lcm(1/3,1/6)
411            1/3
412           
413            In this example:
414            $$
415               1/3 = 2^0*3^{-1}
416            $$
417            $$
418               1/6 = 2^{-1}*3^{-1}
419            $$
420            $$
421               lcm(1/3,1/6)=2^0*3^{-1}=1/3
422            $$
423        """
424        d = self.denom()*other.denom()
425        self_d = self.numer()*other.denom()
426        other_d = other.numer()*self.denom()
427        return self_d.lcm(other_d) / d
428
429    def gcd(self, Rational other):
430        """
431        Return the least common multiple of self and other.
432       
433        One way to define this notion is the following:
434       
435        Note that each rational positive rational number can be written
436        as a product of primes  with integer (positive or negative)
437        powers in a unique way.
438       
439        Then, the GCD of two rational numbers x,y can be defined by
440        specifying that the exponent of every prime p in gcd(x,y) is
441        the infimum of the exponents of p in x,
442        and  the exponent of p  in y
443        (The primes that does not appear in the decomposition of x or y
444        are considered to have exponent zero).
445       
446        This definition is consistent with the definition of the GCD
447        in the rational integers.  Our hopefully interesting notion of GCD
448        for rational numbers is illustrated in the examples below.
449       
450        EXAMPLES:
451           
452            sage: gcd(2/3,1/5)
453            1/15
454           
455            This is consistent with the definition above, since:
456            $$
457               2/3 = 2^1 * 3^{-1}*5^0
458            $$
459            $$
460               1/5 = 2^0 * 3^0   *5^{-1}
461            $$
462            and hence,
463            $$
464               gcd(2/3,1/5)= 2^0*3^{-1}*5^{-1} = 1/15
465            $$
466           
467            sage: gcd(2/3,7/5)
468            1/15
469           
470            In this example:
471            $$
472              2/3 = 2^1*3^{-1}*5^0    * 7^0
473            $$
474            $$
475              7/5 = 2^0*3^0   *5^{-1} * 7^1
476            $$
477            $$
478              gcd(2/3,7/5) = 2^0*3^{-1}*5^{-1}*7^0 = 1/15
479            $$
480           
481            sage: gcd(1/3,1/6)
482            1/6
483           
484            In this example:
485            $$
486              1/3 = 2^0*3^{-1}
487            $$
488            $$
489              1/6 = 2^{-1}*3^{-1}
490            $$
491            $$
492              gcd(1/3,1/6)=2^{-1}*3^{-1}=1/6
493            $$
494           
495            sage: gcd(6/7,9/7)
496            3/7
497
498            In this example:
499            $$
500              6/7 = 2^1*3^1*7^{-1}
501            $$
502            $$
503              9/7 = 2^0*3^2*7^{-1}
504            $$
505            $$
506              gcd(6/7,9/7)=2^0*3^1*7^{-1}=3/7
507            $$
508        """
509        d = self.denom()*other.denom()
510        self_d = self.numer()*other.denom()
511        other_d = other.numer()*self.denom()
512        return self_d.gcd(other_d) / d
513
514    def valuation(self, p):
515        return self.numerator().valuation(p) - self.denominator().valuation(p)
516       
517    def is_square(self):
518        """
519        EXAMPLES:
520            sage: x = 9/4
521            sage: x.is_square()
522            True
523            sage: x = (7/53)^100
524            sage: x.is_square()
525            True
526            sage: x = 4/3
527            sage: x.is_square()
528            False
529            sage: x = -1/4
530            sage: x.is_square()
531            False
532        """
533        return mpq_sgn(self.value) >= 0 and mpz_perfect_square_p(mpq_numref(self.value)) and mpz_perfect_square_p(mpq_denref(self.value))
534
535    def sqrt_approx(self, prec=None, all=False):
536        """
537        EXAMPLES:
538            sage: (5/3).sqrt_approx()
539            1.29099444873581
540            sage: (990829038092384908234098239048230984/4).sqrt_approx()
541            497701978620837137.47374920870362581922510725585130996993055116540856385
542            sage: (5/3).sqrt_approx(prec=200)
543            1.2909944487358056283930884665941332036109739017638636088625
544            sage: (9/4).sqrt_approx()       
545            3/2
546        """
547        try:
548            return self.sqrt(extend=False,all=all)
549        except ValueError:
550            pass
551        if prec is None:
552            prec = max(max(53, 2*(mpz_sizeinbase(mpq_numref(self.value), 2)+2)),
553                   2*(mpz_sizeinbase(mpq_denref(self.value), 2)+2))
554        return self.sqrt(prec=prec, all=all)
555
556    def val_unit(self, p):
557        r"""
558        Returns a pair: the p-adic valuation of self, and the p-adic
559        unit of self, as a Rational.
560
561        We do not require the p be prime, but it must be at least 2.
562        For more documentation see \code{Integer.val_unit}
563
564        AUTHOR:
565            -- David Roe (4/12/07)
566        """
567        return self._val_unit(p)
568
569    cdef _val_unit(Rational self, integer.Integer p):
570        cdef integer.Integer v
571        cdef Rational u
572        if mpz_cmp_ui(p.value, 2) < 0:
573            raise ValueError, "p must be at least 2."
574        if mpq_sgn(self.value) == 0:
575            import sage.rings.infinity
576            u = PY_NEW(Rational)
577            mpq_set_ui(u.value, 1, 1)
578            return (sage.rings.infinity.infinity, u)
579        v = PY_NEW(integer.Integer)
580        u = PY_NEW(Rational)
581        _sig_on
582        mpz_set_ui(v.value, mpz_remove(mpq_numref(u.value), mpq_numref(self.value), p.value))
583        _sig_off
584        if mpz_sgn(v.value) != 0:
585            mpz_set(mpq_denref(u.value), mpq_denref(self.value))
586        else:
587            _sig_on
588            mpz_set_ui(v.value, mpz_remove(mpq_denref(u.value), mpq_denref(self.value), p.value))
589            _sig_off
590            mpz_neg(v.value, v.value)
591        return (v, u)
592
593    def sqrt(self, prec=None, extend=True, all=False):
594        r"""
595        The square root function.
596       
597        INPUT:
598            prec -- integer (default: None): if None, returns an exact
599                 square root; otherwise returns a numerical square
600                 root if necessary, to the given bits of precision.
601            extend -- bool (default: True); if True, return a square
602                 root in an extension ring, if necessary. Otherwise,
603                 raise a ValueError if the square is not in the base
604                 ring.
605            all -- bool (default: False); if True, return all square
606                   roots of self, instead of just one.
607
608        EXAMPLES:
609            sage: x = 25/9
610            sage: x.sqrt()
611            5/3
612            sage: x = 64/4
613            sage: x.sqrt()
614            4
615            sage: x = 100/1
616            sage: x.sqrt()
617            10
618            sage: x.sqrt(all=True)
619            [10, -10]
620            sage: x = 81/5
621            sage: x.sqrt()
622            9/sqrt(5)
623            sage: x = -81/3
624            sage: x.sqrt()
625            3*sqrt(3)*I
626
627            sage: n = 2/3
628            sage: n.sqrt()
629            sqrt(2)/sqrt(3)
630            sage: n.sqrt(prec=10)
631            0.82
632            sage: n.sqrt(prec=100)
633            0.81649658092772603273242802490
634            sage: n.sqrt(prec=100)^2
635            0.66666666666666666666666666667
636            sage: n.sqrt(prec=53, all=True)
637            [0.816496580927726, -0.816496580927726]
638            sage: n.sqrt(extend=False, all=True)
639            Traceback (most recent call last):
640            ...
641            ValueError: square root of 2/3 not a rational number
642            sage: sqrt(-2/3, all=True)
643            [sqrt(2)*I/sqrt(3), -sqrt(2)*I/sqrt(3)]
644            sage: sqrt(-2/3, prec=53)
645            0.816496580927726*I
646            sage: sqrt(-2/3, prec=53, all=True)
647            [0.816496580927726*I, -0.816496580927726*I]
648           
649        AUTHOR:
650            -- Naqi Jaffery (2006-03-05): some examples
651        """
652        if mpq_sgn(self.value) == 0:
653            return [self] if all else self
654           
655        if mpq_sgn(self.value) < 0:
656            if not extend:
657                raise ValueError, "square root of negative number not rational"
658            if prec:
659                from sage.rings.complex_field import ComplexField
660                K = ComplexField(prec)
661                return K(self).sqrt(all=all)
662            from sage.calculus.calculus import sqrt
663            return sqrt(self, all=all)
664
665        cdef Rational z = <Rational> PY_NEW(Rational)           
666        cdef mpz_t tmp
667        cdef int non_square = 0
668
669        _sig_on
670        mpz_init(tmp)
671        mpz_sqrtrem(mpq_numref(z.value), tmp, mpq_numref(self.value))
672        if mpz_sgn(tmp) != 0:
673            non_square = 1
674        else:
675            mpz_sqrtrem(mpq_denref(z.value), tmp, mpq_denref(self.value))
676            if mpz_sgn(tmp) != 0:
677                non_square = 1
678        mpz_clear(tmp)
679        _sig_off
680
681        if non_square:
682            if not extend:
683                raise ValueError, "square root of %s not a rational number"%self
684            if prec:
685                from sage.rings.real_mpfr import RealField
686                K = RealField(prec)
687                return K(self).sqrt(all=all)
688            from sage.calculus.calculus import sqrt
689            return sqrt(self, all=all)
690        if all:
691            return [z, -z]
692        return z
693
694    def period(self):
695        r"""
696        Return the period of the repeating part of the decimal
697        expansion of this rational number.
698
699        ALGORITHM: When a rational number $n/d$ with $(n,d)==1$ is
700        expanded, the period begins after $s$ terms and has length
701        $t$, where $s$ and $t$ are the smallest numbers satisfying
702        $10^s=10^(s+t) (mod d)$.  When $d$ is coprime to 10, this
703        becomes a purely periodic decimal with $10^t=1 (mod d)$.
704        (Lehmer 1941 and Mathworld).
705
706        EXAMPLES:
707            sage: (1/7).period()
708            6
709            sage: RR(1/7)
710            0.142857142857143
711            sage: (1/8).period()
712            1
713            sage: RR(1/8)
714            0.125000000000000
715            sage: RR(1/6)
716            0.166666666666667
717            sage: (1/6).period()
718            1
719            sage: x = 333/106
720            sage: x.period()
721            13
722            sage: RealField(200)(x)
723            3.1415094339622641509433962264150943396226415094339622641509
724        """
725        cdef unsigned int alpha, beta
726        d = self.denominator()
727        alpha = d.valuation(2)
728        beta = d.valuation(5)
729        P = d.parent()
730        if alpha > 0 or beta > 0:
731            d = d//(P(2)**alpha * P(5)**beta)
732        from sage.rings.integer_mod import Mod
733        a = Mod(P(10),d)
734        return a.multiplicative_order()
735
736    def nth_root(self, int n):
737        r"""
738        Computes the nth root of self, or raises a \exception{ValueError}
739        if self is not a perfect nth power.
740
741        INPUT:
742            n -- integer (must fit in C int type)
743
744        AUTHOR:
745           -- David Harvey (2006-09-15)
746
747        EXAMPLES:
748          sage: (25/4).nth_root(2)
749          5/2
750          sage: (125/8).nth_root(3)
751          5/2
752          sage: (-125/8).nth_root(3)
753          -5/2
754          sage: (25/4).nth_root(-2)
755          2/5
756
757          sage: (9/2).nth_root(2)
758          Traceback (most recent call last):
759          ...
760          ValueError: not a perfect nth power
761
762          sage: (-25/4).nth_root(2)
763          Traceback (most recent call last):
764          ...
765          ValueError: cannot take even root of negative number
766
767        """
768        # TODO -- this could be quicker, by using GMP directly.
769        cdef integer.Integer num
770        cdef integer.Integer den
771        cdef int negative
772
773        if n > 0:
774            negative = 0
775        elif n < 0:
776            n = -n
777            negative = 1
778        else:
779            raise ValueError, "n cannot be zero"
780
781        num, exact = self.numerator().nth_root(n, 1)
782        if not exact:
783            raise ValueError, "not a perfect nth power"
784
785        den, exact = self.denominator().nth_root(n, 1)
786        if not exact:
787            raise ValueError, "not a perfect nth power"
788
789        if negative:
790            return den / num
791        else:
792            return num / den
793       
794       
795    def str(self, int base=10):
796        if base < 2 or base > 36:
797            raise ValueError, "base (=%s) must be between 2 and 36"%base
798        cdef size_t n
799        cdef char *s
800
801        n = mpz_sizeinbase (mpq_numref(self.value), base) \
802            + mpz_sizeinbase (mpq_denref(self.value), base) + 3
803        s = <char *>PyMem_Malloc(n)
804        if s == NULL:
805            raise MemoryError, "Unable to allocate enough memory for the string representation of an integer."
806       
807        _sig_on
808        mpq_get_str(s, base, self.value)
809        _sig_off
810        k = <object> PyString_FromString(s)
811        PyMem_Free(s)
812        return k
813
814    def __float__(self):
815        return mpq_get_d(self.value)
816       
817    def __hash__(self):
818        cdef long n, d
819        n = mpz_pythonhash(mpq_numref(self.value))
820        d = mpz_pythonhash(mpq_denref(self.value))
821        if d == 1:
822            return n
823        n = n ^ d
824        if n == -1:
825            return -2
826        return n
827           
828##         cdef char *s
829##         cdef int h
830       
831##         s = mpz_get_str(NULL, 16, mpq_numref(self.value))
832##         h = hash(long(s,16))
833##         free(s)
834##         if mpz_cmp_si(mpq_denref(self.value), 1) == 0:
835##             return h
836##         else:
837##             s = mpz_get_str(NULL, 16, mpq_denref(self.value))
838##             h = h ^ hash(long(s,16))     # xor
839##             if h == -1:    # -1 is not a valid return value
840##                 h = -2
841##             free(s)
842##         return h
843       
844        #cdef int n
845        #n = mpz_get_si(mpq_numref(self.value)) * \
846        #    mpz_get_si(mpq_denref(self.value))
847        #if n == -1:
848        #    return -2     # since -1 is not an allowed Python hash for C ext -- it's an error indicator.
849       
850
851    def __getitem__(self, int n):
852        if n == 0:
853            return self
854        raise IndexError, "index n (=%s) out of range; it must be 0"%n
855   
856    def set_si(self, signed long int n):
857        mpq_set_si(self.value, n, 1)
858
859    def set_str(self, s, base=10):
860        valid = mpq_set_str(self.value, s, base)
861        if valid != 0 or mpz_cmp_si(mpq_denref(self.value), 0) == 0:
862            mpq_set_si(self.value, 0, 1)  # so data is valid -- but don't waste time making backup.
863            raise ValueError, "invalid literal (%s); object set to 0"%s
864        mpq_canonicalize(self.value)
865
866    ################################################################
867    # Optimized arithmetic
868    ################################################################
869    cdef ModuleElement _add_c_impl(self, ModuleElement right):
870        cdef Rational x
871        x = <Rational> PY_NEW(Rational)
872        mpq_add(x.value, self.value, (<Rational>right).value)
873        return x
874
875    cdef ModuleElement _sub_c_impl(self, ModuleElement right):
876        # self and right are guaranteed to be Integers
877        cdef Rational x
878        x = <Rational> PY_NEW(Rational)
879        mpq_sub(x.value, self.value, (<Rational>right).value)
880        return x
881
882    cdef ModuleElement _neg_c_impl(self):
883        cdef Rational x
884        x = <Rational> PY_NEW(Rational)
885        mpq_neg(x.value, self.value)
886        return x
887
888    cdef RingElement _mul_c_impl(self, RingElement right):
889        cdef Rational x
890        x = <Rational> PY_NEW(Rational)
891        if mpz_sizeinbase (mpq_numref(self.value), 2)  > 100000 or \
892             mpz_sizeinbase (mpq_denref(self.value), 2) > 100000:
893            # We only use the signal handler (to enable ctrl-c out) in case
894            # self is huge, so the product might actually take a while to compute.
895            _sig_on
896            mpq_mul(x.value, self.value, (<Rational>right).value)
897            _sig_off
898        else:
899            mpq_mul(x.value, self.value, (<Rational>right).value)
900        return x
901
902    cdef RingElement _div_c_impl(self, RingElement right):
903        """
904        EXAMPLES:
905            sage: 2/3
906            2/3
907            sage: 3/0
908            Traceback (most recent call last):
909            ...
910            ZeroDivisionError: Rational division by zero
911        """
912        if mpq_cmp_si((<Rational> right).value, 0, 1) == 0:
913            raise ZeroDivisionError, "Rational division by zero"
914        cdef Rational x
915        x = <Rational> PY_NEW(Rational)               
916        mpq_div(x.value, self.value, (<Rational>right).value)
917        return x
918       
919    ################################################################
920    # Other arithmetic operations.
921    ################################################################
922
923    def __invert__(self):
924        if self.is_zero():
925            raise ZeroDivisionError, "rational division by zero"       
926        cdef Rational x
927        x = <Rational> PY_NEW(Rational)                       
928        mpq_inv(x.value, self.value)
929        return x
930       
931    def __pow__(self, n, dummy):
932        """
933        Raise self to the integer power n.
934
935        EXAMPLES:
936            sage: (2/3)^5
937            32/243
938            sage: (-1/1)^(1/3)
939            -1
940
941        We raise to some interesting powers:
942            sage: (2/3)^I
943            2^I/3^I
944            sage: (2/3)^sqrt(2)
945            2^sqrt(2)/3^sqrt(2)
946            sage: x,y,z,n = var('x,y,z,n')
947            sage: (2/3)^(x^n + y^n + z^n)
948            3^(-z^n - y^n - x^n)*2^(z^n + y^n + x^n)
949            sage: (-7/11)^(tan(x)+exp(x))
950            11^(-tan(x) - e^x)*-7^(tan(x) + e^x)
951            sage: (2/3)^(3/4)
952            2^(3/4)/3^(3/4)
953            sage: (-1/3)^0
954            1
955            sage: (0/1)^0
956            Traceback (most recent call last):
957            ...
958            ArithmeticError: 0^0 is undefined.
959            sage: s = (1/2)^(2^100)
960            Traceback (most recent call last):
961            ...
962            RuntimeError: exponent must be at most 4294967294  # 32-bit
963            RuntimeError: exponent must be at most 18446744073709551614 # 64-bit
964            sage: s = (1/2)^(-2^100)
965            Traceback (most recent call last):
966            ...
967            RuntimeError: exponent must be at most 4294967294  # 32-bit
968            RuntimeError: exponent must be at most 18446744073709551614 # 64-bit
969        """
970        cdef Rational _self = self
971        cdef unsigned long _n
972
973        if not PY_TYPE_CHECK(self, Rational):  #this is here for no good reason apparent to me... should be removed in the future.
974            assert False, "BUG:  Rational.__pow__ called on a non-Rational"
975            return self.__pow__(float(n)) #whose idea was it to float(n)?
976
977        try:
978            _n = n = PyNumber_Index(n)
979        except TypeError:
980            if PY_TYPE_CHECK(n, Rational):
981                # this is the only sensible answer that avoids rounding and
982                # an infinite recursion.
983                from sage.calculus.calculus import SR
984                return SR(self)**SR(n)
985            if PY_TYPE_CHECK(n, Element):
986                return (<Element>n)._parent(self)**n
987            try:
988                return n.parent()(self)**n
989            except AttributeError:
990                try:
991                    return type(n)(self)**n
992                except:
993                    raise TypeError, "exponent (=%s) must be an integer.\nCoerce your numbers to real or complex numbers first."%n
994
995        if PY_TYPE(n) != <void*>int or n > MAX_UNSIGNED_LONG or -n > MAX_UNSIGNED_LONG:
996            raise RuntimeError, "exponent must be at most %s"%MAX_UNSIGNED_LONG
997       
998        cdef Rational x = <Rational> PY_NEW(Rational)                       
999       
1000        if not _n:
1001            if not mpq_sgn(_self.value):
1002                raise ArithmeticError, "0^0 is undefined."
1003            else:
1004                mpq_set_si(x.value, 1, 1)
1005                return x
1006
1007        cdef mpz_t num, den
1008
1009        _sig_on
1010        mpz_init(num)
1011        mpz_init(den)
1012        if n < 0:  # we used to call (self**(-n)).__invert__()) -- this should be loads faster
1013            _n = PyNumber_Index(-n)
1014            mpz_pow_ui(den, mpq_numref(_self.value), _n) #we switch den and num to invert
1015            mpz_pow_ui(num, mpq_denref(_self.value), _n)
1016        else:
1017            mpz_pow_ui(num, mpq_numref(_self.value), _n)
1018            mpz_pow_ui(den, mpq_denref(_self.value), _n)
1019        mpq_set_num(x.value, num)
1020        mpq_set_den(x.value, den)
1021        mpz_clear(num)
1022        mpz_clear(den)
1023        _sig_off
1024       
1025        return x
1026
1027    def __pos__(self):
1028        return self
1029
1030    def __neg__(self):
1031        cdef Rational x
1032        x = <Rational> PY_NEW(Rational)                       
1033        mpq_neg(x.value, self.value)
1034        return x
1035
1036    def __nonzero__(self):
1037        # A rational number is zero iff its numerator is zero.
1038        return mpq_sgn(self.value) != 0
1039
1040    def __abs__(self):
1041        cdef Rational x
1042        x = <Rational> PY_NEW(Rational)                       
1043        mpq_abs(x.value, self.value)
1044        return x
1045
1046    def mod_ui(Rational self, unsigned long int n):
1047        cdef unsigned int num, den, a
1048
1049        # Documentation from GMP manual:
1050        # "For the ui variants the return value is the remainder, and
1051        # in fact returning the remainder is all the div_ui functions do."
1052        _sig_on
1053        num = mpz_fdiv_ui(mpq_numref(self.value), n)
1054        den = mpz_fdiv_ui(mpq_denref(self.value), n)
1055        _sig_off
1056        return int((num * ai.inverse_mod_int(den, n)) % n)
1057
1058    def __mod__(Rational self, other):
1059        other = integer.Integer(other)
1060        if not other:
1061            raise ZeroDivisionError, "Rational modulo by zero"
1062        n = self.numer() % other
1063        d = self.denom() % other
1064        _sig_on
1065        d = d.inverse_mod(other)
1066        _sig_off
1067        return (n*d)%other
1068
1069    def norm(self):
1070        """
1071        Returns the norm from Q to Q of x (which is just x). This was
1072        added for compatibility with NumberFields.
1073
1074        EXAMPLES:
1075            sage: (1/3).norm()
1076             1/3
1077           
1078        AUTHOR:
1079          -- Craig Citro
1080        """
1081        return self
1082
1083    def trace(self):
1084        """
1085        Returns the trace from Q to Q of x (which is just x). This was
1086        added for compatibility with NumberFields.
1087
1088        EXAMPLES:
1089            sage: (1/3).trace()
1090             1/3
1091                       
1092        AUTHOR:
1093          -- Craig Citro
1094        """
1095        return self
1096       
1097    def charpoly(self, var):
1098        """
1099        Return the characteristic polynomial of this rational number.
1100        This will always be just x - self; this is really here
1101        so that code written for number fields won't crash when
1102        applied to rational numbers.
1103
1104        EXAMPLES:
1105            sage: (1/3).charpoly('x')
1106             x - 1/3
1107           
1108        AUTHOR:
1109          -- Craig Citro           
1110        """
1111        QQ = self.parent()
1112        return QQ[var]([-self,1])
1113
1114    def minpoly(self):
1115        """
1116        Return the minimal polynomial of this rational number.
1117        This will always be just x - self; this is really here
1118        so that code written for number fields won't crash when
1119        applied to rational numbers.
1120
1121        EXAMPLES:
1122            sage: (1/3).minpoly()
1123             x - 1/3
1124                       
1125        AUTHOR:
1126          -- Craig Citro
1127        """
1128        QQ = self.parent()
1129        return QQ['x']([-self,1])
1130
1131    cdef integer.Integer _integer_c(self):
1132        if not mpz_cmp_si(mpq_denref(self.value), 1) == 0:
1133            raise TypeError, "no coercion of this rational to integer"
1134        cdef integer.Integer n
1135        n = PY_NEW(integer.Integer)
1136        n.set_from_mpz(mpq_numref(self.value))
1137        return n
1138
1139    def _integer_(self):
1140        return self._integer_c()
1141
1142    def numer(self):
1143        """
1144        Return the numerator of this rational number.
1145
1146        EXAMPLE:
1147            sage: x = -5/11
1148            sage: x.numer()
1149            -5
1150        """
1151        cdef integer.Integer n
1152        n = PY_NEW(integer.Integer)
1153        n.set_from_mpz(mpq_numref(self.value))
1154        return n
1155   
1156    def numerator(self):
1157        """
1158        Return the numerator of this rational number.
1159
1160        EXAMPLE:
1161            sage: x = 5/11
1162            sage: x.numerator()
1163            5
1164
1165            sage: x = 9/3
1166            sage: x.numerator()
1167            3
1168        """
1169        cdef integer.Integer n
1170        n = PY_NEW(integer.Integer)
1171        n.set_from_mpz(mpq_numref(self.value))
1172        return n
1173
1174       
1175    def __int__(self):
1176        """
1177        Return coercion of self to Python int.
1178
1179        This takes the floor of self if self has a denominator (which
1180        is consistent with Python's long(floats)).
1181
1182        EXAMPLES:
1183            sage: int(7/3)
1184            2
1185            sage: int(-7/3)
1186            -3
1187        """
1188        return int(self.__long__())
1189       
1190    def __long__(self):
1191        """
1192        Return coercion of self to Python long.
1193
1194        This takes the floor of self if self has a denominator (which
1195        is consistent with Python's long(floats)).
1196
1197        EXAMPLES:
1198            sage: long(7/3)
1199            2L
1200            sage: long(-7/3)
1201            -3L
1202        """
1203        cdef mpz_t x
1204        if mpz_cmp_si(mpq_denref(self.value),1) != 0:
1205            mpz_init(x)
1206            mpz_fdiv_q(x, mpq_numref(self.value), mpq_denref(self.value))
1207            n = mpz_get_pylong(x)
1208            mpz_clear(x)
1209            return n
1210        else:
1211            return mpz_get_pylong(mpq_numref(self.value))
1212
1213    def denom(self):
1214        """
1215        self.denom(): Return the denominator of this rational number.
1216       
1217        EXAMPLES:
1218            sage: x = 5/13
1219            sage: x.denom()
1220            13
1221            sage: x = -9/3
1222            sage: x.denom()
1223            1
1224        """
1225        cdef integer.Integer n
1226        n = PY_NEW(integer.Integer)
1227        n.set_from_mpz(mpq_denref(self.value))
1228        return n
1229
1230    def denominator(self):
1231        """
1232        self.denominator(): Return the denominator of this rational number.
1233       
1234        EXAMPLES:
1235            sage: x = -5/11
1236            sage: x.denominator()
1237            11
1238            sage: x = 9/3
1239            sage: x.denominator()
1240            1
1241        """
1242        cdef integer.Integer n
1243        n = PY_NEW(integer.Integer)
1244        n.set_from_mpz(mpq_denref(self.value))
1245        return n
1246
1247    def factor(self):
1248        return sage.rings.rational_field.factor(self)
1249
1250    def floor(self):
1251        """
1252        self.floor(): Return the floor of this rational number as an integer.
1253
1254        EXAMPLES:
1255            sage: n = 5/3; n.floor()
1256            1
1257            sage: n = -17/19; n.floor()
1258            -1
1259            sage: n = -7/2; n.floor()
1260            -4
1261            sage: n = 7/2; n.floor()
1262            3
1263            sage: n = 10/2; n.floor()
1264            5
1265        """
1266        cdef integer.Integer n
1267        n = integer.Integer()
1268        mpz_fdiv_q(n.value, mpq_numref(self.value), mpq_denref(self.value))
1269        return n
1270
1271    def ceil(self):
1272        """
1273        self.ceil(): Return the ceiling of this rational number.
1274
1275        If this rational number is an integer, this returns this
1276        number, otherwise it returns the floor of this number +1.
1277
1278        EXAMPLES:
1279            sage: n = 5/3; n.ceil()
1280            2
1281            sage: n = -17/19; n.ceil()
1282            0
1283            sage: n = -7/2; n.ceil()
1284            -3
1285            sage: n = 7/2; n.ceil()
1286            4
1287            sage: n = 10/2; n.ceil()
1288            5
1289        """
1290        cdef integer.Integer n
1291        n = integer.Integer()
1292        mpz_cdiv_q(n.value, mpq_numref(self.value), mpq_denref(self.value))
1293        return n
1294       
1295    def height(self):
1296        """
1297        The max absolute value of the numerator and denominator of self,
1298        as an Integer.
1299
1300        EXAMPLES:
1301            sage: a = 2/3
1302            sage: a.height()
1303            3
1304            sage: a = 34/3
1305            sage: a.height()
1306            34
1307            sage: a = -97/4
1308            sage: a.height()
1309            97
1310
1311        AUTHOR:
1312            -- Naqi Jaffery (2006-03-05): examples
1313        """
1314        x = abs(self.numer())
1315        if x > self.denom():
1316            return x
1317        return self.denom()
1318
1319    def _lcm(self, Rational other):
1320        """
1321        Returns the least common multiple, in the rational numbers,
1322        of self and other.  This function returns either 0 or 1 (as
1323        a rational number).
1324        """
1325        if mpz_cmp_si(mpq_numref(self.value), 0) == 0 and \
1326               mpz_cmp_si(mpq_numref(other.value), 0) == 0:
1327            return Rational(0)
1328        return Rational(1)
1329       
1330    def _gcd(self, Rational other):
1331        """
1332        Returns the least common multiple, in the rational numbers,
1333        of self and other.  This function returns either 0 or 1 (as
1334        a rational number).
1335        """
1336        if mpz_cmp_si(mpq_numref(self.value), 0) == 0 and \
1337               mpz_cmp_si(mpq_numref(other.value), 0) == 0:
1338            return Rational(0)
1339        return Rational(1)
1340       
1341
1342    def additive_order(self):
1343        """
1344        Return the additive order of self.
1345
1346        EXAMPLES:
1347            sage: QQ(0).additive_order()
1348            1
1349            sage: QQ(1).additive_order()
1350            +Infinity
1351        """
1352        import sage.rings.infinity
1353        if self.is_zero():
1354            return integer.Integer(1)
1355        else:
1356            return sage.rings.infinity.infinity
1357
1358       
1359    def multiplicative_order(self):
1360        """
1361        Return the multiplicative order of self.
1362
1363        EXAMPLES:
1364            sage: QQ(1).multiplicative_order()
1365            1
1366            sage: QQ('1/-1').multiplicative_order()
1367            2
1368            sage: QQ(0).multiplicative_order()
1369            +Infinity
1370            sage: QQ('2/3').multiplicative_order()
1371            +Infinity
1372            sage: QQ('1/2').multiplicative_order()
1373            +Infinity
1374        """
1375        import sage.rings.infinity
1376        if self.is_one():
1377            return integer.Integer(1)
1378        elif mpz_cmpabs(mpq_numref(self.value),mpq_denref(self.value))==0:
1379            # if the numerator and the denominator are equal in absolute value,
1380            # then the rational number is -1
1381            return integer.Integer(2)
1382        else:
1383            return sage.rings.infinity.infinity   
1384   
1385    def is_one(self):
1386        r"""
1387        Determine if a rational number is one.
1388       
1389        EXAMPLES:       
1390            sage: QQ(1/2).is_one()
1391            False
1392            sage: QQ(4/4).is_one()
1393            True
1394        """
1395        # A rational number is equal to 1 iff its numerator and denominator are equal
1396        return mpz_cmp(mpq_numref(self.value),mpq_denref(self.value))==0
1397        r"""Test if a rational number is zero
1398       
1399        EXAMPLES:
1400       
1401        sage: QQ(1/2).is_zero()
1402        False
1403        sage: QQ(0/4).is_zero()
1404        True
1405        """
1406       
1407    def is_integral(self):
1408        r"""
1409        Determine if a rational number is integral (i.e is in $\Z$).
1410       
1411        EXAMPLES:
1412            sage: QQ(1/2).is_integral()
1413            False
1414            sage: QQ(4/4).is_integral()
1415            True
1416        """
1417        return bool(self in ZZ)
1418   
1419    cdef _lshift(self, long int exp):
1420        r"""
1421        Return $self*2^exp$
1422        """
1423        cdef Rational x
1424        x = <Rational> PY_NEW(Rational)
1425        _sig_on
1426        if exp < 0:
1427            mpq_div_2exp(x.value,self.value,-exp)
1428        else:
1429            mpq_mul_2exp(x.value,self.value,exp)       
1430        _sig_off
1431        return x
1432
1433    def __lshift__(x,y):
1434        if isinstance(x, Rational) and isinstance(y, (int, long, integer.Integer)):
1435            return (<Rational>x)._lshift(y)
1436        return bin_op(x, y, operator.lshift)
1437   
1438    cdef _rshift(self, long int exp):
1439        r"""
1440        Return $self/2^exp$
1441        """
1442        cdef Rational x
1443        x = <Rational> PY_NEW(Rational)
1444        _sig_on
1445        if exp < 0:
1446            mpq_mul_2exp(x.value,self.value,-exp)
1447        else:
1448            mpq_div_2exp(x.value,self.value,exp)
1449        _sig_off
1450        return x
1451
1452    def __rshift__(x,y):
1453        if isinstance(x, Rational) and isinstance(y, (int, long, integer.Integer)):
1454            return (<Rational>x)._rshift(y)
1455        return bin_op(x, y, operator.rshift)
1456   
1457    ##################################################
1458    # Support for interfaces
1459    ##################################################
1460   
1461    def _pari_(self):
1462        return self.numerator()._pari_()/self.denominator()._pari_()
1463
1464    def _interface_init_(self):
1465        """
1466        EXAMPLES:
1467            sage: kash(3/1).Type()              # optional
1468            elt-fld^rat
1469            sage: magma(3/1).Type()             # optional
1470            FldRatElt
1471        """
1472        return '%s/%s'%(self.numerator(), self.denominator())
1473
1474
1475
1476def pyrex_rational_reconstruction(integer.Integer a, integer.Integer m):
1477    """
1478    Find the rational reconstruction of a mod m, if it exists.
1479    INPUT:
1480        a -- Integer
1481        m -- Integer
1482    OUTPUT:
1483        x -- rings.rational.Rational
1484    """
1485    cdef Rational x
1486    x = <Rational> PY_NEW(Rational)
1487    mpq_rational_reconstruction(x.value, a.get_value()[0], m.get_value()[0])
1488    return x
1489
1490# def test(n):
1491#     cdef mpz_t a, m
1492#     mpz_init(a); mpz_init(m)
1493#     mpz_set_str(a, "133173434946179363436783721312585228097", 0)
1494#     mpz_set_str(m, "10000000000000000000000000000000000000000", 0)
1495#     cdef int k
1496#     cdef mpq_t x
1497#     mpq_init(x)
1498#     import time
1499#     t = time.clock()
1500#     for k from 0 <= k < n:
1501#         mpq_rational_reconstruction(x, a, m)
1502#     print time.clock() - t
1503
1504
1505def make_rational(s):
1506    r = Rational()
1507    r._reduce_set(s)
1508    return r
1509   
1510cdef class Z_to_Q(Morphism):
1511
1512    def __init__(self):
1513        import integer_ring
1514        import rational_field
1515        import sage.categories.homset
1516        Morphism.__init__(self, sage.categories.homset.Hom(integer_ring.ZZ, rational_field.QQ))
1517       
1518    cdef Element _call_c_impl(self, Element x):
1519        cdef Rational rat
1520        rat = <Rational> PY_NEW(Rational)
1521        mpq_set_z(rat.value, (<integer.Integer>x).value)
1522        return rat
1523
Note: See TracBrowser for help on using the repository browser.