Teorija brojeva u SAGE-u

verzija: SageMath 8.4

In [1]:
%display ascii_art

Skupovi brojeva

In [2]:
[NN,ZZ,QQ,RR,CC]
Out[2]:
[ Non negative integer semiring, Integer Ring, Rational Field,

 Real Field with 53 bits of precision,

 Complex Field with 53 bits of precision ]
In [3]:
[24 in NN, 24 in ZZ, -32 in NN, -32 in ZZ]
Out[3]:
[ True, True, False, True ]
In [4]:
[2/3 in ZZ, 2/3 in QQ, 2.3 in QQ, 2.3 in RR, sqrt(2) in QQ, sqrt(2) in RR]
Out[4]:
[ False, True, True, True, False, True ]
In [5]:
[sqrt(-2) in RR, sqrt(-2) in CC, I in CC]
Out[5]:
[ False, True, True ]
In [6]:
type(NN)
Out[6]:
<type 'sage.misc.lazy_import.LazyImport'>
In [7]:
type(CC)
Out[7]:
<class 'sage.rings.complex_field.ComplexField_class_with_category'>

Cjelobrojne funkcije "pod" i "strop"

In [8]:
floor(5.23)
Out[8]:
5
In [9]:
ceil(5.23)
Out[9]:
6
In [10]:
floor(pi)
Out[10]:
3
In [11]:
ceil(pi)
Out[11]:
4
In [12]:
floor(-5.23)
Out[12]:
-6
In [13]:
ceil(-5.23)
Out[13]:
-5
In [14]:
plot(floor(x),(x,-3,3),thickness=3,figsize=[4,4])
Out[14]:
In [15]:
plot(floor(x),(x,-3,3),thickness=3,exclude=[-3..3],figsize=[4,4])
Out[15]:
In [16]:
plot(ceil(x),(x,-3,3),thickness=3,figsize=[4,4])
Out[16]:
In [17]:
plot(ceil(x),(x,-3,3),thickness=3,exclude=[-3..3],figsize=[4,4])
Out[17]:
In [18]:
slika1=plot(floor(x),(x,-3,3),color="blue",thickness=3,exclude=[-3..3])
slika2=plot(ceil(x),(x,-3,3),color="red",thickness=3,exclude=[-3..3])
show(slika1+slika2,figsize=[4,4])

Teorem o dijeljenju s ostatkom

In [19]:
500//17
Out[19]:
29
In [20]:
500%17
Out[20]:
7
In [21]:
-500//17
Out[21]:
-30
In [22]:
-500%17
Out[22]:
10

Neočekivani rezultati

In [23]:
500//-17
Out[23]:
-30
In [24]:
500%-17
Out[24]:
-10

Definiramo svoje funkcije koje će raditi prema našim definicijama

In [25]:
def kvocijent(a,b):
    if b==0:
        return "Error: dijeljenje s nulom"
    elif b>0:
        return floor(a/b)
    else:
        return ceil(a/b)
In [26]:
def ostatak(a,b):
    if b==0:
        return "Error: dijeljenje s nulom"
    elif b>0:
        return a-floor(a/b)*b
    else:
        return a-ceil(a/b)*b
In [27]:
kvocijent(500,17)
Out[27]:
29
In [28]:
ostatak(500,17)
Out[28]:
7
In [29]:
kvocijent(500,-17)
Out[29]:
-29
In [30]:
ostatak(500,-17)
Out[30]:
7
In [31]:
ostatak(-50,0)
Out[31]:
Error: dijeljenje s nulom

Najveća zajednička mjera

In [32]:
gcd(30,18)
Out[32]:
6
In [33]:
gcd(-30,18)
Out[33]:
6
In [34]:
gcd(-15,-17)
Out[34]:
1
In [35]:
gcd(252,200)
Out[35]:
4
In [36]:
gcd([15,30,48])
Out[36]:
3

Jedno cjelobrojno rješenje jednadžbe 252x+200y=M(252,200)

In [37]:
xgcd(252,200)
Out[37]:
( 4, -23, 29 )

Najmanji zajednički višekratnik

In [38]:
lcm(8,12)
Out[38]:
24
In [39]:
lcm(-8,12)
Out[39]:
24
In [40]:
lcm([12,5,14])
Out[40]:
420
In [41]:
lcm(252,200)
Out[41]:
12600

Prosti brojevi

Ispitivanje da li je broj prost ili složen

In [42]:
is_prime(1023)
Out[42]:
False
In [43]:
is_prime(17)
Out[43]:
True
In [44]:
map(lambda x: is_prime(x),[1023,17])
Out[44]:
[ False, True ]

Da li je prirodni broj potencija prostog broja

In [45]:
is_prime_power(64)
Out[45]:
True
In [46]:
is_prime_power(111934216)
Out[46]:
False

Rastav broja na proste faktore

In [47]:
%display latex
In [48]:
factor(40)
Out[48]:
In [49]:
factor(4063)
Out[49]:
In [50]:
factor(10468)
Out[50]:
In [51]:
factor(10^30+1)
Out[51]:

Djeljitelji prirodnog broja

In [52]:
%display ascii_art
In [53]:
divisors(40)
Out[53]:
[ 1, 2, 4, 5, 8, 10, 20, 40 ]
In [54]:
divisors(4063)
Out[54]:
[ 1, 17, 239, 4063 ]
In [55]:
divisors(10468)
Out[55]:
[ 1, 2, 4, 2617, 5234, 10468 ]

Prosti djeljitelji prirodnog broja

In [56]:
[faktor for faktor in divisors(40) if is_prime(faktor)]
Out[56]:
[ 2, 5 ]
In [57]:
[faktor for faktor in divisors(4063) if is_prime(faktor)]
Out[57]:
[ 17, 239 ]
In [58]:
[faktor for faktor in divisors(10468) if is_prime(faktor)]
Out[58]:
[ 2, 2617 ]
In [59]:
[faktor for faktor in divisors(10^30+1) if is_prime(faktor)]
Out[59]:
[ 61, 101, 3541, 9901, 27961, 4188901, 39526741 ]

Skup prostih brojeva

In [60]:
Primes()
Out[60]:
Set of all prime numbers: 2, 3, 5, 7, ...
In [61]:
23 in Primes()
Out[61]:
True
In [62]:
1 in Primes()
Out[62]:
False

Prvih 30 prostih brojeva

In [63]:
primes_first_n(30)
Out[63]:
[ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,

 73, 79, 83, 89, 97, 101, 103, 107, 109, 113 ]
In [64]:
list_plot(primes_first_n(30),figsize=[5,5])
Out[64]:
In [65]:
list_plot(primes_first_n(30),plotjoined=True,figsize=[5,5])
Out[65]:
In [66]:
list_plot(primes_first_n(100),color="red",figsize=[5,5])
Out[66]:
In [67]:
list_plot(primes_first_n(100),plotjoined=True,color="red",figsize=[5,5])
Out[67]:

Deseti po redu prosti broj

In [68]:
primes_first_n(10)[-1]
Out[68]:
29

Skup prostih brojeva većih ili jednakih od 41 i strogo manjih od 200

In [69]:
prime_range(41,200)
Out[69]:
[ 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,

 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197,

 199 ]

Skup svih brojeva između 41 i 150 koji su potencije prostih brojeva

In [70]:
prime_powers(41,150)
Out[70]:
[ 41, 43, 47, 49, 53, 59, 61, 64, 67, 71, 73, 79, 81, 83, 89, 97, 101, 103,

 107, 109, 113, 121, 125, 127, 128, 131, 137, 139, 149 ]

Najmanji prosti broj veći od n

In [71]:
next_prime(100)
Out[71]:
101
In [72]:
next_prime(31)
Out[72]:
37
In [73]:
next_prime(44)
Out[73]:
47
In [74]:
next_prime(2^100)
Out[74]:
1267650600228229401496703205653

Najveći prosti broj manji od n

In [75]:
previous_prime(100)
Out[75]:
97
In [76]:
previous_prime(31)
Out[76]:
29
In [77]:
previous_prime(44)
Out[77]:
43
In [78]:
previous_prime(2^100)
Out[78]:
1267650600228229401496703205361

Funkcija π(x)

In [79]:
prime_pi(100)
Out[79]:
25
In [80]:
prime_pi(10000)
Out[80]:
1229
In [81]:
prime_pi(34.56)
Out[81]:
11
In [82]:
prime_pi(pi^100)
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-82-8ef177951bc0> in <module>()
----> 1 prime_pi(pi**Integer(100))

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/functions/prime_pi.pyx in sage.functions.prime_pi.PrimePi.__call__ (build/cythonized/sage/functions/prime_pi.c:4011)()
    194         else:
    195             self.__primeBound = 0u if len(args) < 2 else args[1]
--> 196         return super(PrimePi, self).__call__(args[0], coerce=coerce, hold=hold)
    197 
    198     def _eval_(self, x):

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:11949)()
    996             res = self._evalf_try_(*args)
    997             if res is None:
--> 998                 res = super(BuiltinFunction, self).__call__(
    999                         *args, coerce=coerce, hold=hold)
   1000 

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6920)()
    487             res = g_function_evalv(self._serial, vec, hold)
    488         elif self._nargs == 1:
--> 489             res = g_function_eval1(self._serial,
    490                     (<Expression>args[0])._gobj, hold)
    491         elif self._nargs == 2:

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/functions/prime_pi.pyx in sage.functions.prime_pi.PrimePi._eval_ (build/cythonized/sage/functions/prime_pi.c:4146)()
    246             z = arg_to_uint64(x, 'prime_pi', 'x')
    247         except NotImplementedError:
--> 248             raise
    249         except TypeError:
    250             return None

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/functions/prime_pi.pyx in sage.functions.prime_pi.PrimePi._eval_ (build/cythonized/sage/functions/prime_pi.c:4102)()
    244         cdef uint64_t z
    245         try:
--> 246             z = arg_to_uint64(x, 'prime_pi', 'x')
    247         except NotImplementedError:
    248             raise

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/functions/prime_pi.pyx in sage.functions.prime_pi.arg_to_uint64 (build/cythonized/sage/functions/prime_pi.c:3158)()
     48         return 0ull
     49     if mpz_sizeinbase((<Integer>x).value, 2) > 63:
---> 50         raise NotImplementedError("computation of " + s1 + " for x >= "
     51                 + "2^63 is not implemented")
     52     if mpz_sizeinbase((<Integer>x).value, 2) > 43:

NotImplementedError: computation of prime_pi for x >= 2^63 is not implemented
In [83]:
floor(pi^100).ndigits(10)
Out[83]:
50
In [84]:
n(pi^100)
Out[84]:
5.18784831431959e49
In [85]:
prime_pi(pi^20)
Out[85]:
401480918
In [86]:
[prime_pi(n) for n in xrange(1,50)]
Out[86]:
[ 0, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9,

 9, 9, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 14, 14, 14, 14,

 15, 15, 15 ]
In [87]:
[(n,prime_pi(n)) for n in xrange(1,50)]
Out[87]:
[ ( 1, 0 ), ( 2, 1 ), ( 3, 2 ), ( 4, 2 ), ( 5, 3 ), ( 6, 3 ), ( 7, 4 ), ( 8,

 4 ), ( 9, 4 ), ( 10, 4 ), ( 11, 5 ), ( 12, 5 ), ( 13, 6 ), ( 14, 6 ), ( 15,

 6 ), ( 16, 6 ), ( 17, 7 ), ( 18, 7 ), ( 19, 8 ), ( 20, 8 ), ( 21, 8 ), ( 22,

 8 ), ( 23, 9 ), ( 24, 9 ), ( 25, 9 ), ( 26, 9 ), ( 27, 9 ), ( 28, 9 ), ( 29,

 10 ), ( 30, 10 ), ( 31, 11 ), ( 32, 11 ), ( 33, 11 ), ( 34, 11 ), ( 35, 11 ),

 ( 36, 11 ), ( 37, 12 ), ( 38, 12 ), ( 39, 12 ), ( 40, 12 ), ( 41, 13 ), ( 42,

 13 ), ( 43, 14 ), ( 44, 14 ), ( 45, 14 ), ( 46, 14 ), ( 47, 15 ), ( 48, 15 ),

 ( 49, 15 ) ]

Teorem o prostim brojevima

In [88]:
prvi=list_plot([prime_pi(n) for n in xrange(1,101)])
show(prvi,figsize=[7,5])
In [89]:
drugi=plot(x/log(x),(x,1,100),color="red")
show(drugi,figsize=[7,5])
In [90]:
show(prvi+drugi,figsize=[7,5])
In [91]:
def teorem_prosti_brojevi(n,xy=[7,5]):
    slika1=list_plot([prime_pi(k) for k in xrange(1,n)])
    slika2=plot(x/log(x),(x,1,n),color="red")
    slika3=show(slika1+slika2,figsize=xy)
    return slika3
In [92]:
teorem_prosti_brojevi(100)
In [93]:
teorem_prosti_brojevi(1000)
In [94]:
teorem_prosti_brojevi(10000)
In [95]:
teorem_prosti_brojevi(100000)
In [96]:
[2e4,8e4,1e5]
Out[96]:
[ 20000.0000000000, 80000.0000000000, 100000.000000000 ]

Slučajni prosti brojevi

slučajni prosti broj manji od 200

In [97]:
random_prime(200)
Out[97]:
19
In [98]:
random_prime(200)
Out[98]:
89

lista od 10 slučajnih prostih brojeva manjih od 200

In [99]:
[random_prime(200) for n in xrange(10)]
Out[99]:
[ 83, 137, 149, 37, 107, 199, 107, 61, 73, 73 ]

slučajni prosti broj između 200 i 500

In [100]:
random_prime(500,lbound=200)
Out[100]:
239

slučajni prosti broj sa 12 znamenaka

In [101]:
random_prime(10^12,lbound=10^11)
Out[101]:
777845379931

3x4 matrica slučajnih prostih brojeva manjih od 1000

In [102]:
[[random_prime(1000) for n in xrange(3)] for m in xrange(4)]
Out[102]:
[ [ 521, 641, 521 ], [ 307, 307, 101 ], [ 359, 181, 149 ], [ 83, 397, 17 ] ]
In [103]:
matrix([[random_prime(1000) for n in xrange(3)] for m in xrange(4)])
Out[103]:
[823 127 929]
[149 113 743]
[563 859 691]
[911 227 653]

dva slučajna prosta broja za RSA kriptosustav

In [104]:
(p,q)=(random_prime(2^640,lbound=2^639),random_prime(2^640,lbound=2^639))
In [105]:
p
Out[105]:
2698113121867022929099652292218321507741723168272234365317725824976100071758886568647043602143838777575965180311138795041679886677030147861412188290966464334085999598247936954937247159841280187
In [106]:
q
Out[106]:
3808045006935005427134803261202208657005914558716433084503387835626163276728558372558267215273938762979311071738958580965190516684635734789883089959128090046668428868922052718935729546941010829

brojevi p i q u bazi 10 imaju 193 znamenke, a u bazi 2 imaju 640 znamenaka (bitova)

In [107]:
p.ndigits(10)
Out[107]:
193
In [108]:
p.ndigits(2)
Out[108]:
640
In [109]:
[q.ndigits(10),q.ndigits(2)]
Out[109]:
[ 193, 640 ]

lista znamenaka broja p u bazi 2 (pazite, znamenke su navedene u obrnutom redoslijedu)

In [110]:
p.digits(2)
Out[110]:
[ 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,

 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1,

 1, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1,

 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0,

 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,

 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1,

 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0,

 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1,

 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1,

 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1,

 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0,

 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0,

 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1,

 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0,

 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0,

 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1,

 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1,

 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0,

 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0,

 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1,

 0, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1,

 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,

 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1,

 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0,

 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1 ]
In [111]:
broj=4
(broj.digits(2),list(reversed(broj.digits(2))))
Out[111]:
( [ 0, 0, 1 ], [ 1, 0, 0 ] )

opcija "proof=False" daje pseudoprosti broj

In [112]:
time random_prime(2^640,lbound=2^639)
CPU times: user 891 ms, sys: 6.78 ms, total: 898 ms
Wall time: 904 ms
Out[112]:
4120186589216991920349041312839644621495454916407628964323070146857226616247775401360782463721611815307278546074896434637355072665309769426663387113537010011367316988764945055485865495707931337
In [113]:
time random_prime(2^640,lbound=2^639,proof=False)
CPU times: user 40.9 ms, sys: 39 µs, total: 41 ms
Wall time: 40 ms
Out[113]:
3882241706781602693901877966955985167820382265475757699762840555609499657118763747467174785518529823264358561868056441769620353089267229905907445383537337621449104414698619952801346405977064243
In [114]:
time random_prime(2^640,lbound=2^639,proof=True)
CPU times: user 807 ms, sys: 0 ns, total: 807 ms
Wall time: 812 ms
Out[114]:
3229963912907973669965752999930952842691035834917075514683350917792532476062668156159720187569069338640440679533490611805851509863273609050957760212213838720274953344673857484407868124581571937

Linearne kongruencije

Multiplikativni inverz modulo n

$15x\equiv 1\!\pmod{341}$     rješenje: $x=91$

In [115]:
inverse_mod(15,341)
Out[115]:
91

$10x\equiv 1\!\pmod{12}$     nema rješenja

In [116]:
inverse_mod(10,12)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-116-5194dcfc2b0a> in <module>()
----> 1 inverse_mod(Integer(10),Integer(12))

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/arith/misc.pyc in inverse_mod(a, m)
   2041     """
   2042     try:
-> 2043         return a.inverse_mod(m)
   2044     except AttributeError:
   2045         return Integer(a).inverse_mod(m)

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/rings/integer.pyx in sage.rings.integer.Integer.inverse_mod (build/cythonized/sage/rings/integer.c:41205)()
   6579         sig_off()
   6580         if r == 0:
-> 6581             raise ZeroDivisionError(f"inverse of Mod({self}, {m}) does not exist")
   6582         return ans
   6583 

ZeroDivisionError: inverse of Mod(10, 12) does not exist

za rješavanje bilo koje linearne kongruencije možemo lagano napisati svoju funkciju

In [117]:
def lin_kong(a,b,n):
    d=gcd(a,n)
    if b%d!=0:
        return False
    a1=a//d
    b1=b//d
    n1=n//d
    x=(inverse_mod(a1,n1)*b1)%n1
    rjesenja=[x+k*n1 for k in xrange(d)]
    return rjesenja

$79x\equiv 15\!\pmod{722}$     rješenje: $119$

In [118]:
lin_kong(79,15,722)
Out[118]:
[ 119 ]

$155x\equiv 1185\!\pmod{1405}$     rješenja: $198, 479, 760, 1041, 1322$

In [119]:
lin_kong(155,1185,1405)
Out[119]:
[ 198, 479, 760, 1041, 1322 ]

$15x\equiv 1\!\pmod{341}$     rješenje: $91$

In [120]:
lin_kong(15,1,341)
Out[120]:
[ 91 ]

$10x\equiv 1\!\pmod{12}$     nema rješenja

In [121]:
lin_kong(10,1,12)
Out[121]:
False

$65x\equiv 27\!\pmod{169}$     nema rješenja

In [122]:
lin_kong(65,27,169)
Out[122]:
False

Funkcija koja daje rješenja linearne kongruencije $ax\equiv b\!\pmod{n}$ zajedno s tablicom kakvu dobivamo ručnim rješavanjem.

In [123]:
def lin_kong_rucno(a,b,n):
    qi=[]
    yi=[0,1]
    n1,q,r1,d=a,n//a,n%a,a
    while r1!=0:
        qi.append(q)
        yi.append(yi[-2]-q*yi[-1])
        q=n1//r1
        n1,r1=r1,n1%r1
        if r1!=0: d=r1
    if b%d!=0: return False
    a1,b1,n1=a//d,b//d,n//d
    x1=(yi[-1]*b1)%n1
    rjesenja=[x1+k*n1 for k in xrange(d)]
    tablica=table([["$i$"]+range(-1,len(yi)-1)+["rjesenja"],["$q_i$"," "," "]+qi+[rjesenja],["$y_i$"]+yi+
                   ["$(a_1,b_1,n_1)=(%d,%d,%d)$"%(a1,b1,n1)]])
    return tablica

$79x\equiv 15\!\pmod{722}$     rješenje: $119$

In [124]:
lin_kong_rucno(79,15,722)
Out[124]:
rjesenja

$155x\equiv 1185\!\pmod{1405}$     rješenja: $198, 479, 760, 1041, 1322$

In [125]:
lin_kong_rucno(155,1185,1405)
Out[125]:
rjesenja

$10x\equiv 1\!\pmod{12}$     nema rješenja

In [126]:
lin_kong_rucno(10,1,12)
Out[126]:
False

Kineski teorem o ostacima

$x\equiv 1\!\pmod{4},\quad x\equiv 8\!\pmod{9},\quad x\equiv 10\!\pmod{25}$     rješenje: $x=485$

In [127]:
crt([1,8,10],[4,9,25])
Out[127]:
485

$x\equiv 4\!\pmod{5},\quad x\equiv 1\!\pmod{12},\quad x\equiv 7\!\pmod{14}$     rješenje: $x=49$

In [128]:
crt([4,1,7],[5,12,14])
Out[128]:
49

$x\equiv 2\!\pmod{3},\quad x\equiv 1\!\pmod{12},\quad x\equiv 7\!\pmod{14}$     nema rješenja

In [129]:
crt([2,1,7],[3,12,14])
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-129-72356d14d325> in <module>()
----> 1 crt([Integer(2),Integer(1),Integer(7)],[Integer(3),Integer(12),Integer(14)])

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/arith/misc.pyc in crt(a, b, m, n)
   3112     """
   3113     if isinstance(a, list):
-> 3114         return CRT_list(a, b)
   3115 
   3116     try:

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/arith/misc.pyc in CRT_list(v, moduli)
   3217     from sage.arith.functions import lcm
   3218     for i in range(1, len(v)):
-> 3219         x = CRT(x,v[i],m,moduli[i])
   3220         m = lcm(m,moduli[i])
   3221     return x%m

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/arith/misc.pyc in crt(a, b, m, n)
   3126     q, r = f(g)
   3127     if r != 0:
-> 3128         raise ValueError("No solution to crt problem since gcd(%s,%s) does not divide %s-%s" % (m, n, a, b))
   3129     from sage.arith.functions import lcm
   3130 

ValueError: No solution to crt problem since gcd(3,12) does not divide 2-1

Funkcija koja daje rješenje sustava kongruencija pomoću kineskog teorema o ostacima na način kojeg provodimo ručnim rješavanjem. Na izlazu se vraća rječnik u kojemu su dane vrijednosti za $n, k_i, x_i, x_0, x_{0,mod}$.

In [130]:
def KTO_rucno(A,N):
    """Daje rjesenje sustava kongruencija x=a_i(mod n_i) pomocu kineskog teorema o ostacima
    na nacin kojeg provodimo rucnim rjesavanjem.
    
    Moduli n_i moraju biti u parovima relativno prosti. Na izlazu se vraca rjecnik u kojemu su
    dane vrijednosti za n, k_i, x_i, x_0 i x0_mod.
    
    A je lista a_i-ova, a N je lista n_i-ova."""
    for t in Combinations(N,2).list():
        if gcd(t[0],t[1])!=1:
            return "Error: moduli nisu u parovima relativno prosti."
    n=prod(N)
    ki=map(lambda x: n/x,N)
    xi=[]
    for i in xrange(len(N)):
        xi.append(int(solve_mod([ki[i]*x==A[i]],N[i])[0][0]))
    x0=sum(map(lambda x: x[0]*x[1],zip(ki,xi)))
    x0mod=x0%n
    izlaz=dict(zip(['n','k_i','x_i','x0','x0_mod'],[n,ki,xi,x0,x0mod]))
    return izlaz

$x\equiv 1\!\pmod{4},\quad x\equiv 8\!\pmod{9},\quad x\equiv 10\!\pmod{25}$     rješenje: $x=485$

In [131]:
KTO_rucno([1,8,10],[4,9,25])
Out[131]:
{ k_i:[ 225, 100, 36 ], x0:1385, x_i:[ 1, 8, 10 ], x0_mod:485, n:900 }

$x\equiv 4\!\pmod{5},\quad x\equiv 1\!\pmod{12},\quad x\equiv 7\!\pmod{14}$     rješenje: $x=49$

In [132]:
KTO_rucno([4,1,7],[5,12,14])
Out[132]:
Error: moduli nisu u parovima relativno prosti.

Eulerova funkcija

In [133]:
euler_phi(17)
Out[133]:
16
In [134]:
euler_phi(28)
Out[134]:
12
In [135]:
[euler_phi(n) for n in xrange(1,51)]
Out[135]:
[ 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, 12, 10, 22,

 8, 20, 12, 18, 12, 28, 8, 30, 16, 20, 16, 24, 12, 36, 18, 24, 16, 40, 12, 42,

 20, 24, 22, 46, 16, 42, 20 ]
In [136]:
list_plot([euler_phi(n) for n in xrange(1,51)],xmax=55,figsize=[5,4])
Out[136]:

Kako dobiti listu svih pozitivnih brojeva manjih od 100 koji su relativno prosti sa 100?

In [137]:
filter(lambda n: gcd(100,n)==1,range(100))
Out[137]:
[ 1, 3, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49,

 51, 53, 57, 59, 61, 63, 67, 69, 71, 73, 77, 79, 81, 83, 87, 89, 91, 93, 97,

 99 ]
In [138]:
def reducirani_sustav_ostataka(m):
    ostaci=filter(lambda n: gcd(m,n)==1,range(m))
    return ostaci
In [139]:
reducirani_sustav_ostataka(100)
Out[139]:
[ 1, 3, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49,

 51, 53, 57, 59, 61, 63, 67, 69, 71, 73, 77, 79, 81, 83, 87, 89, 91, 93, 97,

 99 ]
In [140]:
reducirani_sustav_ostataka(1000)
Out[140]:
[ 1, 3, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49,

 51, 53, 57, 59, 61, 63, 67, 69, 71, 73, 77, 79, 81, 83, 87, 89, 91, 93, 97,

 99, 101, 103, 107, 109, 111, 113, 117, 119, 121, 123, 127, 129, 131, 133, 137,

 139, 141, 143, 147, 149, 151, 153, 157, 159, 161, 163, 167, 169, 171, 173,

 177, 179, 181, 183, 187, 189, 191, 193, 197, 199, 201, 203, 207, 209, 211,

 213, 217, 219, 221, 223, 227, 229, 231, 233, 237, 239, 241, 243, 247, 249,

 251, 253, 257, 259, 261, 263, 267, 269, 271, 273, 277, 279, 281, 283, 287,

 289, 291, 293, 297, 299, 301, 303, 307, 309, 311, 313, 317, 319, 321, 323,

 327, 329, 331, 333, 337, 339, 341, 343, 347, 349, 351, 353, 357, 359, 361,

 363, 367, 369, 371, 373, 377, 379, 381, 383, 387, 389, 391, 393, 397, 399,

 401, 403, 407, 409, 411, 413, 417, 419, 421, 423, 427, 429, 431, 433, 437,

 439, 441, 443, 447, 449, 451, 453, 457, 459, 461, 463, 467, 469, 471, 473,

 477, 479, 481, 483, 487, 489, 491, 493, 497, 499, 501, 503, 507, 509, 511,

 513, 517, 519, 521, 523, 527, 529, 531, 533, 537, 539, 541, 543, 547, 549,

 551, 553, 557, 559, 561, 563, 567, 569, 571, 573, 577, 579, 581, 583, 587,

 589, 591, 593, 597, 599, 601, 603, 607, 609, 611, 613, 617, 619, 621, 623,

 627, 629, 631, 633, 637, 639, 641, 643, 647, 649, 651, 653, 657, 659, 661,

 663, 667, 669, 671, 673, 677, 679, 681, 683, 687, 689, 691, 693, 697, 699,

 701, 703, 707, 709, 711, 713, 717, 719, 721, 723, 727, 729, 731, 733, 737,

 739, 741, 743, 747, 749, 751, 753, 757, 759, 761, 763, 767, 769, 771, 773,

 777, 779, 781, 783, 787, 789, 791, 793, 797, 799, 801, 803, 807, 809, 811,

 813, 817, 819, 821, 823, 827, 829, 831, 833, 837, 839, 841, 843, 847, 849,

 851, 853, 857, 859, 861, 863, 867, 869, 871, 873, 877, 879, 881, 883, 887,

 889, 891, 893, 897, 899, 901, 903, 907, 909, 911, 913, 917, 919, 921, 923,

 927, 929, 931, 933, 937, 939, 941, 943, 947, 949, 951, 953, 957, 959, 961,

 963, 967, 969, 971, 973, 977, 979, 981, 983, 987, 989, 991, 993, 997, 999 ]

Posljednje dvije znamenke broja $3^{502}\cdot 7^{201}$

In [141]:
3^502*7^201%100
Out[141]:
63

Za velike potencije naredba power_mod puno brže računa ostatak

In [142]:
time 45^33333333%342
CPU times: user 1.15 s, sys: 59.6 ms, total: 1.21 s
Wall time: 1.22 s
Out[142]:
153
In [143]:
time power_mod(45,33333333,342)
CPU times: user 102 µs, sys: 0 ns, total: 102 µs
Wall time: 93.9 µs
Out[143]:
153

Sustavi kongruencija

In [144]:
var('x,y')
Out[144]:
( x, y )

$x+5y\equiv 3\!\pmod{8},\quad 4x+5y\equiv 1\!\pmod{8}$      rješenje: $(2,5)$

In [145]:
solve_mod([x+5*y==3,4*x+5*y==1],8)
Out[145]:
[ ( 2, 5 ) ]

$3x+5y\equiv 14\!\pmod{28},\quad 5x+9y\equiv 6\!\pmod{28}$     rješenja: $(20,2),\,\,(6,16)$

In [146]:
solve_mod([3*x+5*y==14,5*x+9*y==6],28)
Out[146]:
[ ( 20, 2 ), ( 6, 16 ) ]

$2x+3y\equiv 6\!\pmod{7},\quad 3x+y\equiv 2\!\pmod{7}$     rješenja: $(0,2),\,\,(1,6),\,\,(2,3),\,\,(3,0),\,\,(4,4),\,\,(5,1),\,\,(6,5)$

In [147]:
solve_mod([2*x+3*y==6,3*x+y==2],7)
Out[147]:
[ ( 0, 2 ), ( 1, 6 ), ( 2, 3 ), ( 3, 0 ), ( 4, 4 ), ( 5, 1 ), ( 6, 5 ) ]

$2x+3y\equiv 6\!\pmod{5\cdot6\cdot7\cdot11},\quad 3x+y\equiv 2\!\pmod{5\cdot6\cdot7\cdot11}$

In [148]:
solve_mod([2*x+3*y==6,3*x+y==2],5*6*7*11)
Out[148]:
[ ( 0, 2 ), ( 330, 1322 ), ( 660, 332 ), ( 990, 1652 ), ( 1320, 662 ), ( 1650,

 1982 ), ( 1980, 992 ) ]

$3x+5y\equiv 14\!\pmod{100},\quad 5x+9y\equiv 6\!\pmod{100}$     rješenja: $(48,74),\,\,(98,24)$

In [149]:
solve_mod([3*x+5*y==14,5*x+9*y==6],100)
Out[149]:
[ ( 48, 74 ), ( 98, 24 ) ]

$79x\equiv 15\!\pmod{722}$     rješenje: $119$

In [150]:
solve_mod([79*x==15],722)
Out[150]:
[ ( 119 ) ]

Red od a modulo n

red od 2 modulo 7 je jednak 3

In [151]:
Mod(2,7).multiplicative_order()
Out[151]:
3
In [152]:
type(Mod(2,7))
Out[152]:
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>
In [153]:
type(mod(2,7))
Out[153]:
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>

red od 16 modulo 29 je jednak 7

In [154]:
Mod(16,29).multiplicative_order()
Out[154]:
7

red od 12 modulo 28 ne postoji

In [155]:
Mod(12,28).multiplicative_order()
---------------------------------------------------------------------------
ArithmeticError                           Traceback (most recent call last)
<ipython-input-155-5dcad6c834d0> in <module>()
----> 1 Mod(Integer(12),Integer(28)).multiplicative_order()

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/rings/finite_rings/integer_mod.pyx in sage.rings.finite_rings.integer_mod.IntegerMod_abstract.multiplicative_order (build/cythonized/sage/rings/finite_rings/integer_mod.c:22196)()
   1762             return sage.rings.integer.Integer(self.__pari__().znorder())
   1763         except PariError:
-> 1764             raise ArithmeticError("multiplicative order of %s not defined since it is not a unit modulo %s"%(
   1765                 self, self.__modulus.sageInteger))
   1766 

ArithmeticError: multiplicative order of 12 not defined since it is not a unit modulo 28
In [156]:
Mod(2,7) in IntegerModRing(7)
Out[156]:
True
In [157]:
mod(2,7) in IntegerModRing(7)
Out[157]:
True

Primitivni korijen modulo n

primitivni korijen modulo 7

In [158]:
primitive_root(7)
Out[158]:
3

svi primitivni korijeni modulo 7

In [159]:
filter(lambda x: Mod(x,7).multiplicative_order()==euler_phi(7), range(1,7))
Out[159]:
[ 3, 5 ]

svi primitivni korijeni modulo 41

In [160]:
filter(lambda x: Mod(x,41).multiplicative_order()==euler_phi(41), range(1,41))
Out[160]:
[ 6, 7, 11, 12, 13, 15, 17, 19, 22, 24, 26, 28, 29, 30, 34, 35 ]

ne postoji primitivni korijen modulo 12

In [161]:
filter(lambda x: gcd(x,12)==1 and Mod(x,12).multiplicative_order()==euler_phi(12), range(1,12))
Out[161]:
[  ]
In [162]:
primitive_root(12)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-162-dfd011334dbe> in <module>()
----> 1 primitive_root(Integer(12))

/home/damir/programi/sage-8.4/local/lib/python2.7/site-packages/sage/arith/misc.pyc in primitive_root(n, check)
   4027         if m%2 and m.is_prime_power():
   4028             return ZZ(pari(n).znprimroot())
-> 4029     raise ValueError("no primitive root")
   4030 
   4031 def nth_prime(n):

ValueError: no primitive root

primitivni korijen modulo n generira reducirani sustav ostataka modulo n

$n=41$

In [163]:
reducirani_sustav_ostataka(41)
Out[163]:
[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,

 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40 ]
In [164]:
[mod(6^k,41) for k in xrange(euler_phi(41))]
Out[164]:
[ 1, 6, 36, 11, 25, 27, 39, 29, 10, 19, 32, 28, 4, 24, 21, 3, 18, 26, 33, 34,

 40, 35, 5, 30, 16, 14, 2, 12, 31, 22, 9, 13, 37, 17, 20, 38, 23, 15, 8, 7 ]
In [165]:
sorted([mod(6^k,41) for k in xrange(euler_phi(41))])
Out[165]:
[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,

 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40 ]

$n=50$

In [166]:
primitive_root(50)
Out[166]:
27
In [167]:
reducirani_sustav_ostataka(50)
Out[167]:
[ 1, 3, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49 ]
In [168]:
sorted([mod(27^k,50) for k in xrange(euler_phi(50))])
Out[168]:
[ 1, 3, 7, 9, 11, 13, 17, 19, 21, 23, 27, 29, 31, 33, 37, 39, 41, 43, 47, 49 ]

Nelinearne kongruencije

$x^5\equiv 2\!\pmod{7}$

In [169]:
solve_mod([x^5==2],7)
Out[169]:
[ ( 4 ) ]

$x^{12}\equiv 37\!\pmod{41}$

In [170]:
solve_mod([x^12==37],41)
Out[170]:
[ ( 2 ), ( 18 ), ( 23 ), ( 39 ) ]

$7^x\equiv 6\!\pmod{17}$

In [171]:
filter(lambda x: mod(7^x,17)==6, range(17))
Out[171]:
[ 13 ]

$3^x\equiv 2\!\pmod{23}$

In [172]:
filter(lambda x: mod(3^x,23)==2, range(23))
Out[172]:
[ 7, 18 ]

$2x^8\equiv 5\!\pmod{13}$

In [173]:
solve_mod([2*x^8==5],13)
Out[173]:
[ ( 2 ), ( 3 ), ( 10 ), ( 11 ) ]

$x^5\equiv 2\!\pmod{12}$     nema rješenja

In [174]:
solve_mod([x^5==2],12)
Out[174]:
[  ]

$x^5\equiv 4\!\pmod{12}$

In [175]:
solve_mod([x^5==4],12)
Out[175]:
[ ( 4 ), ( 10 ) ]
In [176]:
filter(lambda x: mod(x^5,12)==4, range(12))
Out[176]:
[ 4, 10 ]

$x^7\equiv 17\!\pmod{322122548}$

In [177]:
time solve_mod([x^7==17],322122548)
CPU times: user 4.69 s, sys: 63.1 ms, total: 4.76 s
Wall time: 4.74 s
Out[177]:
[ ( 200490385 ) ]
In [178]:
%display latex
In [179]:
factor(322122548)
Out[179]:

Još jedna napomena

Ovdje želimo naglasiti razliku između naredbi mod (ili Mod) i naredbe % (za dobivanje ostatka pri dijeljenju dva cijela broja). Na prvi pogled se čini da te naredbe rade istu stvar, ali nije baš tako. Doduše, na izlazu daju kao isti rezultate, ali te rezultate SAGE potpuno drukčije doživljava. Pogledajmo na jednom primjeru tu razliku.

In [180]:
%display ascii_art
In [181]:
37%14
Out[181]:
9
In [182]:
type(37%14)
Out[182]:
<type 'sage.rings.integer.Integer'>

Kako na 37%14 SAGE gleda kao na element iz prstena $\mathbb{Z}$, zbrajanje i množenje se odvija u prstenu $\mathbb{Z}$ i dobivamo očekivani rezultat.

In [183]:
37%14+5*14
Out[183]:
79

Da li je tako i s naredbom mod?

In [184]:
mod(37,14)
Out[184]:
9
In [185]:
type(mod(37,14))
Out[185]:
<type 'sage.rings.finite_rings.integer_mod.IntegerMod_int'>

Kako na mod(37,14) SAGE gleda kao na element iz prstena $\mathbb{Z}_{14}$, tada se zbrajanje i množenje odvijaju u prstenu $\mathbb{Z}_{14}$. Stoga rezultat kojeg je vratio SAGE nije bug, već je to stvarno dobar rezultat. Naime, u prstenu $\mathbb{Z}_{14}$ svi višekratnici broja 14 su zapravo jednaki 0, pa je tako i $5\cdot14$ također jednako 0 jer se operacije zbrajanja i množenja u prstenu $\mathbb{Z}_{14}$ rade modulo 14.

In [186]:
mod(37,14)+5*14
Out[186]:
9
In [ ]: