Created
March 7, 2020 17:22
-
-
Save simonlindholm/193be06491e539c734c1ea78782c9125 to your computer and use it in GitHub Desktop.
Pollard rho
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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&÷nd.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_;}} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#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