High-Precision Arithmetic
Design ideas and implementation methods for high-precision arithmetic
Primary author:
(Introductory paragraph) Even the unsigned long long integer type can only represent numbers up to (approximately 18 decimal digits). Obviously, this cannot handle ultra-large numerical computation problems such as those in astrophysics; this is a limitation of computer hardware data processing capability.
We need integer types that occupy more bytes.
We need to implement such high-precision arithmetic through programming.
How to Design It
We need a digital storage and computation method that can change according to actual data length.
Why don't we encounter such problems in daily life? Obviously, when dealing with larger numerical operations in daily life, we often use vertical calculation (even mental arithmetic uses this method).
The fundamental reason computers cannot process large numbers is that they cannot process all digits at once, which is the same reason we have to use vertical calculation for large numbers in daily life — we can't just glance at all digits and get the result.
But using vertical calculation allows us to only consider two digits of the large number at a time (current digit calculation and carry or borrow), and this approach also avoids the problem of computers having to process all digits at once.
Note that the vertical calculation method has strong scalability — as long as the paper is large enough to hold the planned maximum number, basic arithmetic operations of arbitrary length can be performed.
So we can think of using a block of memory composed of many bytes as the "paper" for our vertical calculations. Then arrays and strings, whose bytes can be uniformly managed, are very suitable tools.
We store each digit of the large integers to be calculated separately in sequence in array elements.
C++
std::string a = "1145141919810";
std::string b = "123456789";
// A freely scalable "piece of paper"
// with two lines for writing the two addends respectivelyNext, we can perform basic arithmetic operations using carry or borrow methods.
C++
// Addition
std::string add(std::string a, std::string b)
{
// First declare a function to wrap our functionality
// Function parameters are two strings storing the addends, return value is a string storing the result
// First reserve space on the paper for writing the answer
std::string ans = "";
// Intermediate digit-by-digit calculation
// Finally return the answer string
return ans;
}Now we need to think about how to calculate digit by digit.
C++
// Addition
std::string add(std::string a, std::string b)
{
// Reserve space for the answer
std::string ans = "";
// Calculate digit by digit from least significant to most significant
for (
// Start reverse traversal from position 0 (least significant digit)
int i = 0;
// Continue until traversing the most significant digit of the longer string
i < std::max(a.length(), b.length());
// Traverse one digit at a time
i ++
)
{
// Write the answer starting from one position beyond the least significant digit of the longer string
ans[std::max(a.length(), b.length())-1-i]
// The sum of the current digit is the sum of corresponding digits of both addends
= a[a.length()-1-i] + b[b.length()-1-i]
// Also need to consider character literal issues
- '0';
}
// Return answer
return ans;
}Well, some problems arise: what if the answer is longer than both addends, then an error will occur when carrying over to the most significant digit;
when calculating digits that the longer addend has beyond the shorter one, the shorter addend runs out of digits and needs to be padded with 0.
Moreover, this writing style is too bloated, like forcefully describing the human vertical calculation process in computer language, which makes it difficult to reflect the advantages of programming languages.
My first thought was reverse-order calculation (i.e., the most significant digit of the number is at the front of the string); the advantage of this is convenient carry handling — when carrying, you can just write directly to the back.
If output is needed, traverse the answer string in reverse order.
C++
std::string add(std::string a, std::string b)
// Pass in two reverse-order numbers
{
std::string ans = "";
a.push_back('0'), b.push_back('0');
// Pad with 0 to avoid boundary issues
char temp = '0';
// Temporarily store unit calculation results
for (int i = 0; i < std::max(a.length(), b.length()); i ++)
{
ans.push_back(
(a[std::max(a.length(), i)] + b[std::max(b.length(), i)] - '0') > '9' ?
(a[std::max(a.length(), i)] + b[std::max(b.length(), i)] - '0') % 10 + '0' :
(a[std::max(a.length(), i)] + b[std::max(b.length(), i)] - '0') + '0'
);
}
return ans;
// Return reverse-order answer
}Example Walkthrough
Engineering Applications
Algorithm Templates
Standalone functions, first a basic version that can handle leading zeros, decimal points, positive/negative signs, and memory allocation. (Been wanting to write something like this for ages 😋)
C++
char* add(char* a, char* b) // Addition
{
int sign_a = 1, sign_b = 1;
char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') { sign_a = 1; p_a++; }
char *dot_a = strchr(p_a, '.');
int int_len_a = dot_a ? (int)(dot_a - p_a) : (int)strlen(p_a);
int frac_len_a = dot_a ? (int)strlen(dot_a + 1) : 0;
char *int_start_a = p_a;
char *frac_start_a = dot_a ? dot_a + 1 : NULL;
int eff_int_len_a = (int_len_a == 0) ? 1 : int_len_a;
char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') { sign_b = 1; p_b++; }
char *dot_b = strchr(p_b, '.');
int int_len_b = dot_b ? (int)(dot_b - p_b) : (int)strlen(p_b);
int frac_len_b = dot_b ? (int)strlen(dot_b + 1) : 0;
char *int_start_b = p_b;
char *frac_start_b = dot_b ? dot_b + 1 : NULL;
int eff_int_len_b = (int_len_b == 0) ? 1 : int_len_b;
int max_int_len = (eff_int_len_a > eff_int_len_b) ? eff_int_len_a : eff_int_len_b;
int max_frac_len = (frac_len_a > frac_len_b) ? frac_len_a : frac_len_b;
int total_len = max_int_len + max_frac_len;
char *num1 = (char*)malloc(total_len + 1);
char *num2 = (char*)malloc(total_len + 1);
memset(num1, '0', total_len);
num1[total_len] = '\0';
memset(num2, '0', total_len);
num2[total_len] = '\0';
int offset1 = max_int_len - eff_int_len_a;
if (int_len_a > 0) {
memcpy(num1 + offset1, int_start_a, int_len_a);
}
if (frac_len_a > 0) {
memcpy(num1 + max_int_len, frac_start_a, frac_len_a);
}
int offset2 = max_int_len - eff_int_len_b;
if (int_len_b > 0) {
memcpy(num2 + offset2, int_start_b, int_len_b);
}
if (frac_len_b > 0) {
memcpy(num2 + max_int_len, frac_start_b, frac_len_b);
}
int cmp = strcmp(num1, num2);
int result_sign;
char *big, *small;
int do_add;
if (sign_a == sign_b) {
do_add = 1;
result_sign = sign_a;
big = num1;
small = num2;
} else {
if (cmp == 0) {
free(num1);
free(num2);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
} else if (cmp > 0) {
result_sign = sign_a;
big = num1;
small = num2;
} else {
result_sign = sign_b;
big = num2;
small = num1;
}
do_add = 0;
}
char *res = (char*)malloc(total_len + 2);
res[total_len + 1] = '\0';
if (do_add) {
int carry = 0;
for (int i = total_len - 1; i >= 0; i--) {
int sum = (big[i] - '0') + (small[i] - '0') + carry;
res[i + 1] = (char)((sum % 10) + '0');
carry = sum / 10;
}
res[0] = (char)(carry + '0');
} else {
int borrow = 0;
for (int i = total_len - 1; i >= 0; i--) {
int diff = (big[i] - '0') - (small[i] - '0') - borrow;
if (diff < 0) {
diff += 10;
borrow = 1;
} else {
borrow = 0;
}
res[i + 1] = (char)(diff + '0');
}
res[0] = '0';
}
free(num1);
free(num2);
int int_start = 0;
while (int_start < max_int_len && res[int_start] == '0') {
int_start++;
}
int int_part_len = max_int_len + 1 - int_start;
int frac_start_idx = max_int_len + 1;
int frac_end = total_len;
while (frac_end >= frac_start_idx && res[frac_end] == '0') {
frac_end--;
}
int frac_part_len = frac_end - frac_start_idx + 1;
if (frac_part_len < 0) frac_part_len = 0;
int is_zero = (int_start == max_int_len && res[max_int_len] == '0' && frac_part_len == 0);
if (is_zero) {
free(res);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
int out_len = 0;
if (result_sign == -1) out_len++;
out_len += int_part_len;
if (frac_part_len > 0) out_len += 1 + frac_part_len;
char *out = (char*)malloc(out_len + 1);
char *ptr = out;
if (result_sign == -1) *ptr++ = '-';
memcpy(ptr, res + int_start, int_part_len);
ptr += int_part_len;
if (frac_part_len > 0) {
*ptr++ = '.';
memcpy(ptr, res + frac_start_idx, frac_part_len);
ptr += frac_part_len;
}
*ptr = '\0';
free(res);
return out;
}
char* sub(char* a, char* b) {
char* neg_b;
if (b[0] == '-') {
neg_b = strdup(b + 1);
} else if (b[0] == '+') {
neg_b = (char*)malloc(strlen(b) + 2);
neg_b[0] = '-';
strcpy(neg_b + 1, b + 1);
} else {
neg_b = (char*)malloc(strlen(b) + 2);
neg_b[0] = '-';
strcpy(neg_b + 1, b);
}
char* res = add(a, neg_b); // Reuse
free(neg_b);
return res;
}
char* mul(char* a, char* b) // Multiplication
{
int sign_a = 1, sign_b = 1;
char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') { sign_a = 1; p_a++; }
char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') { sign_b = 1; p_b++; }
int final_sign = sign_a * sign_b;
int frac_a = 0;
int len_a = (int)strlen(p_a);
char *clean_a = (char*)malloc(len_a + 1);
int idx = 0, dot_seen = 0;
for (int i = 0; p_a[i]; i++) {
if (p_a[i] == '.') {
dot_seen = 1;
} else {
clean_a[idx++] = p_a[i];
if (dot_seen) frac_a++;
}
}
clean_a[idx] = '\0';
if (idx == 0) {
clean_a[0] = '0'; clean_a[1] = '\0'; frac_a = 0;
}
int frac_b = 0;
int len_b = (int)strlen(p_b);
char *clean_b = (char*)malloc(len_b + 1);
idx = 0; dot_seen = 0;
for (int i = 0; p_b[i]; i++) {
if (p_b[i] == '.') {
dot_seen = 1;
} else {
clean_b[idx++] = p_b[i];
if (dot_seen) frac_b++;
}
}
clean_b[idx] = '\0';
if (idx == 0) {
clean_b[0] = '0'; clean_b[1] = '\0'; frac_b = 0;
}
int total_frac = frac_a + frac_b;
int len_clean_b = (int)strlen(clean_b);
char *result_int = (char*)malloc(2);
result_int[0] = '0'; result_int[1] = '\0';
for (int i = len_clean_b - 1; i >= 0; i--) {
int digit = clean_b[i] - '0';
if (digit == 0) continue;
char *partial = (char*)malloc(2);
partial[0] = '0'; partial[1] = '\0';
for (int j = 0; j < digit; j++) {
char *temp = add(partial, clean_a);
free(partial);
partial = temp;
}
int shift = len_clean_b - 1 - i;
if (shift > 0) {
int partial_len = (int)strlen(partial);
char *shifted = (char*)malloc(partial_len + shift + 1);
strcpy(shifted, partial);
for (int k = 0; k < shift; k++)
shifted[partial_len + k] = '0';
shifted[partial_len + shift] = '\0';
free(partial);
partial = shifted;
}
char *temp = add(result_int, partial);
free(result_int);
free(partial);
result_int = temp;
}
free(clean_a);
free(clean_b);
if (strcmp(result_int, "0") == 0) {
free(result_int);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
int len_res = (int)strlen(result_int);
char *final_str;
if (total_frac == 0) {
int out_len = len_res + (final_sign == -1 ? 1 : 0);
final_str = (char*)malloc(out_len + 1);
char *ptr = final_str;
if (final_sign == -1) *ptr++ = '-';
strcpy(ptr, result_int);
} else if (len_res <= total_frac) {
int out_len = 2 + total_frac;
if (final_sign == -1) out_len++;
final_str = (char*)malloc(out_len + 1);
char *ptr = final_str;
if (final_sign == -1) *ptr++ = '-';
*ptr++ = '0';
*ptr++ = '.';
int zeros = total_frac - len_res;
for (int k = 0; k < zeros; k++) *ptr++ = '0';
strcpy(ptr, result_int);
} else {
int int_part_len = len_res - total_frac;
int out_len = len_res + 1;
if (final_sign == -1) out_len++;
final_str = (char*)malloc(out_len + 1);
char *ptr = final_str;
if (final_sign == -1) *ptr++ = '-';
memcpy(ptr, result_int, int_part_len);
ptr += int_part_len;
*ptr++ = '.';
strcpy(ptr, result_int + int_part_len);
char *end = final_str + out_len - 1;
while (*end == '0') end--;
if (*end == '.') end--;
*(end + 1) = '\0';
char *check = final_str;
if (*check == '-') check++;
if (strcmp(check, "0") == 0 || strcmp(check, "0.") == 0) {
free(final_str);
final_str = (char*)malloc(2);
final_str[0] = '0'; final_str[1] = '\0';
}
}
free(result_int);
return final_str;
}
char* div(char* a, char* b, int precision) // Division
{
if (!a || !b || precision < 0) return NULL;
int sign_a = 1, sign_b = 1;
char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') { sign_a = 1; p_a++; }
char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') { sign_b = 1; p_b++; }
int final_sign = sign_a * sign_b;
{
char *tmp = p_b;
int has_nonzero = 0;
while (*tmp) {
if (*tmp != '0' && *tmp != '.') { has_nonzero = 1; break; }
tmp++;
}
if (!has_nonzero) return NULL;
}
int dec_b = 0;
int b_len = (int)strlen(p_b);
char *divisor = (char*)malloc(b_len + 1);
int idx = 0, dot_seen = 0;
for (int i = 0; p_b[i]; i++) {
if (p_b[i] == '.') { dot_seen = 1; }
else {
divisor[idx++] = p_b[i];
if (dot_seen) dec_b++;
}
}
divisor[idx] = '\0';
if (idx == 0) { divisor[0] = '0'; divisor[1] = '\0'; }
{
char *s = divisor;
while (*s == '0' && *(s+1)) s++;
if (s != divisor) memmove(divisor, s, strlen(s) + 1);
}
char *int_a = NULL, *frac_a = NULL;
int int_len_a = 0, frac_len_a = 0;
char *dot_a = strchr(p_a, '.');
if (dot_a) {
int_len_a = (int)(dot_a - p_a);
frac_len_a = (int)strlen(dot_a + 1);
int_a = (char*)malloc(int_len_a + 1);
memcpy(int_a, p_a, int_len_a);
int_a[int_len_a] = '\0';
frac_a = (char*)malloc(frac_len_a + 1);
strcpy(frac_a, dot_a + 1);
} else {
int_len_a = (int)strlen(p_a);
int_a = (char*)malloc(int_len_a + 1);
strcpy(int_a, p_a);
frac_a = (char*)malloc(1);
frac_a[0] = '\0';
frac_len_a = 0;
}
int num_len = int_len_a + frac_len_a;
char *num_a = (char*)malloc(num_len + 1);
if (int_len_a > 0) memcpy(num_a, int_a, int_len_a);
if (frac_len_a > 0) memcpy(num_a + int_len_a, frac_a, frac_len_a);
num_a[num_len] = '\0';
char *dividend = NULL;
int new_dot_pos = int_len_a + dec_b;
if (new_dot_pos >= num_len) {
dividend = (char*)malloc(new_dot_pos + 1);
memcpy(dividend, num_a, num_len);
for (int i = num_len; i < new_dot_pos; i++) dividend[i] = '0';
dividend[new_dot_pos] = '\0';
} else {
int int_part_len = new_dot_pos;
int frac_part_len = num_len - new_dot_pos;
dividend = (char*)malloc(int_part_len + 1 + frac_part_len + 1);
memcpy(dividend, num_a, int_part_len);
dividend[int_part_len] = '.';
memcpy(dividend + int_part_len + 1, num_a + int_part_len, frac_part_len);
dividend[int_part_len + 1 + frac_part_len] = '\0';
}
free(num_a);
free(int_a);
free(frac_a);
char *int_part = NULL, *frac_part = NULL;
char *dot_divd = strchr(dividend, '.');
if (dot_divd) {
*dot_divd = '\0';
int_part = strdup(dividend);
frac_part = strdup(dot_divd + 1);
*dot_divd = '.';
} else {
int_part = strdup(dividend);
frac_part = strdup("");
}
{
char *s = int_part;
while (*s == '0' && *(s+1)) s++;
if (s != int_part) memmove(int_part, s, strlen(s) + 1);
}
if (strcmp(int_part, "0") == 0) {
int is_zero = 1;
for (char *p = frac_part; *p; p++) if (*p != '0') { is_zero = 0; break; }
if (is_zero) {
free(dividend); free(int_part); free(frac_part); free(divisor);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
}
char *rem = strdup("0");
int max_int_quot_len = (int)strlen(int_part) + 1;
char *int_quot = (char*)malloc(max_int_quot_len);
int_quot[0] = '\0';
int int_quot_len = 0;
for (int i = 0; int_part[i]; i++) {
char digit = int_part[i];
char *rem10 = (char*)malloc(strlen(rem) + 2);
strcpy(rem10, rem);
strcat(rem10, "0");
char digit_str[2] = {digit, '\0'};
char *cur = add(rem10, digit_str);
free(rem10);
int q = 0;
while (1) {
int len_cur = (int)strlen(cur), len_div = (int)strlen(divisor);
int cmp;
if (len_cur != len_div) cmp = len_cur > len_div ? 1 : -1;
else cmp = strcmp(cur, divisor);
if (cmp < 0) break;
char *tmp = sub(cur, divisor);
free(cur);
cur = tmp;
q++;
}
int_quot[int_quot_len++] = (char)(q + '0');
int_quot[int_quot_len] = '\0';
free(rem);
rem = cur;
}
if (int_quot_len == 0) {
int_quot[0] = '0'; int_quot[1] = '\0'; int_quot_len = 1;
}
int extra = (precision >= 0) ? 1 : 0;
int target_frac_len = precision + extra;
char *frac_quot = (char*)malloc(target_frac_len + 2);
frac_quot[0] = '\0';
int frac_quot_len = 0;
int idx_frac = 0;
int frac_len = (int)strlen(frac_part);
while (frac_quot_len < target_frac_len) {
char digit;
if (idx_frac < frac_len) {
digit = frac_part[idx_frac++];
} else {
if (strcmp(rem, "0") == 0) break;
digit = '0';
}
char *rem10 = (char*)malloc(strlen(rem) + 2);
strcpy(rem10, rem);
strcat(rem10, "0");
char digit_str[2] = {digit, '\0'};
char *cur = add(rem10, digit_str);
free(rem10);
int q = 0;
while (1) {
int len_cur = (int)strlen(cur), len_div = (int)strlen(divisor);
int cmp;
if (len_cur != len_div) cmp = len_cur > len_div ? 1 : -1;
else cmp = strcmp(cur, divisor);
if (cmp < 0) break;
char *tmp = sub(cur, divisor);
free(cur);
cur = tmp;
q++;
}
frac_quot[frac_quot_len++] = (char)(q + '0');
frac_quot[frac_quot_len] = '\0';
free(rem);
rem = cur;
if (strcmp(rem, "0") == 0 && idx_frac >= frac_len) break;
}
if (frac_quot_len > precision) {
int round_up = 0;
if (frac_quot_len >= precision + 1 && frac_quot[precision] >= '5') {
round_up = 1;
}
frac_quot[precision] = '\0';
frac_quot_len = precision;
if (round_up) {
for (int i = frac_quot_len - 1; i >= 0; i--) {
if (frac_quot[i] < '9') {
frac_quot[i]++;
round_up = 0;
break;
} else {
frac_quot[i] = '0';
}
}
if (round_up) {
char *carry = (char*)malloc(3);
carry[0] = '1'; carry[1] = '\0';
char *new_int = add(int_quot, carry);
free(int_quot);
int_quot = new_int;
free(carry);
}
}
}
while (frac_quot_len > 0 && frac_quot[frac_quot_len - 1] == '0') {
frac_quot[--frac_quot_len] = '\0';
}
free(rem);
free(int_part);
free(frac_part);
free(dividend);
free(divisor);
char *trimmed_int = int_quot;
while (*trimmed_int == '0' && *(trimmed_int+1)) trimmed_int++;
int is_zero_result = (strcmp(trimmed_int, "0") == 0 && frac_quot_len == 0);
if (is_zero_result) {
free(int_quot);
free(frac_quot);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
int out_len = (final_sign == -1 ? 1 : 0) + (int)strlen(trimmed_int);
if (frac_quot_len > 0) out_len += 1 + frac_quot_len;
char *out = (char*)malloc(out_len + 1);
char *ptr = out;
if (final_sign == -1) *ptr++ = '-';
strcpy(ptr, trimmed_int);
ptr += strlen(trimmed_int);
if (frac_quot_len > 0) {
*ptr++ = '.';
memcpy(ptr, frac_quot, frac_quot_len);
ptr += frac_quot_len;
}
*ptr = '\0';
free(int_quot);
free(frac_quot);
return out;
}
char* mod(char* a, char* b) // Modulo
{
if (!a || !b) return NULL;
int sign_a = 1, sign_b = 1;
char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') p_a++;
char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') p_b++;
int frac_a = 0, frac_b = 0;
char *clean_a = (char*)malloc(strlen(p_a) + 1);
int idx_a = 0, dot_a = 0;
for (int i = 0; p_a[i]; i++) {
if (p_a[i] == '.') dot_a = 1;
else {
clean_a[idx_a++] = p_a[i];
if (dot_a) frac_a++;
}
}
clean_a[idx_a] = '\0';
char *clean_b = (char*)malloc(strlen(p_b) + 1);
int idx_b = 0, dot_b = 0;
for (int i = 0; p_b[i]; i++) {
if (p_b[i] == '.') dot_b = 1;
else {
clean_b[idx_b++] = p_b[i];
if (dot_b) frac_b++;
}
}
clean_b[idx_b] = '\0';
int max_frac = (frac_a > frac_b) ? frac_a : frac_b;
int len_a_int = idx_a + (max_frac - frac_a);
int len_b_int = idx_b + (max_frac - frac_b);
char *int_a = (char*)malloc(len_a_int + 1);
strcpy(int_a, clean_a);
for (int i = idx_a; i < len_a_int; i++) int_a[i] = '0';
int_a[len_a_int] = '\0';
char *int_b = (char*)malloc(len_b_int + 1);
strcpy(int_b, clean_b);
for (int i = idx_b; i < len_b_int; i++) int_b[i] = '0';
int_b[len_b_int] = '\0';
free(clean_a);
free(clean_b);
char *a_int = int_a;
while (*a_int == '0' && *(a_int+1)) a_int++;
char *b_int = int_b;
while (*b_int == '0' && *(b_int+1)) b_int++;
if (strcmp(b_int, "0") == 0) {
free(int_a);
free(int_b);
return NULL;
}
char *rem = strdup("0");
int a_len = (int)strlen(a_int);
for (int i = 0; i < a_len; i++) {
int rem_len = (int)strlen(rem);
char *cur = (char*)malloc(rem_len + 2);
strcpy(cur, rem);
cur[rem_len] = a_int[i];
cur[rem_len + 1] = '\0';
char *cur_stripped = cur;
while (*cur_stripped == '0' && *(cur_stripped+1)) cur_stripped++;
if (cur_stripped != cur) {
memmove(cur, cur_stripped, strlen(cur_stripped) + 1);
}
while (1) {
int cmp = 0;
int cur_len = (int)strlen(cur);
int b_len = (int)strlen(b_int);
if (cur_len > b_len) cmp = 1;
else if (cur_len < b_len) cmp = -1;
else cmp = strcmp(cur, b_int);
if (cmp < 0) break;
char *tmp = sub(cur, b_int);
free(cur);
cur = tmp;
}
free(rem);
rem = cur;
}
free(int_a);
free(int_b);
char *rem_stripped = rem;
while (*rem_stripped == '0' && *(rem_stripped+1)) rem_stripped++;
if (rem_stripped != rem) {
memmove(rem, rem_stripped, strlen(rem_stripped) + 1);
}
if (strcmp(rem, "0") == 0) {
free(rem);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
char *result;
if (sign_a < 0) {
result = (char*)malloc(strlen(rem) + 2);
result[0] = '-';
strcpy(result + 1, rem);
} else {
result = strdup(rem);
}
free(rem);
return result;
}
char* pow(char* a, char* b, int precision) // Exponentiation
{
if (!a || !b) return NULL;
if (precision < 0) precision = 10;
int sign_a = 1;
const char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') p_a++;
int sign_b = 1;
const char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') p_b++;
int b_int_len = 0, b_frac_len = 0;
const char *dot = strchr(p_b, '.');
if (dot) {
b_int_len = (int)(dot - p_b);
b_frac_len = (int)strlen(dot + 1);
} else {
b_int_len = (int)strlen(p_b);
}
char *abs_b = (char*)malloc(b_int_len + b_frac_len + 2);
int idx = 0;
for (int i = 0; i < b_int_len; i++) abs_b[idx++] = p_b[i];
if (b_frac_len > 0) {
abs_b[idx++] = '.';
for (int i = 0; i < b_frac_len; i++) abs_b[idx++] = dot[i+1];
}
abs_b[idx] = '\0';
int b_is_integer = 1;
if (b_frac_len > 0) {
for (int i = 0; i < b_frac_len; i++)
if (dot[i+1] != '0') { b_is_integer = 0; break; }
}
int base_is_zero = 1;
for (const char *s = p_a; *s; s++) {
if (*s != '0' && *s != '.') { base_is_zero = 0; break; }
}
if (base_is_zero) {
if (strcmp(abs_b, "0") == 0 || (b_is_integer && strcmp(abs_b, "0") == 0)) {
free(abs_b);
return NULL;
}
if (sign_b < 0) {
free(abs_b);
return NULL;
}
free(abs_b);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
if (strcmp(abs_b, "0") == 0) {
free(abs_b);
char *one = (char*)malloc(2);
one[0] = '1'; one[1] = '\0';
return one;
}
if (b_is_integer) {
char *exp_str = strdup(p_b);
if (dot) exp_str[b_int_len] = '\0';
char *strip = exp_str;
while (*strip == '0' && *(strip+1)) strip++;
if (strip != exp_str) memmove(exp_str, strip, strlen(strip)+1);
char str_two[] = "2";
char *parity = mod(exp_str, str_two);
int exp_odd = (strcmp(parity, "1") == 0);
free(parity);
char *base = strdup(p_a);
char *exp_copy = strdup(exp_str);
char *result = strdup("1");
while (strcmp(exp_copy, "0") != 0) {
char *mod2 = mod(exp_copy, str_two);
int odd = (strcmp(mod2, "1") == 0);
free(mod2);
if (odd) {
char *tmp = mul(result, base);
free(result);
result = tmp;
}
char *tmp_base = mul(base, base);
free(base);
base = tmp_base;
char *new_exp;
if (strcmp(exp_copy, "0") == 0) {
new_exp = strdup("0");
} else {
int len = strlen(exp_copy);
char *half = (char*)malloc(len + 1);
int carry = 0, idx_h = 0;
for (int i = 0; i < len; i++) {
int d = exp_copy[i] - '0' + carry * 10;
int q = d / 2;
carry = d % 2;
if (q > 0 || idx_h > 0) half[idx_h++] = q + '0';
}
if (idx_h == 0) half[idx_h++] = '0';
half[idx_h] = '\0';
new_exp = half;
}
free(exp_copy);
exp_copy = new_exp;
}
free(base);
free(exp_copy);
free(exp_str);
if (sign_b < 0) {
char one[] = "1";
char *inv = div(one, result, precision);
free(result);
result = inv;
}
if (sign_a < 0 && exp_odd) {
int is_zero = 1;
for (char *p = result; *p; p++)
if (*p != '0' && *p != '.') { is_zero = 0; break; }
if (!is_zero) {
char *signed_res = (char*)malloc(strlen(result) + 2);
signed_res[0] = '-';
strcpy(signed_res+1, result);
free(result);
result = signed_res;
}
}
free(abs_b);
return result;
}
if (sign_a < 0) {
free(abs_b);
return NULL;
}
auto exp_func = [&](char *x) -> char* {
char *sum = strdup("1");
char *term = strdup("1");
char *x_copy = strdup(x);
char k_buf[20];
for (int k = 1; k <= 150; k++) {
sprintf(k_buf, "%d", k);
char *tmp1 = mul(term, x_copy);
char *tmp2 = div(tmp1, k_buf, precision);
free(term); free(tmp1);
term = tmp2;
char *tmp3 = add(sum, term);
free(sum);
sum = tmp3;
int is_zero = 1;
for (char *p = term; *p; p++) {
if (*p != '0' && *p != '.' && *p != '-') { is_zero = 0; break; }
}
if (is_zero) break;
}
free(term);
free(x_copy);
return sum;
};
auto ln_func = [&](char *x) -> char* {
char *y = strdup("0");
char two[] = "2";
for (int iter = 0; iter < 12; iter++) {
char *ey = exp_func(y);
char *diff = sub(x, ey);
char *sum_xe = add(x, ey);
char *term = mul(two, diff);
char *delta = div(term, sum_xe, precision);
free(term); free(diff); free(sum_xe); free(ey);
char *new_y = add(y, delta);
free(y); free(delta);
y = new_y;
}
return y;
};
char *a_abs = strdup(p_a);
char *ln_a = ln_func(a_abs);
free(a_abs);
char *b_signed;
if (sign_b < 0) {
b_signed = (char*)malloc(strlen(abs_b) + 2);
b_signed[0] = '-';
strcpy(b_signed + 1, abs_b);
} else {
b_signed = strdup(abs_b);
}
free(abs_b);
char *power = mul(b_signed, ln_a);
free(b_signed);
free(ln_a);
char *result = exp_func(power);
free(power);
if (result[0] == '-') {
int is_zero = 1;
for (int i = 1; result[i]; i++)
if (result[i] != '0' && result[i] != '.') { is_zero = 0; break; }
if (is_zero) {
free(result);
result = strdup("0");
}
}
return result;
}
char* sqrt(char* a, char* b, int precision) // Square root / Nth root
{
if (!a || !b || precision < 0) return NULL;
int sign_a = 1;
const char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') p_a++;
int sign_b = 1;
const char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') p_b++;
const char *dot_b = strchr(p_b, '.');
int len_b = dot_b ? (int)(dot_b - p_b) : (int)strlen(p_b);
char *exp_str = (char*)malloc(len_b + 1);
memcpy(exp_str, p_b, len_b);
exp_str[len_b] = '\0';
if (dot_b) {
for (const char *s = dot_b + 1; *s; s++)
if (*s != '0') {
free(exp_str);
return NULL;
}
}
char *strip = exp_str;
while (*strip == '0' && *(strip + 1)) strip++;
if (strip != exp_str) memmove(exp_str, strip, strlen(strip) + 1);
if (strcmp(exp_str, "0") == 0) {
free(exp_str);
return NULL;
}
int a_is_zero = 1;
for (const char *s = p_a; *s; s++)
if (*s != '0' && *s != '.') { a_is_zero = 0; break; }
if (a_is_zero) {
free(exp_str);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
if (sign_a < 0) {
char two_str[] = "2";
char *rem = mod(exp_str, two_str);
int even = (strcmp(rem, "0") == 0);
free(rem);
if (even) {
free(exp_str);
return NULL;
}
}
int guard = 10;
int work_prec = precision + guard;
char *a_abs = strdup(p_a);
char *n_str = strdup(exp_str);
char *n_minus_one = sub(n_str, (char*)"1");
int need_reciprocal = (sign_b < 0);
char *x = strdup(a_abs);
if (strcmp(x, "0") == 0) {
free(x); free(a_abs); free(n_str); free(n_minus_one); free(exp_str);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
for (int iter = 0; iter < 200; iter++) {
char *x_pow = strdup("1");
char *exp_copy = strdup(n_minus_one);
char *base_pow = strdup(x);
while (strcmp(exp_copy, "0") != 0) {
char *mod2 = mod(exp_copy, (char*)"2");
int odd = (strcmp(mod2, "1") == 0);
free(mod2);
if (odd) {
char *tmp = mul(x_pow, base_pow);
free(x_pow);
x_pow = tmp;
}
char *tmp_base = mul(base_pow, base_pow);
free(base_pow);
base_pow = tmp_base;
int len = strlen(exp_copy);
char *half = (char*)malloc(len + 1);
int carry = 0, idx_h = 0;
for (int i = 0; i < len; i++) {
int d = exp_copy[i] - '0' + carry * 10;
int q = d / 2;
carry = d % 2;
if (q > 0 || idx_h > 0) half[idx_h++] = q + '0';
}
if (idx_h == 0) half[idx_h++] = '0';
half[idx_h] = '\0';
free(exp_copy);
exp_copy = half;
}
free(exp_copy);
free(base_pow);
char *term1 = mul(n_minus_one, x);
char *term2 = div(a_abs, x_pow, work_prec);
free(x_pow);
char *numerator = add(term1, term2);
free(term1); free(term2);
char *new_x = div(numerator, n_str, work_prec);
free(numerator);
char *diff = sub(new_x, x);
char *diff_abs = (diff[0] == '-') ? diff + 1 : diff;
int th_len = 2 + (work_prec - 1);
char *threshold = (char*)malloc(th_len + 1);
threshold[0] = '0'; threshold[1] = '.';
for (int i = 0; i < work_prec - 1; i++) threshold[2 + i] = '0';
threshold[th_len - 1] = '1';
threshold[th_len] = '\0';
char *check = sub(threshold, diff_abs);
int converged = (check[0] != '-');
free(check);
free(threshold);
free(diff);
if (converged) {
free(x);
x = new_x;
break;
}
free(x);
x = new_x;
}
char *pow10 = strdup("1");
for (int i = 0; i < precision; i++) {
char *tmp = mul(pow10, (char*)"10");
free(pow10);
pow10 = tmp;
}
char *scaled = mul(x, pow10);
char *half_str = strdup("0.5");
char *scaled_plus = add(scaled, half_str);
free(scaled); free(half_str);
char *rounded_scaled;
{
char *dot = strchr(scaled_plus, '.');
if (dot) *dot = '\0';
rounded_scaled = strdup(scaled_plus);
if (dot) *dot = '.';
}
free(scaled_plus);
char *result = div(rounded_scaled, pow10, precision);
free(rounded_scaled);
free(pow10);
if (sign_a < 0) {
char *neg_result = (char*)malloc(strlen(result) + 2);
neg_result[0] = '-';
strcpy(neg_result + 1, result);
free(result);
result = neg_result;
}
if (need_reciprocal) {
char *one = strdup("1");
char *inv = div(one, result, precision);
free(one);
free(result);
result = inv;
}
free(x);
free(a_abs);
free(n_str);
free(n_minus_one);
free(exp_str);
return result;
}
char* log(char* a, char* b, int precision) // Logarithm
{
if (!a || !b || precision < 0) return NULL;
int sign_a = 1;
const char *p_a = a;
if (*p_a == '-') { sign_a = -1; p_a++; }
else if (*p_a == '+') p_a++;
int a_is_zero = 1;
for (const char *s = p_a; *s; s++)
if (*s != '0' && *s != '.') { a_is_zero = 0; break; }
if (a_is_zero || sign_a < 0) return NULL;
{
char *one = strdup("1");
char *a_copy = strdup(p_a);
int alen = strlen(a_copy);
while (alen > 1 && a_copy[alen-1] == '0') { a_copy[--alen] = '\0'; }
if (alen > 1 && a_copy[alen-1] == '.') a_copy[--alen] = '\0';
if (strcmp(a_copy, "1") == 0) {
free(one); free(a_copy);
return NULL;
}
free(one); free(a_copy);
}
int sign_b = 1;
const char *p_b = b;
if (*p_b == '-') { sign_b = -1; p_b++; }
else if (*p_b == '+') p_b++;
int b_is_zero = 1;
for (const char *s = p_b; *s; s++)
if (*s != '0' && *s != '.') { b_is_zero = 0; break; }
if (b_is_zero || sign_b < 0) return NULL;
{
char *b_copy = strdup(p_b);
int blen = strlen(b_copy);
while (blen > 1 && b_copy[blen-1] == '0') { b_copy[--blen] = '\0'; }
if (blen > 1 && b_copy[blen-1] == '.') b_copy[--blen] = '\0';
if (strcmp(b_copy, "1") == 0) {
free(b_copy);
char *zero = (char*)malloc(2);
zero[0] = '0'; zero[1] = '\0';
return zero;
}
free(b_copy);
}
auto exp_func = [&](char *x) -> char* {
char *sum = strdup("1");
char *term = strdup("1");
char *x_copy = strdup(x);
char k_buf[20];
int max_iters = 150;
for (int k = 1; k <= max_iters; k++) {
sprintf(k_buf, "%d", k);
char *tmp1 = mul(term, x_copy);
char *tmp2 = div(tmp1, k_buf, precision + 10);
free(term); free(tmp1);
term = tmp2;
char *tmp3 = add(sum, term);
free(sum);
sum = tmp3;
int is_zero = 1;
for (char *p = term; *p; p++) {
if (*p != '0' && *p != '.' && *p != '-') { is_zero = 0; break; }
}
if (is_zero) break;
}
free(term);
free(x_copy);
return sum;
};
auto ln_func = [&](char *x) -> char* {
char *y = strdup("0");
char two[] = "2";
int max_iters = 12;
for (int iter = 0; iter < max_iters; iter++) {
char *ey = exp_func(y);
char *diff = sub(x, ey);
char *sum_xe = add(x, ey);
char *term = mul(two, diff);
char *delta = div(term, sum_xe, precision + 10);
free(term); free(diff); free(sum_xe); free(ey);
char *new_y = add(y, delta);
free(y); free(delta);
y = new_y;
}
return y;
};
char *ln_b = ln_func(strdup(p_b));
char *ln_a = ln_func(strdup(p_a));
char *result = div(ln_b, ln_a, precision);
free(ln_b);
free(ln_a);
return result;
}
Wrapped in a class, usable like normal floating-point numbers, supports direct assignment and computation. (Been wanting to write this for ages too 😋)
C++
class HighPrecision
{
private:
public:
};