高精度计算

高精度运算
所谓的高精度运算,是指参与运算的数(加数,减数,因子……）范围大大超出了标准数据类型（整型，实型）能表示的范围的运算。

例如，求两个200位的数的和。

这时，就要用到高精度算法了。

在这里，我们先讨论高精度加法。

高精度运算主要解决以下三个问题：一、加数、减数、运算结果的输入和存储
运算因子超出了整型、实型能表示的范围，肯定不能直接用一个数的形式来表示。

在Pascal中，能表示多个数的数据类型有两种：数组和字符串。

数组：每个数组元素存储1位（在优化时，这里是一个重点！），有多少位就需要多少个数组元素；用数组表示数的优点：每一位都是数的形式，可以直接加减；运算时非常方便。

用数组表示数的缺点：数组不能直接输入；输入时每两位数之间必须有分隔符，不符合数值的输入习惯；
字符串：字符串的最大长度是255，可以表示255位。

用字符串表示数的优点：能直接输入输出，输入时，每两位数之间不必分隔符，符合数值的输入习惯；用字符串表示数的缺点：字符串中的每一位是一个字符，不能直接进行运算，必须先将它转化为数值再进行运算；运算时非常不方便；
综合以上所述，对上面两种数据结构取长补短：用字符串读入数据，用数组存储数据：var s1,s2:string;
a,b,c:array [1..260] of integer;
i,l,k1,k2:integer;
begin
write('input s1:');readln(s1);
write('input s2:');readln(s2);
{————读入两个数s1，s2，都是字符串类型}
l:=length(s1);{求出s1的长度，也即s1的位数；有关字符串的知识。

}
k1:=260;
for i:=l downto 1 do
begin
a[k1]:=ord(s1[i])-48;{将字符转成数值}
k1:=k1-1;
end;
k1:=k1+1;
{————以上将s1中的字符一位一位地转成数值并存在数组a中；低位在后（从第260位开始），高位在前（每存完一位，k1减1），完后，k1指向最高位} 对s2的转化过程和上面一模一样。

二、运算过程
在往下看之前，大家先列竖式计算35+86。

注意的问题：
（1）运算顺序：两个数靠右对齐；从低位向高位运算；先计算低位再计算高位；
（2）运算规则：同一位的两个数相加再加上从低位来的进位，成为该位的和；这个和去掉向高位的进位就成为该位的值；如上例：3+8+1=12，向前一位进1，本位的值是2；可借助MOD、DIV运算完成这一步；
（3）最后一位的进位：如果完成两个数的相加后，进位位值不为0，则应添加一位；（4）如果两个加数位数不一样多，则按位数多的一个进行计算；
if k1>k2 then k:=k1 else k:=k2;
y:=0;
for i:=260 downto k do
begin
x:=a[i]+b[i]+y;
c[i]:=x mod 10;
y:=x div 10;
end;
if y<>0 then begin k:=k-1;c[k]:=y; end;
三、结果的输出（这也是优化的一个重点）
按运算结果的实际位数输出
for i:=k to 260 do write(c[i]);
writeln;
例子：求两个数的加法
program sum;
var s,s1,s2:string;
a,b,c:array [1..260] of integer;
i,l,k1,k2:integer;
begin
write('input s1:');readln(s1);
write('input s2:');readln(s2);
l:=length(s1);
k1:=260;
for i:=l downto 1 do
begin
a[k1]:=ord(s1[i])-48;
k1:=k1-1;
end;
k1:=k1+1;
l:=length(s2);
k2:=260;
for I+:=l downto 1 do
begin
b[k2]:=ord(s2[i])-48;
k2:=k2-1;
end;
k2:=k2+1;
if k1>k2 then k:=k2 else k:=k1;
y:=0;
for i:=260 downto k do
begin
x:=a[i]+b[i]+y;
c[i]:=x mod 10;
y:=x div 10;
end;
if y<>0 then begin k:=k-1;c[k]:=y;
end;
for i:=k to 260 do write(c[i]);
writeln;
end.
四、优化：
以上的方法的有明显的缺点：
（1）浪费空间：一个整型变量（-32768~32767）只存放一位（0~9）；
（2）浪费时间：一次加减只处理一位；
针对以上问题，我们做如下优化：一个数组元素存放四位数；（integer的最大范围是32767，5位的话可能导致出界）。

具体方法：
l:=length(s1);
k1:=260;
repeat {————有关字符串的知识}
s:=copy(s1,l-3,4);
val(s,a[k1],code);
k1:=k1-1;
s1:=copy(s1,1,l-4);
l:=l-4;
until l<=0;
k1:=k1+1;
而因为这个改进，算法要相应改变：
（1）运算时：不再逢十进位，而是逢万进位（mod 10000; div 10000);
（2）输出时：最高位直接输出，其余各位，要判断是否足够4位，不足部分要补0；例如：1，23，2345这样三段的数，输出时，应该是100232345而不是1234567。

改进后的算法：
program sum;
var s1,s2,s:string;
a,b,c:array [1..260] of integer;
i,l,k1,k2,code:integer;
begin
write('input s1:');readln(s1);
write('input s2:');readln(s2);
l:=length(s1);
k1:=260;
repeat {————有关字符串的知识}
s:=copy(s1,l-3,4);
val(s,a[k1],code);
k1:=k1-1;
s1:=copy(s1,1,l-4);
l:=l-4;
until l<=0;
k1:=k1+1;
c[i]:=a[i]-b[i];
end;
▲打印结果：
先找到差的第一个非零数，如果差的所有位数都为零，就直接输出零。

如果不是，就输出符号位和差的第一位。

剩下部分，打印补足零:
因为优化后的高精度减法,是把每四个数位分成一段,而每一段则必须有四个
数,当有一段不足四个数时,就得用"0"补足.(如:第一位是'1',第二位是'34',第三位是'345',第四位是'8', 则应写为'1003403450008').注意:第一位不用补零,(如:第一位为'3',则写成'3').
while (c[k]=0) and (k<=260) do k:=k+1;
if k>260 then write('0')
else begin
write(fh,c[k]);{k是差的第1位；}
for i:=k+1 to 260 do
begin
if c[i]<1000 then write('0');
if c[i]<100 then write('0');
if c[i]<10 then write('0');
write(c[i]);
end;
end;
参考程序：
program Zfjianfa;
const n=25;
var s1,s2,s3,s4,s:string;
a,b,c:array[1..n] of integer;
i,k1,k2,l,code,jw:integer;
cq:char;
begin
readln(s1);
readln(s2);
k1:=length(s1);
k2:=length(s2);
if k1>k2 then for i:=1 to k1-k2 do s2:='0'+s2
else for i:=1 to k2-k1 do s1:='0'+s1;
cq:=' ';
if s1<s2 then begin cq:='-'; s:=s1; s1:=s2; s2:=s; end;
l:=length(s1);
k1:=n;
repeat
s3:=copy(s1,l-3,4);
val(s3,a[k1],code);
k1:=k1-1;
delete(s1,l-3,4);
l:=l-4;
until l<=0;
k1:=k1+1;
i:=length(s2);
k2:=n;
repeat
s4:=copy(s2,i-3,4);
val(s4,b[k2],code);
k2:=k2-1;
delete(s2,i-3,4);
i:=i-4;
until i<=0;
k2:=k2+1;
jw:=0;
for i:=n downto k1 do
begin
a[i]:=a[i]-jw;
jw:=0;
if a[i]<b[i] then begin
jw:=1;
a[i]:=a[i]+10000;
end;
c[i]:=a[i]-b[i];
end;
while (c[k1]=0) and (k1<=n) do
k1:=k1+1;
if k1>n then writeln('0')
else begin
write(cq,c[k1]);
for i:=k1+1 to n do
begin
if c[i]<1000 then write('0');
if c[i]<100 then write('0');
if c[i]<10 then write('0');
write(c[i]);
end;
end;
writeln;
end.
高精度乘法基本思想和加法一样。

其基本流程如下：
i:=n; k:=2*n;
repeat
x:=b[i]; z:=0; m:=k; jw:=0;
for j:=n downto k1 do
begin
y:=a[j];
z:=c[m];
x:=x*y+z+jw;
jw:=x div 10000;
c[m]:=x mod 10000;
m:=m-1;
x:=b[i];
end;
if jw<>0 then c[m]:=jw else m:=m+1;
i:=i-1;
k:=k-1;
until i<k2;
k:=m;
end;
procedure daying;
var i:integer;
begin
write(c[k]);
for i:=k+1 to 2*n do
begin
if c[i]<1000 then write('0');
if c[i]<100 then write('0');
if c[i]<10 then write('0');
write(c[i]);
end;
writeln;
end;
begin
init;
jisuan;
daying;
end.
教材“基础编”P87高精乘法参考程序：
program ex3_1;
var
a,b,c:array[0..1000] of word;
procedure init;
s:string;
ok,i,j:integer;
begin
readln(s);
a[0]:=length(s);
for i:=1 to a[0] do
val(s[a[0]-i+1],a[i],ok);
readln(s);
b[0]:=length(s);
b[0]:=length(s);
for i:=1 to b[0] do
val(s[b[0]-i+1],b[i],ok);
end;
procedure highmul;
var i,j,k:integer;
begin
c[0]:=a[0]+b[0];
for i:=1 to b[0] do
for j:=1 to a[0]+1 do
begin
inc(c[i+j-1],a[j]*b[i] mod 10);
c[i+j]:=c[i+j]+(a[j]*b[i] div 10)+(c[i+j-1] div 10);
c[i+j-1]:=c[i+j-1] mod 10;
end;
end;
procedure print;
var i:integer;
begin
while c[c[0]]=0 do dec(c[0]);
for i:=c[0] downto 1 do
write(c[i]);
writeln;
end;
{---main---}
begin
init;
highmul;
print;
end.
习题练习（一）
Today, facing the rapid development of business, ACM (Association of Calculator Management) recognizes that more powerful calculator should be studied, developed and appeared in future market shortly. ACM now invites you attending such amazing research and development work.
In most business application, the top three useful calculation operators are Addition (+), Subtraction (-) and Multiplication (*) between two given integers. Normally, you may think it is just a piece of cake. However, since some integers for calculation in business application may be very big, such as the GDP of the whole world, the calculator becomes harder to develop.
For example, if we have two integers 20000000000000000 and 4000000000000000, the exact results of the addition, subtraction and multiplication are:
20000000000000000 + 4000000000000000 = 24000000000000000 20000000000000000 - 4000000000000000 = 16000000000000000 20000000000000000 * 4000000000000000 = 80000000000000000000000000000000
Note: ACM prefers the exact format of the results rather than the float format or scientific remark format. For instance, we need "24000000000000000" rather than 2.4*1016.
As a programmer in ACM, your current task is to develop a program to obtain the exact results of the addition (a + b), subtraction (a - b) and multiplication (a * b) between two given integers a and b.
Input
The input consists of several test cases. Each case has two separate lines where the first line gives the integer a and the second gives b (|a| < 10201 and |b| < 10201).
When both a and b are equal to zero, the input ends.
Output
For each test case in the input, output three separate lines showing the exact results of addition (a + b), subtraction (a - b) and multiplication (a * b) of that case, one result per line.
Leaving a blank line between two successive test cases in the output.
Sample Input
20000000000000000
4000000000000000
Sample Output
24000000000000000
16000000000000000 80000000000000000000000000000000
#include<iostream>
using namespace std;
bool flag;
const int cell = 10000,ncell = 4;
const int step[5] = { 0,1000,100,10,1};
void abs StrToHi(int *arr,char *str)//base 10 {
int i,len,j,sum;
len=strlen(str);
arr[0] = 0;
for(i=len-1,j=ncell,sum = 0;i>=0;--i) {
sum += step[j] * (str[i] - 48); j--;
if(j == 0)
{
arr[++arr[0]] = sum;
sum = 0;
j=4;
}
}
if(j < 4 && sum > 0) arr[++arr[0]] = sum;
}
void strToHi(char *str,int * arr)
{
if(str[0]=='-')
{
absStrToHi(arr,str+1);
flag = false;
arr[0] *= -1;
}
else
{
absStrToHi(arr,str);
if(arr[0] > 0)
flag = false;
}
}
void prHi(int *arr)//base 10
{
int i;
if(arr[0] < 0)
printf("-");
if(arr[0]!=0) printf("%d",arr[abs(arr[0])]);
else printf("0");
for(i = abs(arr[0]) - 1 ; i >= 1 ; --i)
printf("%04d",arr[i]);//此处根据cellLen改动 printf("\n");
}
int absbignumCmp(int *a,int *b)
{
int i;
if(abs(a[0]) != abs(b[0]))
return abs(a[0]) - abs(b[0]);
for(i = abs(a[0]); i >= 1 ; --i)//对于 0 成立
{
if(a[i] != b[i])
return a[i] - b[i];
}
return 0;
}
void absAdd(int *a,int *b,int *c)
{
int i,la = abs(a[0]),lb = abs(b[0]); int carry;
for(carry=0,i=1;i<= la&& i<=lb;++i) {
c[i] = carry + a[i] + b[i]; if(c[i] >= cell)
{
c[i] -= cell;
carry = 1;
}
else carry = 0;
}
while(i<=la)
{
c[i] = carry + a[i];
if(c[i] >= cell)
{
c[i] -= cell;
carry = 1;
}
else carry = 0;
i++;
}
while(i<=lb)
{
c[i] = carry + b[i];
if(c[i] >= cell)
{
c[i] -= cell;
carry = 1;
}
else carry = 0;
i++;
}
if(carry)
{
c[0] = i ;
c[i]=carry;
}
else c[0] = i - 1;
void absSub(int *a,int *b,int *c)// a>=b and c can be equal to a or b
{
int i,la = abs(a[0]),lb = abs(b[0]);
int borrow;
for(borrow=0,i=1;i<=lb;++i)
{
c[i] = a[i] - borrow - b[i];
if(c[i] < 0)
{
c[i] += cell;
borrow=1;
}
else borrow=0;
}
while(i<=la)
{
c[i] = a[i] - borrow;
if(c[i] < 0)
{
c[i] += cell;
borrow=1;
}
else borrow=0;
i++;
}
i--;
while(i>0 && c[i] == 0)
i--;
c[0] = i;
}
void add(int *a,int *b,int *c)
{
int cmp,t;
if(a[0] == 0 )
{
memcpy(c,b,sizeof(unsigned int) * (abs(b[0]) +1));
return;
}
if(b[0] == 0)
memcpy(c,a,sizeof(unsigned int) * (abs(a[0]) + 1)); return;
}
t = a[0]*b[0];
if( t > 0 )
{
absAdd(a,b,c);
c[0] *= a[0] / abs(a[0]);
}
else
{
cmp = absbignumCmp(a,b);
if(cmp > 0)// a!=0
{
absSub(a,b,c);
c[0] *= a[0] / abs(a[0]);
}
else if(cmp == 0)
c[0] = 0;
else
{
absSub(b,a,c);
c[0] *= b[0] / abs(b[0]);
}
}
}
void sub(int *a,int *b,int *c)
{
b[0] *= -1;
add(a,b,c);
b[0] *= -1;
}
void mul(int *a,int *b,int *c)
{
int i,j,la = abs(a[0]),lb = abs(b[0]),ai,bj,stc,k;
int carry;
memset(c,0,sizeof(unsigned int)*(la+lb+2));
if(a[0] == 0 || b[0] == 0)
c[0] = 0;
else
stc=1;
ai=1;
while(a[ai]==0)
{
ai++;
stc++;
}
bj=1;
while(b[bj]==0)
{
bj++;
stc++;
}
for(i=ai;i<=la;i++)
{
for(j=bj,carry = 0;j<=lb;j++)
{
k = stc + j - bj;
c[k] += a[i] * b[j] + carry; carry = c[k] / cell;
c[k] %= cell;
}
while(carry > 0)
{
c[++k] += carry;
carry = c[k] / cell;
c[k] %= cell;
}
stc++;
}
i=la+lb+1;
while(i>0&&c[i]==0)
i--;
c[0] = i;
if(a[0] * b[0] < 0)
c[0] *= -1;
}
}
int main()
{
char c[251];
int out[111],a[55],b[55];;
bool ff=false;
while(1)
{
flag = true;
scanf("%s",c);
strToHi(c,a);
scanf("%s",c);
strToHi(c,b);
if(flag)
break;
if(ff)
printf("\n");
ff = true;
add(a,b,out);
prHi(out);
sub(a,b,out);
prHi(out);
mul(a,b,out);
prHi(out);
}
return 0;
}
习题练习（二）
/toj/showp.php?pid=2467
T2469
#include<iostream>
using namespace std;
bool flag;
void absStrToHi(unsigned char *arr,char *str,int &ln)//base 10 {
int i,len;
len=strlen(str);
for(i=len-1;i>=0;--i)
arr[len - i] = str[i] - 48;
ln = len;
}
void strToHi(char *str,unsigned char * arr, int &ln)
{
if(str[0]=='-')
{
absStrToHi(arr,str+1,ln);
flag = false;
if(ln > 1 || arr[1] > 0)
arr[0] = 1;
}
else
{
absStrToHi(arr,str,ln);
arr[0] = 0;
if(ln > 1 || arr[1] > 0)
flag = false;
}
}
void prHi(unsigned char *arr,int &ln)//base 10
{
int i;
if(arr[0] == 1)
printf("-");
for(i = ln;i >= 1; --i)
printf("%c",arr[i] + 48);//此处根据cellLen改动
printf("\n");
}
int absbignumCmp(unsigned char *a,unsigned char *b,int &la,int &lb) {
int i;
if(la != lb)
return la - lb;
for(i = la; i>=1 ; --i)//对于 0 成立
{
if(a[i] != b[i])
return a[i] - b[i];
}
return 0;
}
void absAdd(unsigned char *a,unsigned char *b,unsigned char *c,int &la,int &lb,int &lc)// c can be equal to a or b
{
int i;
unsigned char carry;
for(carry=0,i=1;i<=la&& i<=lb;++i)
{
c[i] = carry + a[i] + b[i];
if(c[i] > 9)
{
c[i] -= 10;
carry = 1;
}
else carry = 0;
}
while(i<=la)
{
c[i] = carry + a[i];
if(c[i] > 9)
{
c[i] -= 10;
carry = 1;
}
else carry = 0;
i++;
}
while(i<=lb)
{
c[i] = carry + b[i];
if(c[i] > 9)
{
c[i] -= 10;
carry = 1;
}
else carry = 0;
i++;
}
if(carry) lc = i , c[i]=carry;
else lc = i - 1;
}
void absSub(unsigned char *a,unsigned char *b,unsigned char *c,int &la,int &lb,int &lc)// a>=b and c can be equal to a or b {
int i;
unsigned char borrow;
for(borrow=0,i=1;i<=lb;++i)
{
c[i] = a[i]-borrow-b[i];
if(c[i] > 127)
{
c[i] += 10;
borrow=1;
}
else borrow=0;
}
while(i<=la)
{
c[i] = a[i] - borrow;
if(c[i] > 127)
{
c[i] += 10;
borrow=1;
}
else borrow=0;
i++;
}
i--;
while(i>1 && c[i] == 0)
i--;
lc = i;
}
void add(unsigned char *a,unsigned char *b,unsigned char *c,int &la,int &lb,int &lc)
{
int cmp;
if(a[0] == b[0])
{
absAdd(a,b,c,la,lb,lc);
c[0]=a[0];
}
else
{
cmp = absbignumCmp(a,b,la,lb);
if(cmp > 0)// a!=0
{
absSub(a,b,c,la,lb,lc);
c[0] = a[0];
else if(cmp == 0)
{
c[0] = 0;
c[1] = 0;
lc = 1;
//printf("Zero\n");
}
else
{
absSub(b,a,c,lb,la,lc);
c[0] = b[0];
}
}
}
void sub(unsigned char *a,unsigned char *b,unsigned char *c,int &la,int &lb,int &lc)
{
b[0] = 1 - b[0];
add(a,b,c,la,lb,lc);
b[0] = 1 - b[0];
}
void mul(unsigned char *a,unsigned char *b,unsigned char *c,int &la,int &lb,int &lc)
{
int i,j,ai,bj,stc,k;
unsigned char carry;
memset(c,0,sizeof(unsigned char)*(la+lb+2));
if((la==1&&a[1]==0)||(lb==1&&b[1]==0))
lc = 1;
else
{
stc=1;
ai=1;
while(a[ai]==0)
{
ai++;
stc++;
}
bj=1;
while(b[bj]==0)
bj++;
stc++;
}
for(i=ai;i<=la;i++)
{
for(j=bj,carry = 0;j<=lb;j++)
{
k = stc + j - bj;
c[k] += a[i] * b[j] + carry; carry = c[k] / 10;
c[k] %= 10;
}
while(carry > 0)
{
c[++k] += carry;
carry = c[k] / 10;
c[k] %= 10;
}
stc++;
}
i=la+lb+1;
while(i>1&&c[i]==0)
i--;
lc = i;
if(a[0] != b[0])
c[0] = 1;
}
}
int main()
{
char c[251];
unsigned char out[500],a[250],b[250];;
int la,lb,lo;
bool ff=false;
while(1)
{
flag = true;
scanf("%s",c);
strToHi(c,a,la);
scanf("%s",c);
strToHi(c,b,lb);
if(flag)
break;
if(ff)
printf("\n");
ff = true;
add(a,b,out,la,lb,lo);
prHi(out,lo);
sub(a,b,out,la,lb,lo);
prHi(out,lo);
mul(a,b,out,la,lb,lo);
prHi(out,lo);
}
return 0;
}
习题练习（三）
“At a time when people find that it is hard to do a large calculation by hand, then mankind invented computer.” Computers are often used to do large calculation. In this problem, you are to write a program to do some awful arithmetic operation.
Input:
The input consist of several cases. Each case will be two integers in a line, A and B. Note the A and B will be very large number. But to simplify this problem, the length of A and B will be no more then 1000. a line of A=0 and B=0 terminate the input.
Output:
For each A and B. print the result of A+B in a line, and then the result of A*B in a line.
Sample Input:
1 2
10 100000000000000000000000000000000
0 0
Sample Output:
3
2
100000000000000000000000000000010
1000000000000000000000000000000000
//Author CC
//Date from 04.10.8
//The implementation of USACO's idea
//
/***
1. 各种运算的操作数必须是正整数，至于符号转换根据不同程序要求自行完成。

2. int *arr arr[0]表示存储大数所需cell数，arr[0]的符号表示正负
当 arr[0] == 0 时，表示大数为 0
***/
//除法时间 Total time is 537172 //upgraded Total time is 84203
//乘法时间 Total time is 435656 //upgraded Total time is 285766
////有些函数参数顺序变了！！！！！！！！！！！！！！！！！！！
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cassert>
#include<ctime>
#include<cctype>
using namespace std;
#ifdef _MSC_VER
#define i64 __int64
#else
#define i64 long long
#endif
#define abs(x) ((x)>0?(x):(-(x)))
#define max(a,b)((a)>(b)?(a):(b))
#define min(a,b)((a)<(b)?(a):(b))
#define sgn(a) ((a)==0?0:(a)/abs(a))
const unsigned int cell=1000000000,cellLen=9;//每个cell保存多少数
//--------------------------------------------------------------------------
/* type transformation and SYSIO for base ten number
* convert string to the standard expression in the following programme
*/
void absStrToHi(int *arr,char *str)//base 10
{
int i,len,sum;
len=strlen(str);
assert(len!=0);
for(sum=0,i=0;i<len;++i)
{
sum*=10;
sum+=str[i]-'0';
if((len-i-1)%cellLen==0)
{
arr[(len-i-1)/cellLen+1]=sum; sum=0;
}
}
arr[0]=(len-1)/cellLen+1;
//processing leading zeros
for(i=arr[0];i>=1;--i)
if(arr[i]!=0)
break;
arr[0]=i;
}
void strToHi(char *str,int *arr)
{
if(str[0]=='+')
{
absStrToHi(arr,str+1);
return ;
}
else if(str[0]=='-')
{
absStrToHi(arr,str+1);
arr[0]*=-1;
return ;
}
else
absStrToHi(arr,str);
}
void prHi(int *arr)//base 10
{
int i;
if(arr[0]<0)
printf("-");
printf("%d",arr[abs(arr[0])]);//对于 0 成立
for(i=abs(arr[0])-1;i>=1;--i)
printf("%09d",arr[i]);//此处根据cellLen改动
}
char *hiToStr(int *arr,char *ans)//base 10
{
int offset=0,i;
if(arr[0]<0)
ans[0]='-',offset=1;
sprintf(ans+offset,"%d",arr[abs(arr[0])]);//对于 0 成立
offset=strlen(ans);
for(i=abs(arr[0])-1;i>=1;--i,offset+=9)
sprintf(ans+offset,"%09d",arr[i]);//此处根据cellLen改动
return ans;
}
//------------------------------------------------------------------ /* two bignum comparison
* if( a < b ) return a negtive number
* if( a > b ) return a positive number
* if( a = b ) return 0
*/
int bignumCmp(int *a,int *b)
{
if(a[0]!=b[0])
return a[0]-b[0];
int i,re;
for(i=abs(a[0]);i>=1;--i)//对于 0 成立
{
if((re=a[i]-b[i])!=0)
return re*(a[0]/abs(a[0]));
}
return 0;
}
//------------------------------------------------------------------ /* two bignum addition and substraction
* c and equal to any of the other parameter
*/
void absAdd(int *a,int b,int *c)
{
int i,carry,lim=max(1,abs(a[0])),sum;
for(carry=0,i=1;i<=lim;++i)
{
sum=carry;
if(i<=a[0])
sum+=a[i];
if(i==1) sum+=b;
carry=sum/cell;
c[i]=sum%cell;
}
if(carry)
c[lim+1]=1,c[0]=lim+1;
else
c[0]=lim;
}
void absAdd(int *a,int *b,int *c)// c can be equal to a or b {
int i,lim=max(abs(a[0]),abs(b[0])),carry,sum;
for(carry=0,i=1;i<=lim;++i)
{
sum=carry;
if(i<=a[0]) sum+=a[i];
if(i<=b[0]) sum+=b[i];
carry=sum/cell;
c[i]=sum%cell;
}
if(carry)
c[0]=lim+1,c[lim+1]=carry;
else
c[0]=lim;
}
void absSub(int *a,int b,int *c)//a>=b
{
int i,lim=max(1,abs(a[0])),borrow;
if(b==0)
{
memcpy(c,a,sizeof(int)*(abs(a[0])+1));
return ;
}
for(borrow=0,i=1;i<=lim;++i)
{
c[i]=a[i]-borrow;
if(i==1)
c[i]-=b;
if(c[i]<0)
c[i]+=cell,borrow=1;
else
borrow=0;
}
for(i=lim;i>=1;--i)
if(c[i]!=0)
break;
c[0]=i;
assert(borrow==0 && "|B| > |A|; can't handle this");
}
void absSub(int *a,int *b,int *c)// a>=b and c can be equal to a or b
{
int i,lim=max(abs(a[0]),abs(b[0])),borrow,blen=abs(b[0]);
for(borrow=0,i=1;i<=lim;++i)
{
if(i<=blen)
c[i]=a[i]-borrow-b[i];
else
c[i]=a[i]-borrow;
if(c[i]<0)
c[i]+=cell,borrow=1;
else
borrow=0;
}
for(i=lim;i>=1;--i)
if(c[i]!=0)
break;
c[0]=i;
assert(borrow==0 && "|B| > |A|; can't handle this");
}
void add(int *a,int *b,int *c)
{
int signa=sgn(a[0]),signb=sgn(b[0]);
if(signa==signb)
{
absAdd(a,b,c);
c[0]*=signa;
}
else
{
a[0]*=signa,b[0]*=signb;
if(bignumCmp(a,b)>0)// a!=0
{
absSub(a,b,c);
c[0]*=signa;
}
else// if b==0 then a==0
{
absSub(b,a,c);
c[0]*=signb;
}
if(a!=c)
a[0]*=signa;
if(b!=c)
b[0]*=signb;
}
}
void sub(int *a,int *b,int *c)
{
b[0]=-b[0];
add(a,b,c);
b[0]=-b[0];
}
//------------------------------------------------------------------ /* multiplication of a bignum and a 32-bits number */ void mulByScalar(int *a,int b,int *c)// a can be equal to c
{
int i,lim=abs(a[0]),carry,signb=sgn(b);
i64 tmp;
for(carry=0,i=1;i<=lim;++i)
{
tmp=a[i],tmp*=b,tmp+=carry;
carry=tmp/cell;
c[i]=tmp%cell;
}
if(carry)
c[0]=(lim+1)*signb*sgn(a[0]),c[lim+1]=carry;
else
c[0]=lim*signb*sgn(a[0]);
}
void mulAndAdd(int *a,int s,int offset,int *c)//c is always nonnetgive {
if(a[0]==0||s==0)
return ;
//a!=0 && s!=0
int i,sizea=abs(a[0]),carry,sizec=abs(c[0]),lim;
lim=max(sizea+offset,sizec);
i64 sum;
for(carry=0,i=min(offset+1,max(sizec+1,1));i<=lim||carry;++i) {
if(i>offset&&i<=sizea+offset)
{
sum=s;
sum*=a[i-offset];
}
else
sum=0;
if(i<=sizec)
sum+=c[i];
sum+=carry;
c[i]=sum%cell;
carry=sum/cell;
}
c[0]=max(c[0],i-1);
}
void mul(int *a,int *b,int *c)
{
int i,sizeb=abs(b[0]);
c[0]=0;
for(i=1;i<=sizeb;++i)
mulAndAdd(a,b[i],i-1,c);
c[0]*=sgn(a[0])*sgn(b[0]);
}
//------------------------------------------------------------------ /* Division by a scalar'
* return the remainder
* s < 1 000 000 000
*/
int divByScalar(int *a,int s,int *c)// s>=0 a can be equal to c {
int i,signa=sgn(a[0]);
i64 remain;
assert(s!=0);
for(remain=0,i=abs(a[0]);i>=1;--i)
{
remain*=cell;
remain+=a[i];
c[i]=remain/s;
remain%=s;
}
if(c[abs(a[0])]!=0)
c[0]=signa*abs(a[0]);
else
c[0]=signa*(abs(a[0])-1);
return remain*signa;
}
/* Do this in other base, just change cell [be aware of the reverse order of a] */
int modByScalar(int *a,int s)
{
int i;
i64 remain;
for(remain=0,i=abs(a[0]);i>=1;--i)
{
remain*=cell;
remain+=a[i];
remain%=s;
}
return remain*sgn(a[0]);
}
/* need bignumCmp(int *a,int *b) function
* if b == 0 then return false
* only support base 10^9
*/
bool divByBignum(int *a,int *b,int *c,int *rem)
{
if(b[0]==0)
return false;
int i,from,to,mid,sizeb=abs(b[0]),signb=sgn(b[0]),k;
for(k=1;;k*=10)
if(k*b[sizeb]>=100000000)
break;
if(k!=1)
mulByScalar(a,k,a), mulByScalar(b,k,b);
b[0]*=signb;
int *d=new int [sizeb+3];
i64 tmp1,tmp2;
rem[0]=0,c[0]=0;// rem>=0
for(i=abs(a[0]);i>=1;--i)
{
memmove(rem+2,rem+1,rem[0]*sizeof(int));
rem[1]=a[i];
if(rem[rem[0]+1]!=0)
++rem[0];
if(bignumCmp(rem,b)<0)
{
c[i]=0;
continue;
}
// rem>=b so rem/b>=1
tmp1=rem[rem[0]],tmp2=b[sizeb];
if(rem[0]>sizeb)
tmp1*=cell,tmp1+=rem[rem[0]-1];
else if(sizeb>1)
tmp1*=cell,tmp1+=rem[rem[0]-1],tmp2*=cell,tmp2+= b[sizeb-1];
from=max(1,tmp1/(tmp2+1));
to=min(cell-1,(tmp1+1)/tmp2);
mid=-1;//bug here fixed
while(from<to)// from <= x <= to, when from == to find it
{
mid=(from+to)/2+1;
mulByScalar(b,mid,d);
if(bignumCmp(rem,d)>=0)
from=mid;
else
to=mid-1;
}
c[i]=from;
if(mid!=from)
mulByScalar(b,from,d);
absSub(rem,d,rem);
if(c[i]>0&&i>c[0])
c[0]=i;
}
//
if(k!=1)
divByScalar(a,k,a), divByScalar(b,k,b),
divByScalar(rem,k,rem);
b[0]*=signb;
c[0]*=sgn(a[0])*signb;
rem[0]*=sgn(a[0]);
delete [] d;
return true;
}
//------------------------------------------------------------------ /* base conversion */
void absStrToHi(int base, char map[],char *str, int *hi)
{
int i;
hi[0]=0;
for(i=0;str[i];++i)
{
mulByScalar(hi,base,hi);
absAdd(hi,map[str[i]],hi);
}
}
void hiToStr(int base, char map[],int *hi,char *str)
{
if(hi[0]==0)
{
str[0]='0',str[1]=0;
return ;
}
int i,r;
for(i=0;hi[0];++i)
{
r=divByScalar(hi,base,hi);
str[i]=map[r];
}
str[i]=0;
reverse(str,str+i);
}
//------------------------------------------------------------------ /* nonStorage factorial calculation
Total time is 0.189s when n == 4000
Total time is 0.007s when n == 1000
assert(n>=0); */
void fact(int n,int *ans)
{
ans[0]=1,ans[1]=1;
int i;
i64 p;
for(p=1,i=2;i<=n;++i)
{
p*=i;
if(p*(i+1)>cell)
{
mulByScalar(ans,p,ans);
p=1;
}
}
if(p>1)
mulByScalar(ans,p,ans);
}
/* nonStorage combination number calculation
Total time is 0.124s when n == 4000 k == 2000
Total time is 0.005s when n == 1000 k == 500
assert(n>=k&&k>=0); */
void C(int n,int k,int *ans)
{
if(n-k<k)
{
C(n,n-k,ans);
return;
}
ans[0]=ans[1]=1;
int i;
i64 p;
for(p=1,i=1;i<=k;++i)
{
p*=n-i+1;
if(p*(n-i)>cell)
{
mulByScalar(ans,p,ans);
p=1;
}
}
if(p>1)
mulByScalar(ans,p,ans);
for(p=1,i=1;i<=k;++i)
{
p*=i;
if(p*(i+1)>cell)
{
divByScalar(ans,p,ans);
p=1;
}
}
if(p>1)
divByScalar(ans,p,ans);
}
//------------------------------------------------------------------ #define MAXSTRL 1000000
char in1[MAXSTRL+10],in2[MAXSTRL+10];
int a[MAXSTRL/9+10],b[MAXSTRL/9+10],c[MAXSTRL/9+10],d[MAXSTRL/9+10];
int main()
{
freopen("d:\\test.in","r",stdin);
freopen("d:\\test.out","w",stdout);
int cnt=0,t1=clock();
/*
while(scanf("%s %s",in1,in2)==2)
{
strToHi(a,in1);
strToHi(b,in2);
/*
if((re=bignumCmp(a,b))>0)
printf("a>b\n");
else if(re<0)
printf("a<b\n");
else
printf("a=b\n");
*/
/*
if(divByBignum(a,b,c,d)==true)
{
prHi(c);
printf(" ");。