11476 - Factorizing Large Integers

[tanmoy]

it is getting TLE . is there any better idea to solve this problem ?
i wanna learn some better approach.
plz healp
[#include<stdio.h>
#include<math.h>
#include<vector>
using namespace std;
int r=0,e;
void prime_factors(long long a)
{
e=0;
int w=0;
int s=0;
if(a%2==0)
{
w=1;
s=1;
while(((a%2)==0))
{
e++;
a/=2;
}
if(e>1)printf(" 2^%d",e);
else printf(" 2");
e=0;
}
long long f,i;

f=(long)sqrt(a);
for(i=3;i<=f;i+=2)
{

if(a%i ==0)
{
s=1;
while(a%i ==0)
{
e++;
a/=i;
}
if(e>1 && w==1)printf(" * %lld^%d",i,e);
else if(e==1 && w==1) printf(" * %lld",i);
else if(e>1 && w==0) {printf(" %lld^%d",i,e);w=2;}
else if(e==1 && w==0) {printf(" %lld",i);w=2;}
else if(e>1 && w==2) printf(" * %lld^%d",i,e);
else if(e==1 && w==2) {printf(" * %lld",i);}


e=0;
}

}
if(a!=1 && s==1)
{
printf(" * %lld",a);
}
else if(a !=1 && s==0)printf(" %lld",a);
}

int main()
{


long long n;
int y;
int r=0;
scanf("%d",&y);
do
{
scanf("%lld",&n);
printf("%lld =",n);
prime_factors(n);
printf("\n");
r++;

}while(r!=y);

return 0;
}










][/code]

Edit : Title Modified by Moderator

[Leonid]

Well, your algorithm is O(T * sqrt(N)). Where T is a number of tests and N is a given number. Since T <= 800 and sqrt(N) <= 10^8. So the worst expected number of iterations is 8 * 10^2 * 10^8 = 8 * 10^10. That is obviously too much. We (Kaunas TU) haven't solved this task during the contest. But I suspect this task may deal with some more complex than the classical Sieve of Primes and naive factoring algorithm. You may want to study the algorithms from wiki: http://en.wikipedia.org/wiki/Prime_factorization
I would first consider this algorithm as an option: Pollard's p ? 1 algorithm. However, it would be very interesting to hear from someone who has actually solved this task.

[Robert Gerbicz]

I would first consider this algorithm as an option: Pollard's p ? 1 algorithm. However, it would be very interesting to hear from someone who has actually solved this task.

In my AC code I've used the so called Brent-Pollard-rho method, this is almost exactly the Pollard's p-1 method, but faster, this code used 1.87 sec., the time limit was 15 seconds, but note that my algorithm is highly optimized, currently this holds the record on all sites where I submitted, 500 lines long and about 14Kbytes. What I want to say that I don't think that by the simple Pollard's p-1 method you can avoid TLE, try to implement the Brent's version.

[Leonid]

What I want to say that I don't think that by the simple Pollard's p-1 method you can avoid TLE, try to implement the Brent's version.

Thanks for the idea, I'll definitely read more about that stuff once I have time

[sunny]

My code doesn't use Brent's improvement just normal Pollard Rho method.
Robert Gerbicz's solution is much faster than mine.

First I wanted to make the test data little more tight so only the Brent versions could pass, but later
stopped from doing that.

[hyperion]

I've used Fermat's Factorization in my code, but unfortunately I'm getting Wrong Answer. I've tested my program with big primes and semiprimes in range and it works well. It's a little confusing too, because my code consumes just about 30 miliseconds on UVa's judge (I'm not sure my code is really this fast).
Could anybody please give me some Special input numbers, so I can test and find out my problem ?
Thanks in Advance..

[baodog]

You can try your Code on the following similar factoring large integer problems on SPOJ, and see if your
runtime is also too small.

https://www.spoj.pl/problems/DIVSUM2/

If you want to check simple accuracy try:

https://www.spoj.pl/problems/DIVSUM

[Vytenis]

I have tried both Pollard and Brent versions for factorization and got TLE. However, the major bottleneck in my case is calculating (a*b)%c, where a, b and c are 64-bit integers, without an overflow. I am using O(log n) algorithm similar to fast exponentiation. Is there a way to speed up this computation?

[crackerwang]

My code doesn't use Brent's improvement just normal Pollard Rho method.
Robert Gerbicz's solution is much faster than mine.

First I wanted to make the test data little more tight so only the Brent versions could pass, but later
stopped from doing that.

I also use the normal Pollard Rho,but i got 20 TLE during the contest ,and till to now.can
you email you code to me?
or somebody else email me a AC code ,maybe i can find the someplace to optimize.
thx in advance
my email: wangyongliang.wyl@gmail.com

[friggstad]

I'm using Brent's improvement to Pollard's Rho method and I'm still getting TLE. Does anyone know of a good web resource that discusses practical considerations for implementing Brent's improvement? Also, what pseudo-random functions have people found successful?

[rio]

I have tried both Pollard and Brent versions for factorization and got TLE. However, the major bottleneck in my case is calculating (a*b)%c, where a, b and c are 64-bit integers, without an overflow. I am using O(log n) algorithm similar to fast exponentiation. Is there a way to speed up this computation?

I used Karatsuba multiplication. Since my code is not so fast, maybe there is a better way.

-----
Rio

[temper_3243]

hi ,
i'm trying to get accept using brents version, i get the factors but it gets TLE, if i increase the NOTRIES it becomes TLE, otherwise it doesn't find a factor , Is x^1 + C the optimal function ? Can someone give me their fast algorithm ? i want to see how to improve upon this. Please post some inputs and outputs for which it takes lots of time, semiprimes etc.

For isprime , i'm using miller rabin

Code: Select all

#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cassert>
#include<set>
#include<map>
#include<cstring>
#include<algorithm>
using namespace std;
#define mymin(a,b) ((a<b)?(a):(b))
#define MAXRETRIES (1<<16)
typedef unsigned long long ull;

unsigned int cntm = 0, pcnt;
unsigned long primes[MAXRETRIES];
unsigned long parr[MAXRETRIES];

ull arr[MAXRETRIES];
set < ull > pmap;

map < ull, ull > pm;
map < ull, ull >::iterator it;

#define BAD (0xffffffffffffffffULL)

int NOTRIES= (1<<15);
int FTRIES=(1<<10);

ull
compmod (ull a, ull b, ull n)
{
  ull r = 0;
  assert(n!=0);
  a %= n;
  b %= n;
  while (b)
    {
      if (b & 1)
	{
	  r += a;
	  if (r >= n)
	    r = r - n;
	}
      a <<= 1;
      if(a>=n)
	      a=a-n;
      b >>= 1;
    }
  return r;
}

 ull
comp (ull x, ull N, ull C)
{
  assert(N!=0);
  x %= N;
  ull r = compmod (x, x, N); 
  return (r + C) % N;
}

 ull
gcd (ull a, ull b)
{

  if (b == 0)
    return a;

  return gcd (b, a % b);
}

 ull
mypow (ull a, ull b, ull c)
{
  assert(c!=0);
  a %= c;
  ull z = 1, prod = a, u, x, y, prod1, mp;
  int sign = b & 1;
  if (b == 1)
    return a;

  if (sign)
    z = a;

  prod = compmod (a, a, c);
  u = mypow (prod, b / 2, c) % c;

  if (sign)
    u = compmod (z, u, c);
  return u;
} 
bool issquare(ull x)
{
	ull z=sqrt(x);
	if(z*z == x)
		return true;
	return false;
}
ull
brent (ull n, ull C, ull m)
{
  long long y, r, q, pp;
  {
    y = r = q = 1;

    long long x, ys, k, z, g;
    register int i, iter = 0;
    while (++iter < NOTRIES)
      {
	x = y;
	for (i = 1; ++iter < NOTRIES && i <= r; i++)
	  y = comp (y, n, C);

	k = 0;
	while (++iter < NOTRIES)
	  {
	    ys = y;
	    pp = mymin (m, r - k);
	    for (i = 1; ++iter < NOTRIES && i <= pp; i++)
	      {
		y = comp (y, n, C);
		z = abs (x - y);
		q = compmod (q, z, n);
	      }
	    g = gcd (q, n);
	    k += m;
	    if ((k >= r) || (g > 1))
	      break;
	  }
	r <<= 1;
	if (g > 1)
	  break;
      }
    iter = 0;
    if (g == n)
      {
	while (++iter < NOTRIES)
	  {
	    ys = comp (ys, n, C);
	    g = gcd (abs (x - ys), n);
	    if (g > 1)
	      break;
	  }
      }
    if (g > 1)
      return g;
  }

  return BAD;
}

 void
genprimes (void)
{
  int i, j;
  for (i = 3; i * i <= MAXRETRIES; i += 2)
    {
      if (parr[i] == 1)
	continue;
      for (j = i; i * j <= MAXRETRIES; j += 2)
	{
	  parr[i * j] = 1;
	}

    }
  primes[pcnt++] = 2;
  for (i = 3; i <= MAXRETRIES; i += 2)
    {
      if (!parr[i])
	primes[pcnt++] = i;
    }

  return;
}

 
int
isprime (ull n)
{
  ull x, s = 0, d, a, r, z, z1, r2, k, p;
  d = x = n - 1;

  if (n <=3)
    return true;

  if ((n & 1) == 0)
    return false;

  if (n <= MAXRETRIES)
       return ((parr[n] == 0) ? true : false);

  if(pmap.count(n))
	  return 1;

  while ((d &1) == 0)
    {
      d >>=1;
      s++;
    }

  int flag=1;
  
  for (k = 25; k < 250; k+=10 )
    {
      a = primes[k];
      assert(a!=0);
      z = mypow (a, d, n);
      flag=0;
      if(z==1)
	flag =1;
      for (r = 0;  r < s; r++)
	{
	  r2 = (1LL << r);
	  r2 *= d;
	  z1 = mypow (a, r2, n);
	  if(z1 == n-1)
	    flag |=1;
	}
      if(!flag)
	      break;
    }
  
  if(flag)
	  pmap.insert(n);

  return flag;
}

ull
factorizegen (ull n)
{
  int i;
  if(isprime(n))
	  return n;

  for (i = 0; (i < pcnt) && (primes[i] * primes[i]) <= n; i++)
    {
     assert(primes[i]!=0);
      while (n && ((n % primes[i]) == 0))
	{
	  arr[cntm++] = primes[i];
	  n /= primes[i];
	}
    }

  return n;
}
ull fermat(ull N)
{
	if(N==1 )
		return N;

	ull A=ceil(sqrt(N));
	ull Bsq = A*A -N,z;
	int iters=0;
	while(++iters<FTRIES && !issquare(Bsq))
	{
		++A;
		Bsq=A*A-N;
	}

	if(!issquare(Bsq))
		return N;
	
	z=A-sqrt(Bsq);

	if(N%z==0)
	{
		arr[cntm++]=N/z;
		N/=z;
		if(isprime(N)){
			arr[cntm++]=N;
			N=1;
		}
	}
	else 
	{
		z=A+sqrt(Bsq);
		if(N%z==0)
		{
			arr[cntm++]=N/z;
			N/=z;

		}
		if(isprime(N)){
			arr[cntm++]=N;
			N=1;
		}

	}

	return N;

	
}
 
int
main ()
{
  genprimes ();
 
   ull k, z, n, sv, a, b;
  register int C, m, i = 1;
  int no,pp=0;
  scanf (" %d", &no);
  while (no--)
    {
      scanf (" %llu", &n);
      n&=BAD;
      sv = n;
      cntm = 0;
      pp++;
      if(n==1)
      {
	      printf("%llu = %llu\n",n,n);
	      continue;
      }

      if (!isprime (n))
	{
	  n = factorizegen (n);

	  if (isprime (n))
	    {
	      arr[cntm++] = n;
	      n = 1;
	    }
          n=fermat(n);
	  for (C = 1, m = 111; C < 5 && n > (1LL<<32) ;++C, m += 29)
	    {
	      z = brent (n, C, m);
	      if (z <= 1 || z >= n)
		continue;

	      while (n > 1 && (n % z == 0))
		{
		  if (!isprime (z))
		    {
		      k = factorizegen (z);
		    }
		  else
		    k = z;
		  if (k > 1)
		    arr[cntm++] = k;
		  n /= z;
		}

	      if (isprime (n))
		{
		  arr[cntm++] = n;
		  n = 1;
		  break;
		}

	    }

	  if (n > 1)
	    {
	      if (!isprime (n))
		{
		  n = factorizegen (n);
		}
	      if (isprime (n))
		arr[cntm++] = n;
	    }
	}
      else
	arr[cntm++] = n;

      pm.clear ();
      for (i = 0; i < cntm; i++)
	{
	if(arr[i]>1)
	  pm[arr[i]]++;
	}

      i = 0;
      printf ("%llu = ", sv);
      for (it = pm.begin (); it != pm.end (); it++)
	{
	  a = it->first;
	  b = it->second;
	  if (i)
	    printf (" * ");

	  i = 1;

	  printf ("%llu", a);
	  if (b > 1)
	    printf ("^%llu", b);
	}
      printf ("\n");

    }

  return 0;
}

Can someone tell me if the outputs are correct

Code: Select all

70
77841374955853
3531696826961
100001299949
1000609937
8223372036854775841
10000000000000333
331226869951
3531696826961
100001299949
1000609937
8223372036854775841
10000000000000333
331226869951
9293124877756151
27039683174594548
3531696826961
9293124877756151
1809679
2569849
2942509
3115367
3293951
3469043
3808201
4566197
6594901
7262741
8426251
8875067
11569507
12431101
14973047
20760119
20840899
22589837
26290571
32623813
37908601
41165611
52294763
60095729
64940959
71726569
103816283
104754649
211832821
314459137
375957443
438424279
480358751
546013739
616821769
843552299
2153527171
9207818651
9535269647
40100030743
129669709111
38503
917154911
917084617
2083372023534909997
2176151660844999583
270396831745945481
2063512844981574047
3232323
727272727
919191919
765404567
987101789

Outputs are

Code: Select all

77841374955853 = 31153 * 35117 * 71153
3531696826961 = 1879781^2
100001299949 = 100003 * 999983
1000609937 = 10007 * 99991
8223372036854775841 = 61 * 2909 * 46342171761209
10000000000000333 = 2383 * 31727 * 132265613
331226869951 = 579947^2
3531696826961 = 1879781^2
100001299949 = 100003 * 999983
1000609937 = 10007 * 99991
8223372036854775841 = 61 * 2909 * 46342171761209
10000000000000333 = 2383 * 31727 * 132265613
331226869951 = 579947^2
9293124877756151 = 92931259 * 99999989
27039683174594548 = 2^2 * 67 * 100894340203711
3531696826961 = 1879781^2
9293124877756151 = 92931259 * 99999989
1809679 = 1217 * 1487
2569849 = 1283 * 2003
2942509 = 1087 * 2707
3115367 = 1579 * 1973
3293951 = 1453 * 2267
3469043 = 1049 * 3307
3808201 = 1249 * 3049
4566197 = 1283 * 3559
6594901 = 1483 * 4447
7262741 = 2027 * 3583
8426251 = 1489 * 5659
8875067 = 2179 * 4073
11569507 = 3229 * 3583
12431101 = 1873 * 6637
14973047 = 2713 * 5519
20760119 = 3169 * 6551
20840899 = 2579 * 8081
22589837 = 4483 * 5039
26290571 = 4001 * 6571
32623813 = 3121 * 10453
37908601 = 5171 * 7331
41165611 = 3719 * 11069
52294763 = 3463 * 15101
60095729 = 2141 * 28069
64940959 = 7607 * 8537
71726569 = 6131 * 11699
103816283 = 4591 * 22613
104754649 = 10211 * 10259
211832821 = 6101 * 34721
314459137 = 14107 * 22291
375957443 = 18301 * 20543
438424279 = 16421 * 26699
480358751 = 2927 * 164113
546013739 = 13879 * 39341
616821769 = 17761 * 34729
843552299 = 27043 * 31193
2153527171 = 34369 * 62659
9207818651 = 30881 * 298171
9535269647 = 61991 * 153817
40100030743 = 120623 * 332441
129669709111 = 323591 * 400721
38503 = 139 * 277
917154911 = 29333 * 31267
917084617 = 29327 * 31271
2083372023534909997 = 1039414451 * 2004370847
2176151660844999583 = 1224351209 * 1777391687
270396831745945481 = 292028911 * 925924871
2063512844981574047 = 1112041493 * 1855607779
3232323 = 3^2 * 359147
727272727 = 727272727
919191919 = 919191919
765404567 = 765404567
987101789 = 987101789

[temper_3243]

6609411 11476 Factorizing Larget Integers Accepted C++ 2.520 2008-08-24 17:28:24

Did anyone get accepted without caching primes ? Can you please send you code to temper3243@gmail.com

i used brent + fermat.

even with only brent i got accepted.

i got couple of TLE's because of not caching primes. Once you start caching primes > 10^6 in set and semiprimes > 10^9 in a map , you will be able to get accepted.

i think input set is not very hard.

i mean if we have 200 semiprimes (primes > 10^6) or 200 numbers with product of 3 primes with primes greater than > 10^5 ,

If these numbers don't share any common primes (Does anyone have this list )

mine would have been TLE in that case.

i think there are duplicates in the testcase and primes repeated for the numbers, exploit that.

Any suggestions to speed it up ?

[tryit1]

if you are using the function x^2 + 1 and only generating primes until 2^16 and using trial divison to factor numbers like 129669709111( 323591 * 400721) will fail.
if x^2 + 1 is the function you will get TLE for input

Code: Select all

335223032701393  
18308593#18309601=335223032701393
n is 1145285232319
C is 3
1145285232319 = 1070171 * 1070189
n is 1320695116189
C is 3
1320695116189 = 1149209 * 1149221
n is 1425263496311
C is 3
1425263496311 = 1193839 * 1193849


n is 1158781202053
C is 1
1158781202053 = 1076461 * 1076473
n is 1159534851473
C is 1
1159534851473 = 1076813 * 1076821

n is 1440360022379
C is 1
1440360022379 = 1200139 * 1200161
n is 1741988184311
C is 1
1741988184311 = 1319839 * 1319849
n is 2092521368009
C is 1
2092521368009 = 1446551 * 1446559
n is 2615571956513
C is 1
2615571956513 = 1617269 * 1617277
n is 2929150098551
C is 1
2929150098551 = 1711471 * 1711481
n is 3000825179377
C is 1
3000825179377 = 1732277 * 1732301
n is 3085366023437
C is 1
3085366023437 = 1756519 * 1756523
n is 3205467705563
C is 1
3205467705563 = 1790363 * 1790401
n is 3240432014399
C is 1
3240432014399 = 1800119 * 1800121
n is 3376277626081
C is 1
3376277626081 = 1837453 * 1837477
n is 3650671312913
C is 1
3650671312913 = 1910669 * 1910677

pollard rho will not give with x=y=2 and x^2+1 in atleast 10^6 iterations. So generate primes until 2^20 should work fine or use x^2 +3 as function.

in brents again x^2+1 may fail use x^2 +3 , and the value of m should be 1 , 39,67,etc i try until 130

[Sedefcho]

Without any caching (of primes or of divisors already found) I was able
to get ACC in 4.750 secs. With some caching my best ACC was in 1.200 secs.
So if Robert Gerbitz is not doing any caching (between the different test cases)
then for sure his solution is quite better.

This problem was a real monster for me, I tried hard to solve it at
least 4-5 times over the last 1 year and I had already lost hope

So I was really amazed I finally got ACC (and within such good time).

My first ACC program was with certain caching and with a run time of 7.790 secs. But after
some optimizations managed to reduce this (caching version) to 1.200 secs (as mentioned).

I used Brent's variant of Pollard rho but a couple of other optimizations too
in the cases (i.e. for numbers) where Brent's algorithm behaved badly/slowly.
I am really interested if Brent's algorithm behaved badly due to my (not perfect of course)
implementation or due to its natural inabilities to deal with certain types of numbers.

According to wikipedia, it seems to me that both versions (Pollard rho and Brent's variant)
are not yet fully evaluated/estimated (in terms of time complexity).

In general - really nice problem as it makes you learn quite some things
/ unless you already know them of course /