267 lines
7.5 KiB
C
267 lines
7.5 KiB
C
/***************************************************************************
|
|
* __________ __ ___.
|
|
* Open \______ \ ____ ____ | | _\_ |__ _______ ___
|
|
* Source | _// _ \_/ ___\| |/ /| __ \ / _ \ \/ /
|
|
* Jukebox | | ( <_> ) \___| < | \_\ ( <_> > < <
|
|
* Firmware |____|_ /\____/ \___ >__|_ \|___ /\____/__/\_ \
|
|
* \/ \/ \/ \/ \/
|
|
* $Id$
|
|
*
|
|
* Copyright (C) 2018 by Michael A. Sevakis
|
|
*
|
|
* This program is free software; you can redistribute it and/or
|
|
* modify it under the terms of the GNU General Public License
|
|
* as published by the Free Software Foundation; either version 2
|
|
* of the License, or (at your option) any later version.
|
|
*
|
|
* This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY
|
|
* KIND, either express or implied.
|
|
*
|
|
****************************************************************************/
|
|
#include "ap_int.h"
|
|
#include "fixedpoint.h"
|
|
|
|
/* Miscellaneous large-sized integer functions */
|
|
|
|
/* round string, base 10 */
|
|
bool round_number_string10(char *p_rdig, long len)
|
|
{
|
|
/*
|
|
* * p should point to the digit that determines if rounding should occur
|
|
* * buffer is updated in reverse
|
|
* * an additional '1' may be added to the beginning: eg. 9.9 => 10.0
|
|
*/
|
|
#if 1 /* nearest */
|
|
bool round = p_rdig[0] >= '5';
|
|
#else /* even */
|
|
bool round = p_rdig[0] >= '5' && (p_rdig[-1] & 1);
|
|
#endif
|
|
|
|
while (round && len-- > 0) {
|
|
int d = *--p_rdig;
|
|
round = ++d > '9';
|
|
*p_rdig = round ? '0' : d;
|
|
}
|
|
|
|
if (round) {
|
|
/* carry to the next place */
|
|
*--p_rdig = '1';
|
|
}
|
|
|
|
return round;
|
|
}
|
|
|
|
/* format arbitrary-precision base 10 integer */
|
|
char * format_ap_int10(struct ap_int *a,
|
|
char *p_end)
|
|
{
|
|
/*
|
|
* * chunks are in least-to-most-significant order
|
|
* * chunk array is used for intermediate division results
|
|
* * digit string buffer is written high-to-low address order
|
|
*/
|
|
long numchunks = a->numchunks;
|
|
char *p = p_end;
|
|
|
|
if (numchunks == 0) {
|
|
/* fast formatting */
|
|
uint64_t val = a->val;
|
|
|
|
do {
|
|
*--p = val % 10 + '0';
|
|
val /= 10;
|
|
} while (val);
|
|
|
|
a->len = p_end - p;
|
|
return p;
|
|
}
|
|
|
|
uint32_t *chunks = a->chunks;
|
|
long topchunk = numchunks - 1;
|
|
|
|
/* if top chunk(s) are zero, ignore */
|
|
while (topchunk >= 0 && chunks[topchunk] == 0) {
|
|
topchunk--;
|
|
}
|
|
|
|
/* optimized to divide number by the biggest 10^x a uint32_t can hold
|
|
so that r_part holds the remainder (x % 1000000000) at the end of
|
|
the division */
|
|
do {
|
|
uint64_t r_part = 0;
|
|
|
|
for (long i = topchunk; i >= 0; i--) {
|
|
/*
|
|
* Testing showed 29 bits as a sweet spot:
|
|
* * Is a 32-bit constant (good for 32-bit hardware)
|
|
* * No more normalization is required than with 30 and 31
|
|
* (32 bits requires the least but also a large constant)
|
|
* * Doesn't need to be reduced before hand by subtracting the
|
|
* divisor in order to keep it 32-bits which obviates the need
|
|
* to correct with another term of the remainder after
|
|
* multiplying
|
|
*
|
|
* 2305843009 = floor(ldexp(1, 29) / 1000000000.0 * ldexp(1, 32))
|
|
*/
|
|
static const unsigned long c = 2305843009; /* .213693952 */
|
|
uint64_t q_part = r_part*c >> 29;
|
|
r_part = (r_part << 32) | chunks[i];
|
|
r_part -= q_part*1000000000;
|
|
|
|
/* if remainder is still out of modular range, normalize it
|
|
and carry over into quotient */
|
|
while (r_part >= 1000000000) {
|
|
r_part -= 1000000000;
|
|
q_part++;
|
|
}
|
|
|
|
chunks[i] = q_part;
|
|
}
|
|
|
|
/* if top chunk(s) became zero, ignore from now on */
|
|
while (topchunk >= 0 && chunks[topchunk] == 0) {
|
|
topchunk--;
|
|
}
|
|
|
|
/* format each digit chunk, padded to width 9 if not the leading one */
|
|
uint32_t val = r_part;
|
|
int len = 8*(topchunk >= 0);
|
|
|
|
while (len-- >= 0 || val) {
|
|
*--p = (val % 10) + '0';
|
|
val /= 10;
|
|
}
|
|
} while (topchunk >= 0);
|
|
|
|
a->len = p_end - p;
|
|
return p;
|
|
}
|
|
|
|
/* format arbitrary-precision base 10 fraction */
|
|
char * format_ap_frac10(struct ap_int *a,
|
|
char *p_start,
|
|
long precision)
|
|
{
|
|
/*
|
|
* * chunks are in least-to-most-significant order
|
|
* * chunk array is used for intermediate multiplication results
|
|
* * digit string buffer is written low-to-high address order
|
|
* * high bit of fraction must be left-justified to a chunk
|
|
* boundary
|
|
*/
|
|
long numchunks = a->numchunks;
|
|
bool trimlz = precision < 0;
|
|
char *p = p_start;
|
|
|
|
if (trimlz) {
|
|
/* trim leading zeros and provide <precision> digits; a->len
|
|
will end up greater than the specified precision unless the
|
|
value is zero */
|
|
precision = -precision;
|
|
}
|
|
|
|
a->len = precision;
|
|
|
|
if (numchunks == 0) {
|
|
/* fast formatting; shift must be <= 60 as top four bits are used
|
|
for digit carryout */
|
|
if (trimlz && !a->val) {
|
|
/* value is zero */
|
|
trimlz = false;
|
|
}
|
|
|
|
uint64_t val = a->val << (60 - a->shift);
|
|
|
|
while (precision > 0) {
|
|
val *= 10;
|
|
uint32_t c_part = val >> 60;
|
|
|
|
if (trimlz) {
|
|
if (!c_part) {
|
|
a->len++;
|
|
continue;
|
|
}
|
|
|
|
trimlz = false;
|
|
}
|
|
|
|
*p++ = c_part + '0';
|
|
val ^= (uint64_t)c_part << 60;
|
|
precision--;
|
|
}
|
|
|
|
return p;
|
|
}
|
|
|
|
uint32_t *chunks = a->chunks;
|
|
long bottomchunk = 0, topchunk = numchunks;
|
|
|
|
while (topchunk > 0 && chunks[topchunk - 1] == 0) {
|
|
topchunk--;
|
|
}
|
|
|
|
/* optimized to multiply number by the biggest 10^x a uint32_t can hold
|
|
so that c_part holds the carryover into the integer part at the end
|
|
of the multiplication */
|
|
while (precision > 0) {
|
|
/* if bottom chunk(s) are or became zero, skip them */
|
|
while (bottomchunk < numchunks && chunks[bottomchunk] == 0) {
|
|
bottomchunk++;
|
|
}
|
|
|
|
uint32_t c_part = 0;
|
|
|
|
for (long i = bottomchunk; i < topchunk; i++) {
|
|
uint64_t p_part = chunks[i];
|
|
|
|
p_part = p_part * 1000000000 + c_part;
|
|
c_part = p_part >> 32;
|
|
|
|
chunks[i] = p_part;
|
|
}
|
|
|
|
if (topchunk < numchunks && c_part) {
|
|
chunks[topchunk++] = c_part;
|
|
c_part = 0;
|
|
}
|
|
|
|
int len = 9;
|
|
|
|
if (trimlz && bottomchunk < numchunks) {
|
|
if (!c_part) {
|
|
a->len += 9;
|
|
continue;
|
|
}
|
|
|
|
/* first non-zero chunk has leading zeros? */
|
|
for (uint32_t val = c_part; val < 100000000; val *= 10) {
|
|
len--;
|
|
}
|
|
|
|
a->len += 9 - len;
|
|
trimlz = false;
|
|
}
|
|
|
|
/* format each digit chunk, padded to width 9 if not exceeding
|
|
precision */
|
|
precision -= len;
|
|
|
|
if (precision < 0) {
|
|
/* remove extra digits */
|
|
c_part /= ipow(10, -precision);
|
|
len += precision;
|
|
}
|
|
|
|
p += len;
|
|
|
|
char *p2 = p;
|
|
|
|
while (len-- > 0) {
|
|
*--p2 = (c_part % 10) + '0';
|
|
c_part /= 10;
|
|
}
|
|
}
|
|
|
|
return p;
|
|
}
|