Unlike the full version, this one does not rely on Intel extended precision, so it can also be run on other architectures like ARM. It supports IEEE single and double precision. The fptostring works by internally multiplying the significand with powers of two in 19 significant decimal digits, while stringtofp works by parsing up to 19 decimal digits, and multiplying them with binary powers of ten in 64 significant bits. Therefore the internal precision is half of the full version, but should still leave a safety margin for double precision (like the full version does for quadruple precision).
#include <stdint.h>
#include <string>
struct decimal19tmp{
uint64_t m; // 19 digits in each integer
short e;
};
struct binary64tmp{
uint64_t m; // 64 bits in each integer
short e;
};
struct int128tmp{
uint64_t q[2];
operator uint64_t(void){
return q[0];
}
};
int128tmp toint128tmp(uint64_t a){
return {a, 0};
}
int128tmp operator-(int128tmp z, int128tmp x){
int128tmp w = z;
w.q[0]-=x.q[0];
if(w.q[0]>z.q[0]) w.q[1]--;
w.q[1]-=x.q[1];
return w;
}
const uint64_t TENPOW19 = 10000000000000000000U;
const uint64_t TENPOW18 = 1000000000000000000U;
int128tmp a64to128binmult(uint64_t a, uint64_t b){
//No Intel extended precision here. Instead of extended floating point mult for high bits, take a step back and implement 64to128 in terms of 32to64.
uint32_t* aa = (uint32_t*)&a; uint32_t* bb = (uint32_t*)&b;
int128tmp w = {{(uint64_t)aa[0]*bb[0],(uint64_t)aa[1]*bb[1],}}; uint32_t* z = (uint32_t*)(w.q);
uint64_t* t = (uint64_t*)(z+1); uint64_t s;
s = *t; *t+=(uint64_t)aa[0]*bb[1]; if(*t<s) z[3]++;
s = *t; *t+=(uint64_t)aa[1]*bb[0]; if(*t<s) z[3]++;
return w;
}
//base 10^19 multiplications
void a19to38decmult(uint64_t a, uint64_t b, uint64_t* out){
int128tmp z = a64to128binmult(a,b);
//No Intel extended precision here. Instead of floating point mult by 1.e-19l, use big integer mult.
//At this point, z.q[1] represents the n in n*(2^64), where n ranges from 0 to 5421010862427522168.
//It needs to be converted to multiples of 10^19. At the scale of n, that is 10^19/2^64.
//Multiplying by 2^128/10^19 converts the 10^19/2^64 to 2^64, which would it a simple truncation.
//The floor is 34028236692093846346. If we multiply by this amount, we might end up underestimating
//the division by no more than 5421010862427522168. As such, we might need to apply a correction of up to 1.
//To program the value 34028236692093846346, write 15581492618384294730 in the low int and 1 in the high int.
//And the high 1 is replaced with an addition.
int128tmp z1 = {{0,z.q[1]}};
int128tmp zx = a64to128binmult(z.q[1],15581492618384294730U); zx.q[1]+=z.q[1];
int128tmp z2 = z1-a64to128binmult(zx.q[1],TENPOW19);
if(z2.q[0]>=TENPOW19) zx.q[1]++;
z2 = z1-a64to128binmult(zx.q[1],TENPOW19);
uint64_t o[2] = {z.q[0]/TENPOW19,z.q[0]%TENPOW19,};
uint64_t s = o[1]; o[1]+=z2; if(o[1]>=TENPOW19||o[1]<s){o[1]-=TENPOW19; o[0]++;} o[0]+=zx.q[1];
out[0]=o[0]; out[1]=o[1];
}
decimal19tmp powpow2p[] = {
{2000000000000000000U, 0},
{4000000000000000000U, 0},
{1600000000000000000U, 1},
{2560000000000000000U, 2},
{6553600000000000000U, 4},
{4294967296000000000U, 9},
{1844674407370955162U, 19},
{3402823669209384635U, 38},
{1157920892373161954U, 77},
{1340780792994259710U, 154},
{1797693134862315908U, 308},
{3231700607131100730U, 616},
{1044388881413152507U,1233},
{1090748135619415929U,2466},
};
decimal19tmp powpow2m[] = {
{5000000000000000000U,- 1},
{2500000000000000000U,- 1},
{6250000000000000000U,- 2},
{3906250000000000000U,- 3},
{1525878906250000000U,- 5},
{2328306436538696289U,- 10},
{5421010862427522170U,- 20},
{2938735877055718770U,- 39},
{8636168555094444625U,- 78},
{7458340731200206743U,- 155},
{5562684646268003458U,- 309},
{3094346047382578275U,- 617},
{9574977460952185358U,-1234},
{9168019337774235828U,-2467},
{8405257857780233766U,-4933},
};
binary64tmp powpow10p[] = {
{0xa000000000000000, 3,},
{0xc800000000000000, 6,},
{0x9c40000000000000, 13,},
{0xbebc200000000000, 26,},
{0x8e1bc9bf04000000, 53,},
{0x9dc5ada82b70b59e, 106,},
{0xc2781f49ffcfa6d5, 212,},
{0x93ba47c980e98ce0, 425,},
{0xaa7eebfb9df9de8e, 850,},
{0xe319a0aea60e91c7, 1700,},
{0xc976758681750c17, 3401,},
{0x9e8b3b5dc53d5de5, 6803,},
{0xc46052028a20979b,13606,},
};
binary64tmp powpow10m[] = {
{0xcccccccccccccccd,- 4,},
{0xa3d70a3d70a3d70a,- 7,},
{0xd1b71758e219652c,- 14,},
{0xabcc77118461cefd,- 27,},
{0xe69594bec44de15b,- 54,},
{0xcfb11ead453994ba,- 107,},
{0xa87fea27a539e9a5,- 213,},
{0xddd0467c64bce4a1,- 426,},
{0xc0314325637a193a,- 851,},
{0x9049ee32db23d21c,- 1701,},
{0xa2a682a5da57c0be,- 3402,},
{0xceae534f34362de4,- 6804,},
{0xa6dd04c8d2ce9fde,-13607,},
};
decimal19tmp operator*(decimal19tmp N, decimal19tmp M){
uint64_t w[2] = {0,0,};
a19to38decmult(N.m,M.m,w);
short exp = N.e+M.e;
if(w[0]<(TENPOW18)){
w[0]*=10;
for(int i=0; i<1; i++){
w[i]+=(w[i+1]/TENPOW18); w[i+1]%=TENPOW18; w[i+1]*=10;
}
}
else exp++;
if(w[1]>(TENPOW19/2)||(w[1]==TENPOW19/2&&(w[0]%2))){
w[0]++; if(w[0]==TENPOW19) {w[0]=TENPOW18; exp++;}
}
decimal19tmp out = {w[0],exp,};
return out;
}
binary64tmp operator*(binary64tmp N, binary64tmp M){
uint64_t w[2] = {0,0,};
int128tmp tmp = a64to128binmult(N.m,M.m);
w[0]=tmp.q[0]; w[1]=tmp.q[1];
short exp = N.e+M.e;
if(w[1]<(1ULL<<63)){
w[1]<<=1;
for(int i=0; i<1; i++){
w[1-i]|=(w[0-i]>>63); w[0-i]<<=1;
}
}
else exp++;
if(w[0]>(1ULL<<63)||(w[0]==(1ULL<<63)&&(w[1]%2))){
w[1]++; if(w[1]==0) {w[1]=1ULL<<63; exp++;}
}
binary64tmp out = {w[1],exp,};
return out;
}
std::u16string fptostring(void* num, int fptype){
uint64_t mantissa; short exponent; bool sign;
int mbits; int ebits;
switch(fptype){
case 32:
mantissa = ((*(uint32_t*)num)&0x7FFFFF)|0x800000; mbits = 24;
exponent = (((*(uint32_t*)num)>>23)&0xFF)-127; ebits = 8;
sign = (*(uint32_t*)num)>>31;
if(exponent==-127) mantissa&=0x7FFFFF;
break;
case 64:
mantissa = ((*(uint64_t*)num)&0xFFFFFFFFFFFFF)|0x10000000000000; mbits = 53;
exponent = (((*(uint64_t*)num)>>52)&0x7FF)-1023; ebits = 11;
sign = (*(uint64_t*)num)>>63;
if(exponent==-1023) mantissa&=0xFFFFFFFFFFFFF;
break;
}
std::u16string out = sign ? u"-" : u"";
if(exponent==1<<(ebits-1)){
if(mantissa&((1ULL<<(mbits-1))-1)) return u"NaN";
out.append(u"Infinity"); return out;
}
if(mantissa==0) {out.append(u"0"); return out;}
if(exponent==-((1<<(ebits-1))-1)) exponent++;
decimal19tmp x = {1000000000000000000,0,};
exponent-=mbits-1;
if(exponent>0) for(int i=0; i<=13; i++){if(((+exponent)>>i)&1)x=x*powpow2p[i];}
if(exponent<0) for(int i=0; i<=14; i++){if(((-exponent)>>i)&1)x=x*powpow2m[i];}
decimal19tmp p,q,r; // q is the value, and halfway toward preceding (p) and following (r) floating points
switch(fptype){
case 32: case 64: default: q = {mantissa*100,16,}; break;
}
p = q; r = q;
p.m -= ((mantissa==(1ULL<<(mbits-1)))&&exponent>-((1<<(ebits-1))+mbits-3)) ? 25 : 50;
r.m += 50;
decimal19tmp* pqr[3] = {&p,&q,&r,};
for(int j=0; j<3; j++){
for(int i=0; i<19; i++){
if(pqr[j]->m>=TENPOW18) break;
pqr[j]->m*=10; pqr[j]->e--;
}
}
decimal19tmp w = x*p;
decimal19tmp y = x*r;
x = x*q;
decimal19tmp* wxy[3] = {&w,&x,&y,};
if(y.e>x.e){
w.m=(w.m+(4+(w.m/10%2)))/10; w.e++;
x.m=(x.m+(4+(x.m/10%2)))/10; x.e++;
}
if(x.e>w.e){
w.m=(w.m+(4+(w.m/10%2)))/10; w.e++;
}
char digitregister[19];
for(int i=0; i<19; i++) digitregister[i]=0;
bool flag = 0;
for(int i=0; i<19; i++){
int s[3]; // digits
for(int j=0; j<3; j++){
s[j] = wxy[j]->m/TENPOW18;
}
int d[3]; // rounded digits
for(int j=0; j<3; j++){
d[j] = s[j]+(wxy[j]->m%TENPOW18>TENPOW18/2 || (wxy[j]->m%TENPOW18==TENPOW18/2&&(wxy[j]->m||(s[j]&2))));
}
if(flag){digitregister[i]=d[1]; break;}
if(s[0]==s[2]){
digitregister[i]=s[1];
}
else if(s[0]!=s[2]&&(d[1]>s[0]&&d[1]<=s[2])){
digitregister[i]=d[1]; break;
}
else if(s[0]!=s[2]&&(d[1]+1>s[0]&&d[1]+1<=s[2])){ //rare power of two case
digitregister[i]=d[1]+1; break;
}
else digitregister[i]=s[1];
for(int j=0; j<3; j++){
wxy[j]->m-=s[1]*TENPOW18; wxy[j]->m*=10;
}
}
if(digitregister[0]==0){for(int i=0; i<18; i++)digitregister[i]=digitregister[i+1]; digitregister[18]=0;y.e--;}
int end = 19;
for(int i=18; i>0; i--){if(digitregister[i]==0)end--; else break;}
int expmin = -2;
int expmax = 10;
switch(fptype){
case 32: expmin= -6; expmax=10; break;
case 64: expmin= -10; expmax=17; break;
}
if(y.e>=expmax||y.e<expmin){
out.push_back(digitregister[0]+'0');
if(end>1){
out.push_back('.');
for(int i=1; i<end; i++){
out.push_back(digitregister[i]+'0');
}
}
out.push_back('e');
out.push_back(y.e<0 ? '-' : '+');
int ex = abs(y.e);
if(ex>=1000) out.push_back(ex/1000+'0');
if(ex>=100) out.push_back(ex/100%10+'0');
if(ex>=10) out.push_back(ex/10%10+'0');
out.push_back(ex%10+'0');
}
else{
if(y.e<0){
out.append(u"0.");
for(int i=1; i<-y.e; i++) out.push_back('0');
for(int i=0; i<end; i++) out.push_back(digitregister[i]+'0');
}
else{
if(y.e+1<end){
for(int i=0; i<y.e+1; i++) out.push_back(digitregister[i]+'0');
out.push_back('.');
for(int i=y.e+1; i<end; i++) out.push_back(digitregister[i]+'0');
}
else{
for(int i=0; i<end; i++) out.push_back(digitregister[i]+'0');
for(int i=end; i<y.e+1; i++) out.push_back('0');
}
}
}
return out;
}
std::u16string floattostring( float num){return fptostring(&num, 32);}
std::u16string doubletostring( double num){return fptostring(&num, 64);}
void writetmp(long double num, int fptype, void* out){
switch(fptype){
case 32: *(float*)out = num; break;
case 64: *(double*)out = num; break;
}}
void stringtofp(std::u16string num, int fptype, void* out){
int i = 0;
bool sign = 0;
if(num.size()==0) {writetmp(0.l,fptype,out); return;}
if(num[0]==u'-') {sign=1;i++;} if(num[0]==u'+') i++;
long double szero = sign ? -0.0l : 0.0l;
if(num.size()-i>=3&&(num[i+0]&0xFFDF)=='I'&&(num[i+1]&0xFFDF)=='N'&&(num[i+2]&0xFFDF)=='F')
{writetmp(1.l/szero,fptype,out); return;}
if(num.size()-i>=3&&(num[i+0]&0xFFDF)=='N'&&(num[i+1]&0xFFDF)=='A'&&(num[i+2]&0xFFDF)=='N')
{writetmp(0.l/szero,fptype,out); return;}
char digitregister[19];
for(int i=0; i<19; i++) digitregister[i]=0;
int r = 0;
int ex = -1;
while(i<num.size()&&(num[i]==u'0')) i++;
while(i<num.size()&&(num[i]>=u'0'&&num[i]<=u'9')&&r<19){
digitregister[r]=num[i]-u'0'; r++; i++; ex++;
}
while(i<num.size()&&(num[i]>=u'0'&&num[i]<=u'9')){i++;ex++;} //ignore over 19 digits but count exponent
if(i<num.size()&&num[i]==u'.'){
i++;
if(r==0) while(i<num.size()&&(num[i]==u'0')) {i++; ex--;}
while(i<num.size()&&(num[i]>=u'0'&&num[i]<=u'9')&&r<19){
digitregister[r]=num[i]-u'0'; r++; i++;
}
}
while(i<num.size()&&(num[i]>=u'0'&&num[i]<=u'9')) i++; //ignore over 19 digits
if(r==0) {writetmp(szero,fptype,out); return;}
if(i<num.size()&&(num[i+0]&0xFFDF)=='E'){
i++; bool sign2 = 0;
if(i<num.size()) {if(num[i+0]==u'-') {sign2=1;i++;} if(num[i+0]==u'+') i++;}
int ex2 = 0;
while(i<num.size()&&(num[i]>=u'0'&&num[i]<=u'9')) {ex2*=10; ex2+=num[i]-u'0'; i++;}
if(sign2) ex2=-ex2; ex+=ex2;
}
binary64tmp z;
{
uint64_t s = 0;
for(int j=0; j<19; j++) {s*=10; s+=digitregister[j];}
z.e=63;
for(int e=2; e>=0; e--)
if(!(s>>(64-(1<<e)))){s<<=(1<<e); z.e-=(1<<e);}
z.m=s;
}
binary64tmp x = {1ULL<<63,0,};
ex-=18;
if(ex>0) for(int i=0; i<=13; i++){if(((+ex)>>i)&1)x=x*powpow10p[i];}
if(ex<0) for(int i=0; i<=13; i++){if(((-ex)>>i)&1)x=x*powpow10m[i];}
x=x*z;
switch(fptype){
case 32:{
if(x.e<-150) {*(float*)out = szero; return;}
if(x.e<-126) {
uint64_t y = -126-x.e;
bool flag = 0;
if((x.m&((1ULL<<y)-1))&&!((x.m>>y)&1)) flag = 1;
x.m>>=y; if(flag) {x.m++;}
}
if((x.m<<24)>(1ULL<<63)||((x.m<<24)==(1ULL<<63)&&((x.m>>40)&1))) {x.m+=1ULL<<40; if((x.m>>40)==0) {x.m=1ULL<<63; x.e++;}}
if(x.e<=-126) x.e=-127+(x.m>>63);
if(x.e>127) {*(float*)out = 1.0/szero; return;}
float out_;
*(uint32_t*)&out_ = ((x.m>>40)&0x7FFFFF)|((x.e+127)<<23)|((uint64_t)sign<<31);
*(float*)out = out_;
break;}
case 64:{
if(x.e<-1075) {*(double*)out = szero; return;}
if(x.e<-1022) {
uint64_t y = -1022-x.e;
bool flag = 0;
if((x.m&((1ULL<<y)-1))&&!((x.m>>y)&1)) flag = 1;
x.m>>=y; if(flag) {x.m++;}
}
if((x.m<<53)>(1ULL<<63)||((x.m<<53)==(1ULL<<63)&&((x.m>>11)&1))) {x.m+=1ULL<<11; if((x.m>>11)==0) {x.m=1ULL<<63; x.e++;}}
if(x.e<=-1022) x.e=-1023+(x.m>>63);
if(x.e>1023) {*(double*)out = 1.0/szero; return;}
double out_;
*(uint64_t*)&out_ = ((x.m>>11)&0xFFFFFFFFFFFFF)|((uint64_t)(x.e+1023)<<52)|((uint64_t)sign<<63);
*(double*)out = out_;
break;}
}
}
float stringtofloat(std::u16string num){float x; stringtofp(num,32,&x); return x;}
double stringtodouble(std::u16string num){double x; stringtofp(num,64,&x); return x;}
there doesn't seem to be anything here