高精度浮点数
转自51nod,除法居然比自己写的Wint还快…
有空慢慢消化吧。
#include<bits/stdc++.h>
using namespace std;
const int SIZE=264000;
complex<double> e[SIZE],A[SIZE],w1[SIZE],w2[SIZE];
void DFT(complex<double> *A, complex<double> *w, int l, int r, int dt)//A系数矩阵 w值
{
int N = r - l;//区间长度
if (N == 1)
{
w[l] = A[l];
return;
}
int N1 = N >> 1;
int mid = (l + r) >> 1;//区间中点
int pos1 = l;
int pos2 = mid;
for (int i = 0; i < N; i++)
{
if (i & 1)
{
w[pos2] = A[i + l];
pos2++;
}
else
{
w[pos1] = A[i + l];
pos1++;
}
}
DFT(w, A, l, mid, dt << 1);
DFT(w, A, mid, r, dt << 1);
for (int i = 0; i < N1; i++)
{
w[l + i] = A[l + i] + e[i*dt] * A[i + mid];
w[i + mid] = A[l + i] - e[i*dt] * A[i + mid];
}
}
struct Wfloat
{
const static int size = SIZE;
int *a;
int length;
int sign;
int pos;
Wfloat()
{
a = (int*)malloc(sizeof(int)* size);
a[0] = 0;
length = 1;
sign = 0;
pos = 0;
}
Wfloat(const Wfloat &ans)
{
a = (int*)malloc(sizeof(int)* size);
length = ans.length;
sign = ans.sign;
pos = ans.pos;
for (int i = 0; i < length; i++)
{
a[i] = ans.a[i];
}
}
Wfloat(char *word)
{
a = (int*)malloc(sizeof(int)* size);
int wordlength = strlen(word);
pos = 0;
length = 0;
sign = 0;
for (int i = wordlength - 1; i >= 0; i--)
{
char c = word[i];
if (c == '.')
{
pos = length + 1;
}
else if (c == '-')
{
sign = 1;
}
else
{
a[length] = c - '0';
length++;
}
}
}
Wfloat operator = (const Wfloat &ans)
{
length = ans.length;
sign = ans.sign;
pos = ans.pos;
for (int i = 0; i < length; i++)
{
a[i] = ans.a[i];
}
return *this;
}
Wfloat operator =(char *word)//不含小数点
{
int wordlength = strlen(word);
pos = 0;
length = 0;
sign = 0;
for (int i = wordlength - 1; i >= 0; i--)
{
char c = word[i];
if (c == '.')
{
pos = length + 1;
}
else if (c == '-')
{
sign = 1;
}
else
{
a[length] = c - '0';
length++;
}
}
return *this;
}
~Wfloat()
{
free(a);
}
};
void Print(const Wfloat &ans)
{
if (ans.sign == 1) printf("-");
if (ans.pos >= ans.length)
{
printf("0.");
for (int i = 0; i < ans.pos - ans.length; i++)
{
printf("0");
}
for (int i = ans.length - 1; i >= 0; i--)
{
printf("%d", ans.a[i]);
}
}
else
{
for (int i = ans.length - 1; i >= ans.pos; i--)
{
printf("%d", ans.a[i]);
}
if (ans.pos != 0)
{
printf(".");
for (int i = ans.pos - 1; i >= 0; i--)
{
printf("%d", ans.a[i]);
}
}
}
printf("\n");
}
Wfloat operator * (const Wfloat &number1, int number2)
{
Wfloat ans = number1;
if (number2 < 0)
{
ans.sign = -number1.sign;
number2 = -number2;
}
for (int i = 0; i < ans.length; i++)
{
ans.a[i] *= number2;
}
for (int i = 0; i < ans.length; i++)
{
if (ans.a[i] >= 10)
{
if (i == ans.length - 1)
{
ans.length++;
ans.a[i + 1] = 0;
}
ans.a[i + 1] += ans.a[i] / 10;
ans.a[i] %= 10;
}
}
return ans;
}
Wfloat operator / (const Wfloat &number1, int number2)
{
Wfloat answer;
int lastnumber = 0;
int length = -1;
for (int i = number1.length - 1; i >= 0; i--)
{
lastnumber *= 10;
lastnumber += number1.a[i];
answer.a[i] = lastnumber / number2;
lastnumber -= answer.a[i] * number2;
if (length == -1 && answer.a[i] != 0)
{
length = i + 1;
}
}
answer.length = length;
answer.pos = number1.pos;
answer.sign = number1.sign;
return answer;
}
Wfloat Multiplication(const Wfloat &number1, const Wfloat &number2, int eps)
{
Wfloat ans;
double pi = 3.1415926535897932384626;
int length1 = min(eps, number1.length);
int length2 = min(eps, number2.length);
if (length1 > 1000 && length2 > 1000)
{
int tempLength = max(length1 / 3, length2 / 3);
int N = 1;
while (N < tempLength)
{
N <<= 1;
}
N <<= 1;
for (int i = 0; i < N; i++)
{
e[i] = complex<double>(cos(2 * pi * i / N), sin(2 * pi * i / N));
}
int lastPos;
for (int i = 0; i < length1; i++)
{
lastPos = i / 3;
if (i % 3 == 0)
{
A[lastPos] = complex<double>(number1.a[number1.length - length1 + i], 0);
}
else if (i % 3 == 1)
{
A[lastPos] += complex<double>(number1.a[number1.length - length1 + i] * 10.0, 0);
}
else if (i % 3 == 2)
{
A[lastPos] += complex<double>(number1.a[number1.length - length1 + i] * 100.0, 0);
}
}
for (int i = lastPos + 1; i < N; i++)
{
A[i] = complex<double>(0, 0);
}
DFT(A, w1, 0, N, 1);
for (int i = 0; i < length2; i++)
{
lastPos = i / 3;
if (i % 3 == 0)
{
A[lastPos] = complex<double>(number2.a[number2.length - length2 + i], 0);
}
else if (i % 3 == 1)
{
A[lastPos] += complex<double>(number2.a[number2.length - length2 + i] * 10.0, 0);
}
else if (i % 3 == 2)
{
A[lastPos] += complex<double>(number2.a[number2.length - length2 + i] * 100.0, 0);
}
}
for (int i = lastPos + 1; i < N; i++)
{
A[i] = complex<double>(0, 0);
}
DFT(A, w2, 0, N, 1);
for (int i = 0; i < N; i++)
{
w2[i] *= w1[i];
}
for (int i = 0; i < N; i++)
{
e[i] = complex<double>(cos(-2 * pi * i / N), sin(-2 * pi * i / N));
}
DFT(w2, w1, 0, N, 1);
long long temp = 0;
long long temp1 = 0;
for (int i = 0; i < N; i++)
{
temp += w1[i].real() / N + 0.001;
temp1 = temp / 1000;
temp %= 1000;
ans.a[i * 3] = temp % 10;
ans.a[i * 3 + 1] = temp % 100 / 10;
ans.a[i * 3 + 2] = temp / 100;
temp = temp1;
if (ans.a[i * 3 + 2] != 0)
{
ans.length = i * 3 + 3;
}
else if (ans.a[i * 3 + 1] != 0)
{
ans.length = i * 3 + 2;
}
else if (ans.a[i * 3] != 0)
{
ans.length = i * 3 + 1;
}
}
ans.pos = number1.pos + number2.pos - (number1.length - length1 + number2.length - length2);
ans.sign = number1.sign^number2.sign;
return ans;
}
else
{
int length3 = length1 + length2 + 4;
for (int i = 0; i < length3; i++)
{
ans.a[i] = 0;
}
ans.length = 0;
for (int i = 0; i < length1; i++)
{
int number = number1.a[number1.length - length1 + i];
for (int j = 0; j < length2; j++)
{
ans.a[i + j] += number2.a[number2.length - length2 + j] * number;
ans.length = max(ans.length, i + j + 1);
}
}
for (int k = 0; k < ans.length; k++)
{
if (ans.a[k] >= 10)
{
if (k == ans.length - 1)
{
ans.length++;
}
ans.a[k + 1] += ans.a[k] / 10;
ans.a[k] %= 10;
}
}
ans.pos = number1.pos + number2.pos - (number1.length - length1 + number2.length - length2);
ans.sign = number1.sign^number2.sign;
return ans;
}
}
Wfloat Multiplication(const Wfloat &number3, const Wfloat &number4)
{
Wfloat number1 = number3;
Wfloat number2 = number4;
Wfloat ans;
int eps = max(number1.length, number2.length);
int length1 = number1.length;
int length2 = number2.length;
if (length1 > length2)
{
int length3 = length1 - length2;
for (int i = 0; i < length3; i++)
{
number2.a[number2.length] = 0;
number2.length++;
}
}
else if (length1 < length2)
{
int length3 = length2 - length1;
for (int i = 0; i < length3; i++)
{
number1.a[number1.length] = 0;
number1.length++;
}
}
/*Print(number1);
Print(number2);
getchar();
getchar();*/
ans = Multiplication(number1, number2,eps);
int temp = ans.length - ans.pos;
for (int i = 0; i < temp; i++)
{
ans.a[i] = ans.a[ans.pos + i];
}
ans.length = temp;
ans.pos = 0;
//Print(ans);
return ans;
}
Wfloat Subtraction(const Wfloat &number3, const Wfloat &number4, int eps)
{
Wfloat number1 = number3;
Wfloat number2 = number4;
/*Print(number1);
Print(number2);*/
Wfloat ans;
int tempLength = eps;
//printf("%d\n", tempLength);
int length1 = number1.pos - number1.length;
int length2 = number2.pos - number2.length;
if (length1 > length2)
{
int length3 = length1 - length2;
for (int i = 0; i < length3; i++)
{
number1.a[number1.length] = 0;
number1.length++;
}
}
else if (length1 < length2)
{
int length3 = length2 - length1;
for (int i = 0; i < length3; i++)
{
number2.a[number2.length] = 0;
number2.length++;
}
}
int pos = 1;
for (int i = tempLength - 1; i >= 0; i--)
{
int b1, b2;
if (number1.length - pos < 0)
{
b1 = 0;
}
else
{
b1 = number1.a[number1.length - pos];
}
if (number2.length - pos < 0)
{
b2 = 0;
}
else
{
b2 = number2.a[number2.length - pos];
}
ans.a[i] = b1 - b2;
pos++;
}
int length = 0;
for (int i = 0; i < tempLength; i++)
{
if (ans.a[i] < 0)
{
ans.a[i + 1]--;
ans.a[i] += 10;
}
if (ans.a[i] != 0) length = i;
}
length++;
ans.length = length;
ans.pos = number1.pos - number1.length + tempLength;
/*Print(ans);
printf("\n");*/
return ans;
}
Wfloat Subtraction(const Wfloat &number1, const Wfloat &number2)
{
Wfloat ans;
int eps = max(number1.length, number2.length);
ans=Subtraction(number1, number2, eps);
int temp = ans.length - ans.pos;
for (int i = 0; i < temp; i++)
{
ans.a[i] = ans.a[ans.pos + i];
}
ans.length = temp;
ans.pos = 0;
return ans;
}
Wfloat Reciprocal(const Wfloat &N, int eps1)
{
Wfloat ans;
if (N.pos == 0)//整数
{
Wfloat x0, x1;
int N1length = N.length;
int Effectbit = 16;
Effectbit = min(Effectbit, N.length);
int lastZero = N1length - Effectbit;
double num = 0;
for (int i = 1; i <= Effectbit; i++)
{
num *= 10;
num += N.a[N1length - i];
}
x0.pos = lastZero;
num = pow(num, -1);
int tempNum;
int count1 = 0;
int flag = 0;
x0.length = 16;
int n1 = 15;
while (count1 < 16)
{
num *= 10;
tempNum = num;
if (tempNum != 0) num -= tempNum;
if (tempNum == 0 && flag == 0) x0.pos++;
else
{
flag = 1;
x0.a[n1] = tempNum;
n1--;
x0.pos++;
}
if (flag) count1++;
}
int lastbiteps = eps1;//最终需要的精度
int biteps = 32;
//Print(x0);
while (1)
{
//x1 = 2 * x0 - N*x0*x0;
x1 = Multiplication(N, x0, biteps);
x1 = Multiplication(x1, x0, biteps);
//Print(x1);
//Print(x0*2);
x0 = Subtraction(x0 * 2, x1, biteps);
if (biteps > lastbiteps) break;
biteps *= 2;
}
ans = x0;
return ans;
}
return ans;
}
Wfloat Division(const Wfloat &number1, const Wfloat &number2, int eps)
{
Wfloat ans;
return ans;
}
Wfloat Division(const Wfloat &number1, const Wfloat &number2)
{
Wfloat ans;
if (number1.length < number2.length)
{
return ans;
}
//Print(number2);
int eps = number1.length - number2.length + 50;
ans = Reciprocal(number2, eps);
//Print(ans);
ans = Multiplication(ans, number1, eps);
int temp = ans.length - ans.pos;
for (int i = 0; i < temp; i++)
{
ans.a[i] = ans.a[ans.pos + i];
}
ans.length = temp;
ans.pos = 0;
return ans;
}
Wfloat sqrt(Wfloat &N)
{
Wfloat ans;
if (N.sign == 1) return ans;
if (N.pos == 0)//整数
{
int eps = 50;
Wfloat x0, x1;
int N1length = N.length;
int Effectbit = 16;
Effectbit = min(Effectbit, N.length);
int lastZero = N1length - Effectbit;
double num = 0;
for (int i = 1; i <= Effectbit; i++)
{
num *= 10;
num += N.a[N1length - i];
}
if (lastZero % 2 != 0)
{
num *= 10;
Effectbit++;
lastZero--;
}
x0.pos = lastZero / 2;
num = pow(num, -0.5);
int tempNum;
int count1 = 0;
int flag = 0;
x0.length = 16;
int n1 = 15;
while (count1 < 16)
{
num *= 10;
tempNum = num;
if (tempNum != 0) num -= tempNum;
if (tempNum == 0 && flag == 0)x0.pos++;
else
{
flag = 1;
x0.a[n1] = tempNum;
n1--;
x0.pos++;
}
if (flag) count1++;
}
int lastbiteps = N.length / 2 + eps;//最终需要的精度
int biteps = 32;
x0.a[1]--;
if (x0.a[1] < 0)
{
for (int i = 1;; i++)
{
if (x0.a[i] < 0)
{
x0.a[i + 1]--;
x0.a[i] += 10;
}
else break;
}
}
while (1)
{
//x1 = x0 * (3 - N*x0 * x0) / 2;
x1 = Multiplication(N, x0, biteps);
x1 = Multiplication(x1, x0, biteps);
x1.a[x1.pos] = 3;
x1.length = x1.pos + 1;
for (int i = 0; i < x1.pos; i++)
{
x1.a[i] = -x1.a[i];
}
for (int i = 0; i < x1.pos; i++)
{
if (x1.a[i] < 0)
{
x1.a[i] += 10;
x1.a[i + 1]--;
}
}
x1.length = x1.pos + 1;
x1 = x1 / 2;
x1 = Multiplication(x0, x1, biteps);
x0 = x1;
if (biteps > lastbiteps) break;
biteps *= 2;
}
ans = Multiplication(x0, N, lastbiteps);
int temp = ans.length - ans.pos;
for (int i = 0; i < temp; i++)
{
ans.a[i] = ans.a[ans.pos + i];
}
ans.length = temp;
ans.pos = 0;
return ans;
}
}
double Reciprocal(double a)
{
double x0 = 0.2;
double x1;
while (1)
{
x1 = 2 * x0 - a*x0*x0;
if (abs(x0 - x1) < 1e-15) break;
printf("%.16lf\n", x1);
getchar();
swap(x0, x1);
}
return x1;
}
char word[100100];
int main()
{
Wfloat A, B;
scanf("%s", word);
A = word;
scanf("%s", word);
B = word;
Wfloat ans = Division(A, B);
Print(ans);
ans = Subtraction(A, Multiplication(ans, B));
Print(ans);
}