From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from sonic304-24.consmr.mail.ne1.yahoo.com (sonic304-24.consmr.mail.ne1.yahoo.com [66.163.191.150]) by sourceware.org (Postfix) with ESMTPS id 57CF43857C76 for ; Sat, 6 Nov 2021 12:06:11 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.4.1 sourceware.org 57CF43857C76 Authentication-Results: sourceware.org; dmarc=none (p=none dis=none) header.from=att.net Authentication-Results: sourceware.org; spf=none smtp.mailfrom=att.net DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=att.net; s=s1024; t=1636200370; bh=bsZWZrZfMBOnsT/6a/eVKSdlduQzCffF3x/EMyHF7bE=; h=Date:From:Subject:To:In-Reply-To:References:From:Subject:Reply-To; b=iLQaQcCk1cF6wzUXUdvhCw0xXUwTelXbMwZoe57+X8KMZlqk8/93LdXdR22LS5i4aYIDrcYbIy0Xn8OpJI+1A97B/OI83EWETF8Y8RC6mTA2yd8W9lo3H3ZGL4VVJxa3ISodKR6ryVzk5IG8LM2raTAekq4WGvsxlQD7HKj43ag= X-SONIC-DKIM-SIGN: v=1; a=rsa-sha256; c=relaxed/relaxed; d=yahoo.com; s=s2048; t=1636200370; bh=LfUgOoS/YZ2BeKQH0Q7pZQBwAUAIMm9Vs8alqSGbfrX=; h=X-Sonic-MF:Date:From:Subject:To:From:Subject; b=tbg5q8ZxLxXuwJTjcqA8GT6T4RJI/i3rUQ/W9S3kjum0OEhbhJoKeHp6c4cO+dFKylZ8gPc4eP28/ESdRx6DLutSWg7mEXWZb6xee9WJaBJQk2BqWu9PSZXWh1potsAQ2hvbV/AVLsiqX9sHDH2Z9GKLrw9pahYol4BGQWiO3sCFeYFuV90/4OVceNmyA2lmpukGkFsJJzEUVErElwsGzV1U5uL0rXO7ufZvzj29vsFQLLXY1W3oq68nYQWVE4DEeJ+Tdiwf0GXVQTEYc+Un6V14gVheiChlqeDzPUXBVWGt2W/kU+VNg6qkYJHFi5aMYpIhGOpld65J89+OfxHl1w== X-YMail-OSG: syJq0KEVM1kj9_oZv0blrTSB7EnE.bMY7odvQ4qhHi5PMea4qTsoxna9eXnjQsw _W..6WRMCJioXFqjb_3siNaNddqoV4IfWHEM6SfnIsyiwmNIi2goFy8oxNmvb3auR8PkF_ueE4yy Wp0ZvxMRCdkGRrGJvGO8J5n6zE8SGNGn9VR41To2HqnuxSjttVru8lMCJnF0gZQ5MthIzp.KHjw6 5jHOzefHtPq3Z5tE3VQqwhdU.5nXFjfVheMo0JGqKLO3zvxBcQ4_aysSl1zvWuu8GqaFzsoUbmGO Ah8ootKNdy3c5SsWC1OeDWnmt1bTW4_gHVBnXFN5cwqLTybfX7zexq3BlzEz7q1S3O0G6pH1zEpy kZujo0JyQVNOwuIHAedIICZV0JDKESLVUMVPQFGbClnU3aVoZYIqIlVqtI.gIHvtVqiFCvi5kyk. FK2znScwV.8heuNkHyB8x3Z_pFFKBykDh_YkRafSKXC8XmjBjDKFWoGvsojeoE4XjdU.K_FN6NzS 7bbqfwztab6lrr03XYBzu9MLEPRT.Lb2UgYWAaeRZwVf8LHxLUBp5DPpy_uBhI3ynd4Lx54eS82F huTDx_3LLKJ2t182ZNPBNah2Uk2vMw25a_oz2bZhBGjkF8ju.hehEhVn1Z3pL2Mm0FWixRJGK_IW 5DX3yp_Vs68t9TZBu2KQtoy4hMCEXf3IzMQg2sd23iROvJA4qvqIp3gl85f.zYpmpPVCFUptLm6M gDOq3UMKfSERw67.HWIPxa.dVXYkqVmUijmYS3Lo5cIFXYRb0Lq_jSQ1DNlBvynRSVADtkEPGE8k fTwbcHLEkI3qSDCzYxg5kIel9RBwa9auv_LfZjgCURZROXfEYg5ReV7h27R4BoPxGw9obmOCPrso LF5PIrQYEjCVUNeB.JwF_vVUrT8diaqivBs0hDxK0y1nIMPON7YY8SsSMPy6VJrxfed6hBF3d._V 417wOYFaFuaBEuftPz08CCJYrB1Ce0HvW6X4em8t3OlQlAGRhFHwr1wGnxTs_w_Q4XaCqiC1GbBN n8NEAd2xGyDhick7Ij2ZeXsCRzHD1hkolPXhWql5gsIsLdtUFdPUKPRLyR90j53bvQfdKnrN3BH7 .YIOYoPCKwAXSy.r.EUo8tI1CUHrWKAako4gwQkpVgIEAPoBN0gX5ZQ0R1EyRZf1eeJ6g6SSuByy vrpadqsMoUoZz_EuJTEsAruvvOSzZig5R9nmRabvV5AxYvZcEL22z1bELS.kYZDT.U__zxEe7_Wv rSI_bmHlYeVWbgob13tBw5HurEDSe1CWsyPI6fOlH48MErqFtqcDJ91GwRXKA1Rl1S51ZT_bh3N3 FbyE9gKfz99HsFYBIOOLvUjgI5NLU5v.Bcf66YqTCRMtZmNEIW.meKQDK7nzn2CcXO_M24f_CbOE brBF7dn0Sd4O2TBfaDBEsY55OprdO.Zp.e57osH4guNslpKybUvSLkTyaf4ZQj5AQ4azuCzIaMLu Ulvzped8uJIH_KWdz6mWPFYTQt6_ojRxBGagKpNt8XkjFQHJIC7YMrXtsDOoBAMQEbWFVXVST9rv f2Ya2QhcCEh9HeFIBTNncWMc.DIKQxmgKq574fdzfI5FuuWKY9yWu0thR4.adZm_seTAs8j3KWuf aVGvS_UEeDBGR8ldq_roNYXHszyTXHxwBkWNllXTu9eQgb9hzigAvqStnO1xARhE3eGoQEMnb3De on9Y6By3EuwrOVzs8mfm0a4a_JJ4OY9FTzUKmECjpClACNvPsmZop6XBOILYDX2MC3SNQ9zInzu9 t4sPPLwl2AaI0SOJALSfB.FsI6A.l7gYWzDVCngzYDE7p1zEldSE7BcrtHTDoSmexq_Q.mVkoQ3S LtvJsZ7qZZTCJIIaS7aTk7O1eCxXoL_2yV04ylZG9UXT_14JJqp2SgEbowRip2C7ikAGS_mrcsc9 n1bWntWW2jAMRL.JG0JaxFg4XweB3BIwWehDSPXxhGn_t9_wR5EhZxdP1PBeL7sPXdv2fajqAlON r2RaXZikQw6gG95ICbxR0uNhrgehTp4rvWHSxuO3M634ICmN4yJGiSarPUp5xAoOoOPi_iLL0svu GtPYZeILSaY6UwO.lXHiKSpPm8Oivzi2TnWH99RZTR4DNjyiuRLisJiF8jQa9CpTtndB4.tVdZ5P PotNc_lhEX.EmaZYulzejUiEMozpOVLGLXScQGTN0KLiNIsDJ5HRnsHYGz1lfWh3PUl4BjZhGe9X cV4EGM6b25A-- X-Sonic-MF: Received: from sonic.gate.mail.ne1.yahoo.com by sonic304.consmr.mail.ne1.yahoo.com with HTTP; Sat, 6 Nov 2021 12:06:10 +0000 Received: by kubenode549.mail-prod1.omega.bf1.yahoo.com (VZM Hermes SMTP Server) with ESMTPA ID 5366d6a8ef4cf73701c28d5196937243; Sat, 06 Nov 2021 12:06:08 +0000 (UTC) Date: Sat, 06 Nov 2021 08:06:05 -0400 From: Steven Abner Subject: Re: gcc 11.1.0: printf("%.43f\n", 0x1.52f8a8e32e982p-140): printed value is incorrectly rounded To: Newlib Message-Id: <1636200365.25270.0@smtp.mail.att.net> In-Reply-To: <1636128049.25340.0@smtp.mail.att.net> References: <87ilx7vpco.fsf@keithp.com> <90598863-5b2c-30a1-ddb4-0d6d7421a835@nadler.com> <874k8rvm49.fsf@keithp.com> <6b742e91-5a56-8a4d-c9a5-b33d2b33bb66@SystematicSw.ab.ca> <87r1bvtipt.fsf@keithp.com> <1636128049.25340.0@smtp.mail.att.net> X-Mailer: pantheon-mail/1.0.8 MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="=-es46Z1o6jxz9SwEeQxX6" X-Spam-Status: No, score=-2.1 required=5.0 tests=BAYES_00, DKIM_SIGNED, DKIM_VALID, DKIM_VALID_AU, DKIM_VALID_EF, FREEMAIL_FROM, RCVD_IN_DNSWL_NONE, SPF_HELO_NONE, SPF_NONE, TXREP autolearn=ham autolearn_force=no version=3.4.4 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on server2.sourceware.org X-BeenThere: newlib@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Newlib mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 06 Nov 2021 12:06:16 -0000 --=-es46Z1o6jxz9SwEeQxX6 Content-Type: text/plain; charset=us-ascii; format=flowed I didn't mean to stop the conversion! Here is a butchered, self contained use. Hopefully this provides better means to evaluate. Low precision code removed, since I assume embedded systems won't use. In case of not obvious, entry from printf() into realto functions expected a 64 byte buffer to write to. realto functions has a static buffer to cover double range. It uses malloc() if it determines that with localization and precision, if it needs to increase past 64 byte return. realto functions includes decimal encoding of up to 4 byte utf8 encoding, which at the time of writing no script came close to, think 2 bytes was max? There is commented out code for 32 bit float, but not sure would ever get used. sample use: steven@steven-ryzen:~/Development/Projects$ gcc newlib_printf_demo.c -o nldemo steven@steven-ryzen:~/Development/Projects$ ./nldemo 0.0000000000000000000000000000000000000000009 3.14159 Steve --=-es46Z1o6jxz9SwEeQxX6 Content-Type: text/x-csrc Content-Disposition: attachment; filename=newlib_printf_demo.c Content-Transfer-Encoding: quoted-printable #include /* must define a data model */ #ifndef __LP64__ #define __LP32__ 1 #endif /* found out 24 Jan 2019 from new computer build */ /* gcc does not define __LITTLE_ENDIAN__ or __BIG_ENDIAN__ */ #ifndef __LITTLE_ENDIAN__ #if (__BYTE_ORDER__ =3D=3D __ORDER_LITTLE_ENDIAN__) #define __LITTLE_ENDIAN__ #else #define __BIG_ENDIAN__ #endif #endif typedef signed char int8_t; typedef unsigned char uint8_t; typedef short int16_t; typedef unsigned short uint16_t; typedef int int32_t; typedef unsigned int uint32_t; typedef long long int64_t; typedef unsigned long long uint64_t; #if __SIZEOF_INT128__ =3D=3D 16 typedef __int128 int128_t; typedef unsigned __int128 uint128_t; #endif typedef union { #if __SIZEOF_INT128__ =3D=3D 16 uint128_t uo; #endif #ifdef __BIG_ENDIAN__ struct { uint64_t ull1, ull0; } ulls; struct { uint32_t ul3, ul2, ul1, ul0; } uls; struct { uint16_t us7, us6, us5, us4, us3, us2, us1, us0; } uss; #else /* __LITTLE_ENDIAN__ */ struct { uint64_t ull0, ull1; } ulls; struct { uint32_t ul0, ul1, ul2, ul3; } uls; struct { uint16_t us0, us1, us2, us3, us4, us5, us6, us7; } uss; #endif } octle; /* defined in as_stdio.h */ #define quadruple octle /* gcc appears to only defines _LP64 __LP64__ (2019) */ /* only models with 64bit ints */ #if defined (__ILP64__) || defined (__SILP64__) typedef int64_t intmax_t; typedef uint64_t uintmax_t; #elif defined (__LP32__) typedef int16_t intmax_t; typedef uint16_t uintmax_t; #else /* __LP64__ has 32bit int */ typedef int32_t intmax_t; typedef uint32_t uintmax_t; #endif /* these are all register size based */ #if defined (__LP32__) || defined (__ILP32__) typedef int32_t intptr_t; typedef uint32_t uintptr_t; #define __SIZEOF_REGISTER_T__ 4 #define __SIZEOF_SIZE_T__ 4 #else typedef int64_t intptr_t; typedef uint64_t uintptr_t; #define __SIZEOF_REGISTER_T__ 8 #define __SIZEOF_SIZE_T__ 8 #endif typedef intptr_t register_t; typedef intptr_t _intptr_t; typedef uintptr_t u_intptr_t; typedef intptr_t ssize_t; /* Used for a count of bytes or an error = indication. */ typedef uintptr_t size_t; /* Used for sizes of objects. */ typedef intptr_t ptrdiff_t; /* Used for subtracting two pointers. */ typedef uintptr_t clock_t; /* Used for system times in clock ticks o= r CLOCKS_PER_SEC . */ typedef intptr_t time_t; /* Used for time in seconds. */ #define NULL ((void *)0) #ifdef __APPLE__ #define FILE FILE_mac #elif defined __linux /* !__APPLE__ */ typedef struct _IO_FILE _IO_FILE; #define FILE _IO_FILE #else /* !__APPLE__ !__CYGWIN__ !__linux */ #error: unknown FILE structure #endif /* FILE * */ extern void *malloc(size_t); extern void free(void *); extern int puts(const char *); size_t fwrite(const void *restrict, size_t, size_t, FILE *restrict); extern FILE* stdout; /* in ldtostr.c */ extern unsigned char * realtostr(uint32_t [5], unsigned char **, int, int, int, unsigned); extern unsigned char * realtohex(uint32_t [5], unsigned char **, int, int, int, unsigned); /* guarentee use of this printf */ int phx_printf(const char *format, ...); /* local defines */ #define ALT_PREFIX 0002 #define PREC_PREFIX 0400 /* converts to w160 type used by realtostr, * current use as 63-63 with implicit '1"-implicit nil */ static void qtow(quadruple qdbl, uint32_t w160[5]) { /* required strict aliasing conversion */ octle *iqdbl; unsigned char *cqdbl =3D (unsigned char*)&qdbl; iqdbl =3D (octle*)cqdbl; /* long double shifted to high bit of uint128_t, with implicit * bit becoming explicit at bit 127 */ w160[4] =3D iqdbl->uss.us7; *(uint64_t*)&w160[2] =3D (iqdbl->ulls.ull1 << 15) | (iqdbl->uls.ul1 >> 17); /* lower 64 set same position as upper, upper gets 'implicit '1' */ iqdbl->uss.us3 &=3D 1; *(uint64_t*)&w160[0] =3D (iqdbl->ulls.ull0 << 14); /* place implicit bit for all but denormal */ w160[3] |=3D (w160[4] !=3D 0) << 31; } /* converts to w160 type used by realtostr */ static void ldtow(long double ldbl, uint32_t w160[5]) { /* required strict aliasing conversion */ octle *ildbl; unsigned char *cldbl =3D (unsigned char*)&ldbl; ildbl =3D (octle*)cldbl; w160[4] =3D ildbl->uss.us4; *(uint64_t*)&w160[2] =3D ildbl->ulls.ull0; *(uint64_t*)w160 =3D 0; } /* converts to w160 type used by realtostr */ static void dtow(double dbl, uint32_t w160[5]) { /* required strict aliasing conversion */ uint64_t *idbl; unsigned char *cdbl =3D (unsigned char*)&dbl; idbl =3D (uint64_t*)cdbl; *(uint64_t*)&w160[2] =3D *idbl << 11; w160[4] =3D (uint32_t)(*idbl >> 52); if ((w160[4] & 0x7FF) =3D=3D 0x7FF) w160[4] =3D (w160[4] << 4) | 0x000F; else if (((w160[4] & 0x7FF) =3D=3D 0) && (*(uint64_t*)&w160[2] =3D=3D 0))= { w160[4] <<=3D 4; } else { w160[4] =3D ((w160[4] & 0x7FF) + 16383 - 1023) | ((w160[4] & 0x800) << = 4); /* place implicit bit for all but denormal */ w160[3] |=3D ((w160[4] & 0x7FFF) !=3D 0x3C00) << 31; } *(uint64_t*)w160 =3D 0; } /* not sure usefulness, gcc's code forces doubles in printf. Not sure of ov= erhead * created by gcc and/or POSIX requirement of double argument. static void ftow(float flt, uint32_t w160[5]) { *(uint64_t*)w160 =3D 0; *(uint64_t*)&w160[2] =3D 0; w160[4] =3D (*(uint32_t*)&flt) >> 23; w160[3] =3D (*(uint32_t*)&flt) << 9; if ((w160[4] & 0xFF) =3D=3D 0xFF) w160[4] =3D (w160[4] << 8) | 0x00FF; else w160[4] =3D ((w160[4] & 0xFF) + 16383 - 127) | ((w160[4] & 0x100) << 7)= ; } */ /* locate dependent utf8 loading of decimal point */ /* for newlib demo return ascii */ struct __params { unsigned char **s; unsigned char *slimit; //locale_t ld; }; static unsigned __decimalLoad(unsigned char *s, struct __params *params) { (void)params; *s =3D '.'; return (1); } #define BUFSTRSIZE 64 unsigned char * __format(const unsigned char *format, struct __params *params, va_list arg)= { unsigned char buffer[BUFSTRSIZE]; /* size of 64bit binary ouput no pre= c, * field nor alt goes here */ unsigned char **s =3D params->s; unsigned char *slimit =3D params->slimit; unsigned char *strBuf =3D buffer, *strPtr, *strEnd, *sPtr =3D *s; uint64_t val_ull; /* for any sizeof() > 4 bytes, if * unsigned is 64bits it will bypass * val_ul to go here */ unsigned val_ul; /* need unsigned instead of unsigned = long * incase unsigned 64bits */ uint32_t arg_flags; /* arg_flags needs to be 32bits??? */ int arg_prec, arg_field, n; unsigned char ch, sign, arg_size; /* arg_size max value is sizeof(long * double) if int goes 128bits still * within char size */ uintptr_t *vPtr =3D NULL; --format; goto L500; do { strEnd =3D strPtr =3D &strBuf[BUFSTRSIZE]; arg_flags =3D 0; arg_size =3D sizeof(int); arg_prec =3D arg_field =3D -1; sign =3D 0; next_ch: ch =3D *++format; EO_EXTENSIONS: if (!ch) break; switch (ch) { /* Removed Flags for demo */ /* Removed some Length modifiers for demo */ case 'l': if ((ch =3D *++format) !=3D 'l') { /* and ll */ arg_size =3D sizeof(long); goto EO_EXTENSIONS; } arg_size =3D sizeof(long long); goto next_ch; case 'q': arg_size =3D sizeof(octle); goto next_ch; /* 128bi= t */ case 'L': arg_size =3D sizeof(long double); goto next_ch; /* long = double */ /* Removed some Field-Width / Precision for demo */ case '.': ch =3D *++format; arg_flags |=3D PREC_PREFIX; case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': n =3D ch - 48; while ((unsigned)((ch =3D *++format) - '0') <=3D 9) n =3D (n * 10) + (ch - '0'); if (ch !=3D '$') { if (arg_flags & PREC_PREFIX) arg_prec =3D n; else arg_field =3D n; arg_flags &=3D ~PREC_PREFIX; goto EO_EXTENSIONS; } goto next_ch; /* Removed all Formats... but real */ case 'a': case 'e': case 'f': case 'g': case 'A': case 'E': case 'F': case 'G': /* gcc doesn't permit floats, it converts to doubles */ { unsigned char *(*toStr)(uint32_t [5], unsigned char **, int, int, int, unsign= ed); uint32_t w160[5]; uint32_t dp; /* 4 byte decimal point param */ int32_t sch; __decimalLoad((unsigned char*)&dp, params); if ((ch | 0x20) =3D=3D 'a') ch -=3D 'A', toStr =3D realtohex; else ch -=3D 'E', toStr =3D realtostr; sch =3D format[-1]; /* force 128bit long double into quad parsing */ if (sch =3D=3D 'q') { #if (__LDBL_MANT_DIG__ !=3D 64) quad_parse: #endif sch =3D (('q' << 8) | ch); qtow(((vPtr =3D=3D NULL) ? va_arg(arg, octle) : *(octle *)vPtr), = w160); } else if (sch =3D=3D 'L') { #if (__LDBL_MANT_DIG__ !=3D 64) goto quad_parse; #endif sch =3D (('L' << 8) | ch); ldtow(((vPtr =3D=3D NULL) ? va_arg(arg, long double) : *(long dou= ble *)vPtr), w160); } else { sch =3D (('l' << 8) | ch); dtow(((vPtr =3D=3D NULL) ? va_arg(arg, double) : *(double *)vPtr)= , w160); } arg_flags |=3D (strBuf !=3D buffer) << 7; strPtr =3D toStr(w160, &strBuf, sch, arg_prec, (int)arg_flags, dp); } if (strPtr =3D=3D NULL) goto error; val_ul =3D (n =3D (int)(strPtr - strBuf)); strPtr =3D strBuf; arg_prec =3D 0; arg_flags &=3D ~ALT_PREFIX; /* grouping code removed for demo */ } /* end switch */ /* main output section */ /* adjust the different output factors */ if ((arg_prec -=3D n) < 0) arg_prec =3D 0; if ((arg_field -=3D (arg_prec + val_ul)) < 0) arg_field =3D 0; /* is there output space? removed for demo */ /* different formatting removed for demo */ if (sign !=3D 0) *sPtr++ =3D sign; if (arg_prec) do *sPtr++ =3D '0'; while ((--arg_prec)); if (n & 1) *sPtr++ =3D *strPtr++; if (n & 2) { *(int16_t *)sPtr =3D *(int16_t *)strPtr; sPtr +=3D sizeof(int16_t); strPtr +=3D sizeof(int16_t); = } if ((n >>=3D 2)) do { *(int32_t *)sPtr =3D *(int32_t *)strPtr; sPtr +=3D sizeof(int32_t); strPtr +=3D sizeof(int32_t); } while ((--n)); if (arg_field) do *sPtr++ =3D ' '; while ((--arg_field)); /* parsing of format string */ do { L500: ch =3D *++format; if ((sPtr >=3D slimit)) { --format; /* removed buffer extreme */ puts("this is demo, don't be cute!"); goto end; } *sPtr =3D ch; if (!ch) goto end; if (ch =3D=3D '%') break; sPtr++; } while (1); } while (1); *sPtr =3D 0; end: if (strBuf !=3D buffer) free(strBuf); return sPtr; error: if (strBuf !=3D buffer) free(strBuf); *--sPtr =3D 0; return (*s - 12); } #define FORMAT_LINE_MAX 4096 static int phx_vfprintf(FILE *stream, const char *format, va_list arg) { //_vfprintf(FILE *stream, const char *format, va_list arg, locale_t ld) { unsigned char *sPtr, *wrPtr, *rdPtr; ssize_t to_write; size_t u8Sz =3D FORMAT_LINE_MAX; struct __params params; /* make sure writable stream, Remove writable for demo */ if ( (stream =3D=3D NULL) //|| (!_io_writable(stream)) || (format =3D=3D NULL) ) return -1; /* create buffer to 'write' into */ if ((wrPtr =3D (unsigned char*)malloc(u8Sz)) =3D=3D NULL) return -1; /* label rdPtr as . If conversion, rdPtr will be allocated * and filled with UTF-8 version of */ rdPtr =3D (unsigned char*)format; /* UTF-8 input */ params.s =3D &wrPtr; params.slimit =3D (wrPtr + u8Sz); //params.ld =3D ld; sPtr =3D __format((const unsigned char *)rdPtr, ¶ms, arg); /* UTF-8 text ouput */ if ((to_write =3D (ssize_t)(sPtr - wrPtr)) > 0) to_write =3D fwrite((const void *)wrPtr, (size_t)1, (size_t)to_write, s= tream); free(wrPtr); return (int)to_write; } int phx_printf(const char *format, ...) { int written; va_list arg; va_start(arg, format); written =3D phx_vfprintf((FILE *)stdout, format, arg); //written =3D _vfprintf((FILE *)stdout, format, arg, NULL); va_end(arg); return written; } /* * realtostr - convert real numbers to formated latn numerical string * realtohex - convert real numbers to formated hexadecimal string * qtold - convert quadruple-precision number to a system defined long dou= ble * ldtod - convert a system defined long double to a double-precision numb= er */ /* C langauge - string literals */ #ifdef __BIG_ENDIAN__ #define cconst16_t(a,b) (uint16_t)(((uint8_t)(a) << 8) | (uint8_t)(b)) #define cconst32_t(a,b,c,d) (uint32_t)(((uint8_t)(a) << 24) | ((uint8_t)(b= ) << 16) | ((uint8_t)(c) << 8) | (uint8_t)(d)) #else #define cconst16_t(a,b) (uint16_t)(((uint8_t)(b) << 8) | (uint8_t)(a)) #define cconst32_t(a,b,c,d) (uint32_t)(((uint8_t)(d) << 24) | ((uint8_t)(c= ) << 16) | ((uint8_t)(b) << 8) | (uint8_t)(a)) #endif /* compiler convience macro */ #ifndef __GNUC_PREREQ #if defined __GNUC__ && defined __GNUC_MINOR__ #define __GNUC_PREREQ(major,minor) \ (((__GNUC__ << 16) + __GNUC_MINOR__) >=3D (((major) << 16) + (minor))) #else #define __GNUC_PREREQ(major,minor) 0 #endif #endif /* __GNUC_PREREQ */ #ifndef __LDBL_MANT_DIG__ #define __LDBL_MANT_DIG__ LDBL_MANT_DIG #endif /* endian break down of uint64_t */ #ifdef __BIG_ENDIAN__ typedef struct { uint16_t w3; uint16_t w2; uint16_t w1; uint16_t w0; } = _ull_ws; typedef struct { uint32_t hi; uint32_t lo; } _ull_ls; #elif defined (__LITTLE_ENDIAN__) typedef struct { uint16_t w0; uint16_t w1; uint16_t w2; uint16_t w3; } = _ull_ws; typedef struct { uint32_t lo; uint32_t hi; } _ull_ls; #else #error: undefined endianness #endif typedef union { uint64_t ull; _ull_ls uls; _ull_ws uss; } _ull_t; #ifdef __ppc__ #define CNTLZW(x,y) \ do __asm__("cntlzw %0,%1" : "=3Dr" (x) : "r" (y)); while(0) #elif defined (__i386__) #define CNTLZW(x,y) \ do { __asm__("bsr %1,%0" : "=3Dr" (x) : "r" (y)); x =3D 31 - x; } wh= ile(0) #elif defined (__x86_64__) #define CNTLZW(x,y) \ do { __asm__("bsr %1,%0" : "=3Dr" (x) : "r" (y)); x =3D 31 - x; } wh= ile(0) #define CNTLZD(x,y) \ do { __asm__("bsr %1,%0" : "=3Dr" (x) : "r" (y)); x =3D 63 - x; } wh= ile(0) #else /* about 50% slower, 3 instructions vs 31 instructions */ #define CNTLZW(x,y) \ do { unsigned t =3D y; \ t |=3D t >> 1; \ t |=3D t >> 2; \ t |=3D t >> 4; \ t |=3D t >> 8; \ t |=3D t >> 16; \ t =3D t - ((t >> 1) & 0x55555555); \ t =3D (t & 0x33333333) + ((t >> 2) & 0x33333333); \ x =3D 32 - (((t + (t >> 4) & 0xF0F0F0F) * 0x1010101) >> 24); \ } while (0); #endif /* used for __fractional128 stack zeroing */ extern void *memset(void *, int, size_t); long double qtold(quadruple); double ldtod(long double); #define _CARRY 1000000000000000000ULL /* 0x 0DE0B6B3 A7640000 */ /* for rounding */ static const uint64_t _fcarry[] =3D { 5, 50, 500, 5000, 50000, 500000, 5000000, 50000000, 500000000, 5000000000= ULL, 50000000000ULL, 500000000000ULL, 5000000000000ULL, 50000000000000ULL, 500= 000000000000ULL, 5000000000000000ULL, 50000000000000000ULL, 500000000000000000ULL, /*18*/ 1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000= , 10000000000ULL, 100000000000ULL, 1000000000000ULL, 10000000000000ULL, 100= 000000000000ULL, 1000000000000000ULL, 10000000000000000ULL, 100000000000000000ULL, _CARRY }; /* for cases: nan, inf, hex */ static const char *_const_outStr =3D "0123456789abcdefnanpinfe0x0.00000123456789ABCDEFNANPINFE0X0.0000"; /*** internal routines ***/ /* generic, div() moved to 64_t to address endianess */ static uint64_t divlu2r(uint64_t) __attribute__ ((noinline)); static void _ulltobuffer_18lr(uint64_t, unsigned char *) __attribute__ ((n= oinline)); static void _ulltobuffer_nlr(uint64_t, unsigned char *, int) __attribute__= ((noinline)); static int _first_digit(uint64_t); /* ldtostr specific */ static unsigned char *__n_char_sequence(uint32_t [5], unsigned char *, int)= ; static uint32_t __fractional(_ull_t *, _ull_t); static uint32_t __fractional_128(_ull_t *, _ull_t); static uint32_t __positive_exponent_fp(_ull_t *, int); static uint32_t __negative_exponent_fp(_ull_t *, int); /* from fprintf... needs to be the same */ #define ZERO_PREFIX 0001 #define ALT_PREFIX 0002 #define LEFT_PREFIX 0004 #define SIGN_PREFIX 0010 #define SPC_PREFIX 0020 #define STAR_PREFIX 0040 #define GRP_PREFIX 0100 #define USGN_PREFIX 0200 #define PREC_PREFIX 0400 #define ALLOC_BUF USGN_PREFIX #define BUFSTRSIZE 64 /* XXX using upper 16 to pass 'unit' max from calcation prior to entry of x_exponent */ #define UPPER_TYPE 0x20 #define IS_UPPER_TYPE(a) (!((a) & UPPER_TYPE)) #define IS_QUAD(a) ((uint8_t)((a) >> 8) =3D=3D (uint8_t)'q') #define IS_LDBL(a) ((uint8_t)((a) >> 8) =3D=3D (uint8_t)'L') #define IS_DBL(a) ((uint8_t)((a) >> 8) =3D=3D (uint8_t)'l') #define IS_FLT(a) ((uint8_t)((a) >> 8) =3D=3D 0) #define TYPE_ZERO 0x40 #define IS_TYPE_ZERO(a) ((a) & TYPE_ZERO) #define TYPE_E 0 #define TYPE_F 1 #define TYPE_G 2 #define IS_TYPE_E(a) (((a) & 3) =3D=3D 0) #define IS_TYPE_F(a) ((a) & TYPE_F) #define IS_TYPE_G(a) ((a) & TYPE_G) #define SURPRESS_NSEQ(a) (((a) & 4) !=3D 0) #pragma mark *** Support Mathematics *** static uint64_t divlu2r(uint64_t c) { _ull_t q; uint32_t mConst =3D 2305843009U; // 2^61 / 10^9 uint32_t tConst =3D 1000000000U; q.uls.hi =3D (uint32_t)(c >> 32); q.uls.lo =3D ((uint32_t)c); if (q.uls.hi !=3D 0) { uint32_t shift; CNTLZW(shift, q.uls.hi); q.uls.hi =3D (uint32_t)((((q.uls.hi << shift) | (q.uls.lo >> (32 - shif= t))) * (uint64_t)mConst) >> = 32); if (shift > 3) q.uls.hi >>=3D (shift - 3); } else { q.uls.hi =3D (uint32_t)(((q.uls.lo | 1) * (uint64_t)mConst) >> (32 + 29= )); } q.uls.lo =3D (uint32_t)(c - ((uint64_t)q.uls.hi * (uint32_t)tConst)); if (q.uls.lo >=3D tConst) q.uls.hi +=3D 1, q.uls.lo -=3D tConst; return q.ull; } static void _ulltobuffer_18lr(uint64_t c, unsigned char *buf) { _ull_t a, b; int n; --buf; /* in 4.60, 2.30 fixed point integers */ b.ull =3D (uint64_t)(((uint32_t)c) << 2) * 2882303762U; a.ull =3D (uint64_t)(((uint32_t)(c >> 32)) << 2) * 2882303762U; n =3D 9; goto L1; do { b.uls.hi =3D (((b.uls.hi << 4) >> 4) * 10) + (b.uls.lo >> 28); a.uls.hi =3D (((a.uls.hi << 4) >> 4) * 10) + (a.uls.lo >> 28); b.uls.lo <<=3D 4; a.uls.lo <<=3D 4; L1: *++buf =3D (a.uls.hi >> 28) + '0'; buf[9] =3D (b.uls.hi >> 28) + '0'; } while ((--n)); } /* went for speed */ static void _ulltobuffer_nlr(uint64_t c, unsigned char *buf, int m) { _ull_t b; int n; --buf; /* in 4.60, 2.30 fixed point integers */ c <<=3D 2; b.ull =3D (uint64_t)((uint32_t)(c >> 32)) * 2882303762U; n =3D (m <=3D 9) ? m : 9; m -=3D 9; goto L1; do { b.uls.hi =3D (((b.uls.hi << 4) >> 4) * 10) + (b.uls.lo >> 28); b.uls.lo <<=3D 4; L1: *++buf =3D (b.uls.hi >> 28) + '0'; } while ((--n)); if (m > 0) { n =3D m; b.ull =3D (uint64_t)((uint32_t)c) * 2882303762U; goto L2; do { b.uls.hi =3D (((b.uls.hi << 4) >> 4) * 10) + (b.uls.lo >> 28); b.uls.lo <<=3D 4; L2: *++buf =3D (b.uls.hi >> 28) + '0'; } while ((--n)); } } /*#pragma GCC diagnostic ignored "-Wmisleading-indentation"*/ #if __GNUC_PREREQ(4,5) #pragma GCC diagnostic ignored "-Wmisleading-indentation" #endif /* uses 2 64s as swap loads, this prevents stalls */ /* returns 1 less than actual digits, the power of 10, or log10 value */ static int _first_digit(uint64_t a) { _ull_t reg; const uint64_t *fc; uint64_t b; reg.ull =3D a; fc =3D _fcarry; if (reg.uls.hi) { a =3D fc[35]; b =3D fc[34]; if (reg.ull >=3D a) return 17; a =3D fc[33]; if (reg.ull >=3D b) return 16; b =3D fc[32]; if (reg.ull >=3D a) return 15; a =3D fc[31]; if (reg.ull >=3D b) return 14; b =3D fc[30]; if (reg.ull >=3D a) return 13; a =3D fc[29]; if (reg.ull >=3D b) return 12; b =3D fc[28]; if (reg.ull >=3D a) return 11; if (reg.ull >=3D b) return 10; return 9; } else { uint32_t c, d; c =3D (uint32_t)fc[27]; d =3D (uint32_t)fc[26]; if (reg.uls.lo >=3D c) return 9; c =3D (uint32_t)fc[25]; if (reg.uls.lo >=3D d) return 8; d =3D (uint32_t)fc[24]; if (reg.uls.lo >=3D c) return 7; c =3D (uint32_t)fc[23]; if (reg.uls.lo >=3D d) return 6; if (reg.uls.lo >=3D c) return 5; if (reg.uls.lo >=3D 10000U) return 4; if (reg.uls.lo >=3D 1000U) return 3; if (reg.uls.lo >=3D 100U) return 2; if (reg.uls.lo >=3D 10U) return 1; return 0; } } /*#pragma GCC diagnostic pop*/ #if __GNUC_PREREQ(4,5) #pragma GCC diagnostic pop #endif #pragma mark *** Support Subroutines *** static unsigned char * __n_char_sequence(uint32_t w160[5], unsigned char *sPtr, int type) { const char *_outvalues; unsigned is_zero; /* e80 long double type, 63rd bit required set reguardless if nan or in= f */ /* if is_zero is , then is infinity, else nan */ /* any info in nan is considered n-char-sequence */ /* bit 62 of nan, when 0, and others [61-0] set, is considered signalli= ng nan. * bit 62 of nan, when 1, if others 0, special quiet nan, others [61-0]= set, * considered quiet signaling nan. */ is_zero =3D ((w160[3] & 0x7FFFFFFF) | w160[2] | w160[1] | w160[0]); _outvalues =3D _const_outStr + 16; if (IS_UPPER_TYPE(type)) _outvalues +=3D 32; if (!is_zero){ *((uint32_t *)sPtr) =3D ((int32_t *)_outvalues)[1]; return (sPtr + 3); } *((uint32_t *)sPtr) =3D ((int32_t *)_outvalues)[0]; sPtr +=3D 3; if (SURPRESS_NSEQ(type)) return sPtr; /* for doubles 53bits exposed, shift off hidden */ if (!IS_LDBL(type)) *(uint64_t*)&w160[2] <<=3D 1; if (IS_QUAD(type)) { w160[2] |=3D (w160[1] & 0x40000000) >> 30; *(uint64_t*)&w160[0] <<=3D 2; } /* now output the n-char-sequence */ /* (0x hexdigits ) * e80,16nibbles, doubles,13nibbles, floats,5.5nibbles, */ *((uint16_t *)sPtr) =3D cconst16_t('(','0'); sPtr +=3D 2; *sPtr =3D _outvalues[9]; /* 'x' or 'X' */ sPtr++; _outvalues -=3D 16; sPtr[ 0] =3D _outvalues[(w160[3] >> 28)]; sPtr[ 1] =3D _outvalues[((w160[3] << 4) >> 28)]; sPtr[ 2] =3D _outvalues[((w160[3] << 8) >> 28)]; sPtr[ 3] =3D _outvalues[((w160[3] << 12) >> 28)]; sPtr[ 4] =3D _outvalues[((w160[3] << 16) >> 28)]; sPtr[ 5] =3D _outvalues[((w160[3] << 20) >> 28)]; if (IS_FLT(type)) { sPtr +=3D 6; goto end_brace; } sPtr[ 6] =3D _outvalues[((w160[3] << 24) >> 28)]; sPtr[ 7] =3D _outvalues[((w160[3] << 28) >> 28)]; sPtr[ 8] =3D _outvalues[(w160[2] >> 28)]; sPtr[ 9] =3D _outvalues[((w160[2] << 4) >> 28)]; sPtr[10] =3D _outvalues[((w160[2] << 8) >> 28)]; sPtr[11] =3D _outvalues[((w160[2] << 12) >> 28)]; sPtr[12] =3D _outvalues[((w160[2] << 16) >> 28)]; if (IS_DBL(type)) { sPtr +=3D 13; goto end_brace; } sPtr[13] =3D _outvalues[((w160[2] << 20) >> 28)]; sPtr[14] =3D _outvalues[((w160[2] << 24) >> 28)]; sPtr[15] =3D _outvalues[((w160[2] << 28) >> 28)]; sPtr +=3D 16; if (IS_LDBL(type)) goto end_brace; /* undo w160 formatting of '[0]' for number */ sPtr[ 0] =3D _outvalues[(w160[1] >> 28)]; sPtr[ 1] =3D _outvalues[((w160[1] << 4) >> 28)]; sPtr[ 2] =3D _outvalues[((w160[1] << 8) >> 28)]; sPtr[ 3] =3D _outvalues[((w160[1] << 12) >> 28)]; sPtr[ 4] =3D _outvalues[((w160[1] << 16) >> 28)]; sPtr[ 5] =3D _outvalues[((w160[1] << 20) >> 28)]; sPtr[ 6] =3D _outvalues[((w160[1] << 24) >> 28)]; sPtr[ 7] =3D _outvalues[((w160[1] << 28) >> 28)]; sPtr[ 8] =3D _outvalues[(w160[0] >> 28)]; sPtr[ 9] =3D _outvalues[((w160[0] << 4) >> 28)]; sPtr[10] =3D _outvalues[((w160[0] << 8) >> 28)]; sPtr[11] =3D _outvalues[((w160[0] << 12) >> 28)]; sPtr +=3D 12; end_brace: *sPtr =3D ')'; return (sPtr + 1); } /* converts fractional portion of IEEE754 value to * base 10 representation, in 10^18 ull format. */ static uint32_t __fractional(_ull_t *uPtr, _ull_t reg) { uint32_t last =3D 1; uint64_t mlConst =3D 3814697265625ULL; /* 5^18 (42bits) 0x378 2DACE9D9 *= / /* 0x7FFFFFFF FFFFFFFF start 0x7FFFE 1FFF FFFFFFFF split */ uPtr[0].ull =3D (reg.uls.hi >> 13) * mlConst; /* 18 + 45 =3D=3D 63 */ /* test added, rather than just '&=3D'. This is guess that small 'stand= ard' * fractions might be used */ if (((reg.uls.hi &=3D 0x00001FFF) | reg.uls.lo) !=3D 0) { /* the left over 13bits from hi with 5bits from lo (45-18=3D>27). 0x00001FFF FFFFFFFF start 0x00001FFFF8 7FFFFFF split */ uint64_t u64; last +=3D 1; uPtr[0].ull +=3D (u64 =3D ((uint32_t)(reg.ull >> 27)) * mlConst) >> 18; uPtr[1].ull =3D ((uint32_t)u64 & 0x3FFFF) * mlConst; if ((reg.uls.lo &=3D 0x7FFFFFF) !=3D 0) { /* 0x00000000 07FFFFFF start 0x 7FFFE 1FF split */ _ull_t scratch; last +=3D 1; uPtr[0].ull +=3D (u64 =3D (reg.uls.lo >> 9) * mlConst) >> 36; scratch.ull =3D (uint32_t)((uint32_t)u64 & 0x3FFFF) * mlConst; uPtr[2].ull =3D (scratch.uls.lo & 0x3FFFF) * mlConst; u64 =3D uPtr[1].ull + (scratch.ull >> 18) + (((uint32_t)(u64 >> 18)) & 0x3FFFF) * mlConst; if (u64 >=3D _CARRY) u64 -=3D _CARRY, uPtr[0].ull +=3D 1ULL; uPtr[1].ull =3D u64; if ((reg.uls.lo &=3D 0x1FF) !=3D 0) { /* 0x00000000 000001FF end of 63bit resolution */ uint32_t u32; last +=3D 1; // 54*10^18 u32 =3D (uint32_t)((scratch.ull =3D reg.uls.lo * mlConst) >> 27); uPtr[0].ull +=3D (uint64_t)(scratch.uls.hi >> 13); // 54*10^18*10^18 /* first 18 of 45 remainder */ u64 =3D ((u32 & 0x3FFFF) * mlConst) + uPtr[1].ull; /* second 18 of 45 remainder, scratch.uls.lo &=3D 0x07FFFFFF; */ u32 =3D scratch.uls.lo & 0x07FFFFFF; scratch.ull =3D (u32 >> 9) * mlConst; u64 +=3D (scratch.ull >> 18); uPtr[2].ull +=3D (scratch.uls.lo & 0x3FFFF) * mlConst; if (uPtr[2].ull >=3D _CARRY) uPtr[2].ull -=3D _CARRY, u64 +=3D 1UL= L; uPtr[1].ull =3D u64; /* last 9 of 45 remainder */ scratch.ull =3D ((u32 & 0x000001FF) << 9) * mlConst; uPtr[1].ull +=3D (uint64_t)(scratch.uls.hi >> 4); if (uPtr[1].ull >=3D _CARRY) uPtr[1].ull -=3D _CARRY, uPtr[0].ull = +=3D 1ULL; /* of scratch2 36 remainder */ scratch.uls.hi &=3D 0x0000000F; uPtr[2].ull +=3D ((uint32_t)(scratch.ull >> 18)) * mlConst; uPtr[2].ull +=3D ((scratch.uls.lo & 0x3FFFF) * mlConst) >> 18; if (uPtr[2].ull >=3D _CARRY) uPtr[2].ull -=3D _CARRY, uPtr[1].ull = +=3D 1ULL; /* of 18 remainder */ scratch.uls.lo =3D (uint32_t)((scratch.uls.lo & 0x3FFFF) * mlConst)= & 0x3FFFF; uPtr[3].ull =3D scratch.uls.lo * mlConst; } } } return last; } /* this is for second register of frational */ /* expected entry use: __fractional(uPtr, reg2); * expected range of reg2: 0x7fffffffffffc000 thru 0x0000000000004000 * reg2 expectation: shifted right 1, matching reg 1, where high bit is bit= 62 * to avoid initial allocation problems, create stack use instead of uPtr. * also avoids re-zeroing of uPtr regs * uPtr parameter for adding to reg1's converted fractional, the final resu= lt */ static uint32_t __fractional_128(_ull_t *uPtr, _ull_t reg) { uint32_t last; _ull_t sreg[15]; /* StackREGisters */ memset(sreg, 0, sizeof(sreg)); /* create regs of multiplly by 5^63 (3 max) */ __fractional(sreg, reg); /* creates sr3-sr6, affects uPtr[4] */ last =3D __fractional(&sreg[ 3], sreg[0]); if ((sreg[1].ull)) { last =3D 5; // creates uPtr[5] value /* creates sr7-sr10 */ __fractional(&sreg[ 7], sreg[1]); sreg[4].ull +=3D sreg[7].ull; sreg[5].ull +=3D sreg[8].ull; sreg[6].ull +=3D sreg[9].ull; sreg[7].ull =3D sreg[10].ull; if ((sreg[2].ull)) { /* creates sr11-sr14*/ __fractional(&sreg[11], sreg[2]); sreg[5].ull +=3D sreg[11].ull; sreg[6].ull +=3D sreg[12].ull; sreg[7].ull +=3D sreg[13].ull; uPtr[6].ull =3D sreg[14].ull; } if (sreg[7].ull >=3D _CARRY) sreg[6].ull +=3D 1ULL, sreg[7].ull -=3D _= CARRY; uPtr[5].ull =3D sreg[7].ull; while (sreg[6].ull >=3D _CARRY) sreg[5].ull +=3D 1ULL, sreg[6].ull -= =3D _CARRY; } /* sreg3-sreg8 set for addin transfer to uPtr */ uPtr[4].ull =3D sreg[6].ull; uPtr[3].ull +=3D sreg[5].ull; while (uPtr[3].ull >=3D _CARRY) uPtr[2].ull +=3D 1ULL, uPtr[3].ull -=3D = _CARRY; uPtr[2].ull +=3D sreg[4].ull; if (uPtr[2].ull >=3D _CARRY) uPtr[1].ull +=3D 1ULL, uPtr[2].ull -=3D _CA= RRY; uPtr[1].ull +=3D sreg[3].ull; if (uPtr[1].ull >=3D _CARRY) uPtr[0].ull +=3D 1ULL, uPtr[1].ull -=3D _CA= RRY; return last; } static uint32_t __positive_lowexp128(_ull_t *uPtr, unsigned pw2, _ull_t reg0, _ull_t reg1) = { uint32_t last, pw10 =3D 0; uint32_t unit_18 =3D 0; if (pw2 < 64) { /* create shifted reg0, reg1 */ /* special case: reg1 has 49bits, a shift of >=3D 49 bits (*2^49+) * means only 1 reg used for __fractional() */ uint32_t shiftb =3D 63 - pw2; if (pw2 <=3D 48) { /* lowest drop of 15 bits, max of 63 */ uPtr[0].ull =3D reg0.ull >> shiftb; /* pw2 =3D=3D 0 clears explicit '1' */ reg0.ull =3D reg0.ull << pw2, reg0.uls.hi &=3D 0x7fffffff; /* pw2 =3D=3D 0 ORs high bit explicit '0' */ reg0.ull |=3D reg1.ull >> shiftb; reg1.ull =3D reg1.ull << pw2, reg1.uls.hi &=3D 0x7fffffff; last =3D __fractional(&uPtr[1], reg0); if (reg1.ull !=3D 0) last =3D __fractional_128(&uPtr[1], reg1); } else { uint32_t hi; /* @49 shiftb =3D 14, @63 shiftb =3D 0 */ if (shiftb !=3D 0) { reg1.ull =3D (reg1.ull >> shiftb); reg1.uls.hi |=3D (reg0.uls.lo << (32 - shiftb)) >> 1; reg0.ull >>=3D shiftb; } else { /* values 8-15 possible, reducing >=3D _CARRY loop * 8 9 10 11 12 13 14 15 * actual 9, 10, 11, 12, 13, 14, 16, 17 e18 starting values */ hi =3D (uint32_t)((reg0.uss.w3 >> 12) + 1); goto hi_adjusted; } /* reg0 holds shifted whole value */ /* default action */ uPtr[0].ull =3D reg0.ull; if ((shiftb < 4) || ((shiftb =3D=3D 4) && (reg0.ull >=3D _CARRY))) { /* values 0-7 possible, reducing >=3D _CARRY loop * 0 1 2 3 4 5 6 7 * actual 1, 1, 2, 3, 4, 5, 6, 8, e18 starting values */ hi =3D (uint32_t)(reg0.uss.w3 >> 12); if ((hi)) hi_adjusted: reg0.ull -=3D hi * _CARRY; while (reg0.ull >=3D _CARRY) hi +=3D 1, reg0.ull -=3D _CARRY; uPtr[0].ull =3D (uint64_t)hi, ++uPtr; uPtr[0].ull =3D reg0.ull; pw10 =3D 18; } last =3D __fractional(&uPtr[1], reg1); unit_18 =3D 1; } } else { uint32_t n28s; _ull_t reg; /* most used 64 bit gpr */ /* 64 <=3D pw2 < 112, 1 fractional reg, multiple whole regs */ uint32_t hi =3D (uint32_t)((reg0.uss.w3 >> 12) + 1); reg0.ull -=3D hi * _CARRY; while (reg0.ull >=3D _CARRY) hi +=3D 1, reg0.ull -=3D _CARRY; uPtr[0].ull =3D (uint64_t)hi; uPtr[1].ull =3D reg0.ull; /* max 50 ((111-63) =3D 49) [0-49] */ pw2 -=3D 63; uint32_t shiftb =3D 63 - pw2; if ((n28s =3D (pw2 >=3D 28))) pw2 -=3D 28; if ((pw2)) { /* shift range [27-1] */ _ull_t scratch; /* shift [18-9] left [27-1] */ uPtr[0].uls.lo <<=3D pw2; /* shift remainder (uPtr[1], size < 10^19) left [27-1] */ reg.ull =3D divlu2r(uPtr[1].ull); scratch.ull =3D divlu2r((uint64_t)reg.uls.hi << pw2); uPtr[0].ull +=3D scratch.uls.hi; reg.ull =3D ((uint64_t)scratch.uls.lo * 1000000000U) + ((uint64_t)reg.uls.lo << pw2); /* do carry, set uPtr[1] */ if (reg.ull >=3D _CARRY) uPtr[0].ull +=3D 1ULL, reg.ull -=3D _CARRY; uPtr[1].ull =3D reg.ull; } unit_18 =3D 1; pw10 =3D 18; if (n28s) { /* shift of 28 multiples */ /* set initial limiter to current regs valid (<=3D unit_max) */ _ull_t *cap =3D &uPtr[unit_18]; _ull_t *input, *output, next; reg.ull =3D (output =3D (input =3D uPtr))->ull; output->ull =3D 0; next.ull =3D input[1].ull; unit_18++; if (reg.ull >=3D 3725290298ULL) goto L1; unit_18--; output->ull =3D reg.ull << 28; output--; do { _ull_t scratch0; uint32_t u32; reg.ull =3D next.ull; next.ull =3D (++input)[1].ull; output++; L1: scratch0.ull =3D divlu2r(reg.ull); u32 =3D (uint32_t)(((uint64_t)scratch0.uls.hi * 2305843009) >> 32) = >> 1; scratch0.ull =3D ((((uint64_t)scratch0.uls.hi << 28) - ((uint64_t)u32 * 1000000000U)) * 1000000000U) + ((uint64_t)scratch0.uls.lo << 28); if (scratch0.ull >=3D _CARRY) u32 +=3D 1, scratch0.ull -=3D _CARRY; output[1].ull =3D scratch0.ull; if ((output->ull +=3D u32) >=3D _CARRY) output->ull -=3D _CARRY, output[-1].ull +=3D 1; } while (input < cap); } pw10 +=3D ((unit_18 - 1) * 18); /* reg0 has addin to shifted whole value */ pw2 =3D 63 - shiftb; reg0.ull =3D (reg1.ull >> shiftb); reg0.ull +=3D uPtr[1].ull; if (reg0.ull >=3D _CARRY) uPtr[0].ull +=3D 1, reg0.ull -=3D _CARRY; uPtr[1].ull =3D reg0.ull; /* put fractional bits into place */ reg1.ull =3D reg1.ull << pw2, reg1.uls.hi &=3D 0x7fffffff; last =3D __fractional(&uPtr[2], reg1); last++; } return ((pw10 << 16) | ((last + unit_18) << 4) | last); } static uint32_t __positive_lowexp(_ull_t *uPtr, unsigned pw2, _ull_t reg0) { /* do less than 63 code */ uint16_t last, whole_incr, pw10; uint64_t up0; pw10 =3D (whole_incr =3D 0); if (pw2 < 31) { up0 =3D (uint64_t)(reg0.uls.hi >> (31 - pw2)); reg0.uls.hi =3D (reg0.uls.hi << (pw2 + 1)) >> (pw2 + 1); } else if (pw2 =3D=3D 31) { up0 =3D (uint64_t)reg0.uls.hi; reg0.uls.hi =3D 0; } else { uint32_t shift =3D pw2 - 32; up0 =3D reg0.ull >> (31 - shift); reg0.uls.lo =3D (reg0.uls.lo << (shift + 1)) >> (shift + 1); if ((shift > 27) || ((shift =3D=3D 27) && (up0 >=3D _CARRY))) { do uPtr[0].ull +=3D 1ULL; while ((up0 -=3D _CARRY) >=3D _CARRY); uPtr++, whole_incr++, pw10 =3D 18; } } uPtr[0].ull =3D up0; last =3D 1; if ((reg0.ull <<=3D pw2) !=3D 0) last +=3D __fractional(&uPtr[1], reg0); return ((pw10 << 16) | ((last + whole_incr) << 4) | last); } /* now grows low address to high address, or left to right numerically */ /* expects implicit/explicit 1 as 63rd bit */ /* multiplies by 2^28, 2^28 used to keep large values within 10^18 ulls */ /* addition of such values create at max 'overflow' of 1 */ /* valid division range on /28 [0-16383] */ static uint32_t __positive_exponent_fp(_ull_t *uPtr, int type) { unsigned pw2, pw10, n28s; /* power of 2, power of 10, shifts of 28 */ _ull_t reg0; /* most used 64 bit gpr */ /* track reg expansion, initial is to read 2 registers (&uPtr[unit_18])= */ uint32_t unit_18 =3D 1; /* load 'passed' data */ pw2 =3D uPtr[0].uss.w0 - 16383; /* note: retain explicit '1' */ reg0.ull =3D uPtr[1].ull; uPtr[1].ull =3D (uPtr[0].ull =3D 0); if (uPtr[2].ull !=3D 0) { _ull_t reg1; uint32_t hi; reg1.ull =3D uPtr[2].ull; uPtr[2].ull =3D 0; if (pw2 < 112) return __positive_lowexp128(uPtr, pw2, reg0, reg1); /* set up with all whole */ hi =3D (uint32_t)((reg0.uss.w3 >> 12) + 1); reg0.ull -=3D hi * _CARRY; while (reg0.ull >=3D _CARRY) hi +=3D 1, reg0.ull -=3D _CARRY; uPtr[0].ull =3D (uint64_t)hi << 49; uPtr[1].ull =3D reg1.ull >> 14; /* reg1 free to use as gpr */ reg0.ull =3D divlu2r(reg0.ull); reg1.ull =3D divlu2r((uint64_t)reg0.uls.hi << 21); uPtr[0].ull +=3D (uint64_t)reg1.uls.hi << 28; reg0.ull =3D ((uint64_t)reg1.uls.lo * 1000000000U) + ((uint64_t)reg0.uls.lo << 21); reg0.ull =3D divlu2r(reg0.ull); reg1.ull =3D divlu2r((uint64_t)reg0.uls.hi << 28); uPtr[0].ull +=3D (uint64_t)reg1.uls.hi; uPtr[1].ull +=3D ((uint64_t)reg1.uls.lo * 1000000000U) + ((uint64_t)reg0.uls.lo << 28); if (uPtr[1].ull >=3D _CARRY) uPtr[0].ull +=3D 1, uPtr[1].ull -=3D _CAR= RY; pw10 =3D 18; pw2 -=3D 112; /* uPtr[0] range [10384593717069655-5192296858534827] */ if ((pw2 -=3D ((n28s =3D (pw2 * 18725U) >> 19) * 28)) !=3D 0) { _ull_t next, *output; uint32_t u32; next.ull =3D uPtr[1].ull; reg0.ull =3D divlu2r(uPtr[0].ull); /* 2^61 / 10^9 based on << 28 forces all into higher register */ u32 =3D (uint32_t)(((uint64_t)reg0.uls.hi * 2305843009U) >> 32) >> (2= 9 - pw2); reg0.ull =3D ((((uint64_t)reg0.uls.hi << pw2) - ((uint64_t)u32 * 1000000000U)) * 1000000000U) + ((uint64_t)reg0.uls.lo << pw2); if (reg0.ull >=3D _CARRY) u32 +=3D 1, reg0.ull -=3D _CARRY; if ((u32)) { uPtr[1].ull =3D reg0.ull; uPtr[0].ull =3D (uint64_t)u32; /* creates new register, increases uPtr[unit_18] limit (3) */ /* adds pw10 +=3D 18 by return math addition of unit_18 */ unit_18++; output =3D &uPtr[1]; } else { /* no new register, still uPtr[unit_18] limit (2) */ uPtr[0].ull =3D reg0.ull; output =3D &uPtr[0]; } reg0.ull =3D divlu2r(next.ull); u32 =3D (uint32_t)(((uint64_t)reg0.uls.hi * 2305843009U) >> 32) >> (2= 9 - pw2); reg0.ull =3D ((((uint64_t)reg0.uls.hi << pw2) - ((uint64_t)u32 * 1000000000U)) * 1000000000U) + ((uint64_t)reg0.uls.lo << pw2); if (reg0.ull >=3D _CARRY) u32 +=3D 1, reg0.ull -=3D _CARRY; output[1].ull =3D reg0.ull; /* case: * w160[4] =3D 0x4092, w160[3] =3D 0x8001a8fb, w160[2] =3D 0xec8ed2= 3f; */ if ((output[0].ull +=3D u32) >=3D _CARRY) output[-1].ull +=3D 1, output[0].ull -=3D _CARRY; } } else { uint32_t hi; if (pw2 < 63) return __positive_lowexp(uPtr, pw2, reg0); /* below, all use min of 2 regs (uPtr) */ /* creates new register, still uPtr[unit_18] limit (2) */ /* values 8-15 possible, reducing >=3D _CARRY loop * 8 9 10 11 12 13 14 15 * actual 9, 10, 11, 12, 13, 14, 16, 17 e18 starting values */ hi =3D (uint32_t)((reg0.uss.w3 >> 12) + 1); reg0.ull -=3D hi * _CARRY; while (reg0.ull >=3D _CARRY) hi +=3D 1, reg0.ull -=3D _CARRY; uPtr[0].ull =3D (uint64_t)hi; uPtr[1].ull =3D reg0.ull; pw10 =3D 18; pw2 -=3D 63; if ((pw2 -=3D ((n28s =3D (pw2 * 18725U) >> 19) * 28)) !=3D 0) { /* shift range [27-1] */ _ull_t scratch; /* shift [18-9] left [27-1] */ uPtr[0].uls.lo <<=3D pw2; /* shift remainder (uPtr[1], size < 10^19) left [27-1] */ reg0.ull =3D divlu2r(uPtr[1].ull); scratch.ull =3D divlu2r((uint64_t)reg0.uls.hi << pw2); uPtr[0].ull +=3D scratch.uls.hi; reg0.ull =3D ((uint64_t)scratch.uls.lo * 1000000000U) + ((uint64_t)reg0.uls.lo << pw2); /* do carry, set uPtr[1] */ if (reg0.ull >=3D _CARRY) uPtr[0].ull +=3D 1ULL, reg0.ull -=3D _CARR= Y; uPtr[1].ull =3D reg0.ull; } } if (n28s !=3D 0) { /* shift of 28 multiples */ /* set max regs use, cap used as pointer compare limiter */ uint32_t unit_max =3D type >> 16; /* set initial limiter to current regs valid (<=3D unit_max) */ _ull_t *cap =3D &uPtr[unit_18]; do { /* start at base */ /* create loop counter, each 28 loop expects increase of new * reg created. A test used to verify if one is created. loop count= er * then used to increase both pw10, and possibly the ending of * last math reg used (cap) */ _ull_t *input, *output, next; reg0.ull =3D (output =3D (input =3D uPtr))->ull; output->ull =3D 0; next.ull =3D input[1].ull; /* critical point where addin to reg0.ull can cause overflow * and propagation into a prior ull */ unit_18++; if (reg0.ull >=3D 3725290298ULL) goto L1; unit_18--; /* shift of reg0.ull with addin remains below 10^19 */ output->ull =3D reg0.ull << 28; /* undo loop increment of first addin location, * and lcnt (new reg increase) initial assumption */ output--; do { _ull_t scratch0; uint32_t u32; reg0.ull =3D next.ull; next.ull =3D (++input)[1].ull; output++; L1: scratch0.ull =3D divlu2r(reg0.ull); u32 =3D (uint32_t)(((uint64_t)scratch0.uls.hi * 2305843009) >> 32) = >> 1; scratch0.ull =3D ((((uint64_t)scratch0.uls.hi << 28) - ((uint64_t)u32 * 1000000000U)) * 1000000000U) + ((uint64_t)scratch0.uls.lo << 28); if (scratch0.ull >=3D _CARRY) u32 +=3D 1, scratch0.ull -=3D _CARRY; output[1].ull =3D scratch0.ull; if ((output->ull +=3D u32) >=3D _CARRY) output->ull -=3D _CARRY, output[-1].ull +=3D 1; } while (input < cap); /* tail pointer update based on new reg add vs needs for precision = */ /* limit allowed to increase above cap for pw10 counter */ cap =3D &uPtr[((unit_18 > unit_max) ? unit_max : unit_18)]; } while ((--n28s)); /* convert regs used into pw10. add to pre-28 value */ } pw10 +=3D ((unit_18 - 1) * 18); return ((pw10 << 16) | ((unit_18 + 1) << 4)); } static uint32_t __negative_exponent_fp(_ull_t *uPtr, int type) { int pw2, pw10; /* power of 2, power of 10 */ uint32_t index; _ull_t reg0; pw10 =3D -18; pw2 =3D 16383 - uPtr[0].uls.lo; pw2 -=3D (index =3D ((uint32_t)pw2 =3D=3D ((IS_DBL(type)) ? 1023U : 16383= U))); /* remove explicit '1'. add back later as 1e18 if not denormal */ reg0.ull =3D uPtr[1].ull; reg0.uls.hi &=3D 0x7FFFFFFF; uPtr[1].ull =3D 0; if (uPtr[2].ull !=3D 0) { _ull_t reg1; reg1.ull =3D uPtr[2].ull; uPtr[2].ull =3D 0; __fractional(uPtr, reg0); __fractional_128(uPtr, reg1); if (index =3D=3D 0) uPtr->ull +=3D _CARRY; index =3D 8; } else if (reg0.ull !=3D 0) { /* convert fractional to base 10 */ __fractional(uPtr, reg0); /* 2 cases: denormal, normal */ if (index =3D=3D 0) uPtr->ull +=3D _CARRY; index =3D 5; } else { if (pw2 < 18) { /* straight power of 2 */ /* first 17, [-1, -17]... shift of 10^18, * 10^18/18446744073709551616 =3D * 0.0542101086242752217003726400434970855712 * 0.0542101086242752217003726400434970855712 * 2^64 =3D * 10^18 =3D 0x232830643 * 2^32 + 2808348672 */ /* no whole, all fraction */ uPtr->uls.lo =3D ((uint32_t)232830643U << (32 - pw2)) | ((uint32_t)2808348672U >> pw2); uPtr->uls.hi =3D (uint32_t)232830643U >> pw2; return (uint32_t)((-(18 << 16)) | (1 << 4) | 8); } /* there are 18 trailing zeroes to _CARRY */ uPtr[0].ull =3D _CARRY >> 18; pw2 -=3D 18; index =3D 5; } /* 5^18 * 2^2 * 2^16 =3D 10^18 */ /* each shift into new ull is multiplied by 10^18. * or 2^-16 * 10^18 =3D shift * 5^18 * 2^2 */ if (pw2 & 0xf) { /* partial 2^-16, use of shorter, quicker multiply */ /* cant use shift on ulls, when partial gcc adds extra coding, test o= f 32 */ /* use this code, smaller, faster than using unsigned long gcc code *= / uint32_t m =3D pw2 & 0xf; uint32_t shift =3D (16 - m); uint32_t shiftb =3D shift + 16; uint64_t fConst =3D 15258789062500ULL; if (index =3D=3D 8) { uPtr[7].ull =3D ((uint16_t)(uPtr[6].uls.lo << shift)) * fConst; uPtr[6].uls.lo =3D (uPtr[6].uls.hi << shiftb) | (uPtr[6].uls.lo >> m)= ; uPtr[6].uls.hi >>=3D m; uPtr[6].ull +=3D ((uint16_t)(uPtr[5].uls.lo << shift)) * fConst; uPtr[5].uls.lo =3D (uPtr[5].uls.hi << shiftb) | (uPtr[5].uls.lo >> m)= ; uPtr[5].uls.hi >>=3D m; uPtr[5].ull +=3D ((uint16_t)(uPtr[4].uls.lo << shift)) * fConst; uPtr[4].uls.lo =3D (uPtr[4].uls.hi << shiftb) | (uPtr[4].uls.lo >> m)= ; uPtr[4].uls.hi >>=3D m; uPtr[4].ull +=3D ((uint16_t)(uPtr[3].uls.lo << shift)) * fConst; } else { uPtr[4].ull =3D ((uint16_t)(uPtr[3].uls.lo << shift)) * fConst; } uPtr[3].uls.lo =3D (uPtr[3].uls.hi << shiftb) | (uPtr[3].uls.lo >> m); uPtr[3].uls.hi >>=3D m; uPtr[3].ull +=3D ((uint16_t)(uPtr[2].uls.lo << shift)) * fConst; uPtr[2].uls.lo =3D (uPtr[2].uls.hi << shiftb) | (uPtr[2].uls.lo >> m); uPtr[2].uls.hi >>=3D m; uPtr[2].ull +=3D ((uint16_t)(uPtr[1].uls.lo << shift)) * fConst; uPtr[1].uls.lo =3D (uPtr[1].uls.hi << shiftb) | (uPtr[1].uls.lo >> m); uPtr[1].uls.hi >>=3D m; uPtr[1].ull +=3D ((uint16_t)(uPtr[0].uls.lo << shift)) * fConst; uPtr[0].uls.lo =3D (uPtr[0].uls.hi << shiftb) | (uPtr[0].uls.lo >> m); uPtr[0].uls.hi >>=3D m; /* address case of denormals, below 'flag' will write pre-regs otherw= ise */ if ((uPtr[0].uls.hi | uPtr[0].uls.lo) =3D=3D 0) { pw10 -=3D 18; uPtr[0].ull =3D uPtr[1].ull, uPtr[1].ull =3D uPtr[2].ull; uPtr[2].ull =3D uPtr[3].ull, uPtr[3].ull =3D uPtr[4].ull; if (index =3D=3D 8) { uPtr[4].ull =3D uPtr[5].ull, uPtr[5].ull =3D uPtr[6].ull; uPtr[6].ull =3D uPtr[7].ull, uPtr[7].ull =3D 0; } else { uPtr[4].ull =3D 0; } } } /* 2^-16, multiples */ if ((pw2 >> 4) !=3D 0) { _ull_t *regPtr, *limit; pw2 >>=3D 4; limit =3D uPtr + (unsigned)(type >> 16); int16_t fcnt =3D 0; /* code approaches best written assembler, addressed most stalls * inotherwords, use caution rewording, huge penalties. */ do { uint64_t work; uint32_t w0; int32_t flag; regPtr =3D uPtr; w0 =3D (uint32_t)regPtr->uss.w0; fcnt +=3D (flag =3D -((regPtr->ull =3D regPtr->ull >> 16) =3D=3D 0)); work =3D (++regPtr)->ull; do { uint64_t m0 =3D w0 * 15258789062500ULL; w0 =3D (uint16_t)work; regPtr[flag].ull =3D (work =3D m0 + (work >> 16)); ++regPtr; } while ((regPtr <=3D limit) && ((work =3D regPtr->ull) !=3D 0)); regPtr[flag].ull =3D w0 * 15258789062500ULL; } while ((--pw2) !=3D 0); index =3D (int)(regPtr - uPtr) + 1; pw10 +=3D fcnt * 18; } /* only pw10 used in function as exponent of power 10 */ return ((pw10 << 16) | (index << 4) | 8); } /* enters with 'sSize =3D strPtr - sPtr;'. * s, sPtr could have been malloc'd on entry, but expected to have been * cleared, output used, prior to entry. Possibly previous real allocation, * output to 'buffer', FILE or supplied, and not free'd. * printf() supplies BUFSTRSIZE on entry, realtostr supplies 45 ulls. * BUFSTRSIZE considered large for normal requests, * 45 ulls covers all needs for double reals. * all real types generally need 8, [=C2=B1][0][.](precision)[e=C2=B1000[0]= ] * IS_TYPE_F adds exponent to count, [=C2=B1][0](exponent)[.](precision)[e= =C2=B1000[0]] * IS_TYPE_F has possiblity of grouping: * GRP_PREFIX assumes thousands, xxx,xxx,xxx,xxx etc. =3D 6 per ull * can't reuse precision 'reg's as double duty... 18 bytes / 8 byte ull * bufS(string) must come first for free'ing of allocation */ #pragma mark *** External Interfaces *** #if (!__GNUC_PREREQ(4,5)) #pragma GCC diagnostic ignored "-Wuninitialized" #endif unsigned char * realtostr(uint32_t w160[5], unsigned char **s, int type, int arg_prec, int arg_flags, unsigned decimal) { /*smash_stack();*/ int exponent, rpos; /* Rounding POSition */ uint32_t (*math_func)(_ull_t *, int); unsigned char *sPtr, dblt_sign; /* set up for percision of 792 digits, 18/ull * 44ull */ /* malloc if exponent and percision dictate, 44ull =3D 352 bytes */ /* reg: static reserved, regs: actual reg use address (stack or allocat= ed) */ _ull_t reg[45], *regs; /* supress warning due to zero jumps */ #if __GNUC_PREREQ(4,5) #pragma GCC diagnostic ignored "-Wuninitialized" #endif /* utPtr: top of stack, ubPtr: bottom of stack, ufPtr: fractional start= */ _ull_t *utPtr, *ubPtr, *ufPtr; int d0; /* Digit0 (first digit) */ #if __GNUC_PREREQ(4,5) #pragma GCC diagnostic pop #endif uint64_t a; unsigned strsz, pw2, packed_return; /*** input defaults ***/ char *decimalPtr =3D (char*)&decimal; if (arg_prec < 0) arg_prec =3D 6; /* CASE 'g': initial zero(or first digit) is significant digit */ if ((arg_prec) && (IS_TYPE_G(type))) --arg_prec; /*** sign ***/ if ((w160[4] & 0x8000) !=3D 0) dblt_sign =3D '-'; else { dblt_sign =3D 0; if ((arg_flags & (SIGN_PREFIX | SPC_PREFIX)) !=3D 0) dblt_sign =3D (arg_flags & SPC_PREFIX) ? ' ' : '+'; } /*** exponent ***/ /* extra math used for limiting register updating in math units. Also * allows allocation size when precision exceeds static buffersize */ if ((exponent =3D w160[4] & 0x7FFF) =3D=3D 0x7FFF) { /* inf or nan, nan(n-char-sequence) */ *(sPtr =3D *s) =3D dblt_sign, sPtr +=3D (dblt_sign !=3D 0); return __n_char_sequence(w160, sPtr, (type | ((arg_flags & ALT_PREFIX) = << 1))); } if (exponent < 0x3FFF) { uint32_t reg_min =3D (IS_QUAD(type)) ? 7 : 4; math_func =3D __negative_exponent_fp; pw2 =3D reg_min; /* these 'need' pw10 pre-calculations for extra decisions */ /* !E can result in '0', no caculation needed, * >> 4 =3D=3D 0 results in use of default, * default pw2 =3D=3D 4 covers precision < 64 */ /* XXX with !e80, does pw2 cover precision < (NEG_REG_BITS + 1)? */ if ( (!IS_TYPE_E(type)) || ( (((uint32_t)(0x3FFF - exponent) >> 4) !=3D 0) && (arg_prec >=3D 6= 3) ) ) { pw2 =3D 0x3FFF - exponent; /* number of prefixing zeroes */ pw2 =3D (((uint32_t)(((uint64_t)((uint32_t)pw2) * (uint32_t)2585827973U) >> 32)) >> 1); /* check for all zeroes on a 'f' request */ if ( (!IS_TYPE_E(type)) && ((unsigned)arg_prec < pw2) ) { if (pw2 > reg_min) { pw2 =3D reg_min; exponent =3D 0; *(uint64_t*)&w160[2] =3D 0; /* XXX do I need? */ *(uint64_t*)&w160[0] =3D 0; if (IS_TYPE_G(type)) type |=3D TYPE_F; type |=3D TYPE_ZERO; /* must not see as denormal for doubles */ /* remove 'double' label, keep or set 'long double' */ type &=3D ~(UPPER_TYPE << 8); } else { pw2 =3D reg_min; } } else { /* subtract zeroes from significant digits */ uint32_t bits =3D (IS_QUAD(type)) ? 112 : 63; pw2 =3D (0x3FFF - exponent) + bits - pw2; if (pw2 > (unsigned)arg_prec) pw2 =3D arg_prec; /* divide by 18 for number ulls needed */ if (pw2 <=3D 38) pw2 =3D reg_min; else pw2 =3D reg_min + (uint32_t)((((uint64_t)(uint32_t)(pw2 - 38)) * (uint32_t)238609295U) >> = 32); } } } else { math_func =3D __positive_exponent_fp; pw2 =3D 2; /* when 'f', need maximum value */ if ( ((exponent >=3D 0x40ce) && (arg_prec >=3D 27)) || (IS_TYPE_F(type)= ) ) { pw2 =3D exponent - 0x3FFF; pw2 =3D (((uint32_t)(((uint64_t)((uint32_t)pw2) * 2585827973U) >> 32)) >> 1= ); /* 'e' doesn't need full precision */ if ((IS_TYPE_E(type)) && ((unsigned)arg_prec < pw2)) pw2 =3D arg_prec; /* divide by 18 for number ulls needed */ pw2 =3D 2 + (uint32_t)(((uint64_t)((uint32_t)pw2) * 252645135U) >> 32= ); } } /*** allocation ***/ /* allocate reg here based on exponent/precision, see above notes */ /* group allocation, fprintf() to free, returning we allocating 's' */ strsz =3D arg_prec + 8; type |=3D pw2 << 16; /* type F has special requirements, if positive exponent */ if ((type << 31) & ~(exponent - 0x3FFF)) { /* grouping assumption of thousands: 1,000,000 */ /* add in exponent: [=C2=B1][0](exponent)[.](precision)[e=C2=B1000[0]= ] */ /* assume all pw2 as 'whole' */ /* more accuracy, for less allocations, is faster */ uint32_t triplets =3D (exponent - 0x3FFF); strsz +=3D 1 + (((uint32_t)(((uint64_t)((uint32_t)triplets) * 2585827973U) >> = 32)) >> 1); if ((arg_flags & GRP_PREFIX) !=3D 0) strsz +=3D ((uint16_t)21845 * (uint32_t)triplets) >> 16; } regs =3D reg; if ((pw2 > 44) || (strsz > BUFSTRSIZE)) { /* if either static exceeds limit, allocate both */ strsz =3D (strsz + ((1 << 6) - 1)) & ~((1 << 6) - 1); /* test for non-static string buffer */ if (arg_flags & ALLOC_BUF) free(*s); /* assignment to sPtr, sPtr is one tested and free'd in fprintf() */ if ((sPtr =3D malloc((size_t)(strsz + ((pw2 + 1) * sizeof(_ull_t))))) = =3D=3D NULL) return sPtr; *s =3D sPtr; regs =3D (_ull_t*)(sPtr + strsz); } /* zero regs, no longer part of math units */ do regs[pw2].uls.hi =3D 0, regs[pw2].uls.lo =3D 0; while ((--pw2)); /*** fractional load ***/ /* assign 'fraction', 'exponent' info for passage to math units */ /* w160[0,1] used on requested quadruple, or long double 128bit (not e8= 0) */ regs[0].uls.lo =3D exponent; regs[1].ull =3D *(uint64_t*)&w160[2]; regs[2].ull =3D *(uint64_t*)&w160[0]; if ((regs[1].ull | regs[2].ull) =3D=3D 0) { exponent =3D 0; if (IS_TYPE_G(type)) type |=3D TYPE_F; type |=3D TYPE_ZERO; goto skip_round; } /*** math ***/ packed_return =3D (*math_func)(regs, type); /* transfer info returned into addresses, exponent */ ubPtr =3D (utPtr =3D regs) + ((packed_return >> 4) & 0x3FF); exponent =3D (d0 =3D _first_digit(utPtr->ull)); exponent +=3D (int32_t)packed_return >> 16; packed_return &=3D 0x0F; ufPtr =3D ((packed_return =3D=3D 0) ? ubPtr : ((packed_return =3D=3D 0x08) ? regs : ubPtr - (packed_return - 1)= )); /*** rounding ***/ /* assign position 1 pass display precision */ rpos =3D arg_prec + 1; if (IS_TYPE_F(type)) { /* this excludes positive exponents, they will not jump to skip */ if (arg_prec < ~exponent) goto skip_round; /* positive exponents must count pre-radix as precision */ if (~exponent >> 31) arg_prec +=3D exponent; /* when negative, remove pre-zero count, positive add in whole portio= n */ rpos +=3D exponent; } /* rounding position <=3D first digit position, has front end zeroes */ if ((rpos -=3D d0) > 0) { _ull_t *aPtr =3D utPtr; /* test rounding for second ull or higher, assume not large arg_prec = */ while (((++aPtr) < ubPtr) && (rpos > 18)) rpos -=3D 18; /* verify position has digits to round */ if (aPtr < ubPtr) { /* set last ull, if rpos =3D=3D 1 then exclude current aPtr */ ubPtr =3D aPtr + (rpos !=3D 1); if (ufPtr > ubPtr) ufPtr =3D ubPtr; const uint64_t *fc =3D _fcarry; if ((aPtr->ull +=3D fc[18 - rpos]) >=3D _CARRY) { uint64_t t =3D aPtr[-1].ull + 1ULL; /* carry forced propagation, makes ull =3D=3D 0 */ (ubPtr =3D aPtr)->ull =3D 0; (--aPtr)->ull =3D t; if ((aPtr =3D=3D utPtr) && (t >=3D fc[19 + d0])) { exponent++; /* type 'f' has increase in precision if positive exponent */ /* condition of =3D=3D 18 never occurs */ arg_prec +=3D (type & ((uint32_t)~exponent >> 31)); ++d0; } } } } else { /* first ull rounding */ const uint64_t *fc =3D _fcarry; /* limit reads of ulls to next register */ ubPtr =3D utPtr + 1; if (ufPtr > ubPtr) ufPtr =3D ubPtr; /* add '5' to rounding location */ a =3D utPtr->ull + fc[-rpos]; utPtr->ull =3D a; /* test for propagation to new first digit */ if (a >=3D fc[19 + d0]) { exponent++; if ((++d0) =3D=3D 18) { /* rare case. inserted here to avoid alter of ulltobuffer for odd= case * and testing all cases for this condition. */ /* occurs @ 1e18, 1e-18 9.9999950000000026e+17 first precision 6 9.9999999999999987e+17 last 9.9999950000000012e-19 first precision 6 9.9999999999999988e-19 last */ *(sPtr =3D *s) =3D dblt_sign, sPtr +=3D (dblt_sign !=3D 0); /* type f and negative exponent */ if (((uint32_t)exponent >> 31) & type) { unsigned j, k; *sPtr++ =3D '0'; *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; j =3D 1 + (decimalPtr[1] !=3D 0); sPtr +=3D j, decimalPtr +=3D j; if ((j =3D -exponent - 1)) { if (j >=3D (unsigned)arg_prec) goto L24; arg_prec -=3D j; k =3D j & 3; if ((j =3D (int)((unsigned)j >> 2))) do *(uint32_t*)sPtr =3D 0x30303030U, sPtr +=3D 4; while ((-= -j)); *(uint32_t*)sPtr =3D 0x30303030U; sPtr +=3D k; } *sPtr++ =3D '1'; if ((arg_prec -=3D 1) > 0) goto L245; return sPtr; } *sPtr++ =3D '1'; if ( (IS_TYPE_G(type)) && (((unsigned)exponent) <=3D (unsigned)arg_= prec) ) type |=3D TYPE_F; if ( (arg_prec) && (!(IS_TYPE_G(type)) || (arg_flags & ALT_PREFIX))= ) { unsigned dflag =3D (decimalPtr[1] !=3D 0); *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; sPtr +=3D dflag + 1, decimalPtr +=3D dflag + 1; goto L245; } goto L25; } } } /*** adjust for entered type... g ***/ /* rounding can adjust exponent, so now can test if conditions met */ if (IS_TYPE_G(type)) { /* this should exclude negatives */ if (((unsigned)exponent) <=3D (unsigned)arg_prec) goto L23; if (((unsigned)(exponent + 4)) <=3D 4) { arg_prec -=3D exponent; L23: type |=3D TYPE_F; } } /*** output section ***/ skip_round: *(sPtr =3D *s) =3D dblt_sign, sPtr +=3D (dblt_sign !=3D 0); /* type f and negative exponent, */ if ((((uint32_t)exponent >> 31) & type) | (IS_TYPE_ZERO(type))) { unsigned j, k; *sPtr++ =3D '0'; *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; j =3D 1 + (decimalPtr[1] !=3D 0); if (!arg_prec) { if (IS_TYPE_E(type)) goto L11; if (arg_flags & ALT_PREFIX) return (sPtr + j); return sPtr; } sPtr +=3D j, decimalPtr +=3D j; if ((j =3D -exponent - 1)) { if (j >=3D (unsigned)arg_prec) goto L24; arg_prec -=3D j; k =3D j & 3; if ((j =3D (int)((unsigned)j >> 2))) do *(uint32_t*)sPtr =3D 0x30303030U, sPtr +=3D 4; while ((--j)); *(uint32_t*)sPtr =3D 0x30303030U; sPtr +=3D k; } } { /* first ull if utPtr =3D=3D ufPtr then all digits are part of precision. d0 does n= ot count digits but returns log10... log10(1-9) =3D 0 if 'f', data starts at sPtr. Otherwise starts sPtr + size of decimal. will write on last char of decimal, which then that char copies to sPt= r then overwritten by decimal. (!IS_TYPE_F(type)) stringPtr =3D sPtr + decimalsz To determine if first digit is used as precision use (decimalPtr =3D=3D (char*)&decimal) Two-byte decimal determined by (((char*)&decimal)[1] !=3D 0) */ /* deal with d0 */ /* _fcarry multiplier based on d0, then can adjust d0 */ const uint64_t *mConst =3D &_fcarry[35 - d0]; if ((d0 +=3D ((utPtr =3D=3D ufPtr) & IS_TYPE_F(type))) > arg_prec) d0 = =3D arg_prec; unsigned char *stringPtr =3D sPtr; if (!IS_TYPE_F(type)) stringPtr +=3D 1 + (((char*)&decimal)[1] !=3D 0); /* traps d0 =3D [0,9] or all ull multipliers */ if (utPtr->ull < 1000000000) { _ull_t b; /* buf output use pre-incr (++) */ unsigned char *buf =3D stringPtr - 1; uint32_t f =3D (uint32_t)((uint64_t)utPtr->uls.lo * (uint32_t)mConst[= (-9)]); uint32_t n =3D d0 + (decimalPtr =3D=3D (char*)&decimal); b.ull =3D (uint64_t)(f << 2) * 2882303762U; goto L2; do { b.uls.hi =3D (((b.uls.hi << 4) >> 4) * 10) + (b.uls.lo >> 28); b.uls.lo <<=3D 4; L2: *++buf =3D (b.uls.hi >> 28) + '0'; } while ((--n)); } else { _ulltobuffer_nlr((divlu2r(utPtr->ull * (uint32_t)*mConst)), stringPtr= , (d0 + (decimalPtr =3D=3D (char*)&de= cimal))); } /* type 'e' requires radix after first digit */ if (!IS_TYPE_F(type)) { /* move first digit to first spot */ *sPtr =3D *stringPtr; /* can not use short write here, * will overwrite any digits following if sz =3D=3D 1 */ *++sPtr =3D *decimalPtr; if (*++decimalPtr !=3D 0) *++sPtr =3D *decimalPtr, decimalPtr++; /* a must for if arg_prec */ sPtr++; if (!arg_prec) { if (arg_flags & ALT_PREFIX) sPtr +=3D (((char*)&decimal)[1] !=3D 0); else sPtr -=3D 1 + (((char*)&decimal)[1] !=3D 0); goto L11; } } sPtr +=3D d0 + (decimalPtr =3D=3D (char*)&decimal); if ((arg_prec -=3D d0) =3D=3D 0) goto L25; } do { /* case where has whole/fraction */ if ((++utPtr) >=3D ufPtr) { if (decimalPtr =3D=3D (char*)&decimal) { /* this also used as delimiter for whole part, stop point for * trailing zero removal */ unsigned dflag =3D (decimalPtr[1] !=3D 0); *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; sPtr +=3D dflag + 1, decimalPtr +=3D dflag + 1; } /* set to last ull, if not already at */ if (ufPtr =3D=3D ubPtr) break; ufPtr =3D ubPtr; } if ((arg_prec - 18) < 0) { _ulltobuffer_nlr(divlu2r(utPtr->ull), sPtr, arg_prec); sPtr +=3D arg_prec; goto L25; } arg_prec -=3D 18; _ulltobuffer_18lr(divlu2r(utPtr->ull), sPtr); sPtr +=3D 18; } while (1); L24: if ( (arg_prec) && (!(IS_TYPE_G(type)) || (arg_flags & ALT_PREFIX)) ) { /* safe because of added room for exponent in allocation */ unsigned k; L245: k =3D arg_prec & 3; arg_prec =3D (int)((unsigned)arg_prec >> 2); *(uint32_t*)sPtr =3D 0x30303030U; if (arg_prec) do sPtr +=3D 4, *(uint32_t*)sPtr =3D 0x30303030U; while ((--arg_pre= c)); sPtr +=3D k; } L25: /* case 'g', remove any trailing s */ if (IS_TYPE_E(type)) goto L11; if (!(arg_flags & ALT_PREFIX)) { /* test only 'last' byte of decimal point. * '=3D=3D 0' only occurs if written. decimal is at max 2 bytes */ if (*decimalPtr =3D=3D 0) { unsigned char ch; ch =3D *--sPtr; if (IS_TYPE_G(type)) while (ch =3D=3D '0') ch =3D *--sPtr; /* decimal point is at max 2 bytes */ if (ch !=3D (unsigned char)decimalPtr[-1]) ++sPtr; else if (&decimalPtr[-1] !=3D (char*)&decimal) --sPtr; } } else if (*decimalPtr !=3D 0) { unsigned dflag =3D (decimalPtr[1] !=3D 0); *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; sPtr +=3D dflag + 1, decimalPtr +=3D dflag + 1; } /*** exponent output ***/ if (IS_TYPE_F(type)) return sPtr; L11: /* converts g,G to e,E and adds sign (portable character set) */ *(uint16_t*)sPtr =3D cconst16_t(('e' ^ ((type & 0x20) ^ UPPER_TYPE)), ('+' - ((exponent >> 31) << 1))); /* absolute value */ exponent =3D (exponent ^ (exponent >> 31)) - (exponent >> 31); if (exponent < 10) { /* exponent output reqires 2 digit display */ *(uint16_t*)(&sPtr[2]) =3D cconst16_t('0',('0' + exponent)); return (sPtr + 4); } /* quadruple, long doubles range [10,4966] */ uint32_t ecd =3D (uint32_t)(((uint64_t)((uint32_t)exponent << 19) * 2199023989U) >> 32)= ; ++sPtr; if (exponent >=3D 1000) *++sPtr =3D (ecd >> 28) + '0'; ecd =3D ((ecd << 4) >> 4) * 10; if (exponent >=3D 100) *++sPtr =3D (ecd >> 28) + '0'; sPtr[1] =3D ((ecd =3D (((ecd << 4) >> 4) * 10)) >> 28) + '0'; sPtr[2] =3D ((((ecd << 4) >> 4) * 10) >> 28) + '0'; return (sPtr + 3); } /* "[=C2=B1]0xh.hhhhp=C2=B1d" * normalized: h. force to high nibble for higher precision * denornalize: .h force highest nibble as above * output precision based on 'g' style */ unsigned char * realtohex(uint32_t w160[5], unsigned char **s, int type, int arg_prec, int arg_flags, unsigned decimal) { unsigned char *sPtr =3D *s; int exponent; uint64_t frac; /*** input defaults ***/ char *decimalPtr =3D (char*)&decimal; if (arg_prec < 0) { arg_prec =3D (IS_DBL(type)) ? 13 : 15; type |=3D TYPE_G; } /* determine upper of lower case output */ const char *_outvalues =3D _const_outStr; if (IS_UPPER_TYPE(type)) _outvalues +=3D 32; /*** sign && exponent ***/ exponent =3D w160[4] & 0x7FFF; if ((w160[4] & 0x8000) !=3D 0) *sPtr++ =3D '-'; else { if ((arg_flags & (SIGN_PREFIX | SPC_PREFIX)) !=3D 0) *sPtr++ =3D (arg_flags & SPC_PREFIX) ? ' ' : '+'; } if (exponent =3D=3D 0x7FFF) return __n_char_sequence(w160, sPtr, (type | ((arg_flags & ALT_PREFIX) << 1)))= ; exponent -=3D 16383; /*** prefix 0x or 0X ***/ *((uint16_t *)sPtr) =3D ((uint16_t *)_outvalues)[12]; sPtr +=3D 2; /*** fractional load ***/ if ((frac =3D *(uint64_t*)&w160[2]) =3D=3D 0) { exponent =3D 0; *sPtr++ =3D '0'; if (arg_flags & ALT_PREFIX) { unsigned dflag =3D (decimalPtr[1] !=3D 0); *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; sPtr +=3D dflag + 1, decimalPtr +=3D dflag + 1; } goto out_exp; } /* alter exponent by 'shift' amount used to normalize */ if ((frac & 0x8000000000000000) !=3D 0) exponent -=3D 3; /*** rounding ***/ if (arg_prec < 15) { uint64_t rmask =3D 0x0FFFFFFFFFFFFFFFULL >> (arg_prec << 2); uint64_t radd =3D (rmask + 1) >> 1; /* not sure need to mask result */ if ((frac & rmask) > radd) frac =3D (frac + radd) & ~rmask; } /* rotate high nibble to low nibble */ uint32_t snfct =3D 0; if ((frac & 0x8000000000000000ULL) =3D=3D 0) { frac <<=3D 1; snfct =3D 1; *sPtr++ =3D '0'; } else { frac =3D (frac >> 60) | (frac << 4); *sPtr++ =3D _outvalues[((uint32_t)frac & 0x0000000F)]; } /* decimal point */ unsigned dflag =3D (decimalPtr[1] !=3D 0); *(uint16_t*)sPtr =3D *(uint16_t*)decimalPtr; sPtr +=3D dflag + 1, decimalPtr +=3D dflag + 1; /* precision based on type, adjust if denormal */ dflag =3D (arg_prec < ((IS_DBL(type)) ? 13 : 15)) ? arg_prec : ((IS_DBL(type)) ? 13 : 15); arg_prec -=3D dflag; if (snfct =3D=3D 1) dflag++; do { /* rotate into position */ frac =3D (frac >> 60) | (frac << 4); *sPtr++ =3D _outvalues[((uint32_t)frac & 0x0000000F)]; } while ((--dflag)); if (IS_TYPE_G(type)) { unsigned char ch =3D *--sPtr; while (ch =3D=3D '0') ch =3D *--sPtr; if (ch !=3D (unsigned char)decimalPtr[-1]) ++sPtr; else if (&decimalPtr[-1] !=3D (char*)&decimal) --sPtr; } else if (arg_prec) { dflag =3D arg_prec & 3; arg_prec =3D (int)((unsigned)arg_prec >> 2); *(uint32_t*)sPtr =3D 0x30303030U; if (arg_prec) do sPtr +=3D 4, *(uint32_t*)sPtr =3D 0x30303030U; whil= e ((--arg_prec)); sPtr +=3D dflag; } /*** exponent output (p=C2=B1d p[-16494,+16383]) ***/ out_exp: *sPtr++ =3D _outvalues[19]; /* based on portable character set */ *sPtr++ =3D '+' - ((exponent >> 31) << 1); /* absolute value */ exponent =3D (exponent ^ (exponent >> 31)) - (exponent >> 31); if (exponent < 10) { *sPtr =3D exponent + '0'; return (sPtr + 1); } /* divide exponent by 10000 3518437209 */ uint32_t ecd =3D (uint32_t)(((uint64_t)(uint32_t)exponent * 3518437209U) = >> 17); ecd +=3D 1 << 14; if (exponent >=3D 10000) *sPtr++ =3D (ecd >> 28) + '0'; ecd =3D ((ecd << 4) >> 4) * 10; if (exponent >=3D 1000) *sPtr++ =3D (ecd >> 28) + '0'; ecd =3D ((ecd << 4) >> 4) * 10; if (exponent >=3D 100) *sPtr++ =3D (ecd >> 28) + '0'; sPtr[0] =3D ((ecd =3D (((ecd << 4) >> 4) * 10)) >> 28) + '0'; sPtr[1] =3D ((((ecd << 4) >> 4) * 10) >> 28) + '0'; return (sPtr + 2); } #pragma mark *** Utility External Interfaces *** /* need different long double definitions 106,64,53 */ long double qtold(quadruple qdbl) { octle iqdbl; unsigned char *cdbl =3D (unsigned char*)&iqdbl; iqdbl =3D qdbl; #if (__LDBL_MANT_DIG__ =3D=3D 64) iqdbl.ulls.ull0 =3D (iqdbl.ulls.ull1 << 15) | (iqdbl.uss.us3 >> 1); iqdbl.ulls.ull1 >>=3D 48; iqdbl.uls.ul1 =3D (((iqdbl.uss.us4 & 0x7fff) !=3D 0) << 31) | (iqdbl.uls.ul1 & 0x7fffffff); #else int16_t sign =3D (idbl.uss.us7 & 0x8000) >> 4; int16_t e16 =3D (idbl.uss.us7 & 0x7fff) - 16383 + 1023; iqdbl.uss.us7 =3D 0; iqdbl.ulls.ull1 <<=3D 4; iqdbl.uss.us4 |=3D iqdbl.uss.us3 >> 12; /* double-double and double */ #if (__LDBL_MANT_DIG__ =3D=3D 106) /* assumes low double exponent different by 2^-52, keeps sign * of course this doesnt make sense on denormals * thus no extra precision, low double has to be zero */ /* mantissas into correct positions */ iqdbl.ulls.ull0 >>=3D 8; iqdbl.uss.us3 |=3D ((e16 - 52) | sign) << 4; if (e16 <=3D 52) { iqdbl.ulls.ull0 =3D 0; /* unknown if benefits hardware when 0 */ iqdbl.uss.us3 =3D sign << 4; } #endif if (e16 <=3D 0) { if (e16 <=3D -52) iqdbl.ulls.ull1 =3D 0; else iqdbl.ulls.ull1 >>=3D (12 - e16); e16 =3D 0; } else if (e16 > 2047) { /* set to inf */ iqdbl.ulls.ull1 =3D 0; e16 =3D 2047; } iqdbl.uss.us7 |=3D (e16 | sign) << 4; #endif return (*(long double*)cdbl); } double ldtod(long double ldbl) { unsigned char *wc =3D (unsigned char *)&ldbl; #if ((__LDBL_MANT_DIG__ =3D=3D 53) || (__LDBL_MANT_DIG__ =3D=3D 106)) /* ppc double-double XXX untested, just drop lower 64, big endian */ /* wc will point to highest byte of 128 bit */ return (*(double*)wc); #else /* since padding occurs, non gpr registers tend towards 128bit use: */ octle ildbl; int16_t e16; int16_t sign; ildbl =3D *(octle*)wc; wc =3D (unsigned char *)&ildbl; /* 3 known, use IEEE-754 16bit sign/exp, biased */ #if (__LDBL_MANT_DIG__ > 65) /* trap binary128 */ e16 =3D ildbl.uss.us7; ildbl.ulls.ull0 =3D ildbl.ulls.ull1 >> 12; #else /* trap binary64, e80 (64/63 mantissas) */ e16 =3D ildbl.uss.us4; #if (__LDBL_MANT_DIG__ =3D=3D 64) ildbl.uss.us3 &=3D 0x7fff; ildbl.ulls.ull0 >>=3D 11; #else ildbl.ulls.ull0 =3D (((ildbl.ulls.ull1 << 4) | (uint64_t)(ildbl.uss.us3 >> 12))) >> 12; #endif #endif sign =3D (e16 & 0x8000) >> 4; e16 =3D (e16 & 0x7FFF) - 16383 + 1023; if (e16 <=3D 0) { if (e16 <=3D -52) ildbl.ulls.ull0 =3D 0; else ildbl.ulls.ull0 >>=3D (12 - e16); e16 =3D 0; } else if (e16 > 2047) { /* set to inf */ ildbl.ulls.ull0 =3D 0; e16 =3D 2047; } ildbl.uss.us3 |=3D (e16 | sign) << 4; return (*(double*)wc); #endif } int main(void) { phx_printf("%.43f\n", 0x1.52f8a8e32e982p-140); phx_printf("%g\n", 355.0/113.0); return 0; } = --=-es46Z1o6jxz9SwEeQxX6--