Line data Source code
1 : #line 2 "../src/kernel/none/invmod.c"
2 : /* Copyright (C) 2003 The PARI group.
3 :
4 : This file is part of the PARI/GP package.
5 :
6 : PARI/GP is free software; you can redistribute it and/or modify it under the
7 : terms of the GNU General Public License as published by the Free Software
8 : Foundation; either version 2 of the License, or (at your option) any later
9 : version. It is distributed in the hope that it will be useful, but WITHOUT
10 : ANY WARRANTY WHATSOEVER.
11 :
12 : Check the License for details. You should have received a copy of it, along
13 : with the package; see the file 'COPYING'. If not, write to the Free Software
14 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
15 :
16 : /*==================================
17 : * invmod(a,b,res)
18 : *==================================
19 : * If a is invertible, return 1, and set res = a^{ -1 }
20 : * Otherwise, return 0, and set res = gcd(a,b)
21 : *
22 : * This is sufficiently different from bezout() to be implemented separately
23 : * instead of having a bunch of extra conditionals in a single function body
24 : * to meet both purposes.
25 : */
26 :
27 : #ifdef INVMOD_PARI
28 : INLINE int
29 12589432 : invmod_pari(GEN a, GEN b, GEN *res)
30 : #else
31 : int
32 10209444 : invmod(GEN a, GEN b, GEN *res)
33 : #endif
34 : {
35 : GEN v,v1,d,d1,q,r;
36 : pari_sp av, av1;
37 : long s;
38 : ulong g;
39 : ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
40 : int lhmres; /* Lehmer stage return value */
41 :
42 22798876 : if (!signe(b)) { *res=absi(a); return 0; }
43 22798876 : av = avma;
44 22798876 : if (lgefint(b) == 3) /* single-word affair */
45 : {
46 17356843 : ulong d1 = umodiu(a, uel(b,2));
47 17356843 : if (d1 == 0)
48 : {
49 299 : if (b[2] == 1L)
50 299 : { *res = gen_0; return 1; }
51 : else
52 0 : { *res = absi(b); return 0; }
53 : }
54 17356544 : g = xgcduu(uel(b,2), d1, 1, &xv, &xv1, &s);
55 17356544 : set_avma(av);
56 17356544 : if (g != 1UL) { *res = utoipos(g); return 0; }
57 17356531 : xv = xv1 % uel(b,2); if (s < 0) xv = uel(b,2) - xv;
58 17356531 : *res = utoipos(xv); return 1;
59 : }
60 :
61 5442033 : (void)new_chunk(lgefint(b));
62 5442033 : d = absi_shallow(b); d1 = modii(a,d);
63 :
64 5442033 : v=gen_0; v1=gen_1; /* general case */
65 5442033 : av1 = avma;
66 :
67 14667704 : while (lgefint(d) > 3 && signe(d1))
68 : {
69 9225671 : lhmres = lgcdii((ulong*)d, (ulong*)d1, &xu, &xu1, &xv, &xv1, ULONG_MAX);
70 9225671 : if (lhmres != 0) /* check progress */
71 : { /* apply matrix */
72 8316309 : if (lhmres == 1 || lhmres == -1)
73 : {
74 193540 : if (xv1 == 1)
75 : {
76 177581 : r = subii(d,d1); d=d1; d1=r;
77 177581 : a = subii(v,v1); v=v1; v1=a;
78 : }
79 : else
80 : {
81 15959 : r = subii(d, mului(xv1,d1)); d=d1; d1=r;
82 15959 : a = subii(v, mului(xv1,v1)); v=v1; v1=a;
83 : }
84 : }
85 : else
86 : {
87 8122769 : r = subii(muliu(d,xu), muliu(d1,xv));
88 8122769 : a = subii(muliu(v,xu), muliu(v1,xv));
89 8122769 : d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
90 8122769 : v1 = subii(muliu(v,xu1), muliu(v1,xv1)); v = a;
91 8122769 : if (lhmres&1) { togglesign(d); togglesign(v); }
92 4204489 : else { togglesign(d1); togglesign(v1); }
93 : }
94 : }
95 :
96 9225671 : if (lhmres <= 0 && signe(d1))
97 : {
98 1128461 : q = dvmdii(d,d1,&r);
99 1128461 : a = subii(v,mulii(q,v1));
100 1128461 : v=v1; v1=a;
101 1128461 : d=d1; d1=r;
102 : }
103 9225671 : if (gc_needed(av,1))
104 : {
105 885 : if(DEBUGMEM>1) pari_warn(warnmem,"invmod");
106 885 : (void)gc_all(av1, 4, &d,&d1,&v,&v1);
107 : }
108 : } /* end while */
109 :
110 : /* Postprocessing - final sprint */
111 5442033 : if (signe(d1))
112 : {
113 : /* Assertions: lgefint(d)==lgefint(d1)==3, and
114 : * gcd(d,d1) is nonzero and fits into one word
115 : */
116 5369005 : g = xxgcduu(uel(d,2), uel(d1,2), 1, &xu, &xu1, &xv, &xv1, &s);
117 5369005 : if (g != 1UL) { set_avma(av); *res = utoipos(g); return 0; }
118 : /* (From the xgcduu() blurb:)
119 : * For finishing the multiword modinv, we now have to multiply the
120 : * returned matrix (with properly adjusted signs) onto the values
121 : * v' and v1' previously obtained from the multiword division steps.
122 : * Actually, it is sufficient to take the scalar product of [v',v1']
123 : * with [u1,-v1], and change the sign if s==1.
124 : */
125 5368876 : v = subii(muliu(v,xu1),muliu(v1,xv1));
126 5368876 : if (s > 0) setsigne(v,-signe(v));
127 5368876 : set_avma(av); *res = modii(v,b); return 1;
128 : }
129 : /* get here when the final sprint was skipped (d1 was zero already) */
130 73028 : set_avma(av);
131 73028 : if (!equalii(d,gen_1)) { *res = icopy(d); return 0; }
132 72927 : *res = modii(v,b); return 1;
133 : }
|
