Skip to content

Instantly share code, notes, and snippets.

@simonlindholm
Created March 7, 2020 17:22
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save simonlindholm/193be06491e539c734c1ea78782c9125 to your computer and use it in GitHub Desktop.
Save simonlindholm/193be06491e539c734c1ea78782c9125 to your computer and use it in GitHub Desktop.
Pollard rho
namespace ttmath{typedef unsigned int U;typedef signed int sint;typedef uint64_t ulint;typedef int64_t slint;}namespace ttmath{enum LibTypeCode{asm_vc_32=0,asm_gcc_32,asm_vc_64,asm_gcc_64,no_asm_32,no_asm_64};enum ErrorCode{err_ok=0,err_nothing_has_read,err_unknown_character,err_unexpected_final_bracket,err_stack_not_clear,err_unknown_variable,err_division_by_zero,err_interrupt,err_overflow,err_unknown_function,err_unknown_operator,err_unexpected_semicolon_operator,err_improper_amount_of_arguments,err_improper_argument,err_unexpected_end,err_internal_error,err_incorrect_name,err_incorrect_value,err_variable_exists,err_variable_loop,err_functions_loop,err_must_be_only_one_value,err_object_exists,err_unknown_object,err_still_calculating,err_in_short_form_used_function,err_percent_from};struct Conv{U base;bool scient;sint scient_from;bool base_round;sint round;bool trim_zeroes;U comma;U comma2;U group;U group_digits;U group_exp;Conv(){base=10;scient=false;scient_from=15;base_round=true;round=-1;trim_zeroes=true;comma='.';comma2=',';group=0;group_digits=3;group_exp=0;}};class StopCalculating{public:virtual bool WasStopSignal()C volatile{R false;}virtual~StopCalculating(){}};class ExceptionInfo{C char*file;int line;public:ExceptionInfo():file(0),line(0){}ExceptionInfo(C char*f,int l):file(f),line(l){}string Where()C{if(!file)R "unknown";ostringstream result;result<<file<<":"<<line;R result.str();}};class ReferenceError:public logic_error,public ExceptionInfo{public:ReferenceError():logic_error("reference error"){}ReferenceError(C char*f,int l):logic_error("reference error"),ExceptionInfo(f,l){}string Where()C{R ExceptionInfo::Where();}};class RuntimeError:public runtime_error,public ExceptionInfo{public:RuntimeError():runtime_error("internal error"){}RuntimeError(C char*f,int l):runtime_error("internal error"),ExceptionInfo(f,l){}string Where()C{R ExceptionInfo::Where();}};}namespace ttmath{class Misc{public:static void AssignString(string&result,C char*str){result=str;}static void AssignString(wstring&result,C char*str){result.clear();for(;*str;++str)result+=*str;}static void AssignString(wstring&result,C string&str){R AssignString(result,str.c_str());}static void AssignString(string&result,C wchar_t*str){result.clear();for(;*str;++str)result+=static_cast<char>(*str);}static void AssignString(string&result,C wstring&str){R AssignString(result,str.c_str());}static void AddString(string&result,C char*str){result+=str;}static void AddString(wstring&result,C char*str){for(;*str;++str)result+=*str;}template<class char_type>static void SkipWhiteCharacters(C char_type*&c){while((*c==' ')||(*c=='\t')||(*c==13)||(*c=='\n'))++c;}static U CharToDigit(U c){if(c>='0'&&c<='9')R c-'0';if(c>='a'&&c<='z')R c-'a'+10;R c-'A'+10;}static sint CharToDigit(U c,U base){if(c>='0'&&c<='9')c=c-'0';else if(c>='a'&&c<='z')c=c-'a'+10;else if(c>='A'&&c<='Z')c=c-'A'+10;else R -1;if(c>=base)R -1;R sint(c);}static U DigitToChar(U digit){if(digit<10)R digit+'0';R digit-10+'A';}};}namespace ttmath{template<U S>class V{public:U T[S];template<class ostream_type>void PrintTable(ostream_type&output)C{C int columns=8;int c=1;for(int i=S-1;i>=0;--i){output<<"0x"<<setfill('0');output<<setw(8);output<<hex<<T[i];if(i>0){output<<",";if(++c>columns){output<<endl;c=1;}}}output<<dec<<endl;}template<class char_type,class ostream_type>static void PrintVectorLog(C char_type*msg,ostream_type&output,C U*vector,U vector_len){output<<msg<<endl;for(U i=0;i<vector_len;++i)output<<"T["<<i<<"]:"<<vector[i]<<endl;}template<class char_type,class ostream_type>static void PrintVectorLog(C char_type*msg,U carry,ostream_type&output,C U*vector,U vector_len){PrintVectorLog(msg,output,vector,vector_len);output<<"carry:"<<carry<<endl;}template<class char_type,class ostream_type>void PrintLog(C char_type*msg,ostream_type&output)C{PrintVectorLog(msg,output,T,S);}template<class char_type,class ostream_type>void PrintLog(C char_type*msg,U carry,ostream_type&output)C{PrintVectorLog(msg,output,T,S);output<<"carry:"<<carry<<endl;}U Size()C{R S;}void SetZero(){for(U i=0;i<S;++i)T[i]=0;}void SetOne(){SetZero();T[0]=1;}void SetMax(){for(U i=0;i<S;++i)T[i]=4294967295u;}void SetMin(){SetZero();}void Swap(V<S>&ss2){for(U i=0;i<S;++i){U temp=T[i];T[i]=ss2.T[i];ss2.T[i]=temp;}}void SetFromTable(C U*temp_table,U temp_table_len){U temp_table_index=0;sint i;for(i=S-1;i>=0&&temp_table_index<temp_table_len;--i,++temp_table_index)T[i]=temp_table[temp_table_index];if(temp_table_index<temp_table_len){if((temp_table[temp_table_index]&2147483648u)!=0){if(T[0]!=4294967295u)++T[0];}}for(;i>=0;--i)T[i]=0;}U AddOne(){R AddInt(1);}U SubOne(){R SubInt(1);}private:void RclMoveAllWords(U&rest_bits,U&last_c,U bits,U c){rest_bits=bits %32u;U all_words=bits/32u;U mask=(c)?4294967295u:0;if(all_words>=S){if(all_words==S&&rest_bits==0)last_c=T[0]&1;for(U i=0;i<S;++i)T[i]=mask;rest_bits=0;}else if(all_words>0){sint first,second;last_c=T[S-all_words]&1;for(first=S-1,second=first-all_words;second>=0;--first,--second)T[first]=T[second];for(;first>=0;--first)T[first]=mask;}}public:U Rcl(U bits,U c=0){U last_c=0;U rest_bits=bits;if(bits==0)R 0;if(bits>=32u)RclMoveAllWords(rest_bits,last_c,bits,c);if(rest_bits==0){R last_c;}if(rest_bits==1){last_c=Rcl2_one(c);}else if(rest_bits==2){Rcl2_one(c);last_c=Rcl2_one(c);}else{last_c=Rcl2(rest_bits,c);}R last_c;}private:void RcrMoveAllWords(U&rest_bits,U&last_c,U bits,U c){rest_bits=bits %32u;U all_words=bits/32u;U mask=(c)?4294967295u:0;if(all_words>=S){if(all_words==S&&rest_bits==0)last_c=(T[S-1]&2147483648u)?1:0;for(U i=0;i<S;++i)T[i]=mask;rest_bits=0;}else if(all_words>0){U first,second;last_c=(T[all_words-1]&2147483648u)?1:0;for(first=0,second=all_words;second<S;++first,++second)T[first]=T[second];for(;first<S;++first)T[first]=mask;}}public:U Rcr(U bits,U c=0){U last_c=0;U rest_bits=bits;if(bits==0)R 0;if(bits>=32u)RcrMoveAllWords(rest_bits,last_c,bits,c);if(rest_bits==0){R last_c;}if(rest_bits==1){last_c=Rcr2_one(c);}else if(rest_bits==2){Rcr2_one(c);last_c=Rcr2_one(c);}else{last_c=Rcr2(rest_bits,c);}R last_c;}U CompensationToLeft(){U moving=0;sint a;for(a=S-1;a>=0&&T[a]==0;--a);if(a<0)R moving;if(a!=S-1){moving+=(S-1-a)*32u;sint i;for(i=S-1;a>=0;--i,--a)T[i]=T[a];for(;i>=0;--i)T[i]=0;}U moving2=FindLeadingBitInWord(T[S-1]);moving2=32u-moving2-1;Rcl(moving2);R moving+moving2;}bool FindLeadingBit(U&table_id,U&index)C{for(table_id=S-1;table_id!=0&&T[table_id]==0;--table_id);if(table_id==0&&T[table_id]==0){index=0;R false;}index=FindLeadingBitInWord(T[table_id]);R true;}bool FindLowestBit(U&table_id,U&index)C{for(table_id=0;table_id<S&&T[table_id]==0;++table_id);if(table_id>=S){index=0;table_id=0;R false;}index=FindLowestBitInWord(T[table_id]);R true;}U GetBit(U bit_index)C{U index=bit_index/32u;U bit=bit_index %32u;U temp=T[index];U res=SetBitInWord(temp,bit);R res;}U SetBit(U bit_index){U index=bit_index/32u;U bit=bit_index %32u;U res=SetBitInWord(T[index],bit);R res;}void BitAnd(C V<S>&ss2){for(U x=0;x<S;++x)T[x]&=ss2.T[x];}void BitOr(C V<S>&ss2){for(U x=0;x<S;++x)T[x]|=ss2.T[x];}void BitXor(C V<S>&ss2){for(U x=0;x<S;++x)T[x]^=ss2.T[x];}void BitNot(){for(U x=0;x<S;++x)T[x]=~T[x];}void BitNot2(){U table_id,index;if(FindLeadingBit(table_id,index)){for(U x=0;x<table_id;++x)T[x]=~T[x];U mask=4294967295u;U shift=32u-index-1;if(shift)mask>>=shift;T[table_id]^=mask;}else T[0]=1;}public:U MulInt(U ss2){U r1,r2,x1;U c=0;V<S>u(*this);SetZero();if(ss2==0){R 0;}for(x1=0;x1<S-1;++x1){MulTwoWords(u.T[x1],ss2,&r2,&r1);c+=AddTwoInts(r2,r1,x1);}MulTwoWords(u.T[x1],ss2,&r2,&r1);c+=(r2!=0)?1:0;c+=AddInt(r1,x1);R (c==0)?0:1;}template<U result_size>void MulInt(U ss2,V<result_size>&result)C{U r2,r1;U x1size=S;U x1start=0;result.SetZero();if(ss2==0){R ;}if(S>2){for(x1size=S;x1size>0&&T[x1size-1]==0;--x1size);if(x1size==0){R ;}for(x1start=0;x1start<x1size&&T[x1start]==0;++x1start);}for(U x1=x1start;x1<x1size;++x1){MulTwoWords(T[x1],ss2,&r2,&r1);result.AddTwoInts(r2,r1,x1);}R ;}U Mul(C V<S>&ss2,U algorithm=100){switch(algorithm){case 1:R Mul1(ss2);case 2:R Mul2(ss2);case 3:R Mul3(ss2);case 100:default:R MulFastest(ss2);}}void MulBig(C V<S>&ss2,V<S*2>&result,U algorithm=100){switch(algorithm){case 1:R Mul1Big(ss2,result);case 2:R Mul2Big(ss2,result);case 3:R Mul3Big(ss2,result);case 100:default:R MulFastestBig(ss2,result);}}private:U Mul1Ref(C V<S>&ss2){V<S>ss1(*this);SetZero();for(U i=0;i<S*32u;++i){if(Add(*this)){R 1;}if(ss1.Rcl(1))if(Add(ss2)){R 1;}}R 0;}public:U Mul1(C V<S>&ss2){if(this==&ss2){V<S>copy_ss2(ss2);R Mul1Ref(copy_ss2);}else{R Mul1Ref(ss2);}}void Mul1Big(C V<S>&ss2_,V<S*2>&result){V<S*2>ss2;U i;for(i=0;i<S;++i){result.T[i]=T[i];ss2.T[i]=ss2_.T[i];}for(;i<S*2;++i){result.T[i]=0;ss2.T[i]=0;}result.Mul1(ss2);}U Mul2(C V<S>&ss2){V<S*2>result;U i,c=0;Mul2Big(ss2,result);for(i=0;i<S;++i)T[i]=result.T[i];for(;i<S*2;++i)if(result.T[i]!=0){c=1;break;}R c;}void Mul2Big(C V<S>&ss2,V<S*2>&result){Mul2Big2<S>(T,ss2.T,result);}private:template<U ss_size>void Mul2Big2(C U*ss1,C U*ss2,V<ss_size*2>&result){U x1size=ss_size,x2size=ss_size;U x1start=0,x2start=0;if(ss_size>2){for(x1size=ss_size;x1size>0&&ss1[x1size-1]==0;--x1size);for(x2size=ss_size;x2size>0&&ss2[x2size-1]==0;--x2size);for(x1start=0;x1start<x1size&&ss1[x1start]==0;++x1start);for(x2start=0;x2start<x2size&&ss2[x2start]==0;++x2start);}Mul2Big3<ss_size>(ss1,ss2,result,x1start,x1size,x2start,x2size);}template<U ss_size>void Mul2Big3(C U*ss1,C U*ss2,V<ss_size*2>&result,U x1start,U x1size,U x2start,U x2size){U r2,r1;result.SetZero();if(x1size==0||x2size==0)R ;for(U x1=x1start;x1<x1size;++x1){for(U x2=x2start;x2<x2size;++x2){MulTwoWords(ss1[x1],ss2[x2],&r2,&r1);result.AddTwoInts(r2,r1,x2+x1);}}}public:U Mul3(C V<S>&ss2){V<S*2>result;U i,c=0;Mul3Big(ss2,result);for(i=0;i<S;++i)T[i]=result.T[i];for(;i<S*2;++i)if(result.T[i]!=0){c=1;break;}R c;}void Mul3Big(C V<S>&ss2,V<S*2>&result){Mul3Big2<S>(T,ss2.T,result.T);}private:template<U ss_size>void Mul3Big2(C U*ss1,C U*ss2,U*result){C U*x1,*x0,*y1,*y0;if(ss_size>1&&ss_size<20){V<ss_size*2>res;Mul2Big2<ss_size>(ss1,ss2,res);for(U i=0;i<ss_size*2;++i)result[i]=res.T[i];R ;}else if(ss_size==1){R MulTwoWords(*ss1,*ss2,&result[1],&result[0]);}if((ss_size&1)==1){x0=ss1;y0=ss2;x1=ss1+ss_size/2+1;y1=ss2+ss_size/2+1;Mul3Big3<ss_size/2+1,ss_size/2,ss_size*2>(x1,x0,y1,y0,result);}else{x0=ss1;y0=ss2;x1=ss1+ss_size/2;y1=ss2+ss_size/2;Mul3Big3<ss_size/2,ss_size/2,ss_size*2>(x1,x0,y1,y0,result);}}template<U first_size,U second_size,U result_size>void Mul3Big3(C U*x1,C U*x0,C U*y1,C U*y0,U*result){U i,c,xc,yc;V<first_size>temp,temp2;V<first_size*3>z1;Mul3Big2<first_size>(x0,y0,result);Mul3Big2<second_size>(x1,y1,result+first_size*2);xc=AddVector(x0,x1,first_size,second_size,temp.T);yc=AddVector(y0,y1,first_size,second_size,temp2.T);Mul3Big2<first_size>(temp.T,temp2.T,z1.T);for(i=first_size*2;i<first_size*3;++i)z1.T[i]=0;if(xc){c=AddVector(z1.T+first_size,temp2.T,first_size*3-first_size,first_size,z1.T+first_size);}if(yc){c=AddVector(z1.T+first_size,temp.T,first_size*3-first_size,first_size,z1.T+first_size);}if(xc&&yc){for(i=first_size*2;i<first_size*3;++i)if(++z1.T[i]!=0)break;}c=SubVector(z1.T,result+first_size*2,first_size*3,second_size*2,z1.T);c=SubVector(z1.T,result,first_size*3,first_size*2,z1.T);if(first_size>second_size){U z1_size=result_size-first_size;for(i=z1_size;i<first_size*3;++i){}c=AddVector(result+first_size,z1.T,result_size-first_size,z1_size,result+first_size);}else{c=AddVector(result+first_size,z1.T,result_size-first_size,first_size*3,result+first_size);}}public:U MulFastest(C V<S>&ss2){V<S*2>result;U i,c=0;MulFastestBig(ss2,result);for(i=0;i<S;++i)T[i]=result.T[i];for(;i<S*2;++i)if(result.T[i]!=0){c=1;break;}R c;}void MulFastestBig(C V<S>&ss2,V<S*2>&result){if(S<20)R Mul2Big(ss2,result);U x1size=S,x2size=S;U x1start=0,x2start=0;for(x1size=S;x1size>0&&T[x1size-1]==0;--x1size);for(x2size=S;x2size>0&&ss2.T[x2size-1]==0;--x2size);if(x1size==0||x2size==0){result.SetZero();R ;}for(x1start=0;x1start<x1size&&T[x1start]==0;++x1start);for(x2start=0;x2start<x2size&&ss2.T[x2start]==0;++x2start);U distancex1=x1size-x1start;U distancex2=x2size-x2start;if(distancex1<3||distancex2<3)R Mul2Big3<S>(T,ss2.T,result,x1start,x1size,x2start,x2size);Mul3Big(ss2,result);}public:U DivInt(U divisor,U*remainder=0){if(divisor==0){if(remainder)*remainder=0;R 1;}if(divisor==1){if(remainder)*remainder=0;R 0;}V<S>dividend(*this);SetZero();sint i;U r=0;for(i=S-1;i>0&&dividend.T[i]==0;--i);for(;i>=0;--i)DivTwoWords(r,dividend.T[i],divisor,&T[i],&r);if(remainder)*remainder=r;R 0;}U DivInt(U divisor,U&remainder){R DivInt(divisor,&remainder);}U Div(C V<S>&divisor,V<S>*remainder=0,U algorithm=3){switch(algorithm){case 1:R Div1(divisor,remainder);case 2:R Div2(divisor,remainder);case 3:default:R Div3(divisor,remainder);}}U Div(C V<S>&divisor,V<S>&remainder,U algorithm=3){R Div(divisor,&remainder,algorithm);}private:U Div_StandardTest(C V<S>&v,U&m,U&n,V<S>*remainder=0){switch(Div_CalculatingSize(v,m,n)){case 4:if(remainder)remainder->SetZero();SetOne();R 0;case 3:if(remainder)*remainder=*this;SetZero();R 0;case 2:if(remainder)remainder->SetZero();SetZero();R 0;case 1:R 1;}R 2;}U Div_CalculatingSize(C V<S>&v,U&m,U&n){m=n=S-1;for(;n!=0&&v.T[n]==0;--n);if(n==0&&v.T[n]==0)R 1;for(;m!=0&&T[m]==0;--m);if(m==0&&T[m]==0)R 2;if(m<n)R 3;else if(m==n){U i;for(i=n;i!=0&&T[i]==v.T[i];--i);if(T[i]<v.T[i])R 3;else if (T[i]==v.T[i])R 4;}R 0;}public:U Div1(C V<S>&divisor,V<S>*remainder=0){U m,n,test;test=Div_StandardTest(divisor,m,n,remainder);if(test<2)R test;if(!remainder){V<S>rem;R Div1_Calculate(divisor,rem);}R Div1_Calculate(divisor,*remainder);}U Div1(C V<S>&divisor,V<S>&remainder){R Div1(divisor,&remainder);}private:U Div1_Calculate(C V<S>&divisor,V<S>&rest){if(this==&divisor){V<S>divisor_copy(divisor);R Div1_CalculateRef(divisor_copy,rest);}else{R Div1_CalculateRef(divisor,rest);}}U Div1_CalculateRef(C V<S>&divisor,V<S>&rest){sint loop;sint c;rest.SetZero();loop=S*32u;c=0;div_a:c=Rcl(1,c);c=rest.Add(rest,c);c=rest.Sub(divisor,c);c=!c;if(!c)goto div_d;div_b:--loop;if(loop)goto div_a;c=Rcl(1,c);R 0;div_c:c=Rcl(1,c);c=rest.Add(rest,c);c=rest.Add(divisor);if(c)goto div_b;div_d:--loop;if(loop)goto div_c;c=Rcl(1,c);c=rest.Add(divisor);R 0;}public:U Div2(C V<S>&divisor,V<S>*remainder=0){if(this==&divisor){V<S>divisor_copy(divisor);R Div2Ref(divisor_copy,remainder);}else{R Div2Ref(divisor,remainder);}}U Div2(C V<S>&divisor,V<S>&remainder){R Div2(divisor,&remainder);}private:U Div2Ref(C V<S>&divisor,V<S>*remainder=0){U bits_diff;U status=Div2_Calculate(divisor,remainder,bits_diff);if(status<2)R status;if(CmpBiggerEqual(divisor)){Div2(divisor,remainder);SetBit(bits_diff);}else{if(remainder)*remainder=*this;SetZero();SetBit(bits_diff);}R 0;}U Div2_Calculate(C V<S>&divisor,V<S>*remainder,U&bits_diff){U table_id,index;U divisor_table_id,divisor_index;U status=Div2_FindLeadingBitsAndCheck(divisor,remainder,table_id,index,divisor_table_id,divisor_index);if(status<2){R status;}bits_diff=index-divisor_index;V<S>divisor_copy(divisor);divisor_copy.Rcl(bits_diff,0);if(CmpSmaller(divisor_copy,table_id)){divisor_copy.Rcr(1);--bits_diff;}Sub(divisor_copy,0);R 2;}U Div2_FindLeadingBitsAndCheck(C V<S>&divisor,V<S>*remainder,U&table_id,U&index,U&divisor_table_id,U&divisor_index){if(!divisor.FindLeadingBit(divisor_table_id,divisor_index)){R 1;}if(!FindLeadingBit(table_id,index)){SetZero();if(remainder)remainder->SetZero();R 0;}divisor_index+=divisor_table_id*32u;index+=table_id*32u;if(divisor_table_id==0){U r;DivInt(divisor.T[0],&r);if(remainder){remainder->SetZero();remainder->T[0]=r;}R 0;}if(Div2_DivisorGreaterOrEqual(divisor,remainder,table_id,index,divisor_index)){R 0;}R 2;}bool Div2_DivisorGreaterOrEqual(C V<S>&divisor,V<S>*remainder,U table_id,U index,U divisor_index){if(divisor_index>index){if(remainder)*remainder=*this;SetZero();R true;}if(divisor_index==index){U i;for(i=table_id;i!=0&&T[i]==divisor.T[i];--i);if(T[i]<divisor.T[i]){if(remainder)*remainder=*this;SetZero();R true;}else if(T[i]==divisor.T[i]){if(remainder)remainder->SetZero();SetOne();R true;}}R false;}public:U Div3(C V<S>&ss2,V<S>*remainder=0){if(this==&ss2){V<S>copy_ss2(ss2);R Div3Ref(copy_ss2,remainder);}else{R Div3Ref(ss2,remainder);}}U Div3(C V<S>&ss2,V<S>&remainder){R Div3(ss2,&remainder);}private:U Div3Ref(C V<S>&v,V<S>*remainder=0){U m,n,test;test=Div_StandardTest(v,m,n,remainder);if(test<2)R test;if(n==0){U r;DivInt(v.T[0],&r);if(remainder){remainder->SetZero();remainder->T[0]=r;}R 0;}++m;++n;m=m-n;Div3_Division(v,remainder,m,n);R 0;}private:void Div3_Division(V<S>v,V<S>*remainder,U m,U n){V<S+1>uu,vv;V<S>q;U d,u_value_size,u0,u1,u2,v1,v0,j=m;u_value_size=Div3_Normalize(v,n,d);if(j+n==S)u2=u_value_size;else u2=T[j+n];Div3_MakeBiggerV(v,vv);for(U i=j+1;i<S;++i)q.T[i]=0;while(true){u1=T[j+n-1];u0=T[j+n-2];v1=v.T[n-1];v0=v.T[n-2];U qp=Div3_Calculate(u2,u1,u0,v1,v0);Div3_MakeNewU(uu,j,n,u2);Div3_MultiplySubtract(uu,vv,qp);Div3_CopyNewU(uu,j,n);q.T[j]=qp;if(j--==0)break;u2=T[j+n];}if(remainder)Div3_Unnormalize(remainder,n,d);*this=q;}void Div3_MakeNewU(V<S+1>&uu,U j,U n,U u_max){U i;for(i=0;i<n;++i,++j)uu.T[i]=T[j];uu.T[i]=u_max;for(++i;i<S+1;++i)uu.T[i]=0;}void Div3_CopyNewU(C V<S+1>&uu,U j,U n){U i;for(i=0;i<n;++i)T[i+j]=uu.T[i];if(i+j<S)T[i+j]=uu.T[i];}void Div3_MakeBiggerV(C V<S>&v,V<S+1>&vv){for(U i=0;i<S;++i)vv.T[i]=v.T[i];vv.T[S]=0;}U Div3_Normalize(V<S>&v,U n,U&d){U bit=(U)FindLeadingBitInWord(v.T[n-1]);U move=(32u-bit-1);U res=T[S-1];d=move;if(move>0){v.Rcl(move,0);Rcl(move,0);res=res>>(bit+1);}else{res=0;}R res;}void Div3_Unnormalize(V<S>*remainder,U n,U d){for(U i=n;i<S;++i)T[i]=0;Rcr(d,0);*remainder=*this;}U Div3_Calculate(U u2,U u1,U u0,U v1,U v0){V<2>u_temp;U rp;bool next_test;u_temp.T[1]=u2;u_temp.T[0]=u1;u_temp.DivInt(v1,&rp);do{bool decrease=false;if(u_temp.T[1]==1)decrease=true;else{V<2>temp1,temp2;V<2>::MulTwoWords(u_temp.T[0],v0,temp1.T+1,temp1.T);temp2.T[1]=rp;temp2.T[0]=u0;if(temp1>temp2)decrease=true;}next_test=false;if(decrease){u_temp.SubOne();rp+=v1;if(rp>=v1)next_test=true;}}while(next_test);R u_temp.T[0];}void Div3_MultiplySubtract(V<S+1>&uu,C V<S+1>&vv,U&qp){V<S+1>vv_temp(vv);vv_temp.MulInt(qp);if(uu.Sub(vv_temp)){--qp;uu.Add(vv);}}public:U Pow(V<S>pow){if(pow.IsZero()&&IsZero())R 2;V<S>start(*this);V<S>result;result.SetOne();U c=0;while(!c){if(pow.T[0]&1)c+=result.Mul(start);pow.Rcr2_one(0);if(pow.IsZero())break;c+=start.Mul(start);}*this=result;R (c==0)?0:1;}void Sqrt(){V<S>bit,temp;if(IsZero())R ;V<S>value(*this);SetZero();bit.SetZero();bit.T[S-1]=(2147483648u>>1);while(bit>value)bit.Rcr(2);while(!bit.IsZero()){temp=*this;temp.Add(bit);if(value>=temp){value.Sub(temp);Rcr(1);Add(bit);}else{Rcr(1);}bit.Rcr(2);}}void ClearFirstBits(U n){if(n>=S*32u){SetZero();R ;}U*p=T;while(n>=32u){*p++=0;n-=32u;}if(n==0){R ;}U mask=4294967295u;mask=mask<<n;(*p)&=mask;}bool IsTheHighestBitSet()C{R (T[S-1]&2147483648u)!=0;}bool IsTheLowestBitSet()C{R (*T&1)!=0;}bool IsOnlyTheHighestBitSet()C{for(U i=0;i<S-1;++i)if(T[i]!=0)R false;if(T[S-1]!=2147483648u)R false;R true;}bool IsOnlyTheLowestBitSet()C{if(T[0]!=1)R false;for(U i=1;i<S;++i)if(T[i]!=0)R false;R true;}bool IsZero()C{for(U i=0;i<S;++i)if(T[i]!=0)R false;R true;}bool AreFirstBitsZero(U bits)C{U index=bits/32u;U rest=bits %32u;U i;for(i=0;i<index;++i)if(T[i]!=0)R false;if(rest==0)R true;U mask=4294967295u>>(32u-rest);R (T[i]&mask)==0;}template<U argument_size>U FromUInt(C V<argument_size>&p){U min_size=(S<argument_size)?S:argument_size;U i;for(i=0;i<min_size;++i)T[i]=p.T[i];if(S>argument_size){for(;i<S;++i)T[i]=0;}else{for(;i<argument_size;++i)if(p.T[i]!=0){R 1;}}R 0;}template<U argument_size>U FromInt(C V<argument_size>&p){R FromUInt(p);}U FromUInt(U value){for(U i=1;i<S;++i)T[i]=0;T[0]=value;R 0;}U FromInt(U value){R FromUInt(value);}U FromInt(sint value){U c=FromUInt(U(value));if(c||value<0)R 1;R 0;}template<U argument_size>V<S>&O=(C V<argument_size>&p){FromUInt(p);R *this;}V<S>&O=(C V<S>&p){for(U i=0;i<S;++i)T[i]=p.T[i];R *this;}V<S>&O=(U i){FromUInt(i);R *this;}V(U i){FromUInt(i);}V<S>&O=(sint i){FromInt(i);R *this;}V(sint i){FromInt(i);}U FromUInt(ulint n){T[0]=(U)n;if(S==1){U c=((n>>32u)==0)?0:1;R c;}T[1]=(U)(n>>32u);for(U i=2;i<S;++i)T[i]=0;R 0;}U FromInt(ulint n){R FromUInt(n);}U FromInt(slint n){U c=FromUInt(ulint(n));if(c||n<0)R 1;R 0;}V<S>&O=(ulint n){FromUInt(n);R *this;}V(ulint n){FromUInt(n);}V<S>&O=(slint n){FromInt(n);R *this;}V(slint n){FromInt(n);}V(C char*s){FromString(s);}V(C string&s){FromString(s.c_str());}V(C wchar_t*s){FromString(s);}V(C wstring&s){FromString(s.c_str());}V(){}V(C V<S>&u){for(U i=0;i<S;++i)T[i]=u.T[i];}template<U argument_size>V(C V<argument_size>&u){FromUInt(u);}~V(){}U ToUInt()C{R T[0];}U ToUInt(U&result)C{result=T[0];for(U i=1;i<S;++i)if(T[i]!=0)R 1;R 0;}U ToInt(U&result)C{R ToUInt(result);}U ToInt(sint&result)C{result=sint(T[0]);if((result&2147483648u)!=0)R 1;for(U i=1;i<S;++i)if(T[i]!=0)R 1;R 0;}U ToUInt(ulint&result)C{if(S==1){result=T[0];}else{U low=T[0];U high=T[1];result=low;result|=(ulint(high)<<32u);for(U i=2;i<S;++i)if(T[i]!=0)R 1;}R 0;}U ToInt(ulint&result)C{R ToUInt(result);}U ToInt(slint&result)C{ulint temp;U c=ToUInt(temp);result=slint(temp);if(c||result<0)R 1;R 0;}protected:double ToStringLog2(U x)C{static double log_tab[]={1.000000000000000000,0.630929753571457437,0.500000000000000000,0.430676558073393050,0.386852807234541586,0.356207187108022176,0.333333333333333333,0.315464876785728718,0.301029995663981195,0.289064826317887859,0.278942945651129843,0.270238154427319741,0.262649535037193547,0.255958024809815489,0.250000000000000000};if(x<2||x>16)R 0;R log_tab[x-2];}public:template<class string_type>void ToStringBase(string_type&result,U b=10,bool negative=false)C{V<S>temp(*this);U rest,table_id,index,digits;double digits_d;char character;result.clear();if(b<2||b>16)R ;if(!FindLeadingBit(table_id,index)){result='0';R ;}if(negative)result='-';digits_d=table_id;digits_d*=32u;digits_d+=index+1;digits_d*=ToStringLog2(b);digits=static_cast<U>(digits_d)+3;if(result.capacity()<digits)result.reserve(digits);do{temp.DivInt(b,&rest);character=static_cast<char>(Misc::DigitToChar(rest));result.insert(result.end(),character);}while(!temp.IsZero());size_t i1=negative?1:0;size_t i2=result.size()-1;for(;i1<i2;++i1,--i2){char tempc=static_cast<char>(result[i1]);result[i1]=result[i2];result[i2]=tempc;}}void ToString(string&result,U b=10)C{R ToStringBase(result,b);}string ToString(U b=10)C{string result;ToStringBase(result,b);R result;}void ToString(wstring&result,U b=10)C{R ToStringBase(result,b);}wstring ToWString(U b=10)C{wstring result;ToStringBase(result,b);R result;}private:template<class char_type>U FromStringBase(C char_type*s,U b=10,C char_type**after_source=0,bool*value_read=0){V<S>base(b);V<S>temp;sint z;U c=0;SetZero();temp.SetZero();Misc::SkipWhiteCharacters(s);if(after_source)*after_source=s;if(value_read)*value_read=false;if(b<2||b>16)R 1;for(;(z=Misc::CharToDigit(*s,b))!=-1;++s){if(value_read)*value_read=true;if(c==0){temp.T[0]=z;c+=Mul(base);c+=Add(temp);}}if(after_source)*after_source=s;R (c==0)?0:1;}public:U FromString(C char*s,U b=10,C char**after_source=0,bool*value_read=0){R FromStringBase(s,b,after_source,value_read);}U FromString(C string&s,U b=10){R FromString(s.c_str(),b);}V<S>&O=(C char*s){FromString(s);R *this;}V<S>&O=(C string&s){FromString(s.c_str());R *this;}U FromString(C wchar_t*s,U b=10,C wchar_t**after_source=0,bool*value_read=0){R FromStringBase(s,b,after_source,value_read);}U FromString(C wstring&s,U b=10){R FromString(s.c_str(),b);}V<S>&O=(C wchar_t*s){FromString(s);R *this;}V<S>&O=(C wstring&s){FromString(s.c_str());R *this;}bool CmpSmaller(C V<S>&l,sint index=-1)C{sint i;if(index==-1||index>=sint(S))i=S-1;else i=index;for(;i>=0;--i){if(T[i]!=l.T[i])R T[i]<l.T[i];}R false;}bool CmpBigger(C V<S>&l,sint index=-1)C{sint i;if(index==-1||index>=sint(S))i=S-1;else i=index;for(;i>=0;--i){if(T[i]!=l.T[i])R T[i]>l.T[i];}R false;}bool CmpEqual(C V<S>&l,sint index=-1)C{sint i;if(index==-1||index>=sint(S))i=S-1;else i=index;for(;i>=0;--i)if(T[i]!=l.T[i])R false;R true;}bool CmpSmallerEqual(C V<S>&l,sint index=-1)C{sint i;if(index==-1||index>=sint(S))i=S-1;else i=index;for(;i>=0;--i){if(T[i]!=l.T[i])R T[i]<l.T[i];}R true;}bool CmpBiggerEqual(C V<S>&l,sint index=-1)C{sint i;if(index==-1||index>=sint(S))i=S-1;else i=index;for(;i>=0;--i){if(T[i]!=l.T[i])R T[i]>l.T[i];}R true;}bool O<(C V<S>&l)C{R CmpSmaller(l);}bool O>(C V<S>&l)C{R CmpBigger(l);}bool O==(C V<S>&l)C{R CmpEqual(l);}bool O!=(C V<S>&l)C{R !O==(l);}bool O<=(C V<S>&l)C{R CmpSmallerEqual(l);}bool O>=(C V<S>&l)C{R CmpBiggerEqual(l);}V<S>O-(C V<S>&p2)C{V<S>temp(*this);temp.Sub(p2);R temp;}V<S>&O-=(C V<S>&p2){Sub(p2);R *this;}V<S>O+(C V<S>&p2)C{V<S>temp(*this);temp.Add(p2);R temp;}V<S>&O+=(C V<S>&p2){Add(p2);R *this;}V<S>O*(C V<S>&p2)C{V<S>temp(*this);temp.Mul(p2);R temp;}V<S>&O*=(C V<S>&p2){Mul(p2);R *this;}V<S>O/(C V<S>&p2)C{V<S>temp(*this);temp.Div(p2);R temp;}V<S>&O/=(C V<S>&p2){Div(p2);R *this;}V<S>O%(C V<S>&p2)C{V<S>temp(*this);V<S>remainder;temp.Div(p2,remainder);R remainder;}V<S>&O%=(C V<S>&p2){V<S>remainder;Div(p2,remainder);O=(remainder);R *this;}V<S>&O++(){AddOne();R *this;}V<S>O++(int){V<S>temp(*this);AddOne();R temp;}V<S>&O--(){SubOne();R *this;}V<S>O--(int){V<S>temp(*this);SubOne();R temp;}V<S>O~()C{V<S>temp(*this);temp.BitNot();R temp;}V<S>O&(C V<S>&p2)C{V<S>temp(*this);temp.BitAnd(p2);R temp;}V<S>&O&=(C V<S>&p2){BitAnd(p2);R *this;}V<S>O|(C V<S>&p2)C{V<S>temp(*this);temp.BitOr(p2);R temp;}V<S>&O|=(C V<S>&p2){BitOr(p2);R *this;}V<S>O^(C V<S>&p2)C{V<S>temp(*this);temp.BitXor(p2);R temp;}V<S>&O^=(C V<S>&p2){BitXor(p2);R *this;}V<S>O>>(int move)C{V<S>temp(*this);temp.Rcr(move);R temp;}V<S>&O>>=(int move){Rcr(move);R *this;}V<S>O<<(int move)C{V<S>temp(*this);temp.Rcl(move);R temp;}V<S>&O<<=(int move){Rcl(move);R *this;}private:template<class ostream_type,class string_type>static ostream_type&OutputToStream(ostream_type&s,C V<S>&l){string_type ss;l.ToString(ss);s<<ss;R s;}public:friend ostream&O<<(ostream&s,C V<S>&l){R OutputToStream<ostream,string>(s,l);}friend wostream&O<<(wostream&s,C V<S>&l){R OutputToStream<wostream,wstring>(s,l);}private:template<class istream_type,class string_type,class char_type>static istream_type&InputFromStream(istream_type&s,V<S>&l){string_type ss;char_type z;s>>z;while(s.good()&&Misc::CharToDigit(z,10)>=0){ss+=z;z=static_cast<char_type>(s.get());}s.unget();l.FromString(ss);R s;}public:friend istream&O>>(istream&s,V<S>&l){R InputFromStream<istream,string,char>(s,l);}friend wistream&O>>(wistream&s,V<S>&l){R InputFromStream<wistream,wstring,wchar_t>(s,l);}private:U Rcl2_one(U c);U Rcr2_one(U c);U Rcl2(U bits,U c);U Rcr2(U bits,U c);public:static C char*LibTypeStr();static LibTypeCode LibType();U Add(C V<S>&ss2,U c=0);U AddInt(U value,U index=0);U AddTwoInts(U x2,U x1,U index);static U AddVector(C U*ss1,C U*ss2,U ss1_size,U ss2_size,U*result);U Sub(C V<S>&ss2,U c=0);U SubInt(U value,U index=0);static U SubVector(C U*ss1,C U*ss2,U ss1_size,U ss2_size,U*result);static sint FindLeadingBitInWord(U x);static sint FindLowestBitInWord(U x);static U SetBitInWord(U&value,U bit);static void MulTwoWords(U a,U b,U*result_high,U*result_low);static void DivTwoWords(U a,U b,U c,U*r,U*rest);};template<>class V<0>{public:U T[1];void Mul2Big(C V<0>&,V<0>&){};void SetZero(){};U AddTwoInts(U,U,U){R 0;};};}namespace ttmath{template<U S>C char*V<S>::LibTypeStr(){static C char info[]="asm_gcc_32";R info;}template<U S>LibTypeCode V<S>::LibType(){LibTypeCode info=asm_gcc_32;R info;}template<U S>U V<S>::Add(C V<S>&ss2,U c){U b=S;U*p1=T;U*p2=const_cast<U*>(ss2.T);U dummy,dummy2;__asm__ __volatile__("xorl %%edx,%%edx \n""negl %%eax \n""1:\n""movl (%%esi,%%edx,4),%%eax \n""adcl %%eax,(%%ebx,%%edx,4)\n""incl %%edx \n""decl %%ecx \n""jnz 1b \n""adc %%ecx,%%ecx \n":"=c"(c),"=a"(dummy),"=d"(dummy2):"0"(b),"1"(c),"b"(p1),"S"(p2):"cc","memory");R c;}template<U S>U V<S>::AddInt(U value,U index){U b=S;U*p1=T;U c;U dummy,dummy2;__asm__ __volatile__("subl %%edx,%%ecx \n""1:\n""addl %%eax,(%%ebx,%%edx,4)\n""jnc 2f \n""movl $1,%%eax \n""incl %%edx \n""decl %%ecx \n""jnz 1b \n""2:\n""setc %%al \n""movzx %%al,%%edx \n":"=d"(c),"=a"(dummy),"=c"(dummy2):"0"(index),"1"(value),"2"(b),"b"(p1):"cc","memory");R c;}template<U S>U V<S>::AddTwoInts(U x2,U x1,U index){U b=S;U*p1=T;U c;U dummy,dummy2;__asm__ __volatile__("subl %%edx,%%ecx \n""addl %%esi,(%%ebx,%%edx,4)\n""incl %%edx \n""decl %%ecx \n""1:\n""adcl %%eax,(%%ebx,%%edx,4)\n""jnc 2f \n""mov $0,%%eax \n""incl %%edx \n""decl %%ecx \n""jnz 1b \n""2:\n""setc %%al \n""movzx %%al,%%eax \n":"=a"(c),"=c"(dummy),"=d"(dummy2):"0"(x2),"1"(b),"2"(index),"b"(p1),"S"(x1):"cc","memory");R c;}template<U S>U V<S>::AddVector(C U*ss1,C U*ss2,U ss1_size,U ss2_size,U*result){U rest=ss1_size-ss2_size;U c;U dummy1,dummy2,dummy3;__asm__ __volatile__("push %%edx \n""xor %%edx,%%edx \n""1:\n""mov (%%esi,%%edx,4),%%eax \n""adc (%%ebx,%%edx,4),%%eax \n""mov %%eax,(%%edi,%%edx,4)\n""inc %%edx \n""dec %%ecx \n""jnz 1b \n""adc %%ecx,%%ecx \n""pop %%eax \n""or %%eax,%%eax \n""jz 3f \n""xor %%ebx,%%ebx \n""neg %%ecx \n""mov %%eax,%%ecx \n""2:\n""mov (%%esi,%%edx,4),%%eax \n""adc %%ebx,%%eax \n""mov %%eax,(%%edi,%%edx,4)\n""inc %%edx \n""dec %%ecx \n""jnz 2b \n""adc %%ecx,%%ecx \n""3:\n":"=a"(dummy1),"=b"(dummy2),"=c"(c),"=d"(dummy3):"1"(ss2),"2"(ss2_size),"3"(rest),"S"(ss1),"D"(result):"cc","memory");R c;}template<U S>U V<S>::Sub(C V<S>&ss2,U c){U b=S;U*p1=T;U*p2=const_cast<U*>(ss2.T);U dummy,dummy2;__asm__ __volatile__("xorl %%edx,%%edx \n""negl %%eax \n""1:\n""movl (%%esi,%%edx,4),%%eax \n""sbbl %%eax,(%%ebx,%%edx,4)\n""incl %%edx \n""decl %%ecx \n""jnz 1b \n""adc %%ecx,%%ecx \n":"=c"(c),"=a"(dummy),"=d"(dummy2):"0"(b),"1"(c),"b"(p1),"S"(p2):"cc","memory");R c;}template<U S>U V<S>::SubInt(U value,U index){U b=S;U*p1=T;U c;U dummy,dummy2;__asm__ __volatile__("subl %%edx,%%ecx \n""1:\n""subl %%eax,(%%ebx,%%edx,4)\n""jnc 2f \n""movl $1,%%eax \n""incl %%edx \n""decl %%ecx \n""jnz 1b \n""2:\n""setc %%al \n""movzx %%al,%%edx \n":"=d"(c),"=a"(dummy),"=c"(dummy2):"0"(index),"1"(value),"2"(b),"b"(p1):"cc","memory");R c;}template<U S>U V<S>::SubVector(C U*ss1,C U*ss2,U ss1_size,U ss2_size,U*result){U rest=ss1_size-ss2_size;U c;U dummy1,dummy2,dummy3;__asm__ __volatile__("push %%edx \n""xor %%edx,%%edx \n""1:\n""mov (%%esi,%%edx,4),%%eax \n""sbb (%%ebx,%%edx,4),%%eax \n""mov %%eax,(%%edi,%%edx,4)\n""inc %%edx \n""dec %%ecx \n""jnz 1b \n""adc %%ecx,%%ecx \n""pop %%eax \n""or %%eax,%%eax \n""jz 3f \n""xor %%ebx,%%ebx \n""neg %%ecx \n""mov %%eax,%%ecx \n""2:\n""mov (%%esi,%%edx,4),%%eax \n""sbb %%ebx,%%eax \n""mov %%eax,(%%edi,%%edx,4)\n""inc %%edx \n""dec %%ecx \n""jnz 2b \n""adc %%ecx,%%ecx \n""3:\n":"=a"(dummy1),"=b"(dummy2),"=c"(c),"=d"(dummy3):"1"(ss2),"2"(ss2_size),"3"(rest),"S"(ss1),"D"(result):"cc","memory");R c;}template<U S>U V<S>::Rcl2_one(U c){U b=S;U*p1=T;U dummy,dummy2;__asm__ __volatile__("xorl %%edx,%%edx \n""negl %%eax \n""1:\n""rcll $1,(%%ebx,%%edx,4)\n""incl %%edx \n""decl %%ecx \n""jnz 1b \n""adcl %%ecx,%%ecx \n":"=c"(c),"=a"(dummy),"=d"(dummy2):"0"(b),"1"(c),"b"(p1):"cc","memory");R c;}template<U S>U V<S>::Rcr2_one(U c){U b=S;U*p1=T;U dummy;__asm__ __volatile__("negl %%eax \n""1:\n""rcrl $1,-4(%%ebx,%%ecx,4)\n""decl %%ecx \n""jnz 1b \n""adcl %%ecx,%%ecx \n":"=c"(c),"=a"(dummy):"0"(b),"1"(c),"b"(p1):"cc","memory");R c;}template<U S>U V<S>::Rcl2(U bits,U c){U b=S;U*p1=T;U dummy,dummy2,dummy3;__asm__ __volatile__("push %%ebp \n""movl %%ecx,%%esi \n""movl $32,%%ecx \n""subl %%esi,%%ecx \n""movl $-1,%%edx \n""shrl %%cl,%%edx \n""movl %%edx,%%ebp \n""movl %%esi,%%ecx \n""xorl %%edx,%%edx \n""movl %%edx,%%esi \n""orl %%eax,%%eax \n""cmovnz %%ebp,%%esi \n""1:\n""roll %%cl,(%%ebx,%%edx,4)\n""movl (%%ebx,%%edx,4),%%eax \n""andl %%ebp,%%eax \n""xorl %%eax,(%%ebx,%%edx,4)\n""orl %%esi,(%%ebx,%%edx,4)\n""movl %%eax,%%esi \n""incl %%edx \n""decl %%edi \n""jnz 1b \n""and $1,%%eax \n""pop %%ebp \n":"=a"(c),"=D"(dummy),"=S"(dummy2),"=d"(dummy3):"0"(c),"1"(b),"b"(p1),"c"(bits):"cc","memory");R c;}template<U S>U V<S>::Rcr2(U bits,U c){U b=S;U*p1=T;U dummy,dummy2,dummy3;__asm__ __volatile__("push %%ebp \n""movl %%ecx,%%esi \n""movl $32,%%ecx \n""subl %%esi,%%ecx \n""movl $-1,%%edx \n""shll %%cl,%%edx \n""movl %%edx,%%ebp \n""movl %%esi,%%ecx \n""xorl %%edx,%%edx \n""movl %%edx,%%esi \n""addl %%edi,%%edx \n""decl %%edx \n""orl %%eax,%%eax \n""cmovnz %%ebp,%%esi \n""1:\n""rorl %%cl,(%%ebx,%%edx,4)\n""movl (%%ebx,%%edx,4),%%eax \n""andl %%ebp,%%eax \n""xorl %%eax,(%%ebx,%%edx,4)\n""orl %%esi,(%%ebx,%%edx,4)\n""movl %%eax,%%esi \n""decl %%edx \n""decl %%edi \n""jnz 1b \n""roll $1,%%eax \n""andl $1,%%eax \n""pop %%ebp \n":"=a"(c),"=D"(dummy),"=S"(dummy2),"=d"(dummy3):"0"(c),"1"(b),"b"(p1),"c"(bits):"cc","memory");R c;}template<U S>sint V<S>::FindLeadingBitInWord(U x){sint result;U dummy;__asm__ ("movl $-1,%1 \n""bsrl %2,%0 \n""cmovz %1,%0 \n":"=r"(result),"=&r"(dummy):"r"(x):"cc");R result;}template<U S>sint V<S>::FindLowestBitInWord(U x){sint result;U dummy;__asm__ ("movl $-1,%1 \n""bsfl %2,%0 \n""cmovz %1,%0 \n":"=r"(result),"=&r"(dummy):"r"(x):"cc");R result;}template<U S>U V<S>::SetBitInWord(U&value,U bit){U old_bit;U v=value;__asm__ ("btsl %%ebx,%%eax \n""setc %%bl \n""movzx %%bl,%%ebx \n":"=a"(v),"=b"(old_bit):"0"(v),"1"(bit):"cc");value=v;R old_bit;}template<U S>void V<S>::MulTwoWords(U a,U b,U*result_high,U*result_low){U result1_;U result2_;__asm__ ("mull %%edx \n":"=a"(result1_),"=d"(result2_):"0"(a),"1"(b):"cc");*result_low=result1_;*result_high=result2_;}template<U S>void V<S>::DivTwoWords(U a,U b,U c,U*r,U*rest){U r_;U rest_;__asm__ ("divl %%ecx \n":"=a"(r_),"=d"(rest_):"0"(b),"1"(a),"c"(c):"cc");*r=r_;*rest=rest_;}}
#include <bits/stdc++.h>
using namespace std;
#define O operator
#define R return
#define C const
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-but-set-variable"
#ifdef __x86_64__
#error ttmath requires -m32
#endif
// ttmath, BSD licensed:
#include "incl.h"
#pragma GCC diagnostic pop
#undef C
#undef R
#undef O
#define MONTGOMERY
// Twice the numbers to factor should fit in 32*SmallSize bits.
enum {SmallSize = 4, LargeSize = 7, SmallBits = 96};
typedef ttmath::V<SmallSize> ll;
typedef ttmath::V<LargeSize> lll;
const static ll smallMask = (ll(1) << SmallBits) - 1;
// We trial-divide away all factors up to and including trial_lim.
const static ll trial_lim = 10000;
// ... thus, if we ask ourselves whether some number <= trial_lim^2 is a prime,
// the answer is yes.
const static ll trial_lim2 = trial_lim * trial_lim;
// The global modulus, with associated constants.
static ll Mod, NPrime;
static lll LMod, R1Mod, R2Mod;
void SetMod(const ll& x) {
lll B = smallMask; ++B;
lll X = x;
assert((x & 1) != 0);
assert(B > X);
lll R = B % X;
ll xinv = 1, bit = 2;
// We want xinv * x == 1 (mod B). We use Hensel lifting and compute xinv bit by bit.
for (int i = 1; i < SmallBits; i++, bit <<= 1) {
// for bit i:
ll y = xinv * x;
if ((y & bit) != 0)
xinv |= bit;
}
assert(((x * xinv) & smallMask) == 1);
::Mod = x;
::LMod = X;
::R1Mod = R;
::R2Mod = (R * R % X);
::NPrime = (ll)(B - xinv);
}
// Arithmetic modulo ::Mod, in Montgomery form.
#ifdef MONTGOMERY
struct N {
ll x;
N(struct undef* = 0) : x(0) {}
static N redc(const lll& a, const lll& b);
static N raw(const ll& x) { N r; r.x = x; return r; }
static N from(const ll& x) { assert (x < ::Mod); return redc(x, ::R2Mod); }
static N one() { return raw(::R1Mod); }
ll getRaw() const { return x; }
N pow(ll e) const;
};
N operator*(const N& a, const N& b) { return N::redc(a.x, b.x); }
N N::redc(const lll& a, const lll& b) {
lll T = a * b;
ll t = (ll)T & smallMask;
ll m = (t * ::NPrime) & smallMask;
T += (lll)m * ::LMod;
T >>= SmallBits;
if (T >= ::LMod)
T -= ::LMod;
return raw((ll)T);
}
#else
struct N {
ll x;
N(struct undef* = 0) : x(0) {}
static N raw(const ll& x) { N r; r.x = x; return r; }
static N from(const ll& x) { return raw(x); }
static N one() { return raw(1); }
ll getRaw() const { return x; }
N pow(ll e) const;
};
N operator*(const N& a, const N& b) { return N::raw(((lll)a.x * (lll)b.x) % ::LMod); }
#endif
N operator+(const N& a, const N& b) { ll x = a.x + b.x; if (x >= Mod) x -= Mod; return N::raw(x); }
N operator-(const N& a, const N& b) { ll x = a.x - b.x; if (a.x < b.x) x += Mod; return N::raw(x); }
inline N operator-(const N& a) { return 0 - a; }
inline void operator*=(N& a, const N& b) { a = a * b; }
inline void operator+=(N& a, const N& b) { a = a + b; }
inline void operator-=(N& a, const N& b) { a = a - b; }
inline bool operator==(const N& a, const N& b) { return a.x == b.x; }
inline bool operator!=(const N& a, const N& b) { return a.x != b.x; }
inline N db(const N& x) { return x + x; }
N N::pow(ll e) const {
N ret = one();
N pw = *this;
while (e != 0) {
if ((e & 1) != 0)
ret *= pw;
e >>= 1;
pw *= pw;
}
return ret;
}
// Factorization algorithm.
// Extract factors up to trial_lim from x.
void trial_division(vector<ll>& factors, ll& x) {
while (x % 2 == 0) { factors.push_back(2); x /= 2; }
while (x % 3 == 0) { factors.push_back(3); x /= 3; }
while (x % 5 == 0) { factors.push_back(5); x /= 5; }
int lut[8] = {6, 4, 2, 4, 2, 4, 6, 2}, it = 1;
ll i = 7;
while (i*i <= x && i < trial_lim) {
while (x % i == 0) {
factors.push_back(i);
x /= i;
}
i += lut[it++ & 7];
}
if (x != 1 && i*i > x) {
// Then we did trial division up to sqrt(x), so the number is a prime.
factors.push_back(x);
x = 1;
}
}
// Convert a ll to a double. Lossy, naturally.
double to_double(const ll& x) {
double base = pow(256, sizeof(ttmath::U)), pw = 1, res = 0;
for (size_t i = 0; i < x.Size(); i++) {
res += pw * x.T[i];
pw *= base;
}
return res;
}
// Compute whether x can be written as y^n for some n > 1. Returns (y, n)
// if so (for the smallest possible n), else (x, 1).
pair<ll, int> perfect_nth(const ll& x) {
// x has at most 29 digits, so its square root, cube root, etc. are exactly
// representable as doubles. Thus round(pow(y, 0.5)), etc. are exactly the
// nth roots, if there are any.
double y = to_double(x);
for (int i = 2;; i++) {
ll z = (long long)round(pow(y, 1.0 / i));
if (z < trial_lim) break;
ll w = z * z;
for (int j = 2; j < i; j++) w *= z;
if (w == x)
return {z, i};
}
return {x, 1};
}
// Check whether a number is a prime. Deterministic up to trial_lim^2;
// probabilistic (using the Miller-Rabin test) for larger numbers.
bool is_prime(const ll& n) {
if (n <= trial_lim2)
return true;
SetMod(n);
ll d = n - 1;
int r = 0;
while ((d & 1) != 0) {
r++;
d >>= 1;
}
N one = N::one(), minusone = 0 - one;
const int tries[] = {3,5,7,11,13,17,19};
for (size_t tr = 0; tr < sizeof(tries) / sizeof(*tries); tr++) {
N x = N::from(tries[tr]).pow(d);
if (x == one || x == minusone)
continue;
for (int i = 0; i < r-1; i++) {
x *= x;
if (x == one) return false;
if (x == minusone) goto next;
}
return false;
next:;
}
return true; // probably!
}
// Calculate the GCD of two numbers by the Euclidean algorithm.
ll gcd(const ll& a, const ll& b) {
if (b == 0)
return a;
return gcd(b, a % b);
}
// Check if all elements of a list are invertible, and return either 1 if so,
// or a factor of N otherwise (N if the first non-invertible element is 0).
// It's okay to use getRaw() here - Montgomery form preserves factors.
ll test_invertible(const vector<N>& l) {
N cur = N::one(), next;
for (const N& x : l) {
next = cur * x;
if (next == 0) {
// A multiplication resulted in 0, so either x = 0 or
// we found a non-trivial factor.
return gcd(cur.getRaw(), ::Mod);
}
cur = next;
}
ll g = gcd(cur.getRaw(), ::Mod);
return g;
}
// Return a random number modulo ::Mod, with a rather bad RNG.
N rand_mod_n() {
ll x = 0;
for (size_t i = 0; i < SmallSize; i++) {
x <<= 32;
x |= rand();
}
rand();
return N::from(x % ::Mod);
}
// Pollard rho
ll factor(ll n) {
SetMod(n);
for (;;) {
N x = rand_mod_n(), y = x;
vector<N> difs;
N o = N::one();
for (;;) {
x = x * x + o;
y = y * y + o;
y = y * y + o;
if (x == y) break;
difs.push_back(x - y);
if (difs.size() > 1000) {
ll f = test_invertible(difs);
if (f != 1 && f != n)
return f;
difs.clear();
}
}
ll f = test_invertible(difs);
if (f != 1 && f != n)
return f;
}
}
// Main program. Factors numbers from stdin and prints them as p_1^a_1 ... p_k^a^k.
// Terminates when there are no more numbers or the number 0 appears.
int main() {
srand(2);
ll x;
while (cin >> x && x != 0) {
vector<ll> factors;
trial_division(factors, x);
vector<ll> rem;
if (x != 1)
rem.push_back(x);
while (!rem.empty()) {
// Process one factor y of the number. We have three cases: either y
// is a prime, or it is a perfect nth power for some n > 1, or it
// has at least two distinct prime divisors.
ll y = rem.back();
rem.pop_back();
if (is_prime(y)) {
factors.push_back(y);
} else {
pair<ll, int> pn = perfect_nth(y);
if (pn.second > 1) {
for (int i = 0; i < pn.second; i++)
rem.push_back(pn.first);
} else {
ll z = factor(y), z2 = y / z;
assert(z * z2 == y);
rem.push_back(z);
rem.push_back(z2);
}
}
}
sort(factors.begin(), factors.end());
for (size_t i = 0; i < factors.size();) {
size_t j = i+1;
while (j < factors.size() && factors[j] == factors[i]) ++j;
cout << factors[i] << "^" << (j-i) << ' ';
i = j;
}
cout << endl;
}
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment