1 // -*- mode:C++ ; compile-command: "g++-3.4 -I.. -g -c ifactor.cc -DHAVE_CONFIG_H -DIN_GIAC" -*- 2 #include "giacPCH.h" 3 #if !defined __MINGW_H && !defined KHICAS 4 #define GIAC_MPQS // define if you want to use giac for sieving 5 #endif 6 7 #ifdef HAVE_LIBECM 8 #include <ecm.h> 9 #endif 10 11 #include "path.h" 12 /* 13 * Copyright (C) 2003,14 R. De Graeve & B. Parisse, 14 * Institut Fourier, 38402 St Martin d'Heres 15 * 16 * This program is free software; you can redistribute it and/or modify 17 * it under the terms of the GNU General Public License as published by 18 * the Free Software Foundation; either version 3 of the License, or 19 * (at your option) any later version. 20 * 21 * This program is distributed in the hope that it will be useful, 22 * but WITHOUT ANY WARRANTY; without even the implied warranty of 23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 24 * GNU General Public License for more details. 25 * 26 * You should have received a copy of the GNU General Public License 27 * along with this program. If not, see <http://www.gnu.org/licenses/>. 28 */ 29 30 using namespace std; 31 #ifdef GIAC_HAS_STO_38 32 //#undef clock 33 //#undef clock_t 34 #else 35 #include <fstream> 36 //#include <unistd.h> // For reading arguments from file 37 #include "ifactor.h" 38 #include "pari.h" 39 #include "usual.h" 40 #include "sym2poly.h" 41 #include "rpn.h" 42 #include "prog.h" 43 #include "misc.h" 44 #include "giacintl.h" 45 #endif 46 47 #ifdef GIAC_HAS_STO_38 48 #define BESTA_OS 49 #endif 50 // Trying to make ifactor(2^128+1) work on ARM 51 #if defined(RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 52 //#define OLD_AFACT 53 #define GIAC_ADDITIONAL_PRIMES 16// if defined, additional primes are used in sieve 54 #else 55 #define GIAC_ADDITIONAL_PRIMES 32// if defined, additional primes are used in sieve 56 #endif 57 58 59 #ifndef NO_NAMESPACE_GIAC 60 namespace giac { 61 #endif // ndef NO_NAMESPACE_GIAC 62 63 #if 0 64 struct int256 { 65 longlong a; 66 ulonglong b,c,d; 67 }; 68 struct int128 { 69 longlong a; 70 ulonglong b; 71 }; 72 73 void sub(int256 A,int b,int256 & C){ 74 bool Apos=A.a>=0; 75 if (Apos){ 76 if (b<0){ 77 add(A,b,C); 78 return; 79 } 80 bool carry=A.d<b; 81 A.d -= b; 82 if (carry){ 83 carry=A.c<1; 84 A.c--; 85 if (carry){ 86 carry=A.b<1; 87 A.b--; 88 if (carry) 89 A.a--; 90 } 91 } 92 C=A; 93 return; 94 } 95 A.a = -A.a; 96 sub(A,-b,C); 97 C.a = -C.a; 98 } 99 100 void add(const int192 & A,int b,int192 & C){ 101 bool Apos=A.a>=0; 102 if (Apos){ 103 if (b<0){ 104 sub(A,-b,C); 105 return; 106 } 107 A.d += b; 108 if (A.d<b){ 109 A.c++; 110 if (A.c==0){ 111 A.b++; 112 if (A.n==0) 113 A.a++; 114 } 115 } 116 C=A; 117 return; 118 } 119 A.a = -A.a; 120 sub(A,-b,C); 121 C.a = -C.a; 122 } 123 124 inline void int128tounsigned(const int128 & A,unsigned & Aa,unsigned & Ab,unsigned & Ac,unsigned & Ad){ 125 Aa = (* ((unsigned *) (&A.a)+1)) & 0x7fff; 126 Ab = *(unsigned *) A.a; 127 Ac = *( (unsigned *) &A.b)+1; 128 Ad = *(unsigned *) A.d; 129 } 130 131 void mul(int128 A,int128 B,int256 & C){ 132 bool neg=(A.a<0) ^ (B.a<0); 133 unsigned A3,A2,A1,A0,B3,B2,B1,B0,C6,C5,C4,C3,C2,C1,C0; 134 int128tounsigned(A,A3,A2,A1,A0); 135 int128tounsigned(B,B3,B2,B1,B0); 136 ulonglong p1,p2,p3,p; 137 p1=A0; p2=B0; 138 p3=p1*p2; 139 C0=p3; 140 C1=p3>>32; 141 p1=A1; 142 p = p1*p2; 143 p += C1; 144 unsigned short carry = (p<C1); 145 p1=A0; p2=B1; 146 p3=p1*p2; 147 p += p3; 148 carry += (p<p3); 149 C1 = p; 150 C2 = p>>32; 151 C2 += carry; 152 carry = (C2<carry); 153 } 154 155 #endif 156 157 const short int giac_primes[]={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,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063,1069,1087,1091,1093,1097,1103,1109,1117,1123,1129,1151,1153,1163,1171,1181,1187,1193,1201,1213,1217,1223,1229,1231,1237,1249,1259,1277,1279,1283,1289,1291,1297,1301,1303,1307,1319,1321,1327,1361,1367,1373,1381,1399,1409,1423,1427,1429,1433,1439,1447,1451,1453,1459,1471,1481,1483,1487,1489,1493,1499,1511,1523,1531,1543,1549,1553,1559,1567,1571,1579,1583,1597,1601,1607,1609,1613,1619,1621,1627,1637,1657,1663,1667,1669,1693,1697,1699,1709,1721,1723,1733,1741,1747,1753,1759,1777,1783,1787,1789,1801,1811,1823,1831,1847,1861,1867,1871,1873,1877,1879,1889,1901,1907,1913,1931,1933,1949,1951,1973,1979,1987,1993,1997,1999,2003,2011,2017,2027,2029,2039,2053,2063,2069,2081,2083,2087,2089,2099,2111,2113,2129,2131,2137,2141,2143,2153,2161,2179,2203,2207,2213,2221,2237,2239,2243,2251,2267,2269,2273,2281,2287,2293,2297,2309,2311,2333,2339,2341,2347,2351,2357,2371,2377,2381,2383,2389,2393,2399,2411,2417,2423,2437,2441,2447,2459,2467,2473,2477,2503,2521,2531,2539,2543,2549,2551,2557,2579,2591,2593,2609,2617,2621,2633,2647,2657,2659,2663,2671,2677,2683,2687,2689,2693,2699,2707,2711,2713,2719,2729,2731,2741,2749,2753,2767,2777,2789,2791,2797,2801,2803,2819,2833,2837,2843,2851,2857,2861,2879,2887,2897,2903,2909,2917,2927,2939,2953,2957,2963,2969,2971,2999,3001,3011,3019,3023,3037,3041,3049,3061,3067,3079,3083,3089,3109,3119,3121,3137,3163,3167,3169,3181,3187,3191,3203,3209,3217,3221,3229,3251,3253,3257,3259,3271,3299,3301,3307,3313,3319,3323,3329,3331,3343,3347,3359,3361,3371,3373,3389,3391,3407,3413,3433,3449,3457,3461,3463,3467,3469,3491,3499,3511,3517,3527,3529,3533,3539,3541,3547,3557,3559,3571,3581,3583,3593,3607,3613,3617,3623,3631,3637,3643,3659,3671,3673,3677,3691,3697,3701,3709,3719,3727,3733,3739,3761,3767,3769,3779,3793,3797,3803,3821,3823,3833,3847,3851,3853,3863,3877,3881,3889,3907,3911,3917,3919,3923,3929,3931,3943,3947,3967,3989,4001,4003,4007,4013,4019,4021,4027,4049,4051,4057,4073,4079,4091,4093,4099,4111,4127,4129,4133,4139,4153,4157,4159,4177,4201,4211,4217,4219,4229,4231,4241,4243,4253,4259,4261,4271,4273,4283,4289,4297,4327,4337,4339,4349,4357,4363,4373,4391,4397,4409,4421,4423,4441,4447,4451,4457,4463,4481,4483,4493,4507,4513,4517,4519,4523,4547,4549,4561,4567,4583,4591,4597,4603,4621,4637,4639,4643,4649,4651,4657,4663,4673,4679,4691,4703,4721,4723,4729,4733,4751,4759,4783,4787,4789,4793,4799,4801,4813,4817,4831,4861,4871,4877,4889,4903,4909,4919,4931,4933,4937,4943,4951,4957,4967,4969,4973,4987,4993,4999,5003,5009,5011,5021,5023,5039,5051,5059,5077,5081,5087,5099,5101,5107,5113,5119,5147,5153,5167,5171,5179,5189,5197,5209,5227,5231,5233,5237,5261,5273,5279,5281,5297,5303,5309,5323,5333,5347,5351,5381,5387,5393,5399,5407,5413,5417,5419,5431,5437,5441,5443,5449,5471,5477,5479,5483,5501,5503,5507,5519,5521,5527,5531,5557,5563,5569,5573,5581,5591,5623,5639,5641,5647,5651,5653,5657,5659,5669,5683,5689,5693,5701,5711,5717,5737,5741,5743,5749,5779,5783,5791,5801,5807,5813,5821,5827,5839,5843,5849,5851,5857,5861,5867,5869,5879,5881,5897,5903,5923,5927,5939,5953,5981,5987,6007,6011,6029,6037,6043,6047,6053,6067,6073,6079,6089,6091,6101,6113,6121,6131,6133,6143,6151,6163,6173,6197,6199,6203,6211,6217,6221,6229,6247,6257,6263,6269,6271,6277,6287,6299,6301,6311,6317,6323,6329,6337,6343,6353,6359,6361,6367,6373,6379,6389,6397,6421,6427,6449,6451,6469,6473,6481,6491,6521,6529,6547,6551,6553,6563,6569,6571,6577,6581,6599,6607,6619,6637,6653,6659,6661,6673,6679,6689,6691,6701,6703,6709,6719,6733,6737,6761,6763,6779,6781,6791,6793,6803,6823,6827,6829,6833,6841,6857,6863,6869,6871,6883,6899,6907,6911,6917,6947,6949,6959,6961,6967,6971,6977,6983,6991,6997,7001,7013,7019,7027,7039,7043,7057,7069,7079,7103,7109,7121,7127,7129,7151,7159,7177,7187,7193,7207,7211,7213,7219,7229,7237,7243,7247,7253,7283,7297,7307,7309,7321,7331,7333,7349,7351,7369,7393,7411,7417,7433,7451,7457,7459,7477,7481,7487,7489,7499,7507,7517,7523,7529,7537,7541,7547,7549,7559,7561,7573,7577,7583,7589,7591,7603,7607,7621,7639,7643,7649,7669,7673,7681,7687,7691,7699,7703,7717,7723,7727,7741,7753,7757,7759,7789,7793,7817,7823,7829,7841,7853,7867,7873,7877,7879,7883,7901,7907,7919,7927,7933,7937,7949,7951,7963,7993,8009,8011,8017,8039,8053,8059,8069,8081,8087,8089,8093,8101,8111,8117,8123,8147,8161,8167,8171,8179,8191,8209,8219,8221,8231,8233,8237,8243,8263,8269,8273,8287,8291,8293,8297,8311,8317,8329,8353,8363,8369,8377,8387,8389,8419,8423,8429,8431,8443,8447,8461,8467,8501,8513,8521,8527,8537,8539,8543,8563,8573,8581,8597,8599,8609,8623,8627,8629,8641,8647,8663,8669,8677,8681,8689,8693,8699,8707,8713,8719,8731,8737,8741,8747,8753,8761,8779,8783,8803,8807,8819,8821,8831,8837,8839,8849,8861,8863,8867,8887,8893,8923,8929,8933,8941,8951,8963,8969,8971,8999,9001,9007,9011,9013,9029,9041,9043,9049,9059,9067,9091,9103,9109,9127,9133,9137,9151,9157,9161,9173,9181,9187,9199,9203,9209,9221,9227,9239,9241,9257,9277,9281,9283,9293,9311,9319,9323,9337,9341,9343,9349,9371,9377,9391,9397,9403,9413,9419,9421,9431,9433,9437,9439,9461,9463,9467,9473,9479,9491,9497,9511,9521,9533,9539,9547,9551,9587,9601,9613,9619,9623,9629,9631,9643,9649,9661,9677,9679,9689,9697,9719,9721,9733,9739,9743,9749,9767,9769,9781,9787,9791,9803,9811,9817,9829,9833,9839,9851,9857,9859,9871,9883,9887,9901,9907,9923,9929,9931,9941,9949,9967,9973}; 158 159 // #define PRIME_SIEVE 160 #ifdef PRIME_SIEVE 161 // fill Erathosthene sieve crible for searching primes up to 2*crible.size()*32+1 162 // crible is a (packed) bit array, crible[i] is true if 2*i+1 is a prime 163 // crible must be set to true at startup fill_crible(vector<unsigned> & crible,unsigned p)164 void fill_crible(vector<unsigned> & crible,unsigned p){ 165 crible.resize((p-1)/64+1); 166 unsigned cs=crible.size(); 167 unsigned lastnum=64*cs; 168 unsigned lastsieve=int(std::sqrt(double(lastnum))); 169 unsigned primesieved=1; 170 crible[0] = 0xfffffffe; // 1 is not prime and not sieved (2 is not sieved) 171 for (unsigned i=1;i<cs;++i) 172 crible[i]=0xffffffff; 173 for (;primesieved<=lastsieve;primesieved+=2){ 174 // find next prime 175 unsigned pos=primesieved/2; 176 for (;pos<cs;pos++){ 177 if (crible[pos/32] & (1 << (pos %32))) 178 break; 179 } 180 // set mutiples of (2*pos+1) to false 181 primesieved=2*pos+1; 182 unsigned n=3*primesieved; 183 for (;n<lastnum;n+=2*primesieved){ 184 pos=(n-1)/2; 185 crible[(pos/32)] &= ~(1<<(pos %32)); 186 } 187 } 188 } nextprime(vector<unsigned> & crible,unsigned p)189 unsigned nextprime(vector<unsigned> & crible,unsigned p){ 190 // assumes crible has been filled 191 ++p; 192 if (p%2==0) 193 ++p; 194 unsigned pos=(p-1)/2,cs=crible.size()*32; 195 if (2*cs+1<=p) 196 return nextprime(int(p)).val; 197 for (;pos<cs;++pos){ 198 if (crible[pos/32] & (1<<(pos%32))){ 199 pos=2*pos+1; 200 // if (pos!=nextprime(int(p)).val) CERR << "error " << p << '\n'; 201 return pos; 202 } 203 } 204 return nextprime(int(p)).val; // not found 205 } 206 #endif 207 208 static const int giac_last_prime=giac_primes[sizeof(giac_primes)/sizeof(short)-1]; 209 #if defined RTOS_THREADX || defined BESTA_OS || defined NSPIRE 210 const unsigned QS_SIZE=65536; // number of slicetype in a sieve slice 211 typedef unsigned char slicetype; // define to unsigned char if not enough 212 #else 213 const unsigned QS_SIZE=65536; // number of slicetype in a sieve slice 214 typedef unsigned char slicetype; // define to unsigned char if not enough 215 #endif 216 217 #ifdef NSPIRE 218 template<class T> printbool(nio::ios_base<T> & os,const vector<unsigned> & v,int C=1)219 static void printbool(nio::ios_base<T> & os,const vector<unsigned> & v,int C=1){ 220 if (C) 221 C=giacmin(C,int(v.size())); 222 else 223 C=v.size(); 224 for (int c=0;c<C;++c){ 225 for (int s=0;s<32;++s){ 226 os << (((v[c] >> s & 1)==1)?1:0) << " "; 227 } 228 } 229 os << '\n'; 230 } 231 232 template<class T> printbool(nio::ios_base<T> & os,const vector<vector<unsigned>> & m,int L=32)233 void printbool(nio::ios_base<T> & os,const vector< vector<unsigned> > & m,int L=32){ 234 if (L) 235 L=giacmin(L,int(m.size())); 236 else 237 L=m.size(); 238 for (int l=0;l<L;++l){ 239 printbool(os,m[l]); 240 } 241 } 242 #else printbool(ostream & os,const vector<unsigned> & v,int C=1)243 static void printbool(ostream & os,const vector<unsigned> & v,int C=1){ 244 if (C) 245 C=giacmin(C,int(v.size())); 246 else 247 C=int(v.size()); 248 for (int c=0;c<C;++c){ 249 for (int s=0;s<32;++s){ 250 os << (((v[c] >> s & 1)==1)?1:0) << " "; 251 } 252 } 253 os << '\n'; 254 } 255 printbool(ostream & os,const vector<vector<unsigned>> & m,int L=32)256 void printbool(ostream & os,const vector< vector<unsigned> > & m,int L=32){ 257 if (L) 258 L=giacmin(L,int(m.size())); 259 else 260 L=int(m.size()); 261 for (int l=0;l<L;++l){ 262 printbool(os,m[l]); 263 } 264 } 265 #endif 266 267 template <class T> swap(T * & ptr1,T * & ptr2)268 inline void swap(T * & ptr1, T * & ptr2){ 269 T * tmp=ptr1; 270 ptr1=ptr2; 271 ptr2=tmp; 272 } 273 274 #ifdef x86_64 275 #define GIAC_RREF_UNROLL 4 276 #else 277 #define GIAC_RREF_UNROLL 4 278 #endif 279 280 // #define RREF_SORT 281 #ifdef RREF_SORT 282 struct line_t { 283 unsigned * tab; 284 unsigned count; 285 }; 286 operator <(const line_t & l1,const line_t & l2)287 bool operator < (const line_t & l1,const line_t & l2){ 288 if (!l1.count) 289 return false; 290 if (!l2.count) 291 return true; 292 return l1.count<l2.count; 293 } 294 count_ones(unsigned * tab,int C32)295 unsigned count_ones(unsigned * tab,int C32){ 296 register unsigned r=0; 297 register unsigned * tabend=tab+C32; 298 for (;tab!=tabend;++tab){ 299 register unsigned u=*tab; 300 while (u){ 301 r += u & 1; 302 u >>= 1; 303 } 304 } 305 return r; 306 } 307 308 #else 309 struct line_t { 310 unsigned * tab; 311 }; 312 #endif 313 314 315 // mode=0: full reduction, 1 subreduction, 2 finish full reduction from subreduction rref(vector<line_t> & m,int L,int C32,int mode)316 void rref(vector< line_t > & m,int L,int C32,int mode){ 317 int i,l=0,c=0,C=C32*32; 318 for (;l<L && c<C;){ 319 // printbool(CERR,m); 320 int c1=c/32,c2=c%32; 321 // find first non-0 pivot in col c starting at row l 322 for (i=l;i<L;++i){ 323 if ((m[i].tab[c1] >> c2) & 1) 324 break; 325 } 326 if (i==L){ // none found in this column 327 ++c; 328 continue; 329 } 330 if (i!=l) 331 swap(m[i].tab,m[l].tab); // don't care about count... 332 int start=mode==1?l+1:0, end=mode==2?l:L; 333 #ifdef x86_64 334 ulonglong * pivend, * pivbeg; 335 pivbeg = (ulonglong *) (m[l].tab+(c1/GIAC_RREF_UNROLL)*GIAC_RREF_UNROLL); 336 pivend = (ulonglong *) (m[l].tab+C32); 337 #else 338 unsigned * pivbeg = m[l].tab+(c1/GIAC_RREF_UNROLL)*GIAC_RREF_UNROLL, * pivend = m[l].tab+C32; 339 #endif 340 for (i=start;i<end;++i){ 341 if (i==l || ( (m[i].tab[c1] >> c2) & 1)!=1) 342 continue; 343 // line combination l and i 344 #ifdef x86_64 345 ulonglong * curptr=(ulonglong *) (m[i].tab+(c1/GIAC_RREF_UNROLL)*GIAC_RREF_UNROLL); 346 for (ulonglong * pivptr=pivbeg;pivptr!=pivend;curptr += GIAC_RREF_UNROLL/2,pivptr += GIAC_RREF_UNROLL/2){ 347 // small optimization (loop unroll), assumes mult of 4(*32) columns 348 // PREFETCH(curptr+8); 349 *curptr ^= *pivptr; 350 curptr[1] ^= pivptr[1]; 351 #if GIAC_RREF_UNROLL==8 352 curptr[2] ^= pivptr[2]; 353 curptr[3] ^= pivptr[3]; 354 #endif 355 } 356 #else 357 unsigned * curptr=m[i].tab+(c1/GIAC_RREF_UNROLL)*GIAC_RREF_UNROLL; 358 for (unsigned * pivptr=pivbeg;pivptr!=pivend;curptr += GIAC_RREF_UNROLL,pivptr += GIAC_RREF_UNROLL){ 359 // small optimization (loop unroll), assumes mult of 4(*32) columns 360 // PREFETCH(curptr+16); 361 *curptr ^= *pivptr; 362 curptr[1] ^= pivptr[1]; 363 curptr[2] ^= pivptr[2]; 364 curptr[3] ^= pivptr[3]; 365 #if GIAC_RREF_UNROLL==8 366 curptr[4] ^= pivptr[4]; 367 curptr[5] ^= pivptr[5]; 368 curptr[6] ^= pivptr[6]; 369 curptr[7] ^= pivptr[7]; 370 #endif 371 } 372 #endif 373 } 374 ++l; 375 ++c; 376 } 377 } 378 379 template <class T> release_memory(vector<T> & slice)380 void release_memory(vector<T> & slice){ 381 // release memory from slice 382 vector<T> tmp; 383 swap(slice,tmp); 384 } 385 386 #ifdef USE_GMP_REPLACEMENTS modulo(const mpz_t & a,unsigned b)387 int modulo(const mpz_t & a,unsigned b){ 388 if (mpz_cmp_ui(a,0)<0){ 389 mpz_neg(*(mpz_t *)&a,a); 390 int res=modulo(a,b); 391 mpz_neg(*(mpz_t *)&a,a); 392 return b-res; 393 } 394 mp_digit C; 395 mp_mod_d((mp_int *)&a,b,&C); 396 return C; 397 } 398 #else modulo(const mpz_t & a,unsigned b)399 int modulo(const mpz_t & a,unsigned b){ 400 return mpz_fdiv_ui(a,b); 401 } 402 #endif 403 404 #if defined RTOS_THREADX || defined BESTA_OS || defined NSPIRE 405 typedef unsigned short pui_t ; 406 typedef unsigned short ushort_t; 407 typedef short short_t; 408 #else 409 typedef unsigned pui_t ; 410 // #ifndef USE_GMP_REPLACEMENTS // uncomment for Aspen debugging 411 #define PRIMES32 412 // #endif 413 #ifdef PRIMES32 414 typedef unsigned ushort_t; 415 typedef int short_t; 416 #else 417 typedef unsigned short ushort_t; 418 typedef unsigned short int short_t; 419 #endif 420 421 #ifdef EMCC 422 #include <map> 423 #endif 424 #if (defined EMCC || defined(HASH_MAP_NAMESPACE)) && defined(PRIMES32) 425 #define ADDITIONAL_PRIMES_HASHMAP 426 #endif 427 #endif // RTOS_THREADX || BESTA_OS 428 429 struct axbinv { 430 #if 0 431 unsigned short aindex; 432 unsigned short bindex; 433 #else 434 unsigned aindex; 435 unsigned bindex; 436 #endif 437 int shiftpos; 438 pui_t first,second; // indexes in the "puissancestab" table axbinvgiac::axbinv439 axbinv(ushort_t a_,int shiftpos_,ushort_t b_,pui_t f_,pui_t s_):aindex(a_),bindex(b_),shiftpos(shiftpos_),first(f_),second(s_) {}; axbinvgiac::axbinv440 axbinv() {}; 441 }; 442 443 #ifdef ADDITIONAL_PRIMES_HASHMAP largep(const axbinv & A,ushort_t * puissancestab)444 unsigned largep(const axbinv & A,ushort_t * puissancestab) { 445 // return A.largeprime; 446 if (A.second-A.first<3) return 0; 447 #ifdef PRIMES32 448 if (*(puissancestab+A.second-2)!=1) 449 return 0; 450 return *(puissancestab+A.second-1); 451 #else 452 if (*(puissancestab+A.second-3)!=1) 453 return 0; 454 return (unsigned(*(puissancestab+A.second-2)) << 16) + *(puissancestab+A.second-1); 455 #endif 456 } 457 #endif 458 459 #ifdef ADDITIONAL_PRIMES_HASHMAP 460 #ifdef EMCC // container does not seem to be important for <= 70 digits 461 typedef map<unsigned,axbinv> additional_map_t; 462 #else 463 typedef HASH_MAP_NAMESPACE::hash_map<unsigned,axbinv,hash_function_unsigned_object > additional_map_t ; 464 #endif 465 #endif 466 467 #if !defined(RTOS_THREADX) && !defined(BESTA_OS) && !defined NSPIRE 468 // #define WITH_INVA 469 #if defined(__APPLE__) || defined(x86_64) 470 #define LP_TAB_SIZE 15 // slice size will be 2^LP_TAB_SIZE 471 // #define LP_SMALL_PRIMES 472 #define LP_TAB_TOGETHER 473 #define USE_MORE_PRIMES 474 #else 475 #define LP_TAB_SIZE 15 // slice size will be 2^LP_TAB_SIZE 476 #endif // APPLE or 64 bits 477 #endif // !defined RTOS_THREADX and BESTA_OS 478 479 #ifdef LP_TAB_SIZE 480 #define LP_MASK ((1<<LP_TAB_SIZE)-1) 481 struct lp_entry_t { 482 ushort_t pos; 483 ushort_t p; lp_entry_tgiac::lp_entry_t484 lp_entry_t():pos(0),p(0) {}; lp_entry_tgiac::lp_entry_t485 lp_entry_t(ushort_t pos_,ushort_t p_):pos(pos_),p(p_) {}; 486 }; 487 typedef vector<lp_entry_t> lp_tab_t; 488 #endif 489 490 #ifdef LP_TAB_SIZE 491 #define LP_BIT_LIMIT 15 492 #else 493 #define LP_BIT_LIMIT 15 494 #endif 495 496 #if GIAC_ADDITIONAL_PRIMES==16 497 typedef unsigned short additional_t; 498 #else 499 typedef int additional_t; 500 #endif 501 _equalposcomp(const std::vector<additional_t> & v,additional_t w)502 inline int _equalposcomp(const std::vector<additional_t> & v, additional_t w){ 503 int n=1; 504 for (std::vector<additional_t>::const_iterator it=v.begin(),itend=v.end();it!=itend;++it){ 505 if ((*it)==w) 506 return n; 507 else 508 n++; 509 } 510 return 0; 511 } 512 513 // #define SQRTMOD_OUTSIDE 514 #define WITH_LOGP // if defined primes should not exceed 2^24 (perhaps 2^25, choice of sqrt) 515 516 struct small_basis_t { 517 unsigned short root1; 518 unsigned short root2; 519 unsigned short p; 520 unsigned short logp; 521 }; 522 523 #ifdef SQRTMOD_OUTSIDE 524 struct basis_t { 525 unsigned root1; // first root position in slice 526 unsigned root2; // second root position 527 ushort_t p:24; // the prime p 528 #ifdef WITH_LOGP 529 unsigned char logp:8; // could be unsigned char 530 #endif basis_tgiac::basis_t531 basis_t():root1(0),root2(0),p(2) { 532 #ifdef WITH_LOGP 533 logp=sizeinbase2(p); 534 #endif 535 } basis_tgiac::basis_t536 basis_t(ushort_t _p):root1(0),root2(0),p(_p) { 537 #ifdef WITH_LOGP 538 logp=sizeinbase2(p); 539 #endif 540 } 541 } ; 542 543 #else // SQRTMOD_OUTSIDE 544 struct basis_t { 545 unsigned root1; // first root position in slice 546 unsigned root2; // second root position 547 ushort_t p; // the prime p 548 unsigned sqrtmod:24; 549 #ifdef WITH_LOGP 550 unsigned char logp:8; // could be unsigned char 551 #endif basis_tgiac::basis_t552 basis_t():root1(0),root2(0),p(2),sqrtmod(0) { 553 #ifdef WITH_LOGP 554 logp=sizeinbase2(p); 555 #endif 556 } basis_tgiac::basis_t557 basis_t(ushort_t _p):root1(0),root2(0),p(_p),sqrtmod(0) { 558 #ifdef WITH_LOGP 559 logp=sizeinbase2(p); 560 #endif 561 } basis_tgiac::basis_t562 basis_t(ushort_t _p,ushort_t _sqrtmod):root1(0),root2(0),p(_p),sqrtmod(_sqrtmod) { 563 #ifdef WITH_LOGP 564 logp=sizeinbase2(p); 565 #endif 566 } 567 } ; 568 #endif // SQRTMOD_OUTSIDE 569 570 #ifdef LP_SMALL_PRIMES core_sieve(slicetype * slice,small_basis_t * bit,small_basis_t * bitend)571 static inline void core_sieve(slicetype * slice,small_basis_t * bit,small_basis_t * bitend) { 572 for (;bit!=bitend;++bit){ 573 // first root is at bit->root1 574 register unsigned p=bit->p; 575 register unsigned char nbits=bit->logp; 576 register unsigned pos=bit->root1,pos2=bit->root2; 577 if (pos==pos2){ 578 for (;pos<32768; pos += p){ 579 slice[pos] -= nbits; 580 } 581 bit->root2=bit->root1 = pos-32768; // save for next slice 582 } 583 else { 584 for (;pos<32768; pos += p){ 585 slice[pos] -= nbits; 586 } 587 bit->root1 = pos-32768; // save for next slice 588 // second root, polynomial has 2 distinct roots 589 for (;pos2<32768;pos2 += p){ 590 slice[pos2] -= nbits; 591 } 592 bit->root2 = pos2-32768; 593 } 594 } 595 } 596 597 #else // LP_SMALL_PRIMES 598 599 #ifdef LP_TAB_SIZE 600 #define SLICEEND (1<<LP_TAB_SIZE) 601 #else 602 #define SLICEEND ss 603 #endif 604 605 // return position of last prime sieved (useful when large prime hashtable is enabled core_sieve(slicetype * slice,int ss,basis_t * bit,basis_t * bitend)606 static inline basis_t * core_sieve(slicetype * slice,int ss,basis_t * bit,basis_t * bitend) { 607 register unsigned char nbits=sizeinbase2(bit->p); 608 // int next=1 << nbits; 609 for (;bit!=bitend;++bit){ 610 // first root is at bit->root1 611 register ushort_t p=bit->p; 612 #ifdef WITH_LOGP 613 nbits=bit->logp; 614 #else 615 if (p>next){ 616 ++nbits; 617 #if !defined(BESTA_OS) && !defined(RTOS_THREADX) && !defined NSPIRE 618 if (nbits==LP_BIT_LIMIT+1) 619 break; 620 #endif 621 next *=2; 622 } 623 #endif 624 register unsigned pos=bit->root1,pos2=bit->root2; 625 if (pos==pos2){ 626 for (;int(pos)<SLICEEND; pos += p){ 627 slice[pos] -= nbits; 628 } 629 bit->root2=bit->root1 = pos-SLICEEND; // save for next slice 630 } 631 else { 632 for (;int(pos)<SLICEEND; pos += p){ 633 slice[pos] -= nbits; 634 } 635 bit->root1 = pos-SLICEEND; // save for next slice 636 // second root, polynomial has 2 distinct roots 637 for (;int(pos2)<SLICEEND;pos2 += p){ 638 slice[pos2] -= nbits; 639 } 640 bit->root2 = pos2-SLICEEND; 641 } 642 } 643 #if !defined(RTOS_THREADX) && !defined(BESTA_OS) && !defined NSPIRE 644 #ifndef LP_TAB_SIZE 645 for (;bit!=bitend;++bit){ 646 // same as above but we are sieving with primes >2^15, no need to check for nbits increase 647 register ushort_t p=bit->p; 648 register unsigned pos=bit->root1; 649 for (;pos<ss; pos += p){ 650 slice[pos] -= LP_BIT_LIMIT+1; 651 } 652 bit->root1 = pos-ss; // save for next slice 653 // if (sameroot) continue; 654 pos=bit->root2; 655 for (;pos<ss;pos += p){ 656 slice[pos] -= LP_BIT_LIMIT+1; 657 } 658 bit->root2 = pos-ss; 659 } 660 #endif 661 #endif 662 return bit; 663 } 664 #endif // LP_SMALL_PRIMES 665 666 // sieve in [sqrtN+shift,sqrtN+shift+slice.size()-1] 667 // return -1 if memory problem, or the number of relations msieve(const gen & a,const vecteur & sqrtavals,const vecteur & bvals,const mpz_t & c,vector<basis_t> & basis,unsigned lp_basis_pos,vector<small_basis_t> & small_basis,unsigned maxadditional,additional_map_t & additional_primes_map,const gen & N,const gen & isqrtN,slicetype * slice,int ss,int shift,ushort_t * puissancesbegin,ushort_t * & puissancesptr,ushort_t * puissancesend,vector<ushort_t> & curpuissances,vector<ushort_t> & recheck,vector<axbinv> & axbmodn,mpz_t & z1,mpz_t & z2,mpz_t & z3,mpz_t & alloc1,mpz_t & alloc2,mpz_t & alloc3,mpz_t & alloc4,mpz_t & alloc5,const lp_tab_t & lp_tab,GIAC_CONTEXT)668 int msieve(const gen & a,const vecteur & sqrtavals, 669 const vecteur &bvals,const mpz_t& c, 670 vector<basis_t> & basis,unsigned lp_basis_pos, 671 #ifdef LP_SMALL_PRIMES 672 vector<small_basis_t> & small_basis, 673 #endif 674 unsigned maxadditional, 675 #ifdef ADDITIONAL_PRIMES_HASHMAP 676 additional_map_t & additional_primes_map, 677 #else 678 vector<additional_t> & additional_primes,vector<bool> & additional_primes_twice, 679 #endif 680 const gen & N,const gen & isqrtN, 681 slicetype * slice,int ss,int shift, 682 ushort_t * puissancesbegin,ushort_t* & puissancesptr,ushort_t * puissancesend, 683 vector<ushort_t> & curpuissances,vector<ushort_t> &recheck, 684 vector<axbinv> & axbmodn, 685 mpz_t & z1,mpz_t & z2,mpz_t & z3,mpz_t & alloc1,mpz_t & alloc2,mpz_t & alloc3,mpz_t & alloc4,mpz_t & alloc5, 686 #ifdef LP_TAB_SIZE 687 const lp_tab_t & lp_tab, 688 #endif 689 GIAC_CONTEXT){ 690 int nrelations=0; 691 // first fill slice with expected number of bits of 692 // (isqrtN+shift)^2-N = 2*shift*isqrtN + negl. 693 // -> log(2*isqrtN)+log(shift) 694 int shiftss=absint(shift+ss),absshift=absint(shift); 695 int nbits=mpz_sizeinbase(*isqrtN._ZINTptr,2)+sizeinbase2(absshift>shiftss?absshift:shiftss); 696 // int nbits1=int(0.5+std::log(evalf_double(isqrtN,1,context0)._DOUBLE_val/2.*(absshift>shiftss?absshift:shiftss))/std::log(2.)); 697 // int curbits=0; 698 int bs=int(basis.size()); 699 double up_to=1.5; 700 if (nbits>70) 701 up_to += (0.8*(nbits-70))/70; 702 if (debug_infolevel>7) 703 *logptr(contextptr) << CLOCK() << gettext("Sieve tolerance factor ") << up_to << '\n'; 704 unsigned char logB=(unsigned char) (nbits-int(up_to*sizeinbase2(basis.back().p)+.5)); 705 // unsigned char logB=(unsigned char) (nbits-int(up_to*std::log(double(basis.back().p))/std::log(2.0)+.5)); 706 if (debug_infolevel>6) 707 *logptr(contextptr) << CLOCK() << gettext(" reset") << '\n'; 708 // assumes slice type is size 1 byte and multiple of 32 709 #ifdef x86_64 710 ulonglong * ptr=(ulonglong *) &slice[0]; 711 ulonglong * ptrend=ptr+ss/8; 712 ulonglong pattern=(logB <<24)|(logB<<16)|(logB<<8) | logB; 713 pattern = (pattern << 32) | pattern; 714 for (;ptr!=ptrend;++ptr){ 715 *ptr=pattern; 716 } 717 #else 718 unsigned * ptr=(unsigned *) &slice[0]; 719 unsigned * ptrend=ptr+ss/4; 720 unsigned pattern=(logB <<24)|(logB<<16)|(logB<<8) | logB; 721 for (;ptr!=ptrend;++ptr){ 722 *ptr=pattern; 723 } 724 #endif 725 if (debug_infolevel>8) 726 *logptr(contextptr) << CLOCK() << gettext(" end reset, nbits ") << nbits << '\n'; 727 // now for all primes p in basis move in slice from p to p 728 // decrease slice[] by number of bits in p 729 // determines the first prime used in basis 730 #if 0 // def WITH_LOGP 731 nbits=2*mpz_sizeinbase(*isqrtN._ZINTptr,2); 732 int next=50; 733 // note that msieve leaves 20 to 22 primes for normal range, and 15 for large 734 nbits = sizeinbase2(next); 735 #else 736 if (nbits>120) 737 nbits = 7; 738 else { 739 if (nbits>90) 740 nbits = 6; 741 else { 742 if (nbits>78) 743 nbits=5; 744 else 745 nbits = 4; 746 } 747 } 748 int next = 1 << (nbits-1); 749 #endif 750 unsigned bstart; 751 for (bstart=0;bstart<basis.size();++bstart){ 752 int p=basis[bstart].p; 753 if (p>next){ 754 if (debug_infolevel>7) 755 *logptr(contextptr) << gettext("Sieve first prime ") << p << " nbits " << nbits << '\n'; 756 break; 757 } 758 #ifdef LP_SMALL_PRIMES 759 int pos=small_basis[bstart].root1; 760 pos=(pos-ss)%p; 761 if (pos<0) 762 pos+=p; 763 small_basis[bstart].root1=pos; 764 pos=small_basis[bstart].root2; 765 pos=(pos-ss)%p; 766 if (pos<0) 767 pos+=p; 768 small_basis[bstart].root2=pos; 769 #else 770 // update pos_root_mod for later check 771 int pos=basis[bstart].root1; 772 pos=(pos-ss)%p; 773 if (pos<0) 774 pos+=p; 775 basis[bstart].root1=pos; 776 pos=basis[bstart].root2; 777 pos=(pos-ss)%p; 778 if (pos<0) 779 pos+=p; 780 basis[bstart].root2=pos; 781 #endif 782 } 783 next *= 2; 784 if (debug_infolevel>8) 785 *logptr(contextptr) << CLOCK() << gettext(" sieve begin ") << '\n'; 786 // bool sameroot; // Should be there to avoid counting twice the same root but it's faster to ignore it..; 787 #ifdef LP_SMALL_PRIMES 788 small_basis_t * bit=&small_basis[bstart], * bitend=&small_basis[0]+small_basis.size(); 789 core_sieve(slice,bit,bitend); 790 #else 791 basis_t * bit=&basis[bstart], * bitend=&basis[0]+bs; 792 #ifdef LP_TAB_SIZE 793 bitend=core_sieve(slice,ss,bit,&basis[0]+lp_basis_pos); 794 #else 795 bitend=core_sieve(slice,ss,bit,bitend); 796 #endif 797 #endif 798 slicetype * st=slice, * stend=slice+ss; 799 #ifdef LP_TAB_SIZE 800 // sieve for large prime using saved position 801 if (!lp_tab.empty()){ 802 const lp_entry_t * lpit=&lp_tab[0],*lpitend=lpit+lp_tab.size(),*lpitend1=lpitend-8; 803 if (lpitend-lpit>8){ 804 for (;lpit<lpitend1;lpit+=8){ 805 //PREFETCH(lpit + 16); 806 slice[lpit->pos] -= 16; 807 slice[lpit[1].pos] -= 16; 808 slice[lpit[2].pos] -= 16; 809 slice[lpit[3].pos] -= 16; 810 slice[lpit[4].pos] -= 16; 811 slice[lpit[5].pos] -= 16; 812 slice[lpit[6].pos] -= 16; 813 slice[lpit[7].pos] -= 16; 814 } 815 } 816 for (;lpit<lpitend;++lpit) 817 slice[lpit->pos] -= 16; 818 } 819 #endif 820 unsigned cl=0; 821 if (debug_infolevel>6) 822 cl=CLOCK(); 823 if (debug_infolevel>8) 824 *logptr(contextptr) << cl << gettext("relations ") << '\n'; 825 // now find relations 826 st=slice; stend=slice+ss; 827 #ifdef x86_64 828 ulonglong * st8=(ulonglong *) &slice[0],*st8end=st8+ss/8; 829 #else 830 unsigned * st4=(unsigned *) &slice[0],*st4end=st4+ss/4; 831 #endif 832 for ( 833 #ifdef x86_64 834 ;st8!=st8end;st8+=4 835 #else 836 ;st4<st4end;st4+=8 837 #endif 838 ){ 839 // compare slice[pos] to boundary 840 #ifdef x86_64 841 if ( !( (*st8 | st8[1] | st8[2] | st8[3] ) & 0x8080808080808080) ) 842 continue; 843 int pos=int(((slicetype*)st8)-slice); 844 #else 845 if ( !( (*st4 | st4[1] | st4[2] | st4[3] | st4[4] | st4[5] | st4[6] | st4[7]) & 0x80808080) ) 846 continue; 847 int pos=((slicetype*)st4)-slice; 848 #endif 849 st = slice+pos; 850 for (int stpos=0;stpos<32;++st,++pos,++stpos){ 851 if (!(*st&0x80)) 852 continue; 853 // factor (isqrtN+shift+pos)^2-N on basis 854 curpuissances.clear(); recheck.clear(); 855 int shiftpos=shift+pos; 856 #if 0 857 gen tmp=shiftpos; 858 tmp=(a*tmp+2*bvals.back())*tmp+c; 859 tmp.uncoerce(); mpz_set(z1,*tmp._ZINTptr); 860 *logptr(contextptr) << tmp << '\n'; 861 #else 862 mpz_set_si(z1,shiftpos); 863 mpz_mul(z2,z1,*a._ZINTptr); 864 mpz_mul_2exp(z3,*bvals.back()._ZINTptr,1); 865 mpz_add(z2,z2,z3); 866 mpz_mul(z3,z1,z2); 867 mpz_add(z1,z3,c); 868 #endif 869 if (mpz_cmp_si(z1,0)<0){ // if (is_positive(-tmp,context0)) 870 curpuissances.push_back(0xffff); 871 mpz_neg(z1,z1); // tmp=-tmp; 872 } 873 bool done=false; 874 #ifdef LP_TAB_SIZE 875 #ifdef LP_SMALL_PRIMES 876 small_basis_t * basisptr=&small_basis[0], * basisend=basisptr+(bitend-bit); 877 #else 878 basis_t * basisptr=&basis[0], * basisend=basisptr+(bitend-bit); 879 #endif 880 #else // LP_TAB_SIZE 881 basis_t * basisptr=&basis[0], * basisend=basisptr+bs; 882 #endif 883 // we have modified pos_root_mod1 and pos_root_mod2 -> posss 884 int posss=ss-pos; // always positive 885 for (;basisptr!=basisend;++basisptr){ 886 register int bi=basisptr->p; 887 // check if we have a root 888 register int check=bi-(posss%bi); 889 if (check!=bi && check!=int(basisptr->root1) && check!=int(basisptr->root2)) 890 continue; 891 if (check==bi && basisptr->root1 && basisptr->root2) 892 continue; 893 recheck.push_back(bi); 894 } // end for on (small) primes 895 #ifdef LP_TAB_SIZE 896 // add primes from large prime hashtable 897 lp_tab_t::const_iterator lpit=lp_tab.begin(),lpend=lp_tab.end(); 898 int hash_pos=recheck.size(); 899 for (;lpit!=lpend;++lpit){ 900 if (pos==int(lpit->pos)){ 901 recheck.push_back(lpit->p); 902 } 903 } 904 if (int(recheck.size())>hash_pos+1) 905 sort(recheck.begin(),recheck.end()); 906 #endif 907 // now divide first by product of elements of recheck 908 double prod=1,nextprod=1; 909 for (unsigned k=0;k<recheck.size();++k){ 910 nextprod=prod*recheck[k]; 911 if (nextprod< 2147483648. ) 912 prod=nextprod; 913 else { 914 // mpz_fdiv_q_ui(z1,z1,prod); 915 mpz_set_si(z2,int(prod)); 916 #ifdef USE_GMP_REPLACEMENTS 917 mp_grow(&alloc1,z1.used+2); 918 mpz_set_ui(alloc1,0); 919 alloc1.used = z1.used +2 ; 920 mpz_set(alloc2,z1); 921 mpz_set(alloc3,z2); 922 alloc_mp_div(&z1,&z2,&z1,&z3,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 923 #else 924 mpz_divexact(z3,z1,z2); 925 mpz_swap(z1,z3); 926 #endif 927 prod=recheck[k]; 928 } 929 } 930 if (prod!=1){ 931 // mpz_fdiv_q_ui(z1,z1,prod); 932 mpz_set_si(z2,int(prod)); 933 #ifdef USE_GMP_REPLACEMENTS 934 mp_grow(&alloc1,z1.used+2); 935 mpz_set_ui(alloc1,0); 936 alloc1.used = z1.used +2 ; 937 mpz_set(alloc2,z1); 938 mpz_set(alloc3,z2); 939 alloc_mp_div(&z1,&z2,&z1,&z3,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 940 #else 941 mpz_divexact(z3,z1,z2); 942 mpz_swap(z1,z3); 943 #endif 944 } 945 // then set curpuissances 946 bool small_=false; 947 int Z1=0; 948 for (unsigned k=0;k<recheck.size();++k){ 949 int j=0; 950 int bi=recheck[k]; 951 #ifdef USE_GMP_REPLACEMENTS 952 if (!small_){ 953 small_=mpz_sizeinbase(z1,2)<32; 954 if (small_) 955 Z1=mpz_get_si(z1); 956 } 957 if (small_){ 958 for (++j;;++j){ 959 #if 0 // def FXCG 960 if (Z1 % bi) 961 break; 962 Z1 /= bi; 963 #else 964 div_t qr; 965 qr=div(Z1,bi); 966 if (qr.rem) 967 break; 968 Z1=qr.quot; 969 #endif 970 } 971 } 972 else { 973 for (++j;;++j){ 974 #if 1 975 mpz_set_ui(z2,bi); 976 mp_grow(&alloc1,z1.used+2); 977 mpz_set_ui(alloc1,0); 978 alloc1.used = z1.used +2 ; 979 mpz_set(alloc2,z1); 980 mpz_set(alloc3,z2); 981 alloc_mp_div(&z1,&z2,&z2,&z3,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 982 #else 983 mpz_fdiv_qr_ui(z2,z3,z1,bi); 984 #endif 985 if (mpz_cmp_si(z3,0)) 986 break; 987 mpz_set(z1,z2); 988 } 989 } 990 #else 991 for (++j;;++j){ 992 mpz_fdiv_qr_ui(z2,z3,z1,bi); 993 if (mpz_cmp_si(z3,0)) 994 break; 995 mpz_set(z1,z2); 996 } 997 #endif 998 /* 999 while (is_zero(smod(tmp,bi))){ 1000 tmp=tmp/bi; 1001 ++j; 1002 } 1003 */ 1004 if (!done && bi>255){ 1005 curpuissances.push_back(0); 1006 done=true; 1007 } 1008 if (done){ 1009 for (;j;--j) 1010 curpuissances.push_back(bi); 1011 } 1012 else { 1013 for (;j>=256;j-=256) 1014 curpuissances.push_back(bi<<8); 1015 if (j) 1016 curpuissances.push_back( (bi << 8) | j); 1017 } 1018 } 1019 if (small_) 1020 mpz_set_si(z1,Z1); 1021 if (mpz_cmp_si(z1,1)==0){ // is_one(tmp)){ 1022 ++nrelations; 1023 if (debug_infolevel>6) 1024 *logptr(contextptr) << CLOCK() << gettext(" true relation ") << '\n'; 1025 axbmodn.push_back(axbinv(int(sqrtavals.size())-1,shiftpos,int(bvals.size())-1,int(puissancesptr-puissancesbegin),int(puissancesptr-puissancesbegin)+int(curpuissances.size()))); 1026 for (unsigned i=0;i<curpuissances.size();++puissancesptr,++i){ 1027 if (puissancesptr>=puissancesend) 1028 return -1; 1029 *puissancesptr=curpuissances[i]; 1030 } 1031 } 1032 else { 1033 unsigned param2; 1034 #if (GIAC_ADDITIONAL_PRIMES==16) 1035 param2=0xffff; 1036 #else 1037 param2=maxadditional; 1038 #endif 1039 if (mpz_cmp_ui(z1,param2)>0){ 1040 if (debug_infolevel>6) 1041 *logptr(contextptr) << gen(z1) << gettext(" Sieve large remainder:") << '\n'; 1042 } 1043 else { 1044 #ifdef GIAC_ADDITIONAL_PRIMES 1045 additional_t P=mpz_get_ui(z1); 1046 // if (int(P)>2*int(basis.back())) continue; 1047 // if (debug_infolevel>5) 1048 if (debug_infolevel>6) 1049 *logptr(contextptr) << CLOCK() << " " << P << " remain " << '\n'; 1050 #ifdef ADDITIONAL_PRIMES_HASHMAP 1051 // add relation 1052 ++nrelations; 1053 curpuissances.push_back(1); // marker 1054 #if (GIAC_ADDITIONAL_PRIMES==32) 1055 #ifndef PRIMES32 1056 curpuissances.push_back(P >> 16); 1057 #endif 1058 #endif 1059 curpuissances.push_back(P); 1060 for (unsigned i=0;i<curpuissances.size();++puissancesptr,++i){ 1061 if (puissancesptr>=puissancesend) 1062 return -1; 1063 *puissancesptr=curpuissances[i]; 1064 } 1065 additional_map_t::iterator it=additional_primes_map.find(P),itend=additional_primes_map.end(); 1066 if (it!=itend) // build a large prime relation (P is the large prime) 1067 axbmodn.push_back(axbinv(sqrtavals.size()-1,shiftpos,bvals.size()-1,(puissancesptr-puissancesbegin)-curpuissances.size(),(puissancesptr-puissancesbegin))); 1068 else // record a partial relation 1069 additional_primes_map[P]=axbinv(sqrtavals.size()-1,shiftpos,bvals.size()-1,(puissancesptr-puissancesbegin)-curpuissances.size(),(puissancesptr-puissancesbegin)); 1070 #else 1071 int Ppos=_equalposcomp(additional_primes,P); // this is in O(additional^2)=o(B^3) 1072 if (Ppos){ 1073 if (debug_infolevel>6) 1074 *logptr(contextptr) << P << gettext(" already additional") << '\n'; 1075 --Ppos; 1076 additional_primes_twice[Ppos]=true; 1077 } else { 1078 // add a prime in additional_primes if <=QS_B_BOUND 1079 if (int(additional_primes.size())>=4*bs 1080 #if defined(RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 1081 || bs+additional_primes.size()>700 1082 #endif 1083 ) 1084 continue; 1085 additional_primes.push_back(P); 1086 additional_primes_twice.push_back(false); 1087 Ppos=int(additional_primes.size())-1; 1088 } 1089 // add relation 1090 curpuissances.push_back(1); // marker 1091 #if GIAC_ADDITIONAL_PRIMES==32 1092 #ifndef PRIMES32 1093 curpuissances.push_back(P >> 16); 1094 #endif 1095 #endif 1096 curpuissances.push_back(P); 1097 axbmodn.push_back(axbinv(int(sqrtavals.size())-1,shiftpos,int(bvals.size())-1,int(puissancesptr-puissancesbegin),int(puissancesptr-puissancesbegin)+int(curpuissances.size()))); 1098 for (unsigned i=0;i<curpuissances.size();++puissancesptr,++i){ 1099 if (puissancesptr>=puissancesend) 1100 return -1; 1101 *puissancesptr=curpuissances[i]; 1102 } 1103 #endif // ADDITIONAL_PRIMES_HASHMAP 1104 #endif // GIAC_ADDITIONAL_PRIMES 1105 } 1106 } 1107 } 1108 } // end for loop on slice array 1109 if (debug_infolevel>6){ 1110 unsigned cl2=CLOCK(); 1111 *logptr(contextptr) << cl2 << gettext(" end relations ") << cl2-cl << '\n'; 1112 } 1113 return nrelations; 1114 } 1115 1116 // #define MP_MODINV_1 1117 #ifdef MP_MODINV_1 mp_modinv_1(unsigned a,unsigned p)1118 static inline unsigned mp_modinv_1(unsigned a, unsigned p) { 1119 1120 unsigned ps1, ps2, dividend, divisor, rem, q, t; 1121 unsigned parity; 1122 1123 q = 1; rem = a; dividend = p; divisor = a; 1124 ps1 = 1; ps2 = 0; parity = 0; 1125 1126 while (divisor > 1) { 1127 rem = dividend - divisor; 1128 t = rem - divisor; 1129 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1130 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1131 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1132 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1133 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1134 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1135 if (rem >= divisor) { q += ps1; rem = t; t -= divisor; 1136 if (rem >= divisor) { q += ps1; rem = t; 1137 if (rem >= divisor) { 1138 q = dividend / divisor; 1139 rem = dividend % divisor; 1140 q *= ps1; 1141 } 1142 } 1143 } 1144 } 1145 } 1146 } 1147 } 1148 } 1149 } 1150 1151 q += ps2; 1152 parity = ~parity; 1153 dividend = divisor; 1154 divisor = rem; 1155 ps2 = ps1; 1156 ps1 = q; 1157 } 1158 1159 if (parity == 0) 1160 return ps1; 1161 else 1162 return p - ps1; 1163 } 1164 #endif 1165 1166 #if (defined __i386__ || defined __x86_64__) && !defined PIC && !defined _I386_ && !defined __APPLE__ && !defined VISUALC && !defined(FIR_LINUX) && !defined(FIR_ANDROID) 1167 #define _I386_ 1168 #endif 1169 1170 #ifdef _I386_ 1171 // a->a+b*c mod m addmultmod(int & a,int b,int c,int m)1172 inline void addmultmod(int & a,int b,int c,int m){ 1173 asm volatile("testl %%ebx,%%ebx\n\t" /* sign bit=1 if negative */ 1174 "jns .Lok%=\n\t" 1175 "addl %%edi,%%ebx\n" /* a+=m*/ 1176 ".Lok%=:\t" 1177 "imull %%ecx; \n\t" /* b*c in edx:eax */ 1178 "addl %%ebx,%%eax; \n\t" /* b*c+a */ 1179 "adcl $0x0,%%edx; \n\t" /* b*c+a carry */ 1180 "idivl %%edi; \n\t" 1181 :"=d"(a) 1182 :"a"(b),"b"(a),"c"(c),"D"(m) 1183 ); 1184 } 1185 #endif 1186 1187 inline modmult(int a,int b,unsigned p)1188 int modmult(int a,int b,unsigned p){ 1189 #ifdef _I386_ 1190 register int res; 1191 asm volatile("imull %%edx\n\t" /* a*b-> edx:eax */ 1192 "idivl %%ecx\n\t" /* edx:eax div p -> quotient=eax, remainder=edx */ 1193 :"=d"(res) 1194 :"a"(a),"d"(b),"c"(p) 1195 : 1196 ); 1197 return res; 1198 #else 1199 return a*longlong(b) % p; 1200 #endif 1201 } 1202 1203 // assumes b>0 and |a|<b invmodnoerr(int a,int b)1204 int invmodnoerr(int a,int b){ 1205 if (a==1 || a==-1 || a==0) 1206 return a; 1207 if (a<0) // insure a>0 so that all remainders below are >=0 1208 a+=b; 1209 #ifdef _I386_ // works only for ushort_t == unsigned short 1210 // int res=mp_modinv_1(a,b),p=b; 1211 /* GDB: si will step in assembly, info registers show register content, x/i $pc show next ins */ 1212 asm volatile("movl $0,%%edi\n\t" 1213 "movl $1,%%ecx\n\t" 1214 "movl $0,%%edx\n\t" 1215 ".Lloop%=:\t" 1216 "movl %%esi,%%eax\n\t" 1217 "andl $0x80000000,%%esi\n\t" 1218 "xorl $0x80000000,%%esi\n\t" /* parity indicator for sign */ 1219 "andl $0x7fffffff,%%eax\n\t" /* clear high bit of ax */ 1220 "divl %%ebx\n\t" /* divide si by bx, ax=quotient, dx=rem */ 1221 "orl %%ebx,%%esi\n\t" /* copy bx in si but keep high bit of si */ 1222 "movl %%edx,%%ebx\n\t" /* si now contains bx and bx the remainder */ 1223 "mull %%ecx\n\t" /* quotient*cx is in ax (dx=0) */ 1224 "addl %%eax,%%edi\n\t" /* di <- di+q*cx*/ 1225 "xchgl %%edi,%%ecx\n\t" /* cx <- origi di+q*cx, di <- orig cx */ 1226 "testl %%ebx,%%ebx\n\t" 1227 "jne .Lloop%=\n\t" 1228 :"=D"(a),"=S"(b) 1229 :"S"(b),"b"(a) 1230 :"%eax","%ecx","%edx" 1231 ); 1232 if (b<0) 1233 b=b&0x7fffffff; 1234 else 1235 a=-a; 1236 a=(b==1)?a:0; 1237 // if ((a-res)%p) 1238 // CERR << "error" << '\n'; 1239 return a; 1240 #else // i386 1241 1242 #ifdef MP_MODINV_1 1243 return mp_modinv_1(a,b); 1244 #endif 1245 // r0=b=ab*a+1*b 1246 // r1=a=aa*a+0*b 1247 int aa(1),ab(0),ar(0); 1248 #if 0 // def FXCG 1249 ushort_t q,r; 1250 while (a){ 1251 q=b/a; 1252 ar=ab-q*aa; 1253 r=b-q*a; 1254 if (!r) 1255 return a==1?aa:0; 1256 q=a/r; 1257 ab=aa-q*ar; 1258 b=a-q*r; 1259 if (!b) 1260 return r==1?ar:0; 1261 q=r/b; 1262 aa=ar-q*ab; 1263 a=r-q*b; 1264 } 1265 return b==1?ab:0; 1266 #else 1267 div_t qr; 1268 while (a){ 1269 qr=div(b,a); 1270 ar=ab-qr.quot*aa; 1271 b=a; 1272 a=qr.rem; 1273 ab=aa; 1274 aa=ar; 1275 } 1276 if (b==1) 1277 return ab; 1278 return 0; 1279 #endif 1280 #endif // i386 1281 } 1282 1283 #if 0 // def PRIMES32 1284 // assumes |a|<b 1285 longlong invmodnoerr(longlong a,longlong b){ 1286 if (a==1 || a==-1 || a==0) 1287 return a; 1288 if (a<0) // insure a>0 so that all remainders below are >=0 1289 a+=b; 1290 // r0=b=ab*a+1*b 1291 // r1=a=aa*a+0*b 1292 longlong aa(1),ab(0),ar(0); 1293 longlong q,r; 1294 lldiv_t qr; 1295 while (a){ 1296 qr=lldiv(b,a); 1297 ar=ab-qr.quot*aa; 1298 b=a; 1299 a=qr.rem; 1300 ab=aa; 1301 aa=ar; 1302 } 1303 if (b==1) 1304 return ab; 1305 return 0; 1306 } 1307 #endif 1308 find_multiplier(const gen & n,double & delta,GIAC_CONTEXT)1309 static int find_multiplier(const gen & n,double & delta,GIAC_CONTEXT){ 1310 delta=0; 1311 if (n.type!=_ZINT) 1312 return 1; 1313 static const unsigned char mult[] = 1314 { 1, 3, 5, 7, 11, 13, 15, 17, 19, 1315 21, 23, 29, 31, 33, 35, 37, 39, 41, 43, 47}; // only odd values for multiplier 1316 unsigned nmult=sizeof(mult)/sizeof(unsigned char); 1317 double scores[50]; 1318 int nmodp=modulo(*n._ZINTptr,8),knmodp; 1319 // init scores and set value for 2 1320 double ln2=std::log(2.0); 1321 for (unsigned i=0;i<nmult;++i){ 1322 knmodp=(mult[i]*nmodp)%8; 1323 scores[i]=0.5*std::log(double(mult[i])); 1324 switch(knmodp){ 1325 case 1: 1326 scores[i] -= 2*ln2; 1327 break; 1328 case 5: 1329 scores[i] -= ln2; 1330 break; 1331 case 3: case 7: 1332 scores[i] -= 0.5 * ln2; 1333 break; 1334 } 1335 } 1336 // now compute contribution for giac_primes[1..300] 1337 for (unsigned i=1;i<=300;++i){ 1338 int p=giac_primes[i]; 1339 double contrib=std::log(double(p))/(p-1); 1340 nmodp=modulo(*n._ZINTptr,p); 1341 for (unsigned j=0;j<nmult;++j){ 1342 knmodp=(nmodp*mult[j])%p; 1343 if (knmodp==0) 1344 scores[j] -= contrib; 1345 else { 1346 if (!is_undef(sqrt_mod(knmodp,p,true,context0))) 1347 scores[j] -= 2*contrib; 1348 } 1349 } 1350 } 1351 // select the smallest scores 1352 int pos=0; 1353 double minscore=scores[0]-0.1; 1354 for (unsigned i=1;i<nmult;++i){ 1355 if (scores[i]<minscore){ 1356 minscore=scores[i]; 1357 pos=i; 1358 } 1359 } 1360 if (debug_infolevel>6){ 1361 for (unsigned i=0;i<nmult;++i){ 1362 *logptr(contextptr) << gettext("multiplier ") << int(mult[i]) << " score " << scores[i] << '\n'; 1363 } 1364 } 1365 if (pos){ 1366 delta=minscore-scores[0]; 1367 if (debug_infolevel) 1368 *logptr(contextptr) << gettext("Using multiplier ") << int(mult[pos]) << " delta-score " << delta << '\n'; 1369 } 1370 return mult[pos]; 1371 } 1372 add_relation(vector<line_t> relations,unsigned j,ushort_t * curpui,ushort_t * curpuiend,const vector<basis_t> & basis,const vector<additional_t> & additional_primes)1373 void add_relation(vector<line_t> relations,unsigned j,ushort_t * curpui,ushort_t * curpuiend,const vector<basis_t> & basis,const vector<additional_t> & additional_primes){ 1374 unsigned curpuisize=unsigned(curpuiend-curpui); 1375 bool done=false; 1376 unsigned i=0; // position in basis 1377 unsigned k=0; // position in curpui 1378 additional_t p=0; // prime 1379 unsigned bs=unsigned(basis.size()); 1380 for (;k<curpuisize;++k){ 1381 p=curpui[k]; 1382 if (p==0xffff){ 1383 relations[0].tab[j/32] ^= (1 << (j%32)); 1384 continue; 1385 } 1386 if (p==0){ 1387 done=true; 1388 continue; 1389 } 1390 if (p==1){ 1391 #ifndef ADDITIONAL_PRIMES_HASHMAP 1392 p=curpui[k+1]; 1393 #if GIAC_ADDITIONAL_PRIMES==32 1394 #ifndef PRIMES32 1395 p <<= 16; 1396 p += curpui[k+2]; 1397 #endif 1398 #endif 1399 // k must be == curpui.size()-1 1400 // find p in additional_primes and position 1401 int Ppos=_equalposcomp(additional_primes,p); 1402 // *logptr(contextptr) << p << " " << Ppos+bs << " " << relations.size() << '\n'; 1403 relations[bs+Ppos].tab[j/32] |= (1 << (j %32)); 1404 #endif 1405 break; 1406 } 1407 if (!done){ 1408 if (p%2==0) // even exponent? 1409 continue; 1410 p >>= 8; 1411 } 1412 else { 1413 int c=1; 1414 for (;k+1<curpuisize;c++){ 1415 if (curpui[k+1]==unsigned(p)) 1416 ++k; 1417 else 1418 break; 1419 } 1420 if (c%2==0) 1421 continue; 1422 } 1423 // advance to next i in basis 1424 for (;i<bs;++i){ 1425 if (basis[i].p==unsigned(p)) 1426 break; 1427 } 1428 if (i<bs){ 1429 ++i; 1430 relations[i].tab[j/32] ^= (1 << (j %32)); 1431 } 1432 else { 1433 // ERROR 1434 } 1435 } // end loop on k in curpui 1436 } 1437 update_xy(axbinv & A,mpz_t & zx,mpz_t & zy,vector<short_t> & p,vector<short_t> & add_p,const gen & N,const vector<basis_t> & basis,const vector<additional_t> & additional_primes,const vecteur & sqrtavals,const vecteur & bvals,ushort_t * puissancestab,mpz_t & zq,mpz_t & zr,mpz_t & alloc1,mpz_t & alloc2,mpz_t & alloc3,mpz_t & alloc4,mpz_t & alloc5)1438 void update_xy(axbinv & A,mpz_t & zx,mpz_t & zy,vector<short_t> & p,vector<short_t> & add_p,const gen & N,const vector<basis_t> & basis,const vector<additional_t> & additional_primes,const vecteur & sqrtavals,const vecteur & bvals,ushort_t * puissancestab,mpz_t & zq,mpz_t & zr,mpz_t & alloc1, mpz_t & alloc2,mpz_t & alloc3,mpz_t & alloc4, mpz_t & alloc5){ 1439 // x=x*(a*shiftpos+b), y =y*sqrta; 1440 mpz_set_si(alloc2,A.shiftpos); 1441 if (sqrtavals[A.aindex].type==_INT_){ 1442 mpz_mul_ui(alloc1,alloc2,sqrtavals[A.aindex].val); 1443 mpz_mul_ui(alloc2,alloc1,sqrtavals[A.aindex].val); 1444 mpz_mul_ui(zy,zy,sqrtavals[A.aindex].val); 1445 } 1446 else { 1447 mpz_mul(alloc1,alloc2,*sqrtavals[A.aindex]._ZINTptr); 1448 mpz_mul(alloc2,alloc1,*sqrtavals[A.aindex]._ZINTptr); 1449 mpz_mul(zy,zy,*sqrtavals[A.aindex]._ZINTptr); 1450 } 1451 mpz_add(alloc1,alloc2,*bvals[A.bindex]._ZINTptr); 1452 // mpz_mul(alloc2,alloc1,*invsqrtamodnvals[A.aindex]._ZINTptr); 1453 mpz_mul(zr,zx,alloc1); 1454 #ifdef USE_GMP_REPLACEMENTS 1455 mp_grow(&alloc1,zr.used+2); 1456 mpz_set_ui(alloc1,0); 1457 alloc1.used = zr.used +2 ; 1458 mpz_set(alloc2,zr); 1459 mpz_set(alloc3,*N._ZINTptr); 1460 // mpz_set_si(alloc4,0); 1461 // mpz_set_si(alloc5,0); 1462 alloc_mp_div(&zr,N._ZINTptr,&zq,&zx,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 1463 mp_grow(&alloc1,zy.used+2); 1464 mpz_set_ui(alloc1,0); 1465 alloc1.used = zy.used +2 ; 1466 mpz_set(alloc2,zy); 1467 mpz_set(alloc3,*N._ZINTptr); 1468 // mpz_set_si(alloc4,0); 1469 // mpz_set_si(alloc5,0); 1470 alloc_mp_div(&zy,N._ZINTptr,&zq,&zy,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 1471 #else 1472 mpz_tdiv_r(zx,zr,*N._ZINTptr); 1473 mpz_tdiv_r(zy,zy,*N._ZINTptr); 1474 #endif 1475 bool done=false; 1476 unsigned bi=0; 1477 ushort_t * it=puissancestab+A.first,* itend=puissancestab+A.second; 1478 for (;it!=itend;++it){ 1479 if (*it==0xffff) 1480 continue; 1481 if (*it==1){ 1482 ++it; 1483 additional_t p=*it; 1484 #if GIAC_ADDITIONAL_PRIMES==32 1485 #ifndef PRIMES32 1486 p <<= 16; 1487 ++it; 1488 p += *it; 1489 #endif 1490 #endif 1491 int pos=_equalposcomp(additional_primes,p); 1492 if (pos) 1493 ++add_p[pos-1]; 1494 else { 1495 // otherwise ERROR!!! 1496 } 1497 break; 1498 } 1499 if (!*it){ 1500 done=true; 1501 continue; 1502 } 1503 if (done){ 1504 while (bi<basis.size() && basis[bi].p!=*it) 1505 ++bi; 1506 if (bi<basis.size()) 1507 p[bi]++; 1508 else { 1509 // ERROR 1510 } 1511 } 1512 else { 1513 while (basis[bi].p!=(*it>>8)) 1514 ++bi; 1515 p[bi]+=(*it&0xff); 1516 } 1517 } 1518 } 1519 find_bv_be(int tmp,int & bv,int & be)1520 void find_bv_be(int tmp,int & bv,int &be){ 1521 bv=1; be=-1; 1522 while (tmp%2==0){ 1523 ++bv; 1524 tmp /= 2; 1525 } 1526 tmp /= 2; 1527 if (tmp%2) 1528 be=1; 1529 else 1530 be=-1; 1531 } 1532 1533 1534 #ifdef PRIMES32 1535 // Change b coeff of polynomial: update roots for small primes 1536 // for large primes do it depending on LP_TAB_TOGETHER 1537 #ifdef LP_SMALL_PRIMES copy(vector<basis_t> & basis,vector<small_basis_t> & small_basis)1538 void copy(vector<basis_t> & basis,vector<small_basis_t> & small_basis){ 1539 small_basis_t * small_basisptr=&small_basis[0], * small_basisend=small_basisptr+small_basis.size(); 1540 basis_t * basisptr=&basis[0]; 1541 unsigned next=2,logp=1; 1542 if (small_basis[0].p==0){ 1543 for (;small_basisptr<small_basisend;++basisptr,++small_basisptr){ 1544 small_basisptr->root1=basisptr->root1; 1545 small_basisptr->root2=basisptr->root2; 1546 register unsigned short p =basisptr->p; 1547 small_basisptr->p = p; 1548 small_basisptr->logp=logp; 1549 if (p>next){ 1550 ++logp; 1551 next *= 2; 1552 } 1553 } 1554 } 1555 else { 1556 for (;small_basisptr<small_basisend;++basisptr,++small_basisptr){ 1557 small_basisptr->root1=basisptr->root1; 1558 small_basisptr->root2=basisptr->root2; 1559 } 1560 } 1561 } 1562 switch_roots(const vector<int> & bainv2,vector<basis_t> & basis,vector<small_basis_t> & small_basis,unsigned lp_basis_pos,unsigned nslices,unsigned slicesize,unsigned bv,int be,int afact,const vector<ushort_t> & pos,gen b,mpz_t & zq,int M)1563 void switch_roots(const vector<int> & bainv2,vector<basis_t> & basis,vector<small_basis_t> & small_basis,unsigned lp_basis_pos,unsigned nslices,unsigned slicesize,unsigned bv,int be,int afact,const vector<ushort_t> & pos,gen b,mpz_t & zq,int M){ 1564 unsigned bs=basis.size(); 1565 const int * bvpos=&bainv2[(bv-1)*bs]; 1566 #ifdef LP_TAB_TOGETHER 1567 const int * bvposend=bvpos+lp_basis_pos; 1568 #else 1569 const int * bvposend=bvpos+bs; 1570 #endif 1571 basis_t * basisptr=&basis[0]; 1572 if (be>0){ 1573 for (;bvpos<bvposend;++basisptr,++bvpos){ 1574 // PREFETCH(basisptr+4); 1575 // PREFETCH(bvpos+4); 1576 register unsigned p=basisptr->p; 1577 register int r=basisptr->root1-(*bvpos); 1578 if (r<0) 1579 r+=p; 1580 basisptr->root1=r; 1581 r=basisptr->root2-(*bvpos); 1582 if (r<0) 1583 r+=p; 1584 basisptr->root2=r; 1585 } 1586 } 1587 else { 1588 for (;bvpos<bvposend;++basisptr,++bvpos){ 1589 // PREFETCH(basisptr+4); 1590 // PREFETCH(bvpos+4); 1591 register unsigned p=basisptr->p; 1592 register int r=basisptr->root1+(*bvpos); 1593 if (r>p) 1594 r-=p; 1595 basisptr->root1=r; 1596 r=basisptr->root2+(*bvpos); 1597 if (r>p) 1598 r-=p; 1599 basisptr->root2=r; 1600 } 1601 } 1602 // adjust sieve position for prime factors of a, 1603 for (int j=0;j<afact;++j){ 1604 int pj=pos[j]; 1605 ushort_t p=basis[pj].p; 1606 int q,bmodp=p-modulo(*b._ZINTptr,p); 1607 int cmodp=modulo(zq,p); 1608 q=(M+longlong(cmodp)*invmodnoerr((2*bmodp)%p,p))%p; 1609 if (q<0) 1610 q+=p; 1611 basis[pj].root1=q; 1612 basis[pj].root2=q; 1613 } 1614 // set small primes position for sieving 1615 copy(basis,small_basis); 1616 } 1617 1618 #else // LP_SMALL_PRIMES 1619 switch_roots(const vector<int> & bainv2,vector<basis_t> & basis,unsigned lp_basis_pos,unsigned nslices,unsigned slicesize,unsigned bv,int be,int afact,const vector<ushort_t> & pos,gen b,mpz_t & zq,int M)1620 void switch_roots(const vector<int> & bainv2,vector<basis_t> & basis,unsigned lp_basis_pos,unsigned nslices,unsigned slicesize,unsigned bv,int be,int afact,const vector<ushort_t> & pos,gen b,mpz_t & zq,int M){ 1621 unsigned bs=basis.size(); 1622 #ifdef LP_TAB_SIZE 1623 const int * bvpos=&bainv2[(bv-1)*bs],* bvposend=bvpos+lp_basis_pos; 1624 #else 1625 const int * bvpos=&bainv2[(bv-1)*bs],* bvposend=bvpos+bs; 1626 #endif 1627 basis_t * basisptr=&basis[0]; 1628 unsigned decal0=nslices*slicesize; 1629 if (decal0>=basis.back().p){ 1630 if (be<0){ 1631 for (;bvpos<bvposend;++basisptr,++bvpos){ 1632 register unsigned p=basisptr->p; 1633 register unsigned decal = (decal0+(*bvpos))% p; 1634 register unsigned r=basisptr->root1+decal; 1635 if (r>p) 1636 r -= p; 1637 basisptr->root1 = r; 1638 r = basisptr->root2+decal; 1639 if (r>p) 1640 r -= p; 1641 basisptr->root2 = r; 1642 } 1643 } 1644 else { 1645 for (;bvpos<bvposend;++basisptr,++bvpos){ 1646 register unsigned p=basisptr->p; 1647 register unsigned decal = (decal0-(*bvpos))% p; 1648 register unsigned r=basisptr->root1+decal; 1649 if (r>p) 1650 r -= p; 1651 basisptr->root1 = r; 1652 r = basisptr->root2+decal; 1653 if (r>p) 1654 r -= p; 1655 basisptr->root2 = r; 1656 } 1657 } 1658 } 1659 else 1660 { // should not be reached since Mtarget is about basis.back() 1661 for (;bvpos<bvposend;++basisptr,++bvpos){ 1662 register unsigned p=basisptr->p; 1663 register unsigned decal = (decal0+p-be*(*bvpos))% p; 1664 register unsigned r=basisptr->root1+decal; 1665 if (r>p) 1666 r -= p; 1667 basisptr->root1 = r; 1668 r = basisptr->root2+decal; 1669 if (r>p) 1670 r -= p; 1671 basisptr->root2 = r; 1672 } 1673 } 1674 // adjust sieve position for prime factors of a, 1675 for (int j=0;j<afact;++j){ 1676 int pj=pos[j]; 1677 ushort_t p=basis[pj].p; 1678 int q,bmodp=p-modulo(*b._ZINTptr,p); 1679 int cmodp=modulo(zq,p); 1680 q=(M+longlong(cmodp)*invmodnoerr((2*bmodp)%p,p))%p; 1681 if (q<0) 1682 q+=p; 1683 basis[pj].root1=q; 1684 basis[pj].root2=q; 1685 } 1686 #if defined(LP_TAB_SIZE) && !defined(LP_TAB_TOGETHER) 1687 bvposend += bs-lp_basis_pos; 1688 if (be>0){ 1689 for (;bvpos<bvposend;++basisptr,++bvpos){ 1690 // PREFETCH(basisptr+4); 1691 // PREFETCH(bvpos+4); 1692 register unsigned p=basisptr->p; 1693 register int r=basisptr->root1-(*bvpos); 1694 if (r<0) 1695 r+=p; 1696 basisptr->root1=r; 1697 r=basisptr->root2-(*bvpos); 1698 if (r<0) 1699 r+=p; 1700 basisptr->root2=r; 1701 } 1702 } 1703 else { 1704 for (;bvpos<bvposend;++basisptr,++bvpos){ 1705 // PREFETCH(basisptr+4); 1706 // PREFETCH(bvpos+4); 1707 register unsigned p=basisptr->p; 1708 register int r=basisptr->root1+(*bvpos); 1709 if (r>int(p)) 1710 r-=p; 1711 basisptr->root1=r; 1712 r=basisptr->root2+(*bvpos); 1713 if (r>int(p)) 1714 r-=p; 1715 basisptr->root2=r; 1716 } 1717 } 1718 #endif 1719 } 1720 #endif // LP_SMALL_PRIMES 1721 #endif // PRIMES32 1722 1723 // Change a, the leading coeff of polynomial: initialize all roots (small and large primes) init_roots(vector<basis_t> & basis,vector<small_basis_t> & small_basis,vector<ushort_t> & Inva,const vector<ushort_t> & sqrtmod,vector<int> & bainv2,int afact,int afact0,const gen & a,const gen & b,const vecteur & bvalues,mpz_t & zq,unsigned M)1724 void init_roots(vector<basis_t> & basis, 1725 #ifdef LP_SMALL_PRIMES 1726 vector<small_basis_t> & small_basis, 1727 #endif 1728 #ifdef WITH_INVA 1729 vector<ushort_t> & Inva, 1730 #endif 1731 #ifdef SQRTMOD_OUTSIDE 1732 const vector<ushort_t> & sqrtmod, 1733 #endif 1734 #ifdef PRIMES32 1735 vector<int> & bainv2,int afact,int afact0, 1736 #else 1737 ulonglong usqrta, 1738 #endif 1739 const gen & a,const gen & b,const vecteur & bvalues,mpz_t & zq,unsigned M){ 1740 unsigned bs=unsigned(basis.size()); 1741 basis_t * basisptr=&basis.front(),*basisend=basisptr+bs; 1742 #ifdef SQRTMOD_OUTSIDE 1743 vector<ushort_t>::const_iterator sqrtmodit=sqrtmod.begin(); 1744 #endif 1745 for (int i=0;basisptr!=basisend;++i,++basisptr){ 1746 ushort_t p=basisptr->p; 1747 // find inverse of a mod p 1748 #ifdef PRIMES32 1749 int j=invmodnoerr(modulo(*a._ZINTptr,p),p); 1750 // deltar[i]=((2*ulonglong(basis[i].sqrtmod))*j)%p; 1751 #else // PRIMES32 1752 unsigned modu=usqrta%p; 1753 modu=(modu*modu)%p; 1754 int j=invmodnoerr(modu,p); 1755 #endif // PRIMES32 1756 if (j<0) 1757 j += p; 1758 unsigned inva=j; 1759 #ifdef WITH_INVA 1760 Inva[i]=inva; 1761 #else 1762 #ifdef PRIMES32 1763 // set roots change values for all b coeffs for this a 1764 if (afact>afact0){ 1765 int * ptr=&bainv2[i]; 1766 for (int j=1;j<afact;ptr+=bs,++j){ 1767 // PREFETCH(ptr+bs); 1768 *ptr=modmult(modulo(*bvalues[j]._ZINTptr,p),2*inva,p); 1769 } 1770 } 1771 #endif // PRIMES32 1772 #endif // WITH_INVA 1773 // compute roots mod p 1774 #ifdef SQRTMOD_OUTSIDE 1775 ushort_t sqrtm=*sqrtmodit; 1776 ++sqrtmodit; 1777 #else 1778 ushort_t sqrtm=basisptr->sqrtmod; 1779 #endif 1780 int bmodp=p-modulo(*b._ZINTptr,p); 1781 if (inva){ 1782 if (p<=37000){ 1783 // sqrtm<=p/2, bmodp<p, inva<p hence (bmodp+sqrtm)*inva<=(3p/2-1)*(p-1) 1784 // this leaves M up to about 203 millions 1785 basisptr->root1=(M+(bmodp+sqrtm)*inva) % p; 1786 basisptr->root2=(M+(bmodp+p-sqrtm)*inva) % p; 1787 continue; 1788 } 1789 #ifdef _I386_ 1790 register int q=M; 1791 addmultmod(q,bmodp+sqrtm,inva,p); 1792 basisptr->root1=q; 1793 q=M; 1794 addmultmod(q,bmodp+p-sqrtm,inva,p); 1795 basisptr->root2=q; 1796 #else 1797 basisptr->root1=(M+longlong(bmodp+sqrtm)*inva) % p; 1798 basisptr->root2=(M+longlong(bmodp+p-sqrtm)*inva) % p; 1799 #endif 1800 continue; 1801 } 1802 int cmodp=modulo(zq,p); 1803 int q=(M+longlong(cmodp)*invmodnoerr((2*bmodp)%p,p))%p; 1804 if (q<0) 1805 q+=p; 1806 basisptr->root2=q; 1807 basisptr->root1=q; 1808 } 1809 #ifdef WITH_INVA 1810 #ifdef PRIMES32 1811 if (afact>afact0){ 1812 int * bainv2ptr=&bainv2.front(); 1813 basis_t * basisptr,*basisend=&basis.front()+bs; 1814 for (int j=1;j<afact;++j){ 1815 if (bvalues[j].type==_INT_){ 1816 int bjj=bvalues[j].val; 1817 vector<ushort_t>::const_iterator invait=Inva.begin(); 1818 for (basisptr=&basis.front();basisptr<basisend;++invait,++bainv2ptr,++basisptr){ 1819 register int r=(bjj*longlong(2*(*invait))) % basisptr->p; 1820 if (r<0) 1821 r += basisptr->p; 1822 *bainv2ptr=r; 1823 } 1824 } 1825 else { 1826 // longlong up1=up1tmp[2*j]; 1827 // longlong tmp=up1tmp[2*j+1]; 1828 // tmp is <= P^2 where P is the largest factor of a 1829 mpz_t & bz=*bvalues[j]._ZINTptr; 1830 vector<ushort_t>::const_iterator invait=Inva.begin(); 1831 for (basisptr=&basis.front();basisptr<basisend;++invait,++bainv2ptr,++basisptr){ 1832 register int p=basisptr->p; 1833 *bainv2ptr=((modulo(bz,p))*longlong(2*(*invait))) % p; 1834 } 1835 } 1836 } 1837 } 1838 #endif // PRIMES32 1839 #endif // WITH_INVA 1840 1841 #ifdef LP_SMALL_PRIMES // copy primes<2^16 into small_basis 1842 copy(basis,small_basis); 1843 #endif 1844 } 1845 1846 // find relations using (a*x+b)^2=a*(a*x^2+b*x+c) mod n where 1847 // we sieve on [-M,M] for as many polynomials as required 1848 // a is a square, approx sqrt(2*n)/M, and n is a square modulo all primes dividing a 1849 // b satisifies b^2=n mod a (b in [0,a[) 1850 // c=(n-b^2)/a msieve(const gen & n_orig,gen & pn,GIAC_CONTEXT)1851 bool msieve(const gen & n_orig,gen & pn,GIAC_CONTEXT){ 1852 if (n_orig.type!=_ZINT) 1853 return false; 1854 // find multiplier 1855 double delta; 1856 int multiplier=find_multiplier(n_orig,delta,contextptr); 1857 gen N(multiplier*n_orig); 1858 double Nd=evalf_double(N,1,contextptr)._DOUBLE_val; 1859 #if defined RTOS_THREADX || defined NSPIRE 1860 if (Nd>1e40) return false; 1861 #endif 1862 #ifdef BESTA_OS 1863 if (Nd>1e40) return false; 1864 #endif 1865 #ifdef PRIMES32 1866 if (Nd>1e76) return false; 1867 #else 1868 if (Nd>1e63) return false; 1869 #endif 1870 int Ndl=int(std::log10(Nd)-std::log10(double(multiplier))+.5); // +2*delta); 1871 #ifdef LP_TAB_SIZE 1872 int slicesize=(1 << LP_TAB_SIZE); 1873 #else 1874 int slicesize=(QS_SIZE>=65536 && Ndl<61)?32768:QS_SIZE; 1875 #endif 1876 double B=std::exp(std::sqrt(2.0)/4*std::sqrt(std::log(Nd)*std::log(std::log(Nd))))*0.45; 1877 if (B<200) B=200; 1878 int pos1=70,pos0=23,afact=2,afixed=0; // pos position in the basis, afact number of factors 1879 // FIXME Will always include the 3 first primes of the basis 1880 // set a larger Mtarget gives less polynomials but also use less memory 1881 #if defined(RTOS_THREADX) || defined(RTOS_THREADX) || defined NSPIRE 1882 double Mtarget=0.95e5; 1883 if (Nd>1e36) 1884 Mtarget=1.2e5; 1885 #else 1886 double Mtarget=0.55e5; 1887 #ifndef USE_MORE_PRIMES // FIXME improve! in fact use more primes on Core, less on Opteron 1888 if (Ndl>=50){ 1889 Ndl-=50; 1890 short int Btab[]={ 1891 // 50 1892 1900,2100,2300,2500,2700,2900,3100,3400,3700,4000, 1893 // 60 1894 4300,4600,4900,5300,5700,6200,6800,7500,8300,9200, 1895 // 70 1896 10000,11000,12000,13000,14000,15000,16000 1897 }; 1898 if (Ndl<sizeof(Btab)/sizeof(short int)) 1899 B=Btab[Ndl]; 1900 Mtarget=0.66e5; 1901 if (Ndl>7) 1902 Mtarget=0.95e5; 1903 if (Ndl>11) 1904 Mtarget=1.3e5; 1905 if (Ndl>15) 1906 Mtarget=1.6e5; 1907 if (Ndl>19) 1908 Mtarget=1.92e5; 1909 } 1910 #else 1911 if (Ndl>=50) 1912 Mtarget=0.85e5; 1913 if (Ndl>65) 1914 Mtarget=1.3e5; 1915 #endif 1916 #endif 1917 if (debug_infolevel) 1918 *logptr(contextptr) << " " << CLOCK() << gettext(" sieve on ") << N << '\n' << gettext("Number of primes ") << B << '\n'; 1919 // first compute the prime basis and sqrt(N) mod p, p in basis 1920 vector<basis_t> basis; 1921 basis.reserve(unsigned(B)); 1922 #ifdef SQRTMOD_OUTSIDE 1923 vector<ushort_t> sqrtmod; 1924 sqrtmod.reserve(basis.capacity()); 1925 basis.push_back(2); 1926 sqrtmod.push_back(1); 1927 #else 1928 basis.push_back(basis_t(2,1)); // I assume that N is odd... hence has sqrt 1 mod 2 1929 #endif 1930 N.uncoerce(); 1931 // vector<ushort_t> N256; 1932 int i; 1933 mpz_t zx,zy,zq,zr; 1934 mpz_init(zx); mpz_init(zy); mpz_init(zq); mpz_init(zr); 1935 // fastsmod_prepare(N,zx,zy,zr,N256); 1936 for (i=1;i<int(sizeof(giac_primes)/sizeof(short));++i){ 1937 #ifdef TIMEOUT 1938 control_c(); 1939 #endif 1940 if (ctrl_c || interrupted) 1941 break; 1942 ushort_t j=giac_primes[i]; 1943 if (debug_infolevel>6 && (i%500==99)) 1944 *logptr(contextptr) << CLOCK() << gettext(" sieve current basis size ") << basis.size() << '\n'; 1945 #if 1 // def USE_GMP_REPLACEMENTS 1946 // int n=fastsmod_compute(N256,j); 1947 int n=modulo(*N._ZINTptr,j),s; 1948 #else 1949 int n=smod(N,j).val,s; 1950 #endif 1951 if (n<0) 1952 n+=j; 1953 if (n==0){ 1954 #ifdef SQRTMOD_OUTSIDE 1955 basis.push_back(j); 1956 sqrtmod.push_back(0); 1957 #else 1958 basis.push_back(basis_t(j,0)); 1959 #endif 1960 } 1961 else { 1962 if (powmod(n,(unsigned long)((j-1)/2),(int)j)==1){ 1963 s=sqrt_mod(n,int(j),true,contextptr).val; 1964 if (s<0) 1965 s+=j; 1966 #ifdef SQRTMOD_OUTSIDE 1967 basis.push_back(j); 1968 sqrtmod.push_back(s); 1969 #else 1970 basis.push_back(basis_t(j,s)); 1971 #endif 1972 } 1973 } 1974 if (basis.size()>=B) 1975 break; 1976 } 1977 vector<unsigned> crible; 1978 int jp=0; 1979 if (basis.size()<B){ 1980 #ifdef PRIME_SIEVE 1981 fill_crible(crible,int(2.5*B*std::log(B))); 1982 jp=nextprime(crible,basis.back().p+1); 1983 #else 1984 jp=nextprime(int(basis.back().p+1)).val; 1985 #endif 1986 } 1987 unsigned lp_basis_pos=0; // position of first prime > 2^16 in the basis 1988 for (;basis.size()<B;++i){ 1989 #ifdef TIMEOUT 1990 control_c(); 1991 #endif 1992 if (ctrl_c || interrupted) 1993 break; 1994 #ifndef PRIMES32 1995 if (jp>65535){ 1996 break; 1997 } 1998 #endif 1999 if (debug_infolevel>6 && (i%500==99)) 2000 *logptr(contextptr) << CLOCK() << gettext(" sieve current basis size ") << basis.size() << '\n'; 2001 #if 1 // def USE_GMP_REPLACEMENTS 2002 // int n=fastsmod_compute(N256,jp); 2003 int n=modulo(*N._ZINTptr,jp),s; 2004 #else 2005 int n=smod(N,jp).val,s; 2006 #endif 2007 if (n<0) 2008 n+=jp; 2009 if (powmod(n,(unsigned long)((jp-1)/2),jp)==1){ 2010 s=sqrt_mod(n,jp,true,contextptr).val; 2011 if (s<0) 2012 s += jp; 2013 #ifdef LP_TAB_SIZE 2014 if (!lp_basis_pos && jp> (1<<LP_BIT_LIMIT)) 2015 lp_basis_pos=basis.size(); 2016 #endif 2017 #ifdef SQRTMOD_OUTSIDE 2018 basis.push_back(jp); 2019 sqrtmod.push_back(s); 2020 #else 2021 basis.push_back(basis_t(jp,s)); 2022 #endif 2023 } 2024 #ifdef PRIME_SIEVE 2025 jp=nextprime(crible,jp+1); 2026 #else 2027 jp=nextprime(int(jp+1)).val; 2028 #endif 2029 } 2030 if (!lp_basis_pos) 2031 lp_basis_pos=unsigned(basis.size()); 2032 #ifdef LP_SMALL_PRIMES 2033 vector<small_basis_t> small_basis(lp_basis_pos); // will be filled by primes<2^16 2034 #endif 2035 #ifdef TIMEOUT 2036 control_c(); 2037 #endif 2038 if (ctrl_c || interrupted){ 2039 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2040 return false; 2041 } 2042 double dtarget=1.1; 2043 if (Mtarget<basis.back().p*dtarget){ 2044 Mtarget=basis.back().p*dtarget; // (int(basis.back().p*1.1)/slicesize)*slicesize; 2045 } 2046 unsigned ps=sizeinbase2(basis.back().p); 2047 #if !defined(RTOS_THREADX) && !defined(BESTA_OS) && !defined NSPIRE // def USE_MORE_PRIMES 2048 unsigned maxadditional=(2+(basis.back().p>>16))*basis.back().p*ps; 2049 #else 2050 unsigned maxadditional=3*basis.back().p*ps; 2051 #endif 2052 if (debug_infolevel) 2053 *logptr(contextptr) << CLOCK() << gettext(" sieve basis OK, size ") << basis.size() << " largest prime in basis " << basis.back().p << " large prime " << maxadditional << " Mtarget " << Mtarget << '\n' ; 2054 int bs=int(basis.size()); 2055 gen isqrtN=isqrt(N); 2056 isqrtN.uncoerce(); 2057 // now compare isqrtN to a^2 for a in the basis 2058 double seuil=1.414*evalf_double(isqrtN,1,contextptr)._DOUBLE_val/Mtarget; // should be a 2059 seuil=std::sqrt(seuil); // should be product of primes of the basis 2060 #ifdef OLD_AFACT 2061 double dfactors=std::log10(seuil)/3; 2062 // fixed primes are choosen at basis[pos0], variables are choosen around 2000 2063 afact=int(dfactors+.5); 2064 if (afact<=1){ 2065 afact=1; 2066 int i=20; 2067 for (;i<3*bs/4;++i){ 2068 if (seuil<basis[i].p){ 2069 pos1=i; 2070 break; 2071 } 2072 } 2073 if (i>=3*bs/4){ 2074 afact=2; 2075 for (;i<3*bs/4;++i){ 2076 if (seuil<basis[i].p){ 2077 pos1=i; 2078 break; 2079 } 2080 } 2081 } 2082 } 2083 else { 2084 if (afact==2){ 2085 seuil=std::sqrt(seuil); 2086 for (int i=20;i<3*bs/4;++i){ 2087 if (seuil<basis[i].p){ 2088 pos1=i; 2089 break; 2090 } 2091 } 2092 } 2093 else { // afact>=3, 2094 if (dfactors>5.4){ 2095 dfactors -= 3; // 3 large primes 2096 afixed = dfactors/.8; // at least 3 fixed 2097 afact = 3 + afixed; 2098 } 2099 else { 2100 dfactors -= 2; // 2 large primes 2101 afixed = dfactors/.8; 2102 if (afixed==0) 2103 afixed=1; 2104 afact = 2 +afixed; 2105 } 2106 for (int i=0;i<afixed;++i){ 2107 seuil=seuil/basis[pos0+i].p; 2108 } 2109 seuil=std::pow(seuil,1./(afact-afixed)); 2110 for (int i=pos0+afixed+10;i<3*bs/4;++i){ 2111 if (seuil<basis[i].p){ 2112 pos1=i; 2113 break; 2114 } 2115 } 2116 } 2117 } 2118 #else // OLD_AFACT 2119 double logprod=std::log10(seuil); 2120 if (logprod<9){ 2121 afixed=0; 2122 afact=int(std::ceil(logprod/3)); 2123 if (logprod-3*afact<-1) 2124 afixed=1; 2125 } 2126 else { 2127 double logfixed=std::log10(double(basis[pos0].p)); 2128 int maxfixed=int(logprod/logfixed)-2; 2129 if (maxfixed==0) maxfixed=1; 2130 double curseuil=1e10; 2131 afact=0; 2132 for (int i=1;i<=maxfixed;++i){ 2133 double variable=(logprod-i*logfixed)/3; // we want variable primes to be around 1000 2134 int ivariable=int(variable); 2135 double seuiltest=variable/ivariable; 2136 if (i+ivariable>afact){ 2137 afixed=i; 2138 afact=i+ivariable; 2139 curseuil=seuiltest; 2140 } 2141 } 2142 } 2143 for (int i=0;i<afixed;++i) 2144 seuil=seuil/basis[pos0+i].p; 2145 seuil=std::pow(seuil,1./(afact-afixed)); 2146 for (int i=pos0+afixed+10;i<3*bs/4;++i){ 2147 if (seuil<basis[i].p){ 2148 pos1=i-1; // -afact/2; 2149 break; 2150 } 2151 } 2152 #endif // OLD_AFACT 2153 if (debug_infolevel){ 2154 *logptr(contextptr) << gettext("Using ") << afact << " square factors per a coefficient in polynomials" << '\n'; 2155 *logptr(contextptr) << afixed << gettext(" fixed begin at ") << basis[pos0].p << " and " << afact-afixed << " variables at " << basis[pos1].p << '\n'; 2156 } 2157 vector<ushort_t> isqrtN256; 2158 // fastsmod_prepare(isqrtN,zx,zy,zr,isqrtN256); 2159 vector<short_t> isqrtNmodp(bs); 2160 for (int i=0;i<bs;++i){ 2161 #if 1 2162 // isqrtNmodp[i]=fastsmod_compute(isqrtN256,basis[i].p); 2163 isqrtNmodp[i]=modulo(*isqrtN._ZINTptr,basis[i].p); 2164 #else 2165 isqrtNmodp[i]=smod(isqrtN,basis[i].p).val; 2166 #endif 2167 } 2168 #if defined(RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 2169 unsigned puissancestablength=10000; 2170 #else 2171 #if 1 // def USE_MORE_PRIMES 2172 unsigned puissancestablength=bs*(80+bs/500); 2173 #else 2174 unsigned puissancestablength=bs*80; 2175 #endif 2176 #endif 2177 ushort_t * puissancestab=new ushort_t[puissancestablength]; 2178 ushort_t * puissancesptr=puissancestab; 2179 ushort_t * puissancesend=puissancestab+puissancestablength; 2180 slicetype * slice=0; 2181 if (puissancestab) 2182 slice=new slicetype[QS_SIZE]; 2183 if (!slice){ 2184 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2185 return false; 2186 } 2187 // relations will be written in column 2188 vector<axbinv> axbmodn; // contains (sqrta,b,x) 2189 vector<additional_t> additional_primes; 2190 #ifndef ADDITIONAL_PRIMES_HASHMAP 2191 vector<bool> additional_primes_twice; 2192 #endif 2193 #ifdef LP_TAB_SIZE 2194 vector<lp_tab_t> lp_map(128); // at most 128 slices in a sieve 2195 #endif 2196 vecteur sqrtavals,bvals; 2197 #ifdef GIAC_ADDITIONAL_PRIMES 2198 #ifdef ADDITIONAL_PRIMES_HASHMAP 2199 #ifdef EMCC 2200 additional_map_t additional_primes_map; 2201 #else 2202 additional_map_t additional_primes_map(8*bs); 2203 #endif 2204 axbmodn.reserve(bs); 2205 #else 2206 #if defined(RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 2207 additional_primes.reserve(bs); 2208 additional_primes_twice.reserve(bs); 2209 axbmodn.reserve(2*bs); 2210 #else 2211 additional_primes.reserve(4*bs); 2212 additional_primes_twice.reserve(4*bs); 2213 axbmodn.reserve(5*bs); 2214 #endif 2215 sqrtavals.reserve(bs/7); 2216 bvals.reserve(2*bs/7); 2217 #endif // ADDITIONAL_PRIMES_HASHMAP 2218 #else // GIAC_ADDITIONAL_PRIMES 2219 axbmodn.reserve(bs+1); 2220 #endif 2221 // now sieve 2222 unsigned todo_rel; 2223 unsigned marge=bs/100; 2224 if (marge<15) 2225 marge=15; 2226 mpz_t alloc1,alloc2,alloc3,alloc4,alloc5; 2227 mpz_init(alloc1); mpz_init(alloc2); mpz_init(alloc3); mpz_init(alloc4); mpz_init(alloc5); 2228 // vector<ushort_t> a256,b256,tmpv; 2229 vector<ushort_t> curpuissances,recheck,pos(afact); 2230 #ifdef WITH_INVA 2231 vector<ushort_t> Inva(bs); 2232 #endif 2233 vecteur bvalues; // will contain values of b if afact<=afact0 or components of b if afact>afact0 2234 // array for efficient polynomial switch (same a change b) when at least afact0 factors/a 2235 #ifdef PRIMES32 2236 const int afact0=3; 2237 vector<int> bainv2((afact-1)*bs); 2238 vector<longlong> up1tmp; 2239 #endif 2240 for (int i=0;i<afixed;++i) 2241 pos[i]=pos0+i; 2242 for (int i=afixed;i<afact;++i) 2243 pos[i]=pos1+i-afixed; // FIXME should be -afixed 2244 double Mval=1; 2245 for (int i=0;i<afact;++i) 2246 Mval=Mval*basis[pos[i]].p; 2247 Mval=std::sqrt(2*Nd)/(Mval*Mval); 2248 if (debug_infolevel) 2249 *logptr(contextptr) << gettext("First M ") << Mval << '\n'; 2250 Mtarget=Mval; 2251 int avar=afact-afixed; 2252 int end_pos1=2*pos1; 2253 if (avar>1) 2254 end_pos1=pos1+100; 2255 if (avar>2) 2256 end_pos1=pos1+30; 2257 if (int(lp_basis_pos)<end_pos1) 2258 end_pos1=lp_basis_pos; 2259 for (;puissancesptr<puissancesend;++pos.back()){ 2260 double bpos2=1; 2261 if (int(pos.back())>=end_pos1 || basis[pos.back()].p>=45000){ 2262 int i=afact-2; 2263 for (;i>afixed;--i){ 2264 if (int(pos[i])<end_pos1-(afact-i)){ 2265 ++pos[i]; 2266 for (int j=i+1;j<afact;++j) 2267 pos[j]=pos[i]+(j-i); 2268 break; 2269 } 2270 } 2271 if (i<=afixed){ 2272 --pos1; 2273 if (pos1<=5){ 2274 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2275 mpz_clear(alloc1); mpz_clear(alloc2); mpz_clear(alloc3); mpz_clear(alloc4); mpz_clear(alloc5); 2276 delete [] puissancestab; 2277 return false; 2278 } 2279 // reset fixed factors 2280 for (i=0;i<afixed;++i) 2281 pos[i]=pos0+i; 2282 for (i=afixed;i<afact;++i) 2283 pos[i]=pos1+i-afixed; 2284 } 2285 } 2286 for (int i=0;i<afact;++i) 2287 bpos2=bpos2*basis[pos[i]].p; 2288 Mval=std::sqrt(2*Nd)/(bpos2*bpos2); 2289 if (afixed){ 2290 // move "fixed" factors so that Mval becomes closer to Mtarget 2291 // NB: for later threads, afixed should be >=3, so that we can move pos[1] by thread 2292 while (Mval>1.1*Mtarget && int(pos[afixed-1])<int(pos1)-10){ 2293 // Mval is too large, hence one factor of a is too small, increase it 2294 if (int(pos[0])<pos0){ 2295 ++pos[0]; 2296 double coeff=basis[pos[0]].p/double(basis[pos[0]-1].p); 2297 Mval=Mval/(coeff*coeff); 2298 } 2299 else { 2300 ++pos[afixed-1]; 2301 double coeff=basis[pos[afixed-1]].p/double(basis[pos[afixed-1]-1].p); 2302 Mval=Mval/(coeff*coeff); 2303 } 2304 } 2305 while (Mval<0.9*Mtarget && pos[0]>10){ 2306 // Mval is too small, decrease one factor of a 2307 --pos[0]; 2308 if (pos[0]<=10){ 2309 Mval=0; 2310 break; 2311 } 2312 double coeff=basis[pos[0]].p/double(basis[pos[0]+1].p); 2313 Mval=Mval/(coeff*coeff); 2314 } 2315 } 2316 if ( Mval <0.7*Mtarget ){ 2317 if (pos1>pos0+afixed+5 || Mval<32768){ 2318 // CERR << pos ; 2319 int i=afact-1; 2320 for (;i>afixed+1;--i){ 2321 if (pos[i]>pos[i-1]+5) 2322 break; 2323 } 2324 if (i<=afixed+1){ 2325 --pos1; 2326 for (i=0;i<afixed;++i) 2327 pos[i]=pos0+i; 2328 for (i=afixed;i<afact;++i) 2329 pos[i]=pos1+i-afixed; 2330 } 2331 else { 2332 ++pos[i-1]; 2333 for (;i<afact;++i) 2334 pos[i]=pos[i-1]+1; 2335 } 2336 // CERR << pos << '\n'; 2337 } 2338 } 2339 // finished? 2340 #ifdef TIMEOUT 2341 control_c(); 2342 #endif 2343 if (ctrl_c || interrupted) 2344 break; 2345 #ifdef ADDITIONAL_PRIMES_HASHMAP 2346 todo_rel=bs+marge; 2347 #else 2348 todo_rel=bs+marge+unsigned(additional_primes.size()); 2349 #endif 2350 if (axbmodn.size()>=todo_rel) 2351 break; 2352 int nrelationsa=0; 2353 // Not finished yet, construct a new value of a around ad=sqrt(2*n)/M 2354 // using a product of afact square of primes that are in the basis 2355 // and construct a vector of 2^(afact-1) corresponding values of b 2356 // and compute the values of inverses of a mod p 2357 ulonglong usqrta(basis[pos[0]].p); 2358 for (int i=1;i<afact;++i) 2359 usqrta=basis[pos[i]].p*usqrta; // works up to about N=1e86 2360 gen sqrta((int) basis[pos[0]].p); 2361 for (int i=1;i<afact;++i) 2362 sqrta=gen(int(basis[pos[i]].p))*sqrta; 2363 sqrtavals.push_back(sqrta); 2364 gen a=sqrta*sqrta; // a should be about sqrt(Nd/2)/M 2365 a.uncoerce(); 2366 int M=int(std::floor(std::sqrt(Nd*2)/evalf_double(a,1,contextptr)._DOUBLE_val)); 2367 if (debug_infolevel>6) 2368 *logptr(contextptr) << CLOCK() << gettext(" initial value for M= ") << M << '\n'; 2369 int nslices=int(std::ceil((2.*M)/slicesize)); 2370 M=(nslices*slicesize)/2; 2371 bvalues.clear(); 2372 gen curprod=1; 2373 for (int i=0;;){ 2374 #ifdef SQRTMOD_OUTSIDE 2375 int s=sqrtmod[pos[i]]; 2376 #else 2377 int s=basis[pos[i]].sqrtmod; 2378 #endif 2379 int p=basis[pos[i]].p; 2380 longlong p2=p*longlong(p); 2381 // Hensel lift s to be a sqrt of n mod p^2: (s+p*r)^2=s^2+2p*r*s=n => r=(n-s^2)/p*inv(2*s mod p) 2382 int r=p<37000?int((modulo(*N._ZINTptr,p2)-s*s)/p):((smod(N,p2)-s*s)/p).val; 2383 r=(r*invmod(2*s,p))%p; 2384 // overflow should not happen because p is a factor of a hence choosen 2385 // in the 1000 range (perhaps up to 10 000, but not much larger) 2386 // if ((longlong(r)*p)!=r*p) CERR << "overflow" << '\n'; 2387 s += p*r; 2388 #ifdef PRIMES32 2389 if (afact>afact0){ 2390 // store s*(inv( product(basis[pos[j]]^2,j!=i) mod p2)) in bvalues[i] 2391 longlong up1=(usqrta/p); 2392 longlong up=up1%p2; 2393 up=(up*up)%p2; 2394 //longlong tmp=(s*invmodnoerr(up,p2))%p2; 2395 //up1tmp.push_back(up1); 2396 //up1tmp.push_back(tmp); 2397 gen tmp=smod(s*invmod(gen(up),gen(p2)),p2); 2398 if (is_greater(0,tmp,contextptr)) tmp=-tmp; 2399 gen gup1(up1); 2400 bvalues.push_back(gup1*gup1*tmp); 2401 bvalues.back().uncoerce(); 2402 } 2403 else 2404 #endif 2405 { 2406 if (bvalues.empty()) 2407 bvalues.push_back(s); 2408 else { 2409 int js=int(bvalues.size()); 2410 for (int j=0;j<js;++j){ 2411 bvalues.push_back(ichinrem(bvalues[j],-s,curprod,p2)); 2412 bvalues[j]=ichinrem(bvalues[j],s,curprod,p2); 2413 } 2414 } 2415 } 2416 ++i; 2417 if (i==afact) 2418 break; 2419 curprod = p*(p*curprod); 2420 } // end for 2421 // compute inverse of a modulo p (will set to 0 if not invertible) 2422 if (debug_infolevel>6) 2423 *logptr(contextptr) << CLOCK() << gettext(" Computing inverses mod p of the basis ") << '\n'; 2424 // fastsmod_prepare(a,zx,zy,zr,a256); 2425 gen b; 2426 for (int i=0;i< (1<<(afact-1));++i){ 2427 #ifdef TIMEOUT 2428 control_c(); 2429 #endif 2430 if (ctrl_c || interrupted) 2431 break; 2432 #ifdef ADDITIONAL_PRIMES_HASHMAP 2433 todo_rel=bs+marge; 2434 #else 2435 todo_rel=bs+marge+unsigned(additional_primes.size()); 2436 #endif 2437 if (axbmodn.size()>=todo_rel) 2438 break; 2439 if (debug_infolevel>6) 2440 *logptr(contextptr) << CLOCK() << gettext(" Computing c ") << '\n'; 2441 #ifdef PRIMES32 2442 int bv=1,be=-1; 2443 if (afact>afact0){ 2444 if (i==0){ 2445 b=0; 2446 for (unsigned j=0;j<bvalues.size();j++) 2447 b += bvalues[j]; 2448 } 2449 else { 2450 find_bv_be(i,bv,be); 2451 b += (2*be)*bvalues[bv]; 2452 } 2453 } 2454 else 2455 b=bvalues[i]; 2456 #else 2457 b=bvalues[i]; 2458 #endif 2459 b.uncoerce(); 2460 mpz_mul(zx,*b._ZINTptr,*b._ZINTptr); 2461 mpz_sub(zy,zx,*N._ZINTptr); 2462 #ifdef USE_GMP_REPLACEMENTS 2463 mp_grow(&alloc1,zy.used+2); 2464 mpz_set_ui(alloc1,0); 2465 alloc1.used = zy.used +2 ; 2466 mpz_set(alloc2,zy); 2467 mpz_set(alloc3,*a._ZINTptr); 2468 alloc_mp_div(&zy,a._ZINTptr,&zq,&zr,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 2469 #else 2470 mpz_divexact(zq,zy,*a._ZINTptr); 2471 #endif 2472 // gen c=zq; // gen c=(b*b-N)/a; 2473 // c.uncoerce(); 2474 #ifdef PRIMES32 2475 if (afact<=afact0) 2476 #endif 2477 { 2478 bool bneg=mpz_cmp_ui(*b._ZINTptr,0)<0; 2479 if (bneg) 2480 mpz_neg(*b._ZINTptr,*b._ZINTptr); 2481 } 2482 bvals.push_back(b); 2483 if (debug_infolevel>6) 2484 *logptr(contextptr) << CLOCK() << gettext(" Computing roots mod the basis ") << '\n'; 2485 // fastsmod_prepare(b,zx,zy,zr,b256); 2486 #ifdef PRIMES32 2487 if (i && afact>afact0) 2488 switch_roots(bainv2,basis, 2489 #ifdef LP_SMALL_PRIMES 2490 small_basis, 2491 #endif 2492 lp_basis_pos,nslices,slicesize,bv,be,afact,pos,b,zq,M); 2493 else { 2494 init_roots(basis, 2495 #ifdef LP_SMALL_PRIMES 2496 small_basis, 2497 #endif 2498 #ifdef WITH_INVA 2499 Inva, 2500 #endif 2501 #ifdef SQRTMOD_OUTSIDE 2502 sqrtmod, 2503 #endif 2504 bainv2,afact,afact0, 2505 a,b,bvalues,zq,M); 2506 #ifdef LP_TAB_TOGETHER 2507 // init all hashtable for large primes at once 2508 unsigned cl; 2509 if (debug_infolevel>3){ 2510 cl=CLOCK(); 2511 *logptr(contextptr) << cl << gettext(" Init large prime hashtables ") << '\n'; 2512 } 2513 int total=(nslices << (afact-1)); 2514 if (int(lp_map.size()) < total) 2515 lp_map.resize(total); 2516 for (int k=0;k< total;++k) 2517 lp_map[k].clear(); 2518 if (lp_basis_pos){ 2519 for (int k=0;;){ 2520 basis_t * bit=&basis[0]+lp_basis_pos, * bitend=&basis[0]+bs; 2521 unsigned endpos=nslices*slicesize; 2522 lp_tab_t * ptr=&lp_map[0]+k*nslices; 2523 for (;bit!=bitend;++bit){ 2524 register ushort_t p=bit->p; 2525 register unsigned pos=bit->root1; 2526 for (;pos<endpos; pos += p){ 2527 (ptr+(pos >> LP_TAB_SIZE))->push_back(lp_entry_t((pos & LP_MASK),p)); 2528 } 2529 pos=bit->root2; 2530 for (;pos<endpos; pos += p){ 2531 (ptr+(pos >> LP_TAB_SIZE))->push_back(lp_entry_t((pos & LP_MASK),p)); 2532 } 2533 } 2534 ++k; 2535 if (k== (1 << (afact-1))){ 2536 if (debug_infolevel>3){ 2537 unsigned cl2=CLOCK(); 2538 *logptr(contextptr) << cl2 << gettext(" End large prime hashtables ") << cl2-cl << '\n'; 2539 } 2540 break; 2541 } 2542 find_bv_be(k,bv,be); 2543 // switch roots to next polynomial 2544 int * bvpos=&bainv2[(bv-1)*bs],* bvposend=bvpos+bs; 2545 bvpos += lp_basis_pos; 2546 basis_t * basisptr=&basis[0]+lp_basis_pos; 2547 if (be>0){ 2548 for (;bvpos<bvposend;++basisptr,++bvpos){ 2549 register unsigned p=basisptr->p; 2550 register int r=basisptr->root1-(*bvpos); 2551 if (r<0) 2552 r+=p; 2553 basisptr->root1=r; 2554 r=basisptr->root2-(*bvpos); 2555 if (r<0) 2556 r+=p; 2557 basisptr->root2=r; 2558 } 2559 } 2560 else { 2561 for (;bvpos<bvposend;++basisptr,++bvpos){ 2562 register unsigned p=basisptr->p; 2563 register int r=basisptr->root1+(*bvpos); 2564 if (r>int(p)) 2565 r-=p; 2566 basisptr->root1=r; 2567 r=basisptr->root2+(*bvpos); 2568 if (r>int(p)) 2569 r-=p; 2570 basisptr->root2=r; 2571 } 2572 } 2573 } 2574 } 2575 #endif // LP_TAB_TOGETHER 2576 } // end else of if i==0 2577 #if defined(LP_TAB_SIZE) && !defined(LP_TAB_TOGETHER) 2578 if (int(lp_map.size()) < nslices) 2579 lp_map.resize(nslices); 2580 for (int k=0;k< nslices;++k) 2581 lp_map[k].clear(); 2582 if (lp_basis_pos){ 2583 basis_t * bit=&basis[0]+lp_basis_pos, * bitend=&basis[0]+bs; 2584 unsigned endpos=nslices*slicesize; 2585 for (;bit!=bitend;++bit){ 2586 register ushort_t p=bit->p; 2587 register unsigned pos=bit->root1; 2588 for (;pos<endpos; pos += p){ 2589 lp_map[pos >> LP_TAB_SIZE].push_back(lp_entry_t((pos & LP_MASK),p)); 2590 } 2591 pos=bit->root2; 2592 for (;pos<endpos; pos += p){ 2593 lp_map[pos >> LP_TAB_SIZE].push_back(lp_entry_t((pos & LP_MASK),p)); 2594 } 2595 } 2596 } 2597 #endif // LP_TAB_SIZE && !LP_TAB_TOGETHER 2598 #else // PRIMES32 2599 init_roots(basis, 2600 #ifdef WITH_INVA 2601 Inva, 2602 #endif 2603 #ifdef SQRTMOD_OUTSIDE 2604 sqrtmod, 2605 #endif 2606 usqrta,a,b,bvalues,zq,M); 2607 #endif // PRIMES32 2608 // we can now sieve in [-M,M[ by slice of size slicesize 2609 #ifndef GIAC_HAS_STO_38 2610 if (debug_infolevel>5){ 2611 *logptr(contextptr) << CLOCK(); 2612 *logptr(contextptr) << gettext(" Polynomial a,b,M=") << a << "," << b << "," << M << " (" << pos << ")" ; 2613 *logptr(contextptr) << CLOCK() << '\n'; 2614 } 2615 #endif 2616 int nrelationsb=0; 2617 #ifdef LP_TAB_SIZE 2618 #endif 2619 for (int l=0;l<nslices;l++){ 2620 #ifdef TIMEOUT 2621 control_c(); 2622 #endif 2623 if (ctrl_c || interrupted) 2624 break; 2625 #ifdef ADDITIONAL_PRIMES_HASHMAP 2626 todo_rel=bs+marge; 2627 #else 2628 todo_rel=bs+marge+unsigned(additional_primes.size()); 2629 #endif 2630 if (axbmodn.size()>=todo_rel) 2631 break; 2632 int shift=-M+l*slicesize; 2633 int slicerelations=msieve(a,sqrtavals, 2634 bvals,zq,basis,lp_basis_pos, 2635 #ifdef LP_SMALL_PRIMES 2636 small_basis, 2637 #endif 2638 maxadditional, 2639 #ifdef ADDITIONAL_PRIMES_HASHMAP 2640 additional_primes_map, 2641 #else 2642 additional_primes,additional_primes_twice, 2643 #endif 2644 N,isqrtN, 2645 slice,slicesize,shift,puissancestab,puissancesptr,puissancesend,curpuissances,recheck, 2646 axbmodn, 2647 zx,zy,zr,alloc1,alloc2,alloc3,alloc4,alloc5, 2648 #ifdef LP_TAB_SIZE 2649 #ifdef LP_TAB_TOGETHER 2650 lp_map[l+nslices*i], 2651 #else 2652 lp_map[l], 2653 #endif 2654 #endif 2655 contextptr); 2656 if (slicerelations==-1){ 2657 *logptr(contextptr) << gettext("Sieve error: Not enough memory ") << '\n'; 2658 break; 2659 } 2660 nrelationsb += slicerelations; 2661 #ifdef ADDITIONAL_PRIMES_HASHMAP 2662 todo_rel=bs+marge; 2663 #else 2664 todo_rel=bs+marge+unsigned(additional_primes.size()); 2665 #endif 2666 } 2667 if (nrelationsb==0) 2668 bvals.pop_back(); 2669 else 2670 nrelationsa += nrelationsb; 2671 } 2672 #if defined( RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 2673 if (debug_infolevel){ 2674 #ifdef NSPIRE 2675 static int count_print=0; 2676 ++count_print; 2677 if (count_print%4==0) 2678 #endif 2679 *logptr(contextptr) << axbmodn.size() << " of " << todo_rel << " (" << 100-100*(todo_rel-axbmodn.size())/double(bs+marge) << "%)" << '\n'; 2680 } 2681 #endif 2682 if (nrelationsa==0){ 2683 sqrtavals.pop_back(); 2684 } 2685 #if !defined(RTOS_THREADX) && !defined(BESTA_OS) && !defined NSPIRE 2686 if (debug_infolevel>1) 2687 *logptr(contextptr) << CLOCK()<< gettext(" sieved : ") << axbmodn.size() << " of " << todo_rel << " (" << 100-100*(todo_rel-axbmodn.size())/double(bs+marge) << "%), M=" << M << '\n'; 2688 #endif 2689 } // end sieve loop 2690 if (debug_infolevel) 2691 *logptr(contextptr) << gettext("Polynomials a,b in use: #a ") << sqrtavals.size() << " and #b " << bvals.size() << '\n'; 2692 delete [] slice; 2693 #ifdef TIMEOUT 2694 control_c(); 2695 #endif 2696 if (ctrl_c || interrupted || puissancesptr==puissancesend){ 2697 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2698 mpz_clear(alloc1); mpz_clear(alloc2); mpz_clear(alloc3); mpz_clear(alloc4); mpz_clear(alloc5); 2699 delete [] puissancestab; 2700 return false; 2701 } 2702 // We have enough relations, make matrix, reduce it then find x^2=y^2 mod n congruences 2703 if (debug_infolevel) 2704 *logptr(contextptr) << CLOCK() << gettext(" sieve done: used ") << (puissancesptr-puissancestab)*0.002 << " K for storing relations (of " << puissancestablength*0.002 << ")" << '\n'; 2705 release_memory(isqrtNmodp); 2706 #ifdef GIAC_ADDITIONAL_PRIMES 2707 #ifdef ADDITIONAL_PRIMES_HASHMAP 2708 additional_primes.reserve(axbmodn.size()); 2709 vector<axbinv>::const_iterator it=axbmodn.begin(),itend=axbmodn.end(); 2710 for (;it!=itend;++it) { 2711 unsigned u=largep(*it,puissancestab); 2712 if (u) 2713 additional_primes.push_back(u); 2714 } 2715 sort(additional_primes.begin(),additional_primes.end()); // for binary search later 2716 #else 2717 if (debug_infolevel) 2718 *logptr(contextptr) << CLOCK() << gettext(" removing additional primes") << '\n'; 2719 // remove relations with additional primes which are used only once 2720 int lastp=int(axbmodn.size())-1,lasta=int(additional_primes.size())-1; 2721 for (int i=0;i<=lastp;++i){ 2722 ushort_t * curbeg=puissancestab+axbmodn[i].first, * curend=puissancestab+axbmodn[i].second; 2723 bool done=false; 2724 for (;curbeg!=curend;++curbeg){ 2725 if (*curbeg==1) 2726 break; 2727 } 2728 if (curbeg==curend) 2729 continue; 2730 ++curbeg; 2731 additional_t u=*curbeg; 2732 #if GIAC_ADDITIONAL_PRIMES==32 && !defined(PRIMES32) 2733 u <<=16 ; 2734 ++curbeg; 2735 u += *curbeg; 2736 #endif 2737 int pos=_equalposcomp(additional_primes,u); 2738 if (!pos) 2739 continue; 2740 if (pos>lasta){ 2741 // *logptr(contextptr) << cur << '\n'; 2742 continue; 2743 } 2744 --pos; 2745 if (additional_primes_twice[pos]) 2746 continue; 2747 axbmodn[i]=axbmodn[lastp]; 2748 --lastp; 2749 additional_primes[pos]=additional_primes[lasta]; 2750 additional_primes_twice[pos]=additional_primes_twice[lasta]; 2751 --lasta; 2752 --i; // recheck at current index 2753 } 2754 axbmodn.resize(lastp+1); 2755 additional_primes.resize(lasta+1); 2756 if (debug_infolevel) 2757 *logptr(contextptr) << CLOCK() << gettext(" end removing additional primes") << '\n'; 2758 #endif // ADDTIONAL_PRIMES_HASHMAP 2759 #endif // GIAC_ADDITIONAL_PRIMES 2760 // Make relations matrix (currently dense, FIXME improve to sparse and Lanczos algorithm) 2761 int C32=int(std::ceil(axbmodn.size()/32./GIAC_RREF_UNROLL))*GIAC_RREF_UNROLL; 2762 unsigned * tab=new unsigned[axbmodn.size()*C32],*tabend=tab+axbmodn.size()*C32; 2763 if (!tab){ 2764 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2765 mpz_clear(alloc1); mpz_clear(alloc2); mpz_clear(alloc3); mpz_clear(alloc4); mpz_clear(alloc5); 2766 delete [] puissancestab; 2767 return false; 2768 } 2769 // init tab 2770 for (unsigned * ptr=tab;ptr!=tabend;++ptr) 2771 *ptr=0; 2772 int l32=C32*32; 2773 vector< line_t > relations(axbmodn.size()); 2774 for (unsigned i=0;i<axbmodn.size();++i){ 2775 relations[i].tab=tab+i*C32; 2776 } 2777 for (unsigned j=0;j<axbmodn.size();j++){ 2778 ushort_t * curpui=puissancestab+axbmodn[j].first, * curpuiend=puissancestab+axbmodn[j].second; 2779 add_relation(relations,j,curpui,curpuiend,basis,additional_primes); 2780 #ifdef ADDITIONAL_PRIMES_HASHMAP 2781 unsigned u=largep(axbmodn[j],puissancestab); 2782 if (u){ 2783 axbinv & A=additional_primes_map[u]; 2784 curpui=puissancestab+A.first; curpuiend=puissancestab+A.second; 2785 add_relation(relations,j,curpui,curpuiend,basis,additional_primes); 2786 } 2787 #endif 2788 } // end loop on j in puissances 2789 #ifdef RREF_SORT // seems slower 2790 unsigned count0=0,count1=0; 2791 for (int i=0;i<relations.size();++i){ 2792 int c=relations[i].count=count_ones(relations[i].tab,C32); 2793 if (c==0) 2794 ++count0; 2795 if (c==1) 2796 ++count1; 2797 if (debug_infolevel>2){ 2798 cout << i << ", p="; 2799 if (i==0) 2800 cout << "-1"; 2801 else { 2802 if (i<=bs) 2803 cout << basis[i-1].p << " " << relations[i].count << '\n'; 2804 else 2805 cout << '\n'; 2806 } 2807 } 2808 } 2809 if (debug_infolevel) 2810 *logptr(contextptr) << CLOCK() << " begin rref size " << relations.size() << "x" << l32 << " K " << 0.004*relations.size()*C32 << ", " << count0 << " null lines, " << count1 << " 1-line" << '\n'; 2811 #if 0 // debug only 2812 for (int i=0;i<relations.size();++i){ 2813 cout << i << ", p="; 2814 if (i==0) 2815 cout << "-1"; 2816 else { 2817 if (i<=bs) 2818 cout << basis[i-1].p << " " << relations[i].count << '\n'; 2819 else 2820 cout << '\n'; 2821 } 2822 } 2823 #endif 2824 sort(relations.begin(),relations.end()); // put 0 lines at end, otherwise asc. sort 2825 #else // RREF_SORT 2826 if (debug_infolevel) 2827 *logptr(contextptr) << CLOCK() << " begin rref size " << relations.size() << "x" << l32 << " K " << 0.004*relations.size()*C32 << '\n'; 2828 reverse(relations.begin(),relations.end()); 2829 #endif // RREF_SORT 2830 // rref(relations,relations.size(),C32,0); 2831 rref(relations,int(relations.size()),C32,1); 2832 rref(relations,int(relations.size()),C32,2); 2833 if (debug_infolevel) 2834 *logptr(contextptr) << CLOCK() << " end rref" << '\n'; 2835 // printbool(*logptr(contextptr),relations); 2836 // move pivots on the diagonal by inserting 0 lines 2837 vector< unsigned * > relations2(l32); 2838 i=0; 2839 int j=0,rs=int(relations.size()); 2840 for (;i<rs && j<l32;++j){ 2841 if (relations[i].tab[j/32] & (1 << j%32)){ 2842 swap(relations2[j],relations[i].tab); 2843 ++i; 2844 } 2845 } 2846 // printbool(*logptr(contextptr),relations2); 2847 // for each element of the kernel compute x and y / x^2=y^2[N] 2848 // then gcd(x-y,n_orig) 2849 for (i=0;i<l32;++i){ 2850 if (relations2[i] && (relations2[i][i/32] & (1<<i%32))) 2851 continue; 2852 // using column i of relations2 which is in the kernel, build x and y 2853 // for x, we can compute the product of the axbmodn mod N 2854 // for y, compute the product of sqrta mod N and multiply later by ax^2+bx+c factored part 2855 // since (a*x+b)^2=a*(a*x^2+2*b*x+c) mod N 2856 gen x=1,y=1,cur; 2857 mpz_set_ui(zx,1); 2858 mpz_set_ui(zy,1); 2859 vector<short_t> p(bs), add_p(additional_primes.size()); 2860 for (int j=0;j<l32;++j){ 2861 if (j<int(axbmodn.size()) && (i==j || (relations2[j] && (relations2[j][i/32] & (1<<(i%32)))))){ 2862 if (axbmodn[j].aindex>=sqrtavals.size() || axbmodn[j].bindex>=bvals.size()) 2863 return false; // check added because ifactor(nextprime(alog10(17))*nextprime(alog10(19))); fails on Prime (and unable to do parallel debug in giac) 2864 update_xy(axbmodn[j],zx,zy,p,add_p,N,basis,additional_primes,sqrtavals,bvals,puissancestab,zq,zr,alloc1,alloc2,alloc3,alloc4,alloc5); 2865 #ifdef ADDITIONAL_PRIMES_HASHMAP 2866 unsigned u=largep(axbmodn[j],puissancestab); 2867 if (u) 2868 update_xy(additional_primes_map[u],zx,zy,p,add_p,N,basis,additional_primes,sqrtavals,bvals,puissancestab,zq,zr,alloc1,alloc2,alloc3,alloc4,alloc5); 2869 #endif 2870 } // end if (j<axbmodn.size() ...) 2871 } // end for unsigned j=0; j<l32 2872 for (int i=0;i<bs;++i){ 2873 if (p[i] % 2) 2874 *logptr(contextptr) << gettext("error, odd exponent for prime ") << basis[i].p << '\n'; 2875 if (p[i]){ 2876 #if 1 2877 mpz_set_ui(alloc1,basis[i].p); 2878 for (int j=0;j<p[i]/2;++j) 2879 mpz_mul(zy,zy,alloc1); 2880 #ifdef USE_GMP_REPLACEMENTS 2881 mp_grow(&alloc1,zy.used+2); 2882 mpz_set_ui(alloc1,0); 2883 alloc1.used = zy.used +2 ; 2884 mpz_set(alloc2,zy); 2885 mpz_set(alloc3,*N._ZINTptr); 2886 // mpz_set_si(alloc4,0); 2887 // mpz_set_si(alloc5,0); 2888 alloc_mp_div(&zy,N._ZINTptr,&zq,&zr,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 2889 #else 2890 mpz_tdiv_r(zr,zy,*N._ZINTptr); 2891 #endif 2892 mpz_set(zy,zr); 2893 #else 2894 y=y*pow(gen(basis[i].p),int(p[i]/2)); 2895 y=smod(y,N); 2896 #endif 2897 } 2898 } 2899 for (unsigned i=0;i<additional_primes.size();++i){ 2900 if (add_p[i] % 2) 2901 *logptr(contextptr) << gettext("error") << i << '\n'; 2902 if (add_p[i]){ 2903 #if 1 2904 mpz_set_ui(alloc1,additional_primes[i]); 2905 for (int j=0;j<add_p[i]/2;++j) 2906 mpz_mul(zy,zy,alloc1); 2907 #ifdef USE_GMP_REPLACEMENTS 2908 mp_grow(&alloc1,zy.used+2); 2909 mpz_set_ui(alloc1,0); 2910 alloc1.used = zy.used +2 ; 2911 mpz_set(alloc2,zy); 2912 mpz_set(alloc3,*N._ZINTptr); 2913 // mpz_set_si(alloc4,0); 2914 // mpz_set_si(alloc5,0); 2915 alloc_mp_div(&zy,N._ZINTptr,&zq,&zr,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 2916 #else 2917 mpz_tdiv_r(zr,zy,*N._ZINTptr); 2918 #endif 2919 mpz_set(zy,zr); 2920 #else 2921 y=y*pow(gen(additional_primes[i]),int(add_p[i]/2)); 2922 y=smod(y,N); 2923 #endif 2924 } 2925 } 2926 #if 1 2927 y=zy; 2928 x=zx; 2929 #endif 2930 cur=gcd(x-y,n_orig); 2931 if (debug_infolevel>6) 2932 *logptr(contextptr) << CLOCK() << gettext("checking gcd") << cur << " " << N << '\n'; 2933 if ( (cur.type==_INT_ && cur.val>7) || 2934 (cur.type==_ZINT && is_strictly_greater(n_orig,cur,contextptr))){ 2935 pn=cur; 2936 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2937 mpz_clear(alloc1); mpz_clear(alloc2); mpz_clear(alloc3); mpz_clear(alloc4); mpz_clear(alloc5); 2938 delete [] puissancestab; 2939 delete [] tab; 2940 return true; 2941 } 2942 } 2943 mpz_clear(zx); mpz_clear(zy); mpz_clear(zq); mpz_clear(zr); 2944 mpz_clear(alloc1); mpz_clear(alloc2); mpz_clear(alloc3); mpz_clear(alloc4); mpz_clear(alloc5); 2945 delete [] puissancestab; 2946 delete [] tab; 2947 return false; 2948 } 2949 2950 // elliptic curve method, 2951 // http://math.univ-lyon1.fr/~roblot/resources/factorisation.pdf 2952 // This is a very naive implementation 2953 // It does not use Montgomery representation and only phase 1 2954 // For professional implementations, cf. 2955 // https://members.loria.fr/PZimmermann/papers/ecm-submitted.pdf 2956 // https://pdfs.semanticscholar.org/e8eb/13b75292b15dd63c3e7e4b1c8dc334d278ba.pdf 2957 // ecm will be used if available 2958 #define ECM_MAXITER 1000 L(double alpha,double beta,double N)2959 static gen L(double alpha,double beta,double N){ 2960 double lnN=std::log(N); 2961 return std::exp(beta*std::pow(lnN,alpha)*std::pow(std::log(lnN),1-alpha)); 2962 } 2963 2964 2965 #ifndef USE_GMP_REPLACEMENTS 2966 // addition in elliptic curve, returns 1 on success or 0 and m=a divisor of n ecm_add(const mpz_t & x1,const mpz_t & y1,const mpz_t & x2,const mpz_t & y2,const mpz_t & a,const mpz_t & n,mpz_t & m,mpz_t & x,mpz_t & y)2967 int ecm_add(const mpz_t &x1,const mpz_t &y1,const mpz_t & x2,const mpz_t &y2,const mpz_t & a,const mpz_t & n,mpz_t & m,mpz_t & x,mpz_t &y){ 2968 if (mpz_cmp(x1,x2)){ 2969 mpz_sub(x,x2,x1); // x=x2-x1 2970 int res=mpz_invert(m,x,n); // m=inv(x2-x1) mod n 2971 if (res==0){ // not invertible 2972 mpz_gcd(m,x,n); 2973 return 0; // m has non trivial gcd with n 2974 } 2975 mpz_sub(y,y2,y1); 2976 mpz_mul(m,m,y); // m=(y2-y1)*invmod(x2-x1,n); 2977 } 2978 else { 2979 mpz_mul_ui(y,y1,2); 2980 int res=mpz_invert(m,y,n); // m=inv(2*y) mod n 2981 if (res==0){ // not invertible 2982 mpz_gcd(m,y,n); 2983 return 0; // m has non trivial gcd with n 2984 } 2985 mpz_mul(x,x1,x1); 2986 mpz_mul_ui(x,x,3); 2987 mpz_add(x,x,a); 2988 mpz_mul(m,m,x); // m=(3*x1*x1+a)*invmod(2*y1,n); 2989 } 2990 mpz_fdiv_r(m,m,n); // m=mod(m,n); 2991 mpz_mul(x,m,m); 2992 mpz_sub(x,x,x1); 2993 mpz_sub(x,x,x2); 2994 mpz_fdiv_r(x,x,n); // x=mod(m*m-x1-x2,n); 2995 mpz_sub(y,x1,x); 2996 mpz_mul(y,m,y); 2997 mpz_sub(y,y,y1); 2998 mpz_fdiv_r(y,y,n); // y=mod(m*(x1-x)-y1,n); 2999 // smod x and y 3000 mpz_add(m,x,x); 3001 if (mpz_cmp(m,n)>0) 3002 mpz_sub(x,x,n); 3003 mpz_add(m,y,y); 3004 if (mpz_cmp(m,n)>0) 3005 mpz_sub(y,y,n); 3006 return 1; 3007 } 3008 #endif 3009 ecm_add(const gen & x1,const gen & y1,const gen & x2,const gen & y2,const gen & a,const gen & n,gen & m,gen & x,gen & y)3010 gen ecm_add(const gen &x1,const gen &y1,const gen & x2,const gen &y2,const gen & a,const gen & n,gen & m,gen & x,gen &y){ 3011 if (is_inf(x1)){ 3012 x=x2; y=y2; return 1; 3013 } 3014 if (is_inf(x2)){ 3015 x=x1; y=y1; return 1; 3016 } 3017 if (y1+y2==0){ 3018 y=x=unsigned_inf; return 1; 3019 } 3020 #ifndef USE_GMP_REPLACEMENTS 3021 if (x1.type==_ZINT && y1.type==_ZINT && x2.type==_ZINT && y2.type==_ZINT && a.type==_ZINT && n.type==_ZINT ){ 3022 m=gen(1LL<<33); 3023 x=gen(1LL<<33); 3024 y=gen(1LL<<33); 3025 if (!ecm_add(*x1._ZINTptr,*y1._ZINTptr,*x2._ZINTptr,*y2._ZINTptr,*a._ZINTptr,*n._ZINTptr,*m._ZINTptr,*x._ZINTptr,*y._ZINTptr)){ 3026 return gen(*m._ZINTptr); 3027 } 3028 return 1; 3029 } 3030 #endif 3031 if (x1!=x2){ 3032 m=gcd(x2-x1,n); 3033 if (m!=1) 3034 return m; 3035 m=(y2-y1)*invmod(x2-x1,n); 3036 } 3037 else { 3038 m=gcd(y1,n); 3039 if (m!=1) 3040 return m; 3041 m=(3*x1*x1+a)*invmod(2*y1,n); 3042 } 3043 m=smod(m,n); 3044 x=smod(m*m-x1-x2,n); 3045 y=smod(m*(x1-x)-y1,n); 3046 return 1; 3047 } 3048 // multiplication in elliptic curve, ecm_mult(const gen & x1,const gen & y1,ulonglong m,const gen & a,const gen & n,gen & x,gen & y)3049 gen ecm_mult(const gen &x1,const gen &y1,ulonglong m,const gen & a,const gen & n,gen & x,gen &y){ 3050 gen x2(x1),y2(y1),xtmp,ytmp,g,M; 3051 y=x=plus_inf; 3052 while (m){ 3053 if (m%2){ 3054 g=ecm_add(x,y,x2,y2,a,n,M,xtmp,ytmp); 3055 if (g!=1) 3056 return g; 3057 swapgen(x,xtmp); swapgen(y,ytmp);// x=xtmp; y=ytmp; 3058 } 3059 m/=2; 3060 g=ecm_add(x2,y2,x2,y2,a,n,M,xtmp,ytmp); // improve: ecmdup 3061 if (g!=1) 3062 return g; 3063 swapgen(x2,xtmp);swapgen(y2,ytmp);// x2=xtmp; y2=ytmp; 3064 } 3065 return 1; 3066 } _ecm_factor(const gen & n_,GIAC_CONTEXT)3067 gen _ecm_factor(const gen &n_,GIAC_CONTEXT){ 3068 gen B,n(n_); 3069 int maxiter(ECM_MAXITER); 3070 if (n.type==_VECT && n._VECTptr->size()>=2){ 3071 const vecteur & v=*n._VECTptr; 3072 B=v[1]; 3073 if (v.size()>=3 && v[2].type==_INT_) 3074 maxiter=giacmax(1,v[2].val); 3075 n=v.front(); 3076 } 3077 if (!is_integer(n) || is_positive(-n,contextptr)) 3078 return gensizeerr(contextptr); 3079 if (_isprime(n,contextptr)!=0) 3080 return n; 3081 double logp=.5*std::log(evalf_double(n,1,contextptr)._DOUBLE_val); 3082 #ifdef HAVE_LIBECM 3083 double epsilon=.02; // to be adjusted 3084 #else 3085 double epsilon=.45; // to be adjusted 3086 #endif 3087 if (logp>80) // research factors of size not exceeding 35 digits 3088 logp=80; 3089 if (B==0) 3090 B=L(.5,0.707+epsilon,std::exp(logp)); 3091 // B=1000; 3092 B=_ceil(B,contextptr); 3093 #ifdef HAVE_LIBECM 3094 *logptr(contextptr) << "ECM-GMP factor n="<< n << " , B=" << B << ", #curves <=" << maxiter << '\n'; 3095 n.uncoerce(); 3096 double B1=evalf_double(B,1,contextptr)._DOUBLE_val; 3097 /* From ECM README, table of optimal values of B1 3098 digits D optimal B1 default B2 expected curves 3099 N(B1,B2,D) 3100 -power 1 default poly 3101 20 11e3 1.9e6 74 74 [x^1] 3102 25 5e4 1.3e7 221 214 [x^2] 3103 30 25e4 1.3e8 453 430 [D(3)] 3104 35 1e6 1.0e9 984 904 [D(6)] 3105 40 3e6 5.7e9 2541 2350 [D(6)] 3106 45 11e6 3.5e10 4949 4480 [D(12)] 3107 50 43e6 2.4e11 8266 7553 [D(12)] 3108 55 11e7 7.8e11 20158 17769 [D(30)] 3109 60 26e7 3.2e12 47173 42017 [D(30)] 3110 65 85e7 1.6e13 77666 69408 [D(30)] 3111 3112 */ 3113 int res; 3114 gen F(1LL<<33); 3115 for (int i=0;i<maxiter;++i){ 3116 res=ecm_factor(*F._ZINTptr, *n._ZINTptr, B1, 0); 3117 if (res!=0) break; 3118 } 3119 if (res==0) return undef; 3120 return F; 3121 #else 3122 *logptr(contextptr) << "ECM naive factor n="<< n << " , B=" << B << ", #curves <=" << maxiter << '\n'; 3123 for (int i=0;i<maxiter;++i){ 3124 gen a,x,y,b,d,g; 3125 a= rand_interval(makevecteur(0,n-1),true,contextptr);// a=1078104638; 3126 x= smod(rand_interval(makevecteur(0,n-1),true,contextptr),n);// 317359960; 3127 y= smod(rand_interval(makevecteur(0,n-1),true,contextptr),n);// 983830906; 3128 b=smod(y*y-x*x*x-a*x,n); 3129 if (debug_infolevel) 3130 COUT << CLOCK()*1e-6 << " Factor "<< n << " ECM curve " << i << ", B="<<B << ", a=" << a << ", x=" << x << ", y=" << y << '\n'; 3131 d=4*a*a*a-27*b*b; 3132 g=gcd(d,n); 3133 if (g==n) 3134 continue; 3135 if (g!=1) 3136 return g; 3137 gen p(2),pe,tmp,xm,ym; 3138 for (;is_greater(B,pe=p*p,contextptr);p=nextprime(p+1)){ 3139 int e=2; 3140 while (is_greater(B,tmp=pe*p,contextptr)){ 3141 ++e; 3142 pe=tmp; 3143 } 3144 tmp=evalf_double(pe,1,contextptr); 3145 if (pe==tmp){ 3146 g=ecm_mult(x,y,ulonglong(tmp._DOUBLE_val),a,n,xm,ym); 3147 if (g!=1){ 3148 if (debug_infolevel) 3149 COUT << CLOCK()*1e-6 << " ECM success p=" << p << ", p^e=" << pe << '\n'; 3150 return g; 3151 } 3152 swapgen(x,xm); swapgen(y,ym); // x=xm;y=ym; 3153 } 3154 } 3155 B=_floor(1.001*B,contextptr)+1; // to be adjusted 3156 } // end maxiter loop 3157 return undef; 3158 #endif 3159 } 3160 static const char _ecm_factor_s []="ecm_factor"; 3161 static define_unary_function_eval (__ecm_factor,&_ecm_factor,_ecm_factor_s); 3162 define_unary_function_ptr5( at_ecm_factor ,alias_at_ecm_factor,&__ecm_factor,0,true); 3163 3164 // Pollard-rho algorithm 3165 const int POLLARD_GCD=64; 3166 #ifdef GIAC_MPQS 3167 #if defined(RTOS_THREADX) // !defined(BESTA_OS) 3168 const int POLLARD_MAXITER=3000; 3169 #else 3170 const int POLLARD_MAXITER=15000; 3171 #endif 3172 #else 3173 const int POLLARD_MAXITER=15000; 3174 #endif 3175 pollard(gen n,gen k,GIAC_CONTEXT)3176 static gen pollard(gen n, gen k,GIAC_CONTEXT){ 3177 k.uncoerce(); 3178 n.uncoerce(); 3179 int maxiter=POLLARD_MAXITER; 3180 double nd=evalf_double(n,1,contextptr)._DOUBLE_val; 3181 #if defined(RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 3182 int nd1=int(2000*(std::log10(nd)-34)); 3183 #else 3184 int nd1=int(1500*std::pow(16.,(std::log10(nd)-40)/10)); 3185 #endif 3186 if (nd1>maxiter) 3187 maxiter=nd1; 3188 int m,m1,a,a1,j; 3189 m1=m=2; 3190 a1=a=1; 3191 int c=0; 3192 mpz_t g,x,x1,x2,x2k,y,y1,p,q,tmpq,alloc1,alloc2,alloc3,alloc4,alloc5; 3193 mpz_init_set_si(g,1); // ? mp_init_size to specify size 3194 mpz_init_set_si(x,2); 3195 mpz_init_set_si(x1,2); 3196 mpz_init_set_si(y,2); 3197 mpz_init(y1); 3198 mpz_init(x2); 3199 mpz_init(x2k); 3200 mpz_init_set_si(p,1); 3201 mpz_init(q); 3202 mpz_init(tmpq); 3203 mpz_init(alloc1); 3204 mpz_init(alloc2); 3205 mpz_init(alloc3); 3206 mpz_init(alloc4); 3207 mpz_init(alloc5); 3208 while (!ctrl_c && !interrupted && mpz_cmp_si(g,1)==0) { 3209 #ifdef TIMEOUT 3210 control_c(); 3211 #endif 3212 a=2*a+1;//a=2^(e+1)-1=2*l(m)-1 3213 while (!ctrl_c && !interrupted && mpz_cmp_si(g,1)==0 && a>m) { // ok 3214 #ifdef TIMEOUT 3215 control_c(); 3216 #endif 3217 // x=f(x,k,n,q); 3218 #ifdef USE_GMP_REPLACEMENTS 3219 mp_sqr(&x,&x2); 3220 mpz_add(x2k,x2,*k._ZINTptr); 3221 if (mpz_cmp(x2k,*n._ZINTptr)>0){ 3222 mp_grow(&alloc1,x2k.used+2); 3223 mpz_set_ui(alloc1,0); 3224 alloc1.used = x2k.used +2 ; 3225 mpz_set(alloc2,x2k); 3226 mpz_set(alloc3,*n._ZINTptr); 3227 // mpz_set_si(alloc4,0); 3228 // mpz_set_si(alloc5,0); 3229 alloc_mp_div(&x2k,n._ZINTptr,&tmpq,&x,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 3230 } 3231 else 3232 mpz_set(x,x2k); 3233 #else 3234 mpz_mul(x2,x,x); 3235 mpz_add(x2k,x2,*k._ZINTptr); 3236 mpz_tdiv_r(x,x2k,*n._ZINTptr); 3237 #endif 3238 m += 1; 3239 if (debug_infolevel && ((m % 3240 #if defined(RTOS_THREADX) || defined(BESTA_OS) || defined NSPIRE 3241 (1<<10) 3242 #else 3243 (1<<18) 3244 #endif 3245 )==0)) 3246 *logptr(contextptr) << CLOCK() << gettext(" Pollard-rho try ") << m << '\n'; 3247 if (m > maxiter ){ 3248 if (debug_infolevel) 3249 *logptr(contextptr) << CLOCK() << gettext(" Pollard-rho failure, ntries ") << m << '\n'; 3250 mpz_clear(alloc5); 3251 mpz_clear(alloc4); 3252 mpz_clear(alloc3); 3253 mpz_clear(alloc2); 3254 mpz_clear(alloc1); 3255 mpz_clear(tmpq); 3256 mpz_clear(x); 3257 mpz_clear(x1); 3258 mpz_clear(x2); 3259 mpz_clear(x2k); 3260 mpz_clear(y); 3261 mpz_clear(y1); 3262 mpz_clear(p); 3263 mpz_clear(q); 3264 return -1; 3265 } 3266 // p=irem(p*(x1-x),n,q); 3267 mpz_sub(q,x1,x); 3268 mpz_mul(x2,p,q); 3269 #if 0 // def USE_GMP_REPLACEMENTS 3270 if (mpz_cmp(x2,*n._ZINTptr)>0){ 3271 mp_grow(&alloc1,x2.used+2); 3272 mpz_set_ui(alloc1,0); 3273 alloc1.used = x2.used +2 ; 3274 mpz_set(alloc2,x2); 3275 mpz_set(alloc3,*n._ZINTptr); 3276 // mpz_set_si(alloc4,0); 3277 // mpz_set_si(alloc5,0); 3278 alloc_mp_div(&x2,n._ZINTptr,&tmpq,&p,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 3279 } 3280 else 3281 mpz_set(p,x2); 3282 #else 3283 mpz_tdiv_r(p,x2,*n._ZINTptr); 3284 #endif 3285 c += 1; 3286 if (c==POLLARD_GCD) { 3287 // g=gcd(abs(p,context0),n); 3288 mpz_abs(q,p); 3289 my_mpz_gcd(g,q,*n._ZINTptr); 3290 if (mpz_cmp_si(g,1)==0) { 3291 mpz_set(y,x); // y=x; 3292 mpz_set(y1,x1); // y1=x1; 3293 mpz_set_si(p,1); // p=1; 3294 a1=a; 3295 m1=m; 3296 c=0; 3297 } 3298 } 3299 }//m=a=2^e-1=l(m) 3300 if (mpz_cmp_si(g,1)==0) { 3301 mpz_set(x1,x); // x1=x;//x1=x_m=x_l(m)-1 3302 j=3*(a+1)/2; // j=3*iquo(a+1,2); 3303 for (long i=m+1;i<=j;i++){ 3304 // x=f(x,k,n,q); 3305 mpz_mul(x2,x,x); 3306 mpz_add(x2k,x2,*k._ZINTptr); 3307 #if 0 // def USE_GMP_REPLACEMENTS 3308 if (mpz_cmp(x2k,*n._ZINTptr)>0){ 3309 mp_grow(&alloc1,x2k.used+2); 3310 mpz_set_ui(alloc1,0); 3311 alloc1.used = x2k.used +2 ; 3312 mpz_set(alloc2,x2k); 3313 mpz_set(alloc3,*n._ZINTptr); 3314 // mpz_set_si(alloc4,0); 3315 // mpz_set_si(alloc5,0); 3316 alloc_mp_div(&x2k,n._ZINTptr,&tmpq,&x,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 3317 } 3318 else 3319 mpz_set(x,x2); 3320 #else 3321 mpz_tdiv_r(x,x2k,*n._ZINTptr); 3322 #endif 3323 } 3324 m=j; 3325 } 3326 } 3327 //g<>1 ds le paquet de POLLARD_GCD 3328 if (debug_infolevel>5) 3329 CERR << CLOCK() << " Pollard-rho nloops " << m << '\n'; 3330 mpz_set(x,y); // x=y; 3331 mpz_set(x1,y1); // x1=y1; 3332 mpz_set_si(g,1); // g=1; 3333 a=(a1-1)/2; // a=iquo(a1-1,2); 3334 m=m1; 3335 while (!ctrl_c && !interrupted && mpz_cmp_si(g,1)==0) { 3336 #ifdef TIMEOUT 3337 control_c(); 3338 #endif 3339 a=2*a+1; 3340 while (!ctrl_c && !interrupted && mpz_cmp_si(g,1)==0 && a>m) { // ok 3341 #ifdef TIMEOUT 3342 control_c(); 3343 #endif 3344 // x=f(x,k,n,q); 3345 mpz_mul(x2,x,x); 3346 mpz_add(x2k,x2,*k._ZINTptr); 3347 #if 0 // def USE_GMP_REPLACEMENTS 3348 if (mpz_cmp(x2k,*n._ZINTptr)>0){ 3349 mp_grow(&alloc1,x2k.used+2); 3350 mpz_set_ui(alloc1,0); 3351 alloc1.used = x2k.used +2 ; 3352 mpz_set(alloc2,x2k); 3353 mpz_set(alloc3,*n._ZINTptr); 3354 alloc_mp_div(&x2k,n._ZINTptr,&tmpq,&x,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 3355 } 3356 else 3357 mpz_set(x,x2k); 3358 #else 3359 mpz_tdiv_r(x,x2k,*n._ZINTptr); 3360 #endif 3361 m += 1; 3362 if (m > maxiter ){ 3363 mpz_clear(alloc5); 3364 mpz_clear(alloc4); 3365 mpz_clear(alloc3); 3366 mpz_clear(alloc2); 3367 mpz_clear(alloc1); 3368 mpz_clear(tmpq); 3369 mpz_clear(x); 3370 mpz_clear(x1); 3371 mpz_clear(x2); 3372 mpz_clear(x2k); 3373 mpz_clear(y); 3374 mpz_clear(y1); 3375 mpz_clear(p); 3376 mpz_clear(q); 3377 return -1; 3378 } 3379 // p=irem(x1-x,n,q); 3380 mpz_sub(q,x1,x); 3381 #if 0 // def USE_GMP_REPLACEMENTS 3382 if (mpz_cmp(q,*n._ZINTptr)>0){ 3383 mp_grow(&alloc1,q.used+2); 3384 mpz_set_ui(alloc1,0); 3385 alloc1.used = q.used +2 ; 3386 mpz_set(alloc2,q); 3387 mpz_set(alloc3,*n._ZINTptr); 3388 // mpz_set_si(alloc4,0); 3389 // mpz_set_si(alloc5,0); 3390 alloc_mp_div(&q,n._ZINTptr,&tmpq,&p,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 3391 } 3392 else 3393 mpz_set(p,q); 3394 #else 3395 mpz_tdiv_r(p,q,*n._ZINTptr); 3396 #endif 3397 // g=gcd(abs(p,context0),n); // ok 3398 mpz_abs(q,p); 3399 my_mpz_gcd(g,q,*n._ZINTptr); 3400 } 3401 if (mpz_cmp_si(g,1)==0) { 3402 mpz_set(x1,x); // x1=x; 3403 j=3*(a+1)/2; // j=3*iquo(a+1,2); 3404 for (long i=m+1;j>=i;i++){ 3405 // x=f(x,k,n,q); 3406 mpz_mul(x2,x,x); 3407 mpz_add(x2k,x2,*k._ZINTptr); 3408 mpz_tdiv_qr(tmpq,x,x2k,*n._ZINTptr); 3409 } 3410 m=j; 3411 } 3412 } 3413 mpz_clear(alloc5); 3414 mpz_clear(alloc4); 3415 mpz_clear(alloc3); 3416 mpz_clear(alloc2); 3417 mpz_clear(alloc1); 3418 mpz_clear(tmpq); 3419 mpz_clear(x); 3420 mpz_clear(x1); 3421 mpz_clear(x2); 3422 mpz_clear(x2k); 3423 mpz_clear(y); 3424 mpz_clear(y1); 3425 mpz_clear(p); 3426 mpz_clear(q); 3427 #ifdef TIMEOUT 3428 control_c(); 3429 #endif 3430 if (ctrl_c || interrupted){ 3431 mpz_clear(g); 3432 return 0; 3433 } 3434 if (mpz_cmp(g,*n._ZINTptr)==0) { 3435 if (k==1) { 3436 mpz_clear(g); 3437 return(pollard(n,-1,contextptr)); 3438 } 3439 else { 3440 if (k*k==1){ 3441 mpz_clear(g); 3442 return(pollard(n,3,contextptr)); 3443 } 3444 else { 3445 if (is_greater(k,50,contextptr)){ 3446 #if 1 3447 return -1; 3448 #else 3449 ref_mpz_t * ptr=new ref_mpz_t; 3450 mpz_init_set(ptr->z,g); 3451 mpz_clear(g); 3452 return ptr; 3453 #endif 3454 } 3455 else { 3456 mpz_clear(g); 3457 return(pollard(n,k+2,contextptr)); 3458 } 3459 } 3460 } 3461 } 3462 ref_mpz_t * ptr=new ref_mpz_t; 3463 mpz_init_set(ptr->z,g); 3464 mpz_clear(g); 3465 return ptr; 3466 } 3467 3468 // const short int giac_primes[]={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,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997}; 3469 eratosthene(double n,vector<bool> * & v)3470 bool eratosthene(double n,vector<bool> * & v){ 3471 static vector<bool> *ptr=0; 3472 if (!ptr) 3473 ptr=new vector<bool>; 3474 vector<bool> &erato=*ptr; 3475 v=ptr; 3476 if (n+1>erato.size()){ 3477 unsigned N=int(n); 3478 ++N; 3479 #if defined BESTA_OS 3480 if (N>2e6) 3481 return false; 3482 #else 3483 if (N>2e9) 3484 return false; 3485 #endif 3486 N = (N*11)/10; 3487 erato=vector<bool>(N+1,true); 3488 // insure that we won't recompute all again from start for ithprime(i+1) 3489 for (unsigned p=2;;++p){ 3490 while (!erato[p]) // find next prime 3491 ++p; 3492 if (p*p>N) // finished 3493 return true; 3494 for (unsigned i=2*p;i<=N;i+=p) 3495 erato[i]=false; // remove p multiples 3496 } 3497 } 3498 return true; 3499 } 3500 eratosthene2(double n,vector<bool> * & v)3501 bool eratosthene2(double n,vector<bool> * & v){ 3502 static vector<bool> *ptr=0; 3503 if (!ptr) 3504 ptr=new vector<bool>; 3505 vector<bool> &erato=*ptr; 3506 v=ptr; 3507 if (n/2>=erato.size()){ 3508 unsigned N=int(n); 3509 ++N; 3510 #if defined BESTA_OS 3511 if (N>4e6) 3512 return false; 3513 #else 3514 if (N>2e9) 3515 return false; 3516 #endif 3517 // 11/20 insures that we won't recompute all again from start for ithprime(i+1) 3518 N = (N*11)/20; // keep only odd numbers in sieve 3519 erato=vector<bool>(N+1,true); //erato[i] stands for 2*i+1 <-> n corresponds to erato[n/2] 3520 for (unsigned p=3;;p+=2){ 3521 while (!erato[p/2]) // find next prime (first one is p==3) 3522 p+=2; 3523 if (p*p>2*N+1) // finished 3524 return true; 3525 // p is prime, set p*p, (p+2)*p, etc. to be non prime 3526 for (unsigned i=(p*p)/2;i<=N;i+=p) 3527 erato[i]=false; // remove p multiples 3528 } 3529 } 3530 return true; 3531 } 3532 3533 // ithprime(n) is approx invli(n)+invli(sqrt(n))/4 where invli is reciproc. 3534 // of Li(x)=Ei(ln(x)) 3535 // For fast code, cf. https://github.com/kimwalisch/primecount 3536 static const char _ithprime_s []="ithprime"; symb_ithprime(const gen & args)3537 static symbolic symb_ithprime(const gen & args){ 3538 return symbolic(at_ithprime,args); 3539 } ithprime(const gen & g_,GIAC_CONTEXT)3540 static gen ithprime(const gen & g_,GIAC_CONTEXT){ 3541 gen g(g_); 3542 if (!is_integral(g)) 3543 return gentypeerr(contextptr); 3544 if (g.type!=_INT_) 3545 return gensizeerr(contextptr); // symb_ithprime(g); 3546 int i=g.val; 3547 if (i<0) 3548 return gensizeerr(contextptr); 3549 if (i==0) 3550 return 1; 3551 if (i<=int(sizeof(giac_primes)/sizeof(short int))) 3552 return giac_primes[i-1]; 3553 vector<bool> * vptr=0; 3554 #if 1 3555 if (!eratosthene2(i*std::log(double(i))*1.1,vptr)) 3556 return gensizeerr(contextptr); 3557 unsigned count=2; 3558 unsigned s=unsigned(vptr->size()); 3559 for (unsigned k=2;k<s;++k){ 3560 if ((*vptr)[k]){ 3561 ++count; 3562 if (i==count) 3563 return int(2*k+1); 3564 } 3565 } 3566 return undef; 3567 #else 3568 if (!eratosthene(i*std::log(double(i))*1.1,vptr)) 3569 return gensizeerr(contextptr); 3570 unsigned count=2; 3571 unsigned s=vptr->size(); 3572 for (unsigned k=4;k<s;++k){ 3573 if ((*vptr)[k]){ 3574 ++count; 3575 if (i==count) 3576 return int(k); 3577 } 3578 } 3579 return undef; 3580 #endif 3581 } _ithprime(const gen & args,GIAC_CONTEXT)3582 gen _ithprime(const gen & args,GIAC_CONTEXT){ 3583 if ( args.type==_STRNG && args.subtype==-1) return args; 3584 if (args.type==_VECT) 3585 return apply(args,_ithprime,contextptr); 3586 return ithprime(args,contextptr); 3587 } 3588 static define_unary_function_eval (__ithprime,&_ithprime,_ithprime_s); 3589 define_unary_function_ptr5( at_ithprime ,alias_at_ithprime,&__ithprime,0,true); 3590 3591 static const char _nprimes_s []="nprimes"; nprimes(const gen & g_,GIAC_CONTEXT)3592 static gen nprimes(const gen & g_,GIAC_CONTEXT){ 3593 gen g(g_); 3594 if (!is_integral(g)) 3595 return gentypeerr(contextptr); 3596 if (g.type!=_INT_) 3597 return gensizeerr(contextptr); // symb_ithprime(g); 3598 int i=g.val; 3599 if (i<0) 3600 return gensizeerr(contextptr); 3601 if (i<2) 3602 return 0; 3603 vector<bool> * vptr=0; 3604 if (!eratosthene2(i+2,vptr)) 3605 return gensizeerr(contextptr); 3606 unsigned count=1; // 2 is prime, then count odd primes 3607 i=(i-1)/2; 3608 for (int k=1;k<=i;++k){ 3609 if ((*vptr)[k]) 3610 ++count; 3611 } 3612 return int(count); 3613 } _nprimes(const gen & args,GIAC_CONTEXT)3614 gen _nprimes(const gen & args,GIAC_CONTEXT){ 3615 if ( args.type==_STRNG && args.subtype==-1) return args; 3616 if (args.type==_VECT) 3617 return apply(args,_nprimes,contextptr); 3618 return nprimes(args,contextptr); 3619 } 3620 static define_unary_function_eval (__nprimes,&_nprimes,_nprimes_s); 3621 define_unary_function_ptr5( at_nprimes ,alias_at_nprimes,&__nprimes,0,true); 3622 is_divisible_by(const gen & n,unsigned long a)3623 bool is_divisible_by(const gen & n,unsigned long a){ 3624 if (n.type==_ZINT){ 3625 #ifdef USE_GMP_REPLACEMENTS 3626 mp_digit c; 3627 mp_mod_d(n._ZINTptr, a, &c); 3628 return c==0; 3629 #else 3630 return mpz_divisible_ui_p(*n._ZINTptr,a); 3631 #endif 3632 } 3633 return n.val%a==0; 3634 } 3635 3636 // find trivial factors of n, 3637 // if add_last is true the remainder is put in the vecteur, 3638 // otherwise n contains the remainder pfacprem(gen & n,bool add_last,GIAC_CONTEXT)3639 vecteur pfacprem(gen & n,bool add_last,GIAC_CONTEXT){ 3640 gen a; 3641 gen q; 3642 int p,i,prime; 3643 vecteur v(2); 3644 vecteur u; 3645 if (is_zero(n)) 3646 return u; 3647 if (n.type==_ZINT){ 3648 ref_mpz_t * cur = new ref_mpz_t; 3649 mpz_t div,q,r,alloc1,alloc2,alloc3,alloc4,alloc5; 3650 mpz_set(cur->z,*n._ZINTptr); 3651 mpz_init_set(q,*n._ZINTptr); 3652 mpz_init(r); 3653 mpz_init(div); 3654 mpz_init(alloc1); 3655 mpz_init(alloc2); 3656 mpz_init(alloc3); 3657 mpz_init(alloc4); 3658 mpz_init(alloc5); 3659 for (i=0;i<int(sizeof(giac_primes)/sizeof(short int));++i){ 3660 if (mpz_cmp_si(cur->z,1)==0) 3661 break; 3662 prime=giac_primes[i]; 3663 mpz_set_ui(div,prime); 3664 #ifdef USE_GMP_REPLACEMENTS 3665 for (p=0;;p++){ 3666 mp_grow(&alloc1,cur->z.used+2); 3667 mpz_set_ui(alloc1,0); 3668 alloc1.used = cur->z.used +2 ; 3669 mpz_set(alloc2,cur->z); 3670 mpz_set(alloc3,div); 3671 alloc_mp_div(&cur->z,&div,&q,&r,&alloc1,&alloc2,&alloc3,&alloc4,&alloc5); 3672 // mpz_tdiv_qr(q,r,cur->z,div); 3673 if (mpz_cmp_si(r,0)) 3674 break; 3675 mp_exch(&cur->z,&q); 3676 } 3677 // *logptr(contextptr) << "Factor " << prime << " " << p << '\n'; 3678 if (p){ 3679 u.push_back(prime); 3680 u.push_back(p); 3681 } 3682 #else 3683 if (mpz_divisible_ui_p(cur->z,prime)){ 3684 mpz_set_ui(div,prime); 3685 for (p=0;;p++){ 3686 mpz_tdiv_qr(q,r,cur->z,div); 3687 if (mpz_cmp_si(r,0)) 3688 break; 3689 mpz_swap(cur->z,q); 3690 } 3691 // *logptr(contextptr) << "Factor " << prime << " " << p << '\n'; 3692 u.push_back(prime); 3693 u.push_back(p); 3694 } 3695 #endif 3696 } // end for on smal primes 3697 mpz_clear(alloc5); 3698 mpz_clear(alloc4); 3699 mpz_clear(alloc3); 3700 mpz_clear(alloc2); 3701 mpz_clear(alloc1); 3702 mpz_clear(div); mpz_clear(r); mpz_clear(q); 3703 n=cur; 3704 } 3705 else { 3706 for (i=0;i<int(sizeof(giac_primes)/sizeof(short int));++i){ 3707 if (n==1) 3708 break; 3709 a.val=giac_primes[i]; 3710 p=0; 3711 while (is_divisible_by(n,a.val)){ // while (irem(n,a,q)==0){ 3712 n=iquo(n,a); 3713 p=p+1; 3714 } 3715 if (p!=0){ 3716 // *logptr(contextptr) << "Factor " << a << " " << p << '\n'; 3717 u.push_back(a); 3718 u.push_back(p); 3719 } 3720 } 3721 } 3722 if (add_last && i==1229 && !is_one(n)){ 3723 // hack: check if n is a perfect square 3724 double nf=evalf_double(n,1,contextptr)._DOUBLE_val; 3725 nf=std::sqrt(nf); 3726 gen n2=_round(nf,contextptr); 3727 if (n2*n2==n){ 3728 u.push_back(n2); 3729 u.push_back(2); 3730 } 3731 else { 3732 u.push_back(n); 3733 u.push_back(1); 3734 } 3735 n=1; 3736 } 3737 //v[0]=n; 3738 //v[1]=u; 3739 3740 return(u); 3741 } 3742 3743 #ifdef USE_GMP_REPLACEMENTS inpollardsieve(const gen & a,gen k,bool & do_pollard,GIAC_CONTEXT)3744 static gen inpollardsieve(const gen &a,gen k,bool & do_pollard,GIAC_CONTEXT){ 3745 gen b=do_pollard?pollard(a,k,contextptr):-1; 3746 #ifdef TIMEOUT 3747 control_c(); 3748 #endif 3749 #ifdef GIAC_MPQS 3750 if (b==-1 && !ctrl_c && !interrupted){ 3751 do_pollard=false; 3752 if (msieve(a,b,contextptr)) return b; else return -1; } 3753 #endif 3754 return b; 3755 } pollardsieve(const gen & a,gen k,bool & do_pollard,GIAC_CONTEXT)3756 static gen pollardsieve(const gen &a,gen k,bool & do_pollard,GIAC_CONTEXT){ 3757 #if defined( GIAC_HAS_STO_38) || defined(EMCC) || defined NSPIRE 3758 int debug_infolevel_=debug_infolevel; 3759 #if defined RTOS_THREADX || defined NSPIRE 3760 debug_infolevel=2; 3761 if (do_pollard) 3762 *logptr(contextptr) << gettext(" Pollard-rho on ") << a << '\n'; 3763 #else 3764 debug_infolevel=0; 3765 #endif 3766 #endif 3767 gen res=inpollardsieve(a,k,do_pollard,contextptr); 3768 #if defined( GIAC_HAS_STO_38) || defined(EMCC) || defined NSPIRE 3769 debug_infolevel=debug_infolevel_; 3770 #ifdef GIAC_HAS_STO_38 3771 Calc->Terminal.MakeUnvisible(); 3772 #endif 3773 #endif 3774 return res; 3775 } 3776 #else // USE_GMP_REPLACEMENTS pollardsieve(const gen & a,gen k,bool & do_pollard,GIAC_CONTEXT)3777 static gen pollardsieve(const gen &a,gen k,bool & do_pollard,GIAC_CONTEXT){ 3778 gen b=do_pollard?pollard(a,k,contextptr):-1; 3779 #ifdef TIMEOUT 3780 control_c(); 3781 #endif 3782 #ifdef HAVE_LIBECM 3783 if (is_greater(a,1e60,context0) && b==-1 && !ctrl_c && !interrupted && _isprime(a,contextptr)==0){ 3784 int res; 3785 gen F(1LL<<33); 3786 for (int i=0;i<200;++i){ // searching factors of size about 20 digits 3787 res=ecm_factor(*F._ZINTptr, *a._ZINTptr, 11e3, 0); 3788 if (res!=0) break; 3789 } 3790 if (res!=0) 3791 b=F; 3792 } 3793 #endif 3794 #ifdef GIAC_MPQS 3795 if (b==-1 && !ctrl_c && !interrupted){ 3796 do_pollard=false; 3797 if (msieve(a,b,contextptr)) return b; else return -1; } 3798 #else 3799 if (b==-1) 3800 *logptr(contextptr) << "Integer too large for factorization algorithm\n"; 3801 #endif 3802 if (b==-1) 3803 b=a; 3804 return b; 3805 } 3806 #endif // USE_GMP_REPLACEMENTS 3807 ifactor2(const gen & n,vecteur & v,bool & do_pollard,GIAC_CONTEXT)3808 static gen ifactor2(const gen & n,vecteur & v,bool & do_pollard,GIAC_CONTEXT){ 3809 if (is_greater(giac_last_prime*giac_last_prime,n,contextptr) || is_probab_prime_p(n) ){ 3810 v.push_back(n); 3811 return 1; 3812 } 3813 // Check for power of integer: arg must be > 1e4, n*ln(arg)=d => n<d/ln(1e4) 3814 double d=evalf_double(n,1,contextptr)._DOUBLE_val; 3815 int maxpow=int(std::ceil(std::log(d)/std::log(1e4))); 3816 for (int i=2;i<=maxpow;++i){ 3817 if ( (i>2 && i%2==0) || 3818 (i>3 && i%3==0) || 3819 (i>5 && i%5==0) || 3820 (i>7 && i%7==0) ) 3821 continue; 3822 gen u; 3823 if (i==2) 3824 u=isqrt(n); 3825 else { 3826 double x=std::pow(d,1./i); 3827 u=longlong(x); 3828 } 3829 if (pow(u,i,contextptr)==n){ 3830 vecteur w; 3831 do_pollard=true; 3832 ifactor2(u,w,do_pollard,contextptr); 3833 for (int j=0;j<i;j++) 3834 v=mergevecteur(v,w); 3835 return v; 3836 } 3837 } 3838 gen a=pollardsieve(n,1,do_pollard,contextptr); 3839 if (a==-1) 3840 return a; 3841 #ifdef TIMEOUT 3842 control_c(); 3843 #endif 3844 if (ctrl_c || interrupted) 3845 return gensizeerr("Interrupted"); 3846 gen ba=n/a; 3847 if (a!=n) 3848 a=ifactor2(a,v,do_pollard,contextptr); 3849 else { 3850 a=1; 3851 v.push_back(n); 3852 } 3853 if (is_strictly_greater(ba,1,contextptr)) 3854 a=ifactor2(ba,v,do_pollard,contextptr); 3855 return a; 3856 } 3857 facprem(gen & n,GIAC_CONTEXT)3858 static vecteur facprem(gen & n,GIAC_CONTEXT){ 3859 vecteur v; 3860 if (n==1) { return v; } 3861 if ( (n.type==_INT_ && n.val<giac_last_prime*giac_last_prime) || is_probab_prime_p(n)) { 3862 v.push_back(n); 3863 n=1; 3864 return v; 3865 } 3866 if (debug_infolevel>5) 3867 CERR << "Pollard begin " << CLOCK() << '\n'; 3868 bool do_pollard=true; 3869 gen a=ifactor2(n,v,do_pollard,contextptr); 3870 if (a==-1) 3871 return makevecteur(gensizeerr(gettext("Quadratic sieve failure, perhaps number too large"))); 3872 if (is_zero(a)) 3873 return makevecteur(gensizeerr(gettext("Stopped by user interruption"))); 3874 n=1; 3875 return v; 3876 } 3877 mergeifactors(const vecteur & f,const vecteur & g,vecteur & h)3878 void mergeifactors(const vecteur & f,const vecteur &g,vecteur & h){ 3879 h=f; 3880 for (unsigned i=0;i<g.size();i+=2){ 3881 unsigned j=0; 3882 for (;j<f.size();j+=2){ 3883 if (f[j]==g[i]) 3884 break; 3885 } 3886 if (j<f.size()) 3887 h[j+1] += g[i+1]; 3888 else { 3889 h.push_back(g[i]); 3890 h.push_back(g[i+1]); 3891 } 3892 } 3893 } 3894 giac_ifactors(const gen & n0,GIAC_CONTEXT)3895 static vecteur giac_ifactors(const gen & n0,GIAC_CONTEXT){ 3896 if (!is_integer(n0) || is_zero(n0)) 3897 return vecteur(1,gensizeerr(gettext("ifactors"))); 3898 if (is_one(n0)) 3899 return vecteur(0); 3900 gen n(n0); 3901 vecteur f; 3902 vecteur g; 3903 vecteur u; 3904 // First find if |n-k^d|<=1 for d = 2, 3, 5 or 7 3905 double nd=evalf_double(n,1,contextptr)._DOUBLE_val; 3906 double nd2=std::floor(std::sqrt(nd)+.5); 3907 if (absdouble(1-nd2*nd2/nd)<1e-10){ 3908 gen n2=isqrt(n+1); 3909 if (n==n2*n2){ 3910 f=ifactors(n2,contextptr); 3911 iterateur it=f.begin(),itend=f.end(); 3912 for (;it!=itend;++it){ 3913 ++it; 3914 *it = 2 * *it; 3915 } 3916 return f; 3917 } 3918 if (n==n2*n2-1){ 3919 f=ifactors(n2-1,contextptr); 3920 g=ifactors(n2+1,contextptr); 3921 mergeifactors(f,g,u); 3922 return u; 3923 } 3924 } 3925 for (int k=3;;){ 3926 nd2=std::floor(std::pow(nd,1./k)+.5); 3927 if (absdouble(1-std::pow(nd2,k)/nd)<1e-10){ 3928 gen n2=_floor(nd2,contextptr),nf=n2*n2; 3929 for (int j=2;j<k;j++) 3930 nf=nf*n2; 3931 if (n==nf){ 3932 f=ifactors(n2,contextptr); 3933 iterateur it=f.begin(),itend=f.end(); 3934 for (;it!=itend;++it){ 3935 ++it; 3936 *it = k * *it; 3937 } 3938 return f; 3939 } 3940 if (n==nf-1){ // n2^k-1 3941 f=ifactors(n2-1,contextptr); 3942 g=ifactors(n/(n2-1),contextptr); 3943 mergeifactors(f,g,u); 3944 return u; 3945 } 3946 if (n==nf+1){ // n2^k+1 3947 f=ifactors(n2+1,contextptr); 3948 g=ifactors(n/(n2+1),contextptr); 3949 mergeifactors(f,g,u); 3950 return u; 3951 } 3952 } 3953 if (k==11) break; 3954 if (k==7) break; // k=11; 3955 if (k==5) k=7; 3956 if (k==3) k=5; 3957 } 3958 //f=pfacprem(n,false,contextptr); 3959 //cout<<n<<" "<<f<<'\n'; 3960 while (n!=1) { 3961 g=facprem(n,contextptr); 3962 if (is_undef(g)) 3963 return g; 3964 islesscomplexthanf_sort(g.begin(),g.end()); 3965 gen last=0; int p=0; 3966 for (unsigned i=0;i<g.size();++i){ 3967 if (g[i]==last) 3968 ++p; 3969 else { 3970 if (last!=0){ 3971 u.push_back(last); 3972 u.push_back(p); 3973 } 3974 last=g[i]; 3975 p=1; 3976 } 3977 } 3978 u.push_back(last); 3979 u.push_back(p); 3980 } 3981 g=mergevecteur(f,u); 3982 return g; 3983 } 3984 ifactors1(const gen & n0,GIAC_CONTEXT)3985 static vecteur ifactors1(const gen & n0,GIAC_CONTEXT){ 3986 if (is_greater(1e71,n0,contextptr)) 3987 return giac_ifactors(n0,contextptr); 3988 if (n0.type==_VECT && !n0._VECTptr->empty()) 3989 return giac_ifactors(n0._VECTptr->front(),contextptr); 3990 if (!is_integer(n0) || is_zero(n0)) 3991 return vecteur(1,gensizeerr(gettext("ifactors"))); 3992 if (is_one(n0)) 3993 return vecteur(0); 3994 if (_isprime(n0,contextptr)!=0) 3995 return makevecteur(n0,1); 3996 #ifdef HAVE_LIBECM 3997 int res; 3998 gen F(1LL<<33); 3999 double B1=1e6; 4000 for (int i=0;i<ECM_MAXITER;++i){ 4001 res=ecm_factor(*F._ZINTptr, *n0._ZINTptr, B1, 0); 4002 if (res!=0) break; 4003 } 4004 if (res!=0){ 4005 vecteur tmp=ifactors1(n0/F,contextptr); 4006 tmp.push_back(F); 4007 tmp.push_back(1); 4008 return tmp; 4009 } 4010 #endif 4011 #ifdef HAVE_LIBPARI 4012 #ifdef __APPLE__ 4013 return vecteur(1,gensizeerr(gettext("(Mac OS) Large number, you can try pari(); pari_factor(")+n0.print(contextptr)+")")); 4014 #endif 4015 gen g(pari_ifactor(n0),contextptr); 4016 if (g.type==_VECT){ 4017 matrice m(mtran(*g._VECTptr)); 4018 vecteur res; 4019 const_iterateur it=m.begin(),itend=m.end(); 4020 for (;it!=itend;++it){ 4021 if (it->type!=_VECT) return vecteur(1,gensizeerr(gettext("ifactor.cc/ifactors"))); 4022 res.push_back(it->_VECTptr->front()); 4023 res.push_back(it->_VECTptr->back()); 4024 } 4025 return res; 4026 } 4027 #endif // LIBPARI 4028 return giac_ifactors(n0,contextptr); 4029 } 4030 ifactors(const gen & n0,GIAC_CONTEXT)4031 vecteur ifactors(const gen & n0,GIAC_CONTEXT){ 4032 gen n(n0); 4033 vecteur f=pfacprem(n,false,contextptr); 4034 if (is_undef(f)) 4035 return f; 4036 vecteur g=ifactors1(n,contextptr); 4037 if (is_undef(g)) 4038 return g; 4039 return mergevecteur(f,g); 4040 } 4041 ifactors(const gen & r,const gen & i,const gen & ri,GIAC_CONTEXT)4042 vecteur ifactors(const gen & r,const gen & i,const gen & ri,GIAC_CONTEXT){ 4043 gen norm=r*r+i*i; 4044 gen reste(ri); 4045 const vecteur & facto = ifactors(norm,contextptr); 4046 if (is_undef(facto)) 4047 return facto; 4048 int l=int(facto.size())/2; 4049 vecteur res; 4050 for (int i=0;i<l;++i){ 4051 gen prime=facto[2*i]; 4052 int mult=facto[2*i+1].val,multp=0; 4053 int n=smod(prime,4).val; 4054 if (n==2){ 4055 res.push_back(1+cst_i); 4056 res.push_back(mult); 4057 reste=reste/pow(1+cst_i,mult,contextptr); 4058 continue; 4059 } 4060 if (n==-1){ 4061 res.push_back(prime); 4062 res.push_back(mult/2); 4063 reste=reste/pow(prime,mult/2,contextptr); 4064 continue; 4065 } 4066 prime=pa2b2(prime,contextptr); 4067 prime=gen(prime[0],prime[1]); 4068 for (;mult>0;--mult,++multp){ 4069 if (!is_zero(reste % prime)) 4070 break; 4071 reste=reste/prime; 4072 } 4073 if (multp){ 4074 res.push_back(prime); 4075 res.push_back(multp); 4076 } 4077 if (mult){ 4078 prime=conj(prime,contextptr); 4079 res.push_back(prime); 4080 res.push_back(mult); 4081 reste=reste/pow(prime,mult,contextptr); 4082 } 4083 } 4084 if (!is_one(reste)){ 4085 res.insert(res.begin(),1); 4086 res.insert(res.begin(),reste); 4087 } 4088 return res; 4089 } 4090 ifactors(const gen & args,int maplemode,GIAC_CONTEXT)4091 gen ifactors(const gen & args,int maplemode,GIAC_CONTEXT){ 4092 if ( (args.type==_INT_) || (args.type==_ZINT)){ 4093 if (is_zero(args)){ 4094 if (maplemode==1) 4095 return makevecteur(args,vecteur(0)); 4096 else 4097 return makevecteur(args); 4098 } 4099 vecteur v(ifactors(abs(args,contextptr),contextptr)); // ok 4100 if (!v.empty() && is_undef(v.front())) 4101 return v.front(); 4102 if (maplemode!=1){ 4103 if (is_positive(args,context0)) 4104 return v; 4105 return mergevecteur(makevecteur(minus_one,plus_one),v); 4106 } 4107 vecteur res; 4108 const_iterateur it=v.begin(),itend=v.end(); 4109 for (;it!=itend;it+=2){ 4110 res.push_back(makevecteur(*it,*(it+1))); 4111 } 4112 if (is_positive(args,context0)) 4113 return makevecteur(plus_one,res); 4114 else 4115 return makevecteur(minus_one,res); 4116 } 4117 if (args.type==_CPLX && is_integer(*args._CPLXptr) && is_integer(*(args._CPLXptr+1))) 4118 return ifactors(*args._CPLXptr,*(args._CPLXptr+1),args,contextptr); 4119 return gentypeerr(gettext("ifactors")); 4120 } 4121 _ifactors(const gen & args,GIAC_CONTEXT)4122 gen _ifactors(const gen & args,GIAC_CONTEXT){ 4123 if ( args.type==_STRNG && args.subtype==-1) return args; 4124 if (args.type==_VECT && args.subtype==_SEQ__VECT && args._VECTptr->size()==2 && args._VECTptr->back()==at_matrix){ 4125 gen g=args._VECTptr->front(); 4126 g=_ifactors(g,contextptr); 4127 if (g.type!=_VECT || g._VECTptr->size()%2) 4128 return g; 4129 return _matrix(makesequence(g._VECTptr->size()/2,2,g),contextptr); 4130 } 4131 if (args.type==_VECT) 4132 return apply(args,_ifactors,contextptr); 4133 gen g(args); 4134 if (!is_integral(g)) 4135 return gensizeerr(contextptr); 4136 if (calc_mode(contextptr)==1){ // ggb returns factors repeted instead of multiplicites 4137 vecteur res; 4138 gen in=ifactors(g,0,contextptr); 4139 if (in.type==_VECT){ 4140 for (unsigned i=0;i<in._VECTptr->size();i+=2){ 4141 gen f=in[i],m=in[i+1]; 4142 if (m.type==_INT_){ 4143 for (int j=0;j<m.val;++j) 4144 res.push_back(f); 4145 } 4146 } 4147 return res; 4148 } 4149 } 4150 return ifactors(g,0,contextptr); 4151 } 4152 static const char _ifactors_s []="ifactors"; 4153 static define_unary_function_eval (__ifactors,&_ifactors,_ifactors_s); 4154 define_unary_function_ptr5( at_ifactors ,alias_at_ifactors,&__ifactors,0,true); 4155 4156 static const char _facteurs_premiers_s []="facteurs_premiers"; 4157 static define_unary_function_eval (__facteurs_premiers,&_ifactors,_facteurs_premiers_s); 4158 define_unary_function_ptr5( at_facteurs_premiers ,alias_at_facteurs_premiers,&__facteurs_premiers,0,true); 4159 4160 static const char _maple_ifactors_s []="maple_ifactors"; _maple_ifactors(const gen & args,GIAC_CONTEXT)4161 gen _maple_ifactors(const gen & args,GIAC_CONTEXT){ 4162 if ( args.type==_STRNG && args.subtype==-1) return args; 4163 if (args.type==_VECT) 4164 return apply(args,_maple_ifactors,contextptr); 4165 return ifactors(args,1,contextptr); 4166 } 4167 static define_unary_function_eval (__maple_ifactors,&_maple_ifactors,_maple_ifactors_s); 4168 define_unary_function_ptr5( at_maple_ifactors ,alias_at_maple_ifactors,&__maple_ifactors,0,true); 4169 in_factors(const gen & gf,GIAC_CONTEXT)4170 static vecteur in_factors(const gen & gf,GIAC_CONTEXT){ 4171 if (gf.type!=_SYMB) 4172 return makevecteur(gf,plus_one); 4173 unary_function_ptr & u=gf._SYMBptr->sommet; 4174 if (u==at_inv){ 4175 vecteur v=in_factors(gf._SYMBptr->feuille,contextptr); 4176 iterateur it=v.begin(),itend=v.end(); 4177 for (;it!=itend;it+=2) 4178 *(it+1)=-*(it+1); 4179 return v; 4180 } 4181 if (u==at_neg){ 4182 vecteur v=in_factors(gf._SYMBptr->feuille,contextptr); 4183 v.push_back(minus_one); 4184 v.push_back(plus_one); 4185 return v; 4186 } 4187 if ( (u==at_pow) && (gf._SYMBptr->feuille._VECTptr->back().type==_INT_) ){ 4188 vecteur v=in_factors(gf._SYMBptr->feuille._VECTptr->front(),contextptr); 4189 gen k=gf._SYMBptr->feuille._VECTptr->back(); 4190 iterateur it=v.begin(),itend=v.end(); 4191 for (;it!=itend;it+=2) 4192 *(it+1)=k* *(it+1); 4193 return v; 4194 } 4195 if (u!=at_prod) 4196 return makevecteur(gf,plus_one); 4197 vecteur res; 4198 const_iterateur it=gf._SYMBptr->feuille._VECTptr->begin(),itend=gf._SYMBptr->feuille._VECTptr->end(); 4199 for (;it!=itend;++it){ 4200 res=mergevecteur(res,in_factors(*it,contextptr)); 4201 } 4202 return res; 4203 } in_factors1(const vecteur & res,GIAC_CONTEXT)4204 static vecteur in_factors1(const vecteur & res,GIAC_CONTEXT){ 4205 gen coeff(1); 4206 vecteur v; 4207 const_iterateur it=res.begin(),itend=res.end(); 4208 for (;it!=itend;it+=2){ 4209 if (lidnt(*it).empty()) 4210 coeff=coeff*(pow(*it,*(it+1),contextptr)); 4211 else 4212 v.push_back(makevecteur(*it,*(it+1))); 4213 } 4214 return makevecteur(coeff,v); 4215 } factors(const gen & g,const gen & x,GIAC_CONTEXT)4216 vecteur factors(const gen & g,const gen & x,GIAC_CONTEXT){ 4217 gen gf=factor(g,x,false,contextptr); 4218 vecteur res=in_factors(gf,contextptr); 4219 if (xcas_mode(contextptr)!=1) 4220 return res; 4221 return in_factors1(res,contextptr); 4222 } sqff_factors(const gen & g,GIAC_CONTEXT)4223 vecteur sqff_factors(const gen & g,GIAC_CONTEXT){ 4224 gen gf=_sqrfree(g,contextptr); 4225 return in_factors(gf,contextptr); 4226 } 4227 static const char _factors_s []="factors"; _factors(const gen & args,GIAC_CONTEXT)4228 gen _factors(const gen & args,GIAC_CONTEXT){ 4229 if ( args.type==_STRNG && args.subtype==-1) return args; 4230 if (args.type==_VECT && args.subtype==_SEQ__VECT && args._VECTptr->size()>=2 && args._VECTptr->back()==at_matrix){ 4231 gen g; 4232 if (args._VECTptr->size()==2) 4233 g=args._VECTptr->front(); 4234 else 4235 g=gen(vecteur(args._VECTptr->begin(),args._VECTptr->end()-1),_SEQ__VECT); 4236 g=_factors(g,contextptr); 4237 if (g.type!=_VECT || g._VECTptr->size()%2) 4238 return g; 4239 return _matrix(makesequence(g._VECTptr->size()/2,2,g),contextptr); 4240 } 4241 if (args.type==_VECT && args.subtype==_POLY1__VECT){ 4242 gen x(identificateur("xfactors")); 4243 gen res=_poly2symb(makesequence(args,x),contextptr); 4244 res=_factors(res,contextptr); 4245 if (res.type==_VECT && res._VECTptr->size()==2){ 4246 vecteur v(*res._VECTptr); 4247 for (size_t i=0;i<v.size();i+=2){ 4248 v[i]=_symb2poly(makesequence(v[i],x),contextptr); 4249 } 4250 return v; 4251 } 4252 return res; 4253 } 4254 if (args.type==_VECT && args.subtype==_SEQ__VECT && args._VECTptr->size()==2){ 4255 gen j=args._VECTptr->back(); 4256 gen res=_factors(args._VECTptr->front()*j,contextptr); 4257 if (res.type==_VECT && xcas_mode(contextptr)!=1) 4258 res=in_factors1(*res._VECTptr,contextptr); 4259 if (res.type==_VECT && res._VECTptr->size()==2){ 4260 res._VECTptr->front()=recursive_normal(res._VECTptr->front()/j,contextptr); 4261 if (xcas_mode(contextptr)!=1){ 4262 if (is_one(res._VECTptr->front())) 4263 res=res._VECTptr->back(); 4264 else { 4265 j=res._VECTptr->front(); 4266 res=res._VECTptr->back(); 4267 if (res.type==_VECT) 4268 res=mergevecteur(makevecteur(j,1),*res._VECTptr); 4269 } 4270 vecteur v; 4271 aplatir(*res._VECTptr,v,contextptr); 4272 res=v; 4273 } 4274 } 4275 return res; 4276 } 4277 if (args.type==_VECT) 4278 return apply(args,_factors,contextptr); 4279 return factors(args,vx_var,contextptr); 4280 } 4281 static define_unary_function_eval (__factors,&_factors,_factors_s); 4282 define_unary_function_ptr5( at_factors ,alias_at_factors,&__factors,0,true); 4283 ifactors2ifactor(const vecteur & l,bool quote)4284 static gen ifactors2ifactor(const vecteur & l,bool quote){ 4285 int s; 4286 s=int(l.size()); 4287 gen r; 4288 vecteur v(s/2); 4289 for (int j=0;j<s;j=j+2){ 4290 if (!is_one(l[j+1])) 4291 v[j/2]=symbolic(at_pow,gen(makevecteur(l[j],l[j+1]),_SEQ__VECT)); 4292 else 4293 v[j/2]=l[j]; 4294 } 4295 if (v.size()==1){ 4296 #if defined(GIAC_HAS_STO_38) && defined(CAS38_DISABLED) 4297 return symb_quote(v.front()); 4298 #else 4299 if (quote) 4300 return symb_quote(v.front()); 4301 return v.front(); 4302 #endif 4303 } 4304 r=symbolic(at_prod,gen(v,_SEQ__VECT)); 4305 #if defined(GIAC_HAS_STO_38) && defined(CAS38_DISABLED) 4306 r=symb_quote(r); 4307 #endif 4308 if (quote) 4309 return symb_quote(r); 4310 return r; 4311 } ifactor(const gen & n,GIAC_CONTEXT)4312 gen ifactor(const gen & n,GIAC_CONTEXT){ 4313 vecteur l; 4314 l=ifactors(n,contextptr); 4315 if (!l.empty() && is_undef(l.front())) return l.front(); 4316 return ifactors2ifactor(l,calc_mode(contextptr)==1); 4317 } _ifactor(const gen & args,GIAC_CONTEXT)4318 gen _ifactor(const gen & args,GIAC_CONTEXT){ 4319 if ( args.type==_STRNG && args.subtype==-1) return args; 4320 if (args.type==_CPLX && is_integer(*args._CPLXptr) && is_integer(*(args._CPLXptr+1))){ 4321 const vecteur & v=ifactors(*args._CPLXptr,*(args._CPLXptr+1),args,contextptr); 4322 return ifactors2ifactor(v,calc_mode(contextptr)==1); 4323 } 4324 gen n=args; 4325 if (n.type==_VECT && n._VECTptr->size()==1 && is_integer(n._VECTptr->front())) 4326 return ifactor(n,contextptr); 4327 if (n.type==_VECT) 4328 return apply(n,_ifactor,contextptr); 4329 if (!is_integral(n)) 4330 return gensizeerr(contextptr); 4331 if (is_strictly_positive(-n,0)) 4332 return -_ifactor(-n,contextptr); 4333 if (n.type==_INT_ && n.val<=3) 4334 return n; 4335 return ifactor(n,contextptr); 4336 } 4337 static const char _ifactor_s []="ifactor"; 4338 static define_unary_function_eval (__ifactor,&_ifactor,_ifactor_s); 4339 define_unary_function_ptr5( at_ifactor ,alias_at_ifactor,&__ifactor,0,true); 4340 4341 static const char _factoriser_entier_s []="factoriser_entier"; 4342 static define_unary_function_eval (__factoriser_entier,&_ifactor,_factoriser_entier_s); 4343 define_unary_function_ptr5( at_factoriser_entier ,alias_at_factoriser_entier,&__factoriser_entier,0,true); 4344 divis(const vecteur & l3,GIAC_CONTEXT)4345 static vecteur divis(const vecteur & l3,GIAC_CONTEXT){ 4346 vecteur l1(1); 4347 gen d,e; 4348 int s=int(l3.size()); 4349 gen taille=1; 4350 for (int k=0;k<s;k+=2){ 4351 taille=taille*(l3[k+1]+1); 4352 } 4353 if (taille.type!=_INT_ || taille.val>LIST_SIZE_LIMIT) 4354 return vecteur(1,gendimerr(contextptr)); 4355 l1.reserve(taille.val); 4356 l1[0]=1;//l3.push_back(..); 4357 for (int k=0;k<s;k=k+2) { 4358 vecteur l2; 4359 l2.reserve(taille.val); 4360 int s1; 4361 s1=int(l1.size()); 4362 vecteur l4(s1); 4363 d=l3[k]; 4364 e=l3[k+1]; 4365 int ei; 4366 if (e.type==_INT_){ 4367 ei=e.val; 4368 } 4369 else 4370 return vecteur(1,gensizeerr(gettext("Integer too large"))); 4371 for (int j=1;j<=ei;j++){ 4372 gen dj=pow(d,j); 4373 for (int l=0;l<s1;l++){ 4374 l4[l]=l1[l]*dj; 4375 } 4376 // l2=mergevecteur(l2,l4); 4377 iterateur it=l4.begin(),itend=l4.end(); 4378 for (;it!=itend;++it) 4379 l2.push_back(*it); 4380 } 4381 // l1=mergevecteur(l1,l2); 4382 iterateur it=l2.begin(),itend=l2.end(); 4383 for (;it!=itend;++it) 4384 l1.push_back(*it); 4385 } 4386 return(l1); 4387 } idivis(const gen & n,GIAC_CONTEXT)4388 gen idivis(const gen & n,GIAC_CONTEXT){ 4389 vecteur l3(ifactors(n,contextptr)); 4390 if (!l3.empty() && is_undef(l3.front())) return l3.front(); 4391 return divis(l3,contextptr); 4392 } _idivis(const gen & args,GIAC_CONTEXT)4393 gen _idivis(const gen & args,GIAC_CONTEXT){ 4394 if ( args.type==_STRNG && args.subtype==-1) return args; 4395 if (args.type==_VECT) 4396 return apply(args,_idivis,contextptr); 4397 gen n=args; 4398 if (is_zero(n) || (!is_integral(n) && !is_integer(n)) || n.type==_CPLX) 4399 return gentypeerr(contextptr); 4400 return _sort(idivis(abs(n,contextptr),contextptr),contextptr); 4401 } 4402 static const char _idivis_s []="idivis"; 4403 static define_unary_function_eval (__idivis,&_idivis,_idivis_s); 4404 define_unary_function_ptr5( at_idivis ,alias_at_idivis,&__idivis,0,true); 4405 _divis(const gen & args,GIAC_CONTEXT)4406 gen _divis(const gen & args,GIAC_CONTEXT){ 4407 if ( args.type==_STRNG && args.subtype==-1) return args; 4408 if (args.type==_VECT) 4409 return apply(args,_divis,contextptr); 4410 return divis(factors(args,vx_var,contextptr),contextptr); 4411 } 4412 static const char _divis_s []="divis"; 4413 static define_unary_function_eval (__divis,&_divis,_divis_s); 4414 define_unary_function_ptr5( at_divis ,alias_at_divis,&__divis,0,true); 4415 4416 /* 4417 gen ichinreme(const vecteur & a,const vecteur & b){ 4418 vecteur r(2); 4419 gen p=a[1],q=b[1],u,v,d; 4420 egcd(p,q,u,v,d); 4421 if (d!=1) return gensizeerr(contextptr); 4422 r[0]=(u*p*b[0]+v*q*a[0]%p*q); 4423 r[1]=p*q; 4424 return(r); 4425 } 4426 gen _ichinreme(const gen & args){ 4427 if ( args){ 4428 if ( (args.type!=_VECT) || (args._VECTptr->size()!=4) ) 4429 return gensizeerr(contextptr); 4430 vecteur a(2).type==_STRNG && args.subtype==-1{ 4431 if ( (args.type!=_VECT) || (args._VECTptr->size()!=4) ) 4432 return gensizeerr(contextptr); 4433 vecteur a(2))) return args){ 4434 if ( (args.type!=_VECT) || (args._VECTptr->size()!=4) ) 4435 return gensizeerr(contextptr); 4436 vecteur a(2); 4437 if ( (args.type!=_VECT) || (args._VECTptr->size()!=4) ) 4438 return gensizeerr(contextptr); 4439 vecteur a(2),b(2); 4440 a[0]=args[0]; 4441 a[1]=args[1]; 4442 b[0]=args[2]; 4443 b[1]=args[3]; 4444 //gen a=args[0],p=args[1], b=args[2],q=args[3]; 4445 return ichinreme(a,b); 4446 } 4447 static const char _ichinreme_s []="ichinreme"; 4448 static define_unary_function_eval (__ichinreme,&_ichinreme,_ichinreme_s); 4449 define_unary_function_ptr5( at_ichinreme ,alias_at_ichinreme,&__ichinreme,0,true); 4450 */ 4451 euler(const gen & e,GIAC_CONTEXT)4452 gen euler(const gen & e,GIAC_CONTEXT){ 4453 if (e==0) 4454 return e; 4455 vecteur v(ifactors(e,contextptr)); 4456 if (!v.empty() && is_undef(v.front())) return v.front(); 4457 const_iterateur it=v.begin(),itend=v.end(); 4458 for (gen res(plus_one);;){ 4459 if (it==itend) 4460 return res; 4461 gen p=*it; 4462 ++it; 4463 int n=it->val; 4464 res = res * (p-plus_one)*pow(p,n-1); 4465 ++it; 4466 } 4467 } 4468 static const char _euler_s []="euler"; _euler(const gen & args,GIAC_CONTEXT)4469 gen _euler(const gen & args,GIAC_CONTEXT){ 4470 if ( args.type==_STRNG && args.subtype==-1) return args; 4471 if (args.type==_VECT) 4472 return apply(args,_euler,contextptr); 4473 if ( is_integer(args) && is_positive(args,contextptr)) 4474 return euler(args,contextptr); 4475 return gentypeerr(contextptr); 4476 } 4477 static define_unary_function_eval (__euler,&_euler,_euler_s); 4478 define_unary_function_ptr5( at_euler ,alias_at_euler,&__euler,0,true); 4479 pa2b2(const gen & p,GIAC_CONTEXT)4480 gen pa2b2(const gen & p,GIAC_CONTEXT){ 4481 if (p==2) 4482 return makevecteur(1,1); 4483 if (!is_integer(p) || (p%4)!=1 || is_greater(1,p,contextptr)) return gensizeerr(contextptr);// car p!=1 mod 4 4484 gen q=(p-1)/4; 4485 gen a=2; 4486 gen ra; 4487 ra=powmod(a,q,p); 4488 //on cherche ra^2=-1 mod p avec ra!=1 et ra !=p-1 4489 while ((a!=p-1) && ((ra==1)|| (ra==p-1))){ 4490 a=a+1; 4491 ra=powmod(a,q,p); 4492 } 4493 if ((ra==1)||(ra==p-1)) return gensizeerr(contextptr);//car p n'est pas premier 4494 gen ux=1,uy=ra,vx=0,vy=p,wx,wy; 4495 gen m=1; 4496 while(m!=0){ 4497 if (is_positive(vx*vx+vy*vy-ux*ux+uy*uy,0)){ 4498 //on echange u et v 4499 wx=vx; 4500 wy=vy; 4501 vx=ux; 4502 vy=uy; 4503 ux=wx; 4504 uy=wy; 4505 } 4506 gen alpha=inv(2,contextptr)-(ux*vx+uy*vy)*inv(vx*vx+vy*vy,contextptr); 4507 //m=partie entiere de alpha (-v.v/2<(u+mv).v<=v.v/2) 4508 m=_floor(alpha,contextptr); 4509 ux=ux+m*vx; 4510 uy=uy+m*vy; 4511 } 4512 vecteur v(2); 4513 //v repond a la question 4514 v[0]=abs(vx,contextptr); // ok 4515 v[1]=abs(vy,contextptr); // ok 4516 if (vx*vx+vy*vy!=p) 4517 return gensizeerr(contextptr); 4518 return v; 4519 } _pa2b2(const gen & args,GIAC_CONTEXT)4520 gen _pa2b2(const gen & args,GIAC_CONTEXT){ 4521 if ( args.type==_STRNG && args.subtype==-1) return args; 4522 if (!is_integer(args)) 4523 return gensizeerr(contextptr); 4524 gen n=args; 4525 return pa2b2(n,contextptr); 4526 } 4527 static const char _pa2b2_s []="pa2b2"; 4528 static define_unary_function_eval (__pa2b2,&_pa2b2,_pa2b2_s); 4529 define_unary_function_ptr5( at_pa2b2 ,alias_at_pa2b2,&__pa2b2,0,true); 4530 ipropfrac(const gen & a,const gen & b,GIAC_CONTEXT)4531 static gen ipropfrac(const gen & a,const gen & b,GIAC_CONTEXT){ 4532 if (!is_integer(a) || !is_integer(b)) 4533 return gensizeerr(contextptr); 4534 gen r=a%b; 4535 gen q=(a-r)/b; 4536 gen d=gcd(r,b); 4537 r=r/d; 4538 gen b1=b/d; 4539 if (r==0) 4540 return q; 4541 gen v; 4542 v=symbolic(at_division,gen(makevecteur(r,b1),_SEQ__VECT)); 4543 gen w; 4544 w=symbolic(at_plus,gen(makevecteur(q,v),_SEQ__VECT)); 4545 if (calc_mode(contextptr)==1) 4546 return symbolic(at_quote,w); 4547 return w; 4548 } _propfrac(const gen & arg,GIAC_CONTEXT)4549 gen _propfrac(const gen & arg,GIAC_CONTEXT){ 4550 if ( arg.type==_STRNG && arg.subtype==-1) return arg; 4551 gen args(arg); 4552 vecteur v; 4553 if (arg.type==_VECT && arg._VECTptr->size()==2){ 4554 v=vecteur(1,arg._VECTptr->back()); 4555 args=arg._VECTptr->front(); 4556 lvar(args,v); 4557 } 4558 else 4559 v=lvar(arg); 4560 gen g=e2r(args,v,contextptr); 4561 gen a,b; 4562 fxnd(g,a,b); 4563 if (v.empty()) 4564 return ipropfrac(a,b,contextptr); 4565 else { 4566 gen d=r2e(b,v,contextptr); 4567 g=_quorem(makesequence(r2e(a,v,contextptr),d,v.front()),contextptr); 4568 if (is_undef(g)) return g; 4569 vecteur &v=*g._VECTptr; 4570 return v[0]+rdiv(v[1],d,contextptr); 4571 } 4572 } 4573 static const char _propfrac_s []="propfrac"; 4574 static define_unary_function_eval (__propfrac,&_propfrac,_propfrac_s); 4575 define_unary_function_ptr5( at_propfrac ,alias_at_propfrac,&__propfrac,0,true); 4576 step_egcd(int a,int b,GIAC_CONTEXT)4577 void step_egcd(int a,int b,GIAC_CONTEXT){ 4578 gprintf("===============",vecteur(0),1,contextptr); 4579 gprintf("Extended Euclide algorithm for a=%gen and b=%gen",makevecteur(a,b),1,contextptr); 4580 gprintf("L%gen: 1*a+0*b=%gen",makevecteur(1,a),1,contextptr); 4581 gprintf("L%gen: 0*a+1*b=%gen",makevecteur(2,b),1,contextptr); 4582 int i=3; 4583 int u0=1,v0=0,u1=0,v1=1,u2,v2; 4584 for (;b;++i){ 4585 int q=a/b; 4586 u2=u0-q*u1; 4587 v2=v0-q*v1; 4588 int r=a-q*b; 4589 gprintf("iquo(%gen,%gen)=%gen",makevecteur(a,b,q),1,contextptr); 4590 gprintf("L%gen=L%gen-%gen*L%gen: %gen*a+%gen*b=%gen",makevecteur(i,i-2,q,i-1,u2,v2,r),1,contextptr); 4591 u0=u1; 4592 u1=u2; 4593 v0=v1; 4594 v1=v2; 4595 a=b; 4596 b=r; 4597 } 4598 gprintf("Bezout identity %gen*a+%gen*b=%gen",makevecteur(u0,v0,a),1,contextptr); 4599 } 4600 iabcuv(const gen & a,const gen & b,const gen & c,GIAC_CONTEXT)4601 gen iabcuv(const gen & a,const gen & b,const gen & c,GIAC_CONTEXT){ 4602 gen d=gcd(a,b); 4603 if (c%d!=0) return gensizeerr(gettext("No solution in ring")); 4604 gen a1=a/d,b1=b/d,c1=c/d; 4605 gen u,v,w; 4606 if (a1.type==_INT_ && b1.type==_INT_ && step_infolevel(contextptr)) 4607 step_egcd(a1.val,b1.val,contextptr); 4608 egcd(a1,b1,u,v,w); 4609 vecteur r(2); 4610 r[0]=smod(u*c1,b); 4611 r[1]=iquo(c-r[0]*a,b); 4612 return r; 4613 } _iabcuv(const gen & args,GIAC_CONTEXT)4614 gen _iabcuv(const gen & args,GIAC_CONTEXT){ 4615 if ( args.type==_STRNG && args.subtype==-1) return args; 4616 if ( (args.type!=_VECT) || (args._VECTptr->size()!=3) ) 4617 return gensizeerr(contextptr); 4618 gen a=args[0],b=args[1],c=args[2]; 4619 return iabcuv(a,b,c,contextptr); 4620 } 4621 static const char _iabcuv_s []="iabcuv"; 4622 static define_unary_function_eval (__iabcuv,&_iabcuv,_iabcuv_s); 4623 define_unary_function_ptr5( at_iabcuv ,alias_at_iabcuv,&__iabcuv,0,true); 4624 abcuv(const gen & a,const gen & b,const gen & c,const gen & x,GIAC_CONTEXT)4625 gen abcuv(const gen & a,const gen & b,const gen & c,const gen & x,GIAC_CONTEXT){ 4626 gen g=_egcd(makesequence(a,b,x),contextptr); 4627 if (is_undef(g)) return g; 4628 vecteur & v=*g._VECTptr; 4629 gen h=_quorem(makesequence(c,v[2],x),contextptr); 4630 if (is_undef(h)) return h; 4631 vecteur & w=*h._VECTptr; 4632 if (!is_zero(w[1])) 4633 return gensizeerr(gettext("No solution in ring")); 4634 gen U=v[0]*w[0],V=v[1]*w[0]; 4635 if (_degree(makesequence(c,x),contextptr).val<_degree(makesequence(a,x),contextptr).val+_degree(makesequence(b,x),contextptr).val ){ 4636 U=_rem(makesequence(U,b,x),contextptr); 4637 V=_rem(makesequence(V,a,x),contextptr); 4638 } 4639 return makevecteur(U,V); 4640 } _abcuv(const gen & args,GIAC_CONTEXT)4641 gen _abcuv(const gen & args,GIAC_CONTEXT){ 4642 if ( args.type==_STRNG && args.subtype==-1) return args; 4643 if ( (args.type!=_VECT) || (args._VECTptr->size()<3) ) 4644 return gensizeerr(contextptr); 4645 vecteur & v =*args._VECTptr; 4646 if (v.size()>3) 4647 return abcuv(v[0],v[1],v[2],v[3],contextptr); 4648 return abcuv(v[0],v[1],v[2],vx_var,contextptr); 4649 } 4650 static const char _abcuv_s []="abcuv"; 4651 static define_unary_function_eval (__abcuv,&_abcuv,_abcuv_s); 4652 define_unary_function_ptr5( at_abcuv ,alias_at_abcuv,&__abcuv,0,true); 4653 simp2(const gen & a,const gen & b,GIAC_CONTEXT)4654 gen simp2(const gen & a,const gen & b,GIAC_CONTEXT){ 4655 vecteur r(2); 4656 gen d=gcd(a,b); 4657 r[0]=normal(a/d,contextptr); 4658 r[1]=normal(b/d,contextptr); 4659 return r; 4660 } _simp2(const gen & args,GIAC_CONTEXT)4661 gen _simp2(const gen & args,GIAC_CONTEXT){ 4662 if ( args.type==_STRNG && args.subtype==-1) return args; 4663 if ( (args.type!=_VECT) || (args._VECTptr->size()!=2) ) 4664 return gensizeerr(contextptr); 4665 gen a=args[0],b=args[1]; 4666 if ( (a.type==_VECT) || (b.type==_VECT) ) 4667 return gensizeerr(contextptr); 4668 return simp2(a,b,contextptr); 4669 } 4670 static const char _simp2_s []="simp2"; 4671 static define_unary_function_eval (__simp2,&_simp2,_simp2_s); 4672 define_unary_function_ptr5( at_simp2 ,alias_at_simp2,&__simp2,0,true); 4673 fxnd(const gen & a)4674 gen fxnd(const gen & a){ 4675 vecteur v(lvar(a)); 4676 gen g=e2r(a,v,context0); // ok 4677 gen n,d; 4678 fxnd(g,n,d); 4679 return makevecteur(r2e(n,v,context0),r2e(d,v,context0)); // ok 4680 } _fxnd(const gen & args,GIAC_CONTEXT)4681 gen _fxnd(const gen & args,GIAC_CONTEXT){ 4682 if ( args.type==_STRNG && args.subtype==-1) return args; 4683 if (args.type==_VECT) 4684 return apply(args,fxnd); 4685 return fxnd(args); 4686 } 4687 static const char _fxnd_s []="fxnd"; 4688 static define_unary_function_eval (__fxnd,&_fxnd,_fxnd_s); 4689 define_unary_function_ptr5( at_fxnd ,alias_at_fxnd,&__fxnd,0,true); 4690 4691 #ifndef NO_NAMESPACE_GIAC 4692 } // namespace giac 4693 #endif // ndef NO_NAMESPACE_GIAC 4694 4695