/*! @file
* @brief
* contains factoring algorithms which are related to Fermat's method
*/
void fermat_like_method()
{
/* Fermat-ähnliche Methode:
Sei N=p*q (hier ausnahmsweise p,q nichtnotwendigerweise prim)
und sei M=m1*m2 ein Multiplikator, so dass m1*q ungefähr m2*p entspricht.
Dann ist (m1*q)*(m2*p)=(m1*m2)*(q*p)=M*N.
Wir definieren R=floor(sqrt(M*N)).
Der Ansatz ist nun, ein delta zu finden, so dass Q := (R+delta)^2 - M*N
eine Quadratzahl ergibt, denn dann folgt (R+delta)^2 - Q = 0 (mod N),
was gleichbedeutend ist zu (R+delta+sqrt(Q))*(R+delta-sqrt(Q)) = 0 (mod N).
Mit gcd(R+delta+sqrt(Q),N) wäre dann ein Teiler von N gefunden...
Diese Methode ist zwar nicht besonders effizient, da sie eher einer
Probedivision entspricht; sie hilft aber, solche Zahlen besonders
schnell zu faktorisieren, die auf "allzu naive" Weise als
"besonders schwer" zu faktorisierende Zahlen konstruiert sind.
So können zum Beispiel Produkte der Form
N:= p*nextprime(10*p+"kleiner Offset")
effizient zerlegt werden.
*/
#ifdef VERBOSE_NOTICE
cout << "starting Fermat factorization method" << endl;
#endif
if (mpz_sizeinbase(n,10)<50)
{
cout << "number too small, it's not worthwhile do proceed; skipping fermat..." << endl;
return;
}
const unsigned int Witnesses = 4;
const unsigned int Witness[Witnesses] = { 23*41*59, 35*29*37, 31*43*47, 61*67 };
const unsigned int Maxbits = 65536;
unsigned char squareflags[Maxbits];
// initialize bitfield for witnessing non-squares via quadratic residues
for (unsigned int j=0; j<Maxbits; ++j) squareflags[j]=0;
for (unsigned int i=0; i<Witnesses; ++i)
for (unsigned int j=0; j<Witness[i]; ++j) squareflags[(j*j)%Witness[i]] |= 1<<i;
mpz_t c,r,h;
mpz_init(c); mpz_init(r); mpz_init(h);
const unsigned int k_max = 0x100; // or 0x040 (to be less exhaustive)
for (unsigned int k=0; k<k_max; ++k)
{
#ifdef REACT_ON_SIGUSR
if (USRSignalHandler.got_SIGUSR1()) break;
if (USRSignalHandler.got_SIGUSR2()) break;
#endif
mpz_set_ui(c,2*2*3*5*7*11*13*17*19); mpz_pow_ui(c,c,4);
unsigned int m = 1;
if (k&0x001) m*=2;
if (k&0x002) m*=3;
if (k&0x004) m*=5;
if (k&0x008) m*=7;
if (k&0x010) m*=11;
if (k&0x020) m*=13;
if (k&0x040) m*=17;
if (k&0x080) m*=19;
mpz_mul_ui(c,c,m);
// Idea behind these multipliers given by example:
// Let t be a (small) number, then
//
// t^0 t^1 t^2 t^3 t^4
// t^4 covers --- , --- , --- , --- , ---
// t^4 t^3 t^2 t^1 t^0
//
// which is identical to the sequence
//
// 1 1
// --- , --- , 1 , t^2 , t^4
// t^4 t^2
//
// because p*q*t^4 = p*(q*t^4) = (p*t)*(q*t^3) = (p*t^2)*(q*t^2) = ...
//
// 1 1 1
// Similarly t^5 covers --- , --- , - , t , t^3, t^5
// t^5 t^3 t
#if defined(VERBOSE_INFO)
cout << "using Fermat multiplier " << c << " \r" << flush;
#endif
mpz_mul(c,c,n); mpz_sqrt(r,c); mpz_mul(h,r,r); mpz_sub(h,h,c);
mpz_mul_ui(r,r,2);
struct { unsigned int h,r; } ring[Witnesses];
for (unsigned int i=0; i<Witnesses; ++i)
{
ring[i].h=mpz_remainder_ui(h,Witness[i]);
ring[i].r=mpz_remainder_ui(r,Witness[i]);
//cout << "Witness " << i << ": " << ring[i].h << ", " << ring[i].r << endl;
}
for (unsigned int delta=1; delta<50000; ++delta)
{
//if (delta%10000==0) cout << "delta=" << delta << "\r" << flush;
mpz_add(h,h,r); mpz_add_ui(h,h,1); mpz_add_ui(r,r,2); // compute next square: h=(x+1)^2=x^2+2x+1=h+r+1 (using r=2x)
for (unsigned int i=0; i<Witnesses; ++i)
{
ring[i].h+=ring[i].r+1; if (ring[i].h>=Witness[i]) ring[i].h-=Witness[i];
ring[i].r+=2; if (ring[i].r>=Witness[i]) ring[i].r-=Witness[i];
}
if (!(squareflags[ring[0].h] & 0x01)) continue;
if (!(squareflags[ring[1].h] & 0x02)) continue;
if (!(squareflags[ring[2].h] & 0x04)) continue;
if (!(squareflags[ring[3].h] & 0x08)) continue;
//cout << "delta=" << delta << endl;
if (mpz_perfect_square_p(h))
{
#if 1 || defined(VERBOSE_INFO)
cout << endl << "delta=" << delta << ": Square detected..." << endl;
#endif
mpz_div_ui(r,r,2); mpz_sqrt(h,h); mpz_add(h,h,r);
mpz_gcd(h,h,n);
if (mpz_cmp_ui(h,1)!=0)
{
mpz_divexact(n,n,h);
//cout << "t1=" << h << endl;
//cout << "t2=" << n << endl;
if (mpz_probab_prime_p(h,probab_prime_checks))
{
#if 0
// only for debugging...
{
mpz_t sum,N,x;
mpz_init(sum); mpz_init(N); mpz_init(x);
mpz_add(sum,n,h); mpz_div_ui(sum,sum,2); mpz_mul(N,n,h);
mpz_sqrt(x,N); mpz_sub(sum,sum,x);
cout << "phimat-value --> " << sum << endl;
mpz_clear(x); mpz_clear(N); mpz_clear(sum);
}
#endif
cout << h << " is factor." << endl;
std::ostringstream comment;
comment << " [fer]";
Factorization_to_file << MAL(h,comment) << flush;
}
else
{
fermat_like_method(); // try again...
mpz_swap(n,h);
if (mpz_probab_prime_p(h,probab_prime_checks))
{
cout << h << " is factor." << endl;
std::ostringstream comment;
comment << " [fer]";
Factorization_to_file << MAL(h,comment) << flush;
}
else
{
cout << h << " is a composite factor." << endl;
if (Potenztest(h)) Factorization_to_file << " [fer]" << flush;
else
{
std::ostringstream comment;
comment << " [fer] [composite]";
Factorization_to_file << MAL(h,comment) << flush;
}
}
}
if (!mpz_probab_prime_p(h,probab_prime_checks)) fermat_like_method(); // and try again once more...
k=k_max; break;
}
}
}
}
#ifdef VERBOSE_INFO
cout << endl;
#endif
mpz_clear(h); mpz_clear(r); mpz_clear(c);
}
#include "utils.H"
#include <vector>
#include <algorithm>
class entry
{
private:
mp_limb_t hashval;
friend class phimahashvecs;
static unsigned int iL;
static mpz_t L;
public:
static mpz_t Base;
entry(const unsigned int index)
{
if (index==iL+1)
{
// shortcut, inc case that items are inserted "in order"
++iL; mpz_mul(L,L,Base); mpz_mod(L,L,n);
}
else
{
mpz_powm_ui(L,Base,index,n);
}
hashval=mpz_getlimbn(L,0);
}
~entry() { }
const mpz_t& get_mpz(const unsigned int index) const
{
if (index==iL) return L;
else
{
mpz_powm_ui(L,Base,index,n);
return L;
}
}
mp_limb_t get_hash() const { return hashval; }
bool operator< (const entry &rhs) const
{
// important: we do only compare the hash values, not the multiple
// precision numbers!
return (hashval<rhs.hashval);
}
bool operator== (const entry &rhs) const
{
return (hashval==rhs.hashval);
}
signed int compare(const mp_limb_t rhs) const
{
// compare based on hash value
return (hashval<rhs) ? -1 : (hashval>rhs) ? 1 : 0;
}
};
class phimavec : private std::vector<entry>
{
private:
friend class phimahashvecs;
double d_linear_approx_multiplicator;
public:
phimavec() : d_linear_approx_multiplicator(0.0) { }
~phimavec() { }
void sort()
{
//cout << "sorting " << size() << " items." << endl;
std::sort(begin(),end());
//size_t s=size();
erase(std::unique(begin(),end()),end()); // remove duplicate hash values
//if (size()!=s) cout << s-size() << " duplicates removed." << endl;
d_linear_approx_multiplicator = static_cast<double>(size())/(at(size()-1).get_hash()-at(0).get_hash());
}
bool found(const mpz_t &val) const
{
#ifdef VERBOSE_INFO
static unsigned int cs = 0;
#endif
// binary search
signed int L=0, R=size()-1;
const mp_limb_t hash = mpz_getlimbn(val,0);
#if 1
// linear approximation of the index
// (uniform distribution -> much better than binary search!)
const signed int d = static_cast<signed int>(static_cast<double>(hash) * d_linear_approx_multiplicator);
//if ( d<L || d>R ) return false; // should be outside the interval (but there are rounding errors...)
const signed int dd = 127;
__builtin_prefetch(&operator[](d),0,1);
if (d-dd>L && at(d-dd).get_hash()<hash) L=d-dd;
if (d+dd<R && at(d+dd).get_hash()>hash) R=d+dd;
//cout << "-->L: " << L << " R: " << R << " d: " << d << endl;
#endif
while (R>=L)
{
//cout << "searching " << L << " - " << R << endl;
const signed int pos = (L+R)>>1;
__builtin_prefetch(&operator[]((L+pos-1)>>1),0,0);
__builtin_prefetch(&operator[]((R+pos+1)>>1),0,0);
const signed int vgl = operator[](pos).compare(hash);
#ifdef VERBOSE_INFO
++cs;
#endif
if (vgl==0)
{
#ifdef VERBOSE_INFO
cout << cs << " comparisons performed so far" << endl;
#endif
return true;
}
if (vgl<0) L=pos+1; else R=pos-1;
}
//cout << "L: " << L << " R: " << R << " d: " << d << " size:" << size() << endl;
return false;
}
};
class phimahashvecs
{
private:
static const mp_size_t hashlimb = 1; // for mpz_getlimbn(<number>,hashlimb)
static const unsigned int lookup_bins = 0x1000;
static const mp_limb_t limbmask = lookup_bins-1;
phimavec lookup[lookup_bins];
size_t myIntervalsize;
const static unsigned int redundancy = 8;
public:
phimahashvecs() : myIntervalsize(0)
{
entry::iL=0; mpz_init_set_ui(entry::L,1);
}
~phimahashvecs() { }
size_t size() const
{
size_t s = 0;
for (unsigned int i=0; i<lookup_bins; ++i) s+=lookup[i].size();
return s;
}
size_t max_size() const { return lookup[0].max_size(); }
void insert(const unsigned int i)
{
const entry e(i);
const mp_limb_t h=mpz_getlimbn(e.get_mpz(i),hashlimb) & limbmask;
lookup[h].push_back(e);
}
void prepare(const size_t tablesize)
{
// for optimal RAM requirements we precalculate the number of items
// in each lookup bin.
cout << "optimizing memory layout" << endl;
entry::iL=0; mpz_set_ui(entry::L,1);
size_t bin_size[lookup_bins] = { 0 };
for (unsigned int i=0; i<tablesize; ++i)
{
const mp_limb_t h=mpz_getlimbn(entry::L,hashlimb) & limbmask;
++bin_size[h];
mpz_mul(entry::L,entry::L,entry::Base); mpz_mod(entry::L,entry::L,n);
}
entry::iL=0; mpz_set_ui(entry::L,1);
for (unsigned int i=0; i<lookup_bins; ++i)
{
//cout << "binsize[" << i << "]=" << bin_size[i] << endl;
lookup[i].reserve(bin_size[i]);
}
// and now actually insert the items
cout << "inserting items" << endl;
for (size_t i=0; i<tablesize; ++i) insert(i); // implicit constructor call
myIntervalsize = size()-redundancy; // important! call it before duplicate hash values are erased!
#ifdef VERBOSE_INFO
size_t s = 0;
for (unsigned int i=0; i<lookup_bins; ++i)
{
lookup[i].sort();
s+=lookup[i].size();
cout << s << " items prepared.\r" << flush;
}
cout << endl;
#else
for (unsigned int i=0; i<lookup_bins; ++i) lookup[i].sort();
#endif
}
size_t Intervalsize() const { return myIntervalsize; }
bool found(const mpz_t &val) const
{
// first hash into lookup, then do binary search inside this bin
const mp_limb_t hash = mpz_getlimbn(val,hashlimb) & limbmask;
if (!lookup[hash].found(val)) return false;
// now it gets more complicated:
// a possible candidate has been found based on equal hash values;
// perform up to <redundancy> additional checks to reduce the
// risk of false positives
//cout << "candidate found; checking..." << endl;
bool retval=true;
mpz_t x;
mpz_init_set(x,val);
for (unsigned int i=0; i<redundancy; ++i)
{
mpz_mul(x,x,entry::Base); mpz_mod(x,x,n);
const mp_limb_t hash = mpz_getlimbn(x,hashlimb) & limbmask;
retval=lookup[hash].found(x);
//cout << "candidate " << ((retval) ? "passed" : "failed") << " test " << i << endl;
if (!retval) break;
}
mpz_clear(x);
return retval;
}
bool found(const mpz_t &val, unsigned int &the_index) const
{
// first hash into lookup, then do binary search inside this bin
if (!found(val)) return false;
// a candidate has been found with high probability
// to determine "the_index", we perform a binary search!
mpz_t x;
mpz_init(x);
size_t d = Intervalsize()-1;
{
register unsigned int i = 0;
while (d>>=1) ++i;
d=1<<i;
}
cout << "candidate found; performing binary search for the_index" << endl;
size_t hhh=0;
while (d)
{
mpz_powm_ui(x,entry::Base,hhh+d,n);
mpz_invert(x,x,n); mpz_mul(x,x,val); mpz_mod(x,x,n);
if (found(x)) hhh+=d;
d>>=1;
}
mpz_powm_ui(x,entry::Base,hhh,n);
if (mpz_cmp(x,val)!=0)
{
cout << "binary search failed? try slowsearch..." << endl;
mpz_t y; mpz_init(y); mpz_invert(y,entry::Base,n);
while (hhh>0)
{
mpz_mul(x,x,y); mpz_mod(x,x,n); --hhh;
if (mpz_cmp(x,val)==0) break;
}
mpz_clear(y);
}
the_index=hhh;
cout << "the_index = " << the_index << endl;
mpz_clear(x);
return true;
}
};
unsigned int entry::iL = 0;
mpz_t entry::L, entry::Base;
extern unsigned long int allowed_memory_usage_KB(void);
class lambda_delta
{
private:
mpf_t SQRTN, DELTA;
mpf_t fres, f1res;
static const unsigned int bits = 2048;
public:
static void Delta_by_ratio(mpz_t &Delta, const mpz_t n, const double ratio_p=1.0, const double ratio_q=1.0)
{
mpf_set_default_prec(bits);
mpf_t x,y;
mpf_init(x); mpf_init(y);
mpf_set_d(x,ratio_p); mpf_set_d(y,ratio_q);
mpf_div(x,x,y); // full precision division (important!)
mpf_sqrt(x,x);
mpf_ui_div(y,1,x); mpf_add(x,x,y); mpf_sub_ui(x,x,2);
mpf_set_z(y,n); mpf_sqrt(y,y); mpf_mul(x,x,y);
mpz_set_f(Delta,x);
}
lambda_delta(const mpz_t n, const mpz_t delta)
{
mpf_set_default_prec(bits);
mpf_init(SQRTN); mpf_set_z(SQRTN,n); mpf_sqrt(SQRTN,SQRTN); // mpf_floor(SQRTN,SQRTN);
mpf_init(DELTA); mpf_set_z(DELTA,delta);
mpf_init(fres); mpf_init(f1res);
}
~lambda_delta()
{
mpf_clear(f1res); mpf_clear(fres);
mpf_clear(DELTA); mpf_clear(SQRTN);
}
const mpf_t& f(const mpf_t &x)
{
mpf_t h;
mpf_init(h);
mpf_sqrt(fres,x); mpf_ui_div(h,1,fres); mpf_add(fres,fres,h);
mpf_sub_ui(fres,fres,2);
mpf_mul(fres,fres,SQRTN);
mpf_sub(fres,fres,DELTA);
mpf_clear(h);
return fres;
}
const mpf_t& f1(const mpf_t &x)
{
mpf_t h;
mpf_init(h);
mpf_sqrt(f1res,x); mpf_ui_div(f1res,1,f1res);
mpf_pow_ui(h,f1res,3); mpf_sub(f1res,f1res,h);
mpf_div_ui(f1res,f1res,2);
mpf_mul(f1res,f1res,SQRTN);
mpf_clear(h);
return f1res;
}
#if 0
void newton()
{
mpf_t x0,x1;
mpf_init(x0);
mpf_init_set_d(x1,1.1);
do
{
mpf_set(x0,x1);
mpf_div(x1,f(x0),f1(x0)); mpf_sub(x1,x0,x1);
} while (mpf_eq(x1,x0,bits-64)==0);
cout << "newton: "; mpf_out_str(NULL,10,50,x1); cout << endl;
//mpf_sqrt(x0,x1); mpf_div(x0,SQRTN,x0);
//mpz_t p_bound; mpz_init(p_bound); mpz_set_f(p_bound,x0);
//cout << "N=p*q, p<" << p_bound << endl;
//int i; std::cin >> i;
mpf_clear(x1); mpf_clear(x0);
}
#endif
};
void phimat2()
{
cout << "Starting experimental Phimat factoring method." << endl;
cout << "Warning: This method runs very long except for special constructed numbers" << endl;
cout << "of the form N=p*q; p,q prime and q=p+x, x small compared to p." << endl;
mpz_init_set_ui(entry::Base,2);
mpz_t V,x;
mpz_init(x); mpz_sqrt(x,n); mpz_mul_ui(x,x,2);
unsigned int step=2;
mpz_t Delta,initial_Delta,stop_Delta;
mpz_init_set_ui(Delta,0);
mpz_init(stop_Delta);
// if this small block is omitted, then default ratio 1:1 is used.
// Otherwise calculate Delta based on the assumption, that
// n=p*q and p_ratio*p = q_ratio*q (approximately)
lambda_delta::Delta_by_ratio(Delta,n,1.0,1.0);
lambda_delta::Delta_by_ratio(stop_Delta,n,2.0,1.0);
mpz_sub_ui(Delta,Delta,mpz_remainder_ui(Delta,3*32));
mpz_init_set(initial_Delta,Delta);
{
unsigned int r = mpz_remainder_ui(n,32);
unsigned int disp=0; // displacement
if (r%4==1) { disp=2; step=4; }
else if (r%8==3) { disp=4; step=8; }
else if (r%8==7) { disp=0; step=8; }
else { cerr << "unknown state " << r << endl; exit(1); }
if (mpz_remainder_ui(n,6)==5) { step*=3; disp*=3; }
cout << "phimat multiplier is " << step << endl;
step/=2;
r=mpz_remainder_ui(x,step);
mpz_sub_ui(Delta,Delta,r); mpz_add_ui(Delta,Delta,disp);
mpz_sub_ui(x,x,r); mpz_add_ui(x,x,disp);
}
mpz_init_set(V,n); mpz_add_ui(V,V,1); mpz_sub(V,V,x);
mpz_powm(V,entry::Base,V,n);
mpz_powm_ui(entry::Base,entry::Base,step,n);
// remark:
// the "Delta" value for phimat should be exactly twice the value
// of the "Fermat"-Delta value for the same number.
// This is because "phimat" tries to determine p+q, whereas
// the "Fermat" procedure calculates implicitly (p+q)/2.
// Both procedures terminate, when "Delta" has reached this value.
// (Well, maybe phimat terminates earlier with a very low probability,
// because the "congruency" 2^(N+1) == 2^(p+q) (mod N) is not necessarily
// a (perfect) "identity" (and we do not know for certain, that p,q are the
// only prime factors of N).
phimahashvecs entries; // structure of indices (sorted by powers)
{
//const unsigned int tablesize = 0x400000;
const unsigned int t1 = entries.max_size();
unsigned int t2 = static_cast<size_t>((allowed_memory_usage_KB()/sizeof(entry)));
if (t2>65536) t2-=32768;
t2*=1024;
const unsigned int tablesize = MIN(t1,t2);
cout << "Creating table of " << tablesize << " items." << endl;
cout << "This may take a while..." << endl;
entries.prepare(tablesize);
cout << "Okay, items are prepared. Now searching for factors..." << endl;
}
{
// correction for initial Delta (in case it's not zero)
mpz_t x;
mpz_init(x);
mpz_div_ui(x,initial_Delta,step);
mpz_powm(x,entry::Base,x,n);
mpz_invert(x,x,n); // we walk downwards...
mpz_mod(x,x,n);
mpz_mul(V,V,x); mpz_mod(V,V,n);
mpz_clear(x);
}
mpz_t PowIntervalStep;
mpz_init_set(PowIntervalStep,entry::Base);
mpz_powm_ui(PowIntervalStep,PowIntervalStep,entries.Intervalsize(),n);
mpz_invert(PowIntervalStep,PowIntervalStep,n); // we walk downwards...
mpz_mod(PowIntervalStep,PowIntervalStep,n);
mpz_t TableStep;
mpz_init_set_ui(TableStep,entries.Intervalsize()); mpz_mul_ui(TableStep,TableStep,step);
for (unsigned int i=0; /* endlessly */ ; ++i)
{
unsigned int index;
if (entries.found(V,index))
{
mpz_t x,y,z,myDelta;
mpz_init(x); mpz_init(y); mpz_init(z); mpz_init_set(myDelta,Delta);
mpz_set_ui(x,index); mpz_mul_ui(x,x,step); mpz_add(myDelta,myDelta,x);
cout << "Found a candidate at Delta=" << myDelta << endl;
mpz_sqrt(x,n); mpz_mul_ui(x,x,2); mpz_add(x,x,myDelta);
mpz_mul(y,x,x); mpz_mul_ui(z,n,4); mpz_sub(y,y,z);
mpz_sqrt(y,y); mpz_add(y,y,x);
mpz_gcd(z,y,n); mpz_mod(z,z,n);
if (mpz_cmp_ui(z,1)>0)
{
cout << "looks good :-)" << endl;
mpz_set(Delta,myDelta); // accepted... :-)
while (mpz_cmp_ui(n,1)!=0)
{
if (mpz_probab_prime_p(z,probab_prime_checks))
{
cout << z << " is factor." << endl;
std::ostringstream comment;
comment << " [phimat]";
Factorization_to_file << MAL(z,comment) << flush;
}
else
{
cout << z << " is a composite factor." << endl;
if (Potenztest(z)) Factorization_to_file << " [phimat]" << flush;
else
{
std::ostringstream comment;
comment << " [phimat] [composite]";
Factorization_to_file << MAL(z,comment) << flush;
}
}
mpz_divexact(n,n,z);
mpz_set(z,n);
}
}
mpz_clear(myDelta); mpz_clear(z); mpz_clear(y); mpz_clear(x);
if (mpz_cmp_ui(n,1)==0) break;
}
mpz_add(Delta,Delta,TableStep);
mpz_mul(V,V,PowIntervalStep); mpz_mod(V,V,n);
if ((i&0x3ffff)==0)
{
if ((i&0x1ffffff)==0) // status report
{
mpz_t x,y,z;
mpz_init(x); mpz_init(y); mpz_init(z);
mpz_sqrt(x,n); mpz_mul_ui(x,x,2); mpz_add(x,x,Delta);
mpz_mul(y,x,x); mpz_mul_ui(z,n,4); mpz_sub(y,y,z);
mpz_sqrt(y,y); mpz_add(y,y,x);
mpz_div_ui(y,y,2);
mpz_div(x,n,y);
cout << "new bound p<" << x << endl;
mpz_clear(z); mpz_clear(y); mpz_clear(x);
//cout << "stop_Delta: " << stop_Delta << endl;
if (mpz_cmp(Delta,stop_Delta)>0) break;
}
cout << "Delta=" << Delta << '\r' << flush;
}
}
// for validation reasons give the Delta for Fermat,
// where the square could be found...
mpz_div_ui(Delta,Delta,2);
cout << "Fermat Delta would be " << Delta << endl;
mpz_clear(TableStep); mpz_clear(PowIntervalStep);
mpz_clear(initial_Delta); mpz_clear(stop_Delta); mpz_clear(Delta);
mpz_clear(x); mpz_clear(V);
mpz_clear(entry::Base);
}
void phimat(mpz_t n)
{
mpz_t x,y;
mpz_init(x); mpz_init(y);
mpz_mul_ui(x,n,4); mpz_sqrt(x,x); // x=trunc(2*sqrt(n));
mpz_add_ui(y,n,1); mpz_sub(y,y,x); // y=n+1-x
mpz_set_ui(x,2); mpz_powm(x,x,y,n); // x=2^y mod n;
cout << "trying phimat..." << x << endl;
//mpz_add_ui(y,n,1); mpz_div_ui(y,y,2);
unsigned long int z=0;
while (z<1000000000)
{
int p=mpz_scan1(x,0); z+=p;
mpz_tdiv_q_2exp(x,x,p);
//cout << z << ": " << x << endl;
if (mpz_cmp_ui(x,1)==0) break;
mpz_add(x,x,n);
}
cout << z << ":" << x << endl;
if (mpz_cmp_ui(x,1)==0)
{
cout << "x=1 ... " << endl;
mpz_mul_ui(x,n,4); mpz_sqrt(x,x); mpz_add_ui(x,x,z); // x=p+q
mpz_mul(y,x,x);
mpz_sub(y,y,n); mpz_sub(y,y,n); mpz_sub(y,y,n); mpz_sub(y,y,n);
mpz_sqrt(y,y); // y=p+q
mpz_add(y,x,y); mpz_div_ui(y,y,2); // y=p
mpz_div(x,n,y); // x=q
mpz_mul(x,x,y);
cout << x << endl;
if (mpz_cmp(x,n)==0) // wirklich echte Teiler ermittelt?
{
mpz_div(x,n,y); // x=q
if (mpz_probab_prime_p(x,probab_prime_checks))
{
cout << x << " is factor." << endl;
Factorization_to_file << MAL(x) << flush;
}
else
{
cout << x << " is a composite factor." << endl;
Factorization_to_file << MAL(x) << " [composite]" << flush;
}
if (mpz_probab_prime_p(y,probab_prime_checks))
{
cout << y << " is factor." << endl;
Factorization_to_file << MAL(y) << flush;
}
else
{
cout << y << " is a composite factor." << endl;
Factorization_to_file << MAL(y) << " [composite]" << flush;
}
mpz_set_ui(n,1);
}
}
mpz_clear(x); mpz_clear(y);
}