1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993 |
- /* Definitions for the floating point types.
- This file is part of khipu.
- khipu is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 3 of the License, or
- (at your option) any later version.
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU Lesser General Public License for more details.
- You should have received a copy of the GNU Lesser General Public License
- along with this program. If not, see <https://www.gnu.org/licenses/>. */
- #include <cmath>
- #include <climits>
- #include <cstdio>
- #include <limits>
- #include <cstdlib>
- #include "floatp.hpp"
- #include "memory.hpp"
- #include "stream.hpp"
- #include "integer.hpp"
- #include "bvector.hpp"
- #include "io.hpp"
- KP_DECLS_BEGIN
- bool fnan_p (double val)
- {
- return (std::fpclassify (val) == FP_NAN);
- }
- int finf_p (double val)
- {
- return (std::fpclassify (val) != FP_INFINITE ? 0 : (val < 0 ? -1 : 1));
- }
- static int
- fclass_aux (double val, int& sign)
- {
- int ret = std::fpclassify (val);
- switch (ret)
- {
- case FP_INFINITE:
- sign = val < 0 ? -1 : 1;
- break;
- case FP_NAN:
- sign = 0;
- break; // We don't consider -NaN
- default:
- // Assume it's a normalised number.
- ret = FP_NORMAL;
- sign = std::signbit (val);
- }
- return (ret);
- }
- static double
- get_dbl (object obj, int& sign, int& cls)
- {
- double val = as_fltobj(obj)->val;
- if (kp_unlikely (obj == FLT_PINF || obj == FLT_NINF))
- {
- cls = FP_INFINITE;
- #ifdef KP_ARCH_WIDE
- sign = varobj_sign (obj);
- #else
- sign = obj == FLT_NINF;
- #endif
- }
- else if (kp_unlikely (obj == FLT_QNAN))
- {
- cls = FP_NAN;
- sign = 0;
- }
- else
- {
- sign = std::signbit (val);
- #ifdef KP_ARCH_WIDE
- if (varobj_sign (obj))
- {
- sign ^= 1;
- val = -val;
- }
- #endif
- cls = FP_NORMAL;
- }
- return (val);
- }
- #ifdef KP_ARCH_WIDE
- # define make_bigfloat(ptr, sign) \
- ptrtype ((ptr), typecode::BIGFLOAT) | (sign ? SIGN_BIT : 0)
- # define F_ABS(val) val
- #else
- double as_float (object obj)
- {
- int sgn, cls;
- return (get_dbl (obj, sgn, cls));
- }
- static inline object
- make_bigfloat (bigfloat *fp, int sign)
- {
- if (sign)
- fp->len = -fp->len;
- return (fp->as_obj ());
- }
- # define F_ABS abs
- #endif
- object fltobj::make (interpreter *interp, double val)
- {
- fltobj *retp = alloch<fltobj> (TYPE_SHIFT + 1);
- retp->val = val;
- retp->vo_full |= FLAGS_CONST;
- interp->alval = retp->as_obj ();
- gc_register (interp, retp);
- return (interp->alval);
- }
- struct fake_bigfloat : public local_varobj<bigfloat>
- {
- limb_t data[LNDBL_SIZE];
- fake_bigfloat (double x)
- {
- this->local_init (this->data, KP_NELEM (this->data), 0);
- this->expo = dbltoui (fabs (x), this->data, this->len);
- }
- object get_obj (int sign)
- {
- #ifdef KP_ARCH_WIDE
- return (this->as_obj () | (sign ? SIGN_BIT : 0));
- #else
- if (sign)
- this->len = -this->len;
- return (this->as_obj ());
- #endif
- }
- };
- template <typename F1, typename F2>
- static result<object> ff_op (interpreter *interp, object x, object y,
- int y_sign, F1 fn, F2 alt)
- {
- int s1, s2, c1, c2, rs;
- double d1 = get_dbl (x, s1, c1);
- double d2 = get_dbl (y, s2, c2);
- double ret = KP_TRY (fn (interp, d1, d2));
- int rc = fclass_aux (ret, rs);
- if (kp_likely (rc == FP_NORMAL))
- kp_return (fltobj::make (interp, ret));
- else if ((c1 | c2) != FP_NORMAL)
- {
- if (c1 == FP_NORMAL)
- kp_return (y_sign ? neg_f (interp, y) : y);
- else if (c2 == FP_NORMAL)
- kp_return (x);
- kp_return (c1 != c2 || s1 != (s2 ^ y_sign) ? FLT_QNAN : x);
- }
- fake_bigfloat tx (d1), ty (d2);
- return (alt (interp, make_bigfloat (&tx, s1), make_bigfloat (&ty, s2)));
- }
- object add_ff (interpreter *interp, object x, object y)
- {
- return (deref (ff_op (interp, x, y, 0,
- [](void *, double a, double b)
- { return (a + b); }, add_FF)));
- }
- object sub_ff (interpreter *interp, object x, object y)
- {
- return (deref (ff_op (interp, x, y, 1,
- [](void *, double a, double b)
- { return (a - b); }, sub_FF)));
- }
- object mul_ff (interpreter *interp, object x, object y)
- {
- return (deref (ff_op (interp, x, y, 0,
- [](void *, double a, double b)
- { return (a * b); }, mul_FF)));
- }
- result<object> div_ff (interpreter *interp, object x, object y)
- {
- return (ff_op (interp, x, y, 0,
- [](interpreter *ip, double a, double b) -> result<double>
- {
- if (b == 0)
- return (ip->raise ("arith-error",
- "division by zero"));
- return (a / b);
- }, div_FF));
- }
- object neg_f (interpreter *interp, object obj)
- {
- int sign, cls;
- double val = get_dbl (obj, sign, cls);
- if (cls == FP_NAN || val == 0)
- kp_return (obj);
- else if (cls == FP_INFINITE)
- kp_return (obj == FLT_PINF ? FLT_NINF : FLT_PINF);
- #ifdef KP_ARCH_WIDE
- kp_return (obj ^ SIGN_BIT);
- #else
- kp_return (fltobj::make (interp, -val));
- #endif
- }
- object abs_f (interpreter *interp, object obj)
- {
- if (obj == FLT_QNAN)
- kp_return (obj);
- int sign, cls;
- double val = get_dbl (obj, sign, cls);
- if (sign != 0)
- #ifdef KP_ARCH_WIDE
- kp_return (obj ^ SIGN_BIT);
- #else
- kp_return (fltobj::make (interp, -val));
- #endif
- (void)val;
- kp_return (obj);
- }
- result<object> mod_ff (interpreter *interp, object x, object y)
- {
- int s1, c1, s2, c2;
- double v1 = get_dbl (x, s1, c1);
- double v2 = get_dbl (y, s2, c2);
- if (v2 == 0)
- return (interp->raise ("arith-error", "modulo by zero"));
- else if (c1 == FP_NAN || c2 == FP_NAN)
- kp_return (FLT_QNAN);
- kp_return (fltobj::make (interp, fmod (v1, v2)));
- }
- int cmp_ff (interpreter *interp, object x, object y)
- {
- int s1, c1, s2, c2;
- double v1 = get_dbl (x, s1, c1);
- double v2 = get_dbl (y, s2, c2);
- if (c1 == FP_NAN || c2 == FP_NAN)
- return (1); // Any non-zero value will do.
- return (v1 < v2 ? -1 : (v1 > v2 ? +1 : 0));
- }
- bool eq_ff (interpreter *interp, object x, object y)
- {
- int s1, c1, s2, c2;
- double v1 = get_dbl (x, s1, c1);
- double v2 = get_dbl (y, s2, c2);
- if ((c1 | c2) != FP_NORMAL)
- return (false);
- return (v1 == v2);
- }
- uint32_t hash_f (interpreter *interp, object obj)
- {
- int sign, cls;
- double i, dv = get_dbl (obj, sign, cls);
- if (modf (dv, &i) != 0)
- {
- uint32_t ret = hashbuf (&dv, sizeof (dv));
- return (sign ? ~ret : ret);
- }
- int64_t iv = (int64_t)dv;
- if (fitsfixint_p (iv))
- return (hash_i (interp, fixint (iv)));
- uival uv;
- uv.qv = abs (iv);
- limb_t data[] = { uv.limbs.lo, uv.limbs.hi };
- local_varobj<bigint> bi;
- bi.local_init (data, 1 + (data[1] != 0));
- #ifdef KP_ARCH_WIDE
- return (hash_I (interp, bi.as_obj () | (iv < 0 ? SIGN_BIT : 0)));
- #else
- if (iv < 0)
- bi.len = -bi.len;
- return (hash_I (interp, bi.as_obj ()));
- #endif
- }
- static result<int>
- pad_writef (interpreter *interp, stream *strm, int ch,
- int width, int lv, int flags)
- {
- if ((flags & (io_info::FLG_LJUST | io_info::FLG_ZERO)) || lv >= width)
- return (0);
- int rv = 0;
- char buf[256];
- lv = width - lv;
- memset (buf, ch, min (lv, (int)sizeof (buf)));
- for (; lv >= (int)sizeof (buf); lv -= sizeof (buf))
- { rv += KP_TRY (strm->write (interp, buf, sizeof (buf))); }
- rv += KP_TRY (strm->write (interp, buf, lv));
- return (rv);
- }
- static result<int64_t>
- write_flt (interpreter *interp, stream *strm, double dbl, int width,
- int prec, int flags, int radix)
- {
- uint32_t space[(DBL_MANT_DIG + 28) / 29 + 1 +
- (DBL_MAX_EXP + DBL_MANT_DIG + 28 + 8) / 9];
- uint32_t *ap, *dp, *rp, *zp;
- int e2, expo = 0, ix, jx, lx, px = 1, rv = 0;
- char buf[(DBL_MANT_DIG >> 2) + 9];
- char ebuf_st[sizeof (int) * 3];
- char *ebuf = ebuf_st + sizeof (ebuf_st), *estr;
- const char *prefix = "-0X+0X 0X-0x+0x 0x";
- if (dbl < 0)
- dbl = -dbl;
- else if (flags & io_info::FLG_SIGN)
- prefix += 3;
- else if (flags & io_info::FLG_SPACE)
- prefix += 6;
- else
- px = 0, ++prefix;
- dbl = frexpl (dbl, &e2);
- if (abs (radix) == 16)
- { // Hex-float.
- int convdig = radix < 0 ? 'A' : 'a';
- if ((dbl += dbl) != 0)
- --e2;
- if (prec != 0 && prec < (int)((DBL_DIG + 1) * .831) + 1)
- { // Round the mantissa.
- double tail = dbl;
- int q;
- for (q = prec ; ; --q)
- {
- int digit = (int)tail;
- tail -= digit;
- if (q == 0)
- {
- if (digit % 2 != 0 ? tail >= .5 : tail > .5)
- tail = 1 - tail;
- else
- tail = -tail;
- break;
- }
- tail *= 16.l;
- }
- if (tail != 0)
- for (q = prec; q > 0; --q)
- tail *= 0.0625;
- dbl += tail;
- }
- ix = 0; // Return value.
- if ((flags & io_info::FLG_SIGN) || px != 0)
- { ix += KP_TRY (strm->putb (interp, *prefix)); }
- int digit = (int)dbl;
- dbl -= digit;
- {
- char tmp[] = { '0', (char)(convdig - 'A' + 'X'), (char)('0' + digit) };
- ix += KP_TRY (strm->write (interp, tmp, sizeof (tmp)));
- }
- if ((flags & io_info::FLG_ALT) || dbl > 0 || prec > 0)
- {
- ix += KP_TRY (strm->putb (interp, '.')); // XXX: Localization.
- /* The following loop must terminate, based on the above
- * assertion that FLT_RADIX is a power of 2. */
- for (; dbl > 0; --prec)
- {
- dbl *= 16.;
- digit = (int)dbl;
- dbl -= digit;
- ix += KP_TRY (strm->putb (interp, digit < 0 ?
- '0' : convdig - 10 + digit));
- }
- while (prec-- > 0)
- { ix += KP_TRY (strm->putb (interp, '0')); }
- }
- ix += KP_TRY (strm->putb (interp, convdig - 'A' + 'P'));
- ix += KP_TRY (strm->putb (interp, expo < 0 ? '-' : '+'));
- estr = uitostr1 (ebuf, e2, 10);
- ix += KP_TRY (strm->write (interp, estr, ebuf - estr));
- return (ix);
- }
- if ((dbl += dbl) != 0)
- {
- dbl *= 0x1p28;
- e2 -= 29;
- }
- if (prec < 0)
- prec = 6;
- {
- uint32_t *p = e2 < 0 ? space : space + KP_NELEM (space) - DBL_MANT_DIG - 1;
- ap = rp = zp = p;
- }
- do
- {
- *zp = (uint32_t)dbl;
- dbl = 1000000000 * (dbl - *zp++);
- }
- while (dbl != 0);
- while (e2 > 0)
- {
- uint32_t cy = 0;
- int sh = min (29, e2);
- for (dp = zp - 1; dp >= ap; --dp)
- {
- uint64_t twx = ((uint64_t)*dp << sh) + cy;
- *dp = (uint32_t)(twx % 1000000000);
- cy = (uint32_t)(twx / 1000000000);
- }
- if (cy != 0)
- *--ap = cy;
- for (; zp > ap && zp[-1] == 0; --zp) ;
- e2 -= sh;
- }
- while (e2 < 0)
- {
- uint32_t cy = 0;
- int sh = min (9, -e2);
- int need = 1 + (prec + DBL_MANT_DIG / 3 + 8) / 9;
- for (dp = ap; dp < zp; ++dp)
- {
- uint32_t mask = *dp & ((1u << sh) - 1);
- *dp = (*dp >> sh) + cy;
- cy = (1000000000 >> sh) * mask;
- }
- if (*ap == 0)
- ++ap;
- if (cy != 0)
- *zp++ = cy;
- uint32_t *bp = (flags & io_info::FLG_FIXED) ? rp : ap;
- if ((int)(zp - bp) > need)
- zp = bp + need;
- e2 += sh;
- }
- if (ap < zp)
- for (ix = 10, expo = 9 * (int)(rp - ap);
- *ap >= (uint32_t)ix; ix *= 10, ++expo) ;
- jx = prec - expo * !(flags & io_info::FLG_FIXED) -
- (prec != 0 && !(flags & (io_info::FLG_FIXED | io_info::FLG_EXPO)));
- if (jx < 9 * (int)(zp - rp - 1))
- {
- dp = rp + 1 + ((jx + 9 * DBL_MAX_EXP) / 9 - DBL_MAX_EXP);
- jx = (jx + 9 * DBL_MAX_EXP) % 9;
- for (ix = 10, ++jx; jx < 9; ix *= 10, ++jx) ;
- uint32_t tmp = *dp % ix;
- if (tmp != 0 || dp + 1 != zp)
- {
- double round = 2 / DBL_EPSILON;
- if (*dp % 2 != 0 || (ix == 1000000000 &&
- dp > ap && dp[-1] % 2 != 0))
- round += 2;
- double small = tmp < (uint32_t)(ix >> 1) ? 0x0.8p0 :
- tmp == (uint32_t)(ix >> 1) && dp + 1 == zp ? 0x1.0p0 : 0x1.8p0;
- if (px != 0 && *prefix == '-')
- {
- round = -round;
- small = -small;
- }
- *dp -= tmp;
- if (round + small != round)
- {
- *dp += ix;
- while (*dp > 1000000000 - 1)
- {
- *dp-- = 0;
- if (dp < ap)
- *--ap = 0;
- ++*dp;
- }
- for (ix = 10, expo = 9 * (int)(rp - ap);
- *ap >= (uint32_t)ix; ix *= 10, ++expo) ;
- }
- }
- if (zp > dp + 1)
- zp = dp + 1;
- }
- for (; zp > ap && zp[-1] == 0; --zp) ;
- if (!(flags & (io_info::FLG_EXPO | io_info::FLG_FIXED)))
- {
- if (prec == 0)
- prec = 1;
- if (prec > expo && expo >= -4)
- {
- flags |= io_info::FLG_FIXED;
- prec -= expo + 1;
- }
- else
- {
- flags |= io_info::FLG_EXPO;
- --prec;
- }
- if (!(flags & io_info::FLG_ALT))
- {
- if (zp > ap && zp[-1] != 0)
- for (ix = 10, jx = 0; zp[-1] % ix == 0; ix *= 10, ++jx) ;
- else
- jx = 9;
- int uoff = 9 * ((int)(zp - rp) - 1) - jx;
- if (flags & io_info::FLG_EXPO)
- uoff += expo;
- prec = min (prec, max (0, uoff));
- }
- }
- lx = prec || (flags & io_info::FLG_ALT);
- if (prec > INT_MAX - 1 - lx)
- return (-1);
- lx += prec + 1;
- if (flags & io_info::FLG_FIXED)
- {
- if (expo > INT_MAX - lx)
- return (-1);
- else if (expo > 0)
- lx += expo;
- }
- else
- {
- estr = uitostr1 (ebuf, abs (expo), 10);
- while (ebuf - estr < 2)
- *--estr = '0';
- *--estr = expo < 0 ? '-' : '+';
- *--estr = radix < 0 ? 'E' : 'e';
- lx += ebuf - estr;
- }
- if (lx > INT_MAX - px)
- return (-1);
- rv += KP_TRY (pad_writef (interp, strm, ' ', width, px + lx, flags));
- rv += KP_TRY (strm->write (interp, prefix, px));
- rv += KP_TRY (pad_writef (interp, strm, '0', width,
- px + lx, flags ^ io_info::FLG_ZERO));
- if (flags & io_info::FLG_FIXED)
- {
- if (ap > rp)
- ap = rp;
- for (dp = ap; dp <= rp; ++dp)
- {
- char *s = uitostr1 (buf + 9, *dp, 10);
- if (dp != ap)
- while (s > buf)
- *--s = '0';
- else if (s == buf + 9)
- *--s = '0';
- rv += KP_TRY (strm->write (interp, s, buf + 9 - s));
- }
- if (prec || (flags & io_info::FLG_ALT))
- { rv += KP_TRY (strm->putb (interp, '.')); } // XXX: Localization.
- for (; dp < zp && prec > 0; ++dp, prec -= 9)
- {
- char *s = uitostr1 (buf + 9, *dp, 10);
- while (s > buf)
- *--s = '0';
- rv += KP_TRY (strm->write (interp, s, min (9, prec)));
- }
- rv += KP_TRY (pad_writef (interp, strm, '0', prec + 9, 9, 0));
- }
- else
- {
- if (zp <= ap)
- zp = ap + 1;
- for (dp = ap; dp < zp && prec >= 0; ++dp)
- {
- char *s = uitostr1 (buf + 9, *dp, 10);
- if (s == buf + 9)
- *--s = '0';
- if (dp != ap)
- while (s > buf)
- *--s = '0';
- else
- {
- rv += KP_TRY (strm->putb (interp, *s++));
- if (prec > 0 || (flags & io_info::FLG_ALT))
- { rv += KP_TRY (strm->putb (interp, '.')); }
- }
- rv += KP_TRY (strm->write (interp, s,
- min ((int)(buf + 9 - s), prec)));
- prec -= (int)(buf + 9 - s);
- }
- rv += KP_TRY (pad_writef (interp, strm, '0', prec + 18, 18, 0));
- rv += KP_TRY (strm->write (interp, estr, ebuf - estr));
- }
- rv += KP_TRY (pad_writef (interp, strm, ' ', width,
- px + lx, flags ^ io_info::FLG_LJUST));
- return (rv);
- }
- result<int64_t> write_f (interpreter *interp, stream *strm,
- object obj, io_info& info)
- {
- int sign, cls;
- double v = get_dbl (obj, sign, cls);
- if (cls == FP_INFINITE)
- {
- char buf[] = { sign ? '-' : '+', 'I', 'N', 'F' };
- return (strm->write (interp, buf, 4));
- }
- else if (cls != FP_NORMAL)
- return (strm->write (interp, "NaN", 3));
- return (write_flt (interp, strm, v, info.width, info.prec,
- info.flags, info.radix));
- }
- result<int64_t> pack_f (interpreter *interp, stream *strm,
- object obj, pack_info&)
- {
- int sign, cls;
- double v = get_dbl (obj, sign, cls);
- return (strm->write (interp, &v));
- }
- result<object> unpack_f (interpreter *interp, stream *strm,
- pack_info& info, bool)
- {
- double v;
- bool rv = KP_TRY (strm->sread (interp, &v));
- if (rv)
- kp_return (fltobj::make (interp, v));
- return (info.error ("failed to read floating point value"));
- }
- // Big float implementation.
- static_assert ((alignof (bigfloat) % alignof (limb_t)) == 0,
- "invalid alignment for big floats");
- bigfloat* bigfloat::alloc_raw (int len)
- {
- bigfloat *retp = (bigfloat *)alloch (sizeof (*retp) + sizeof (limb_t) * len,
- typecode::BIGFLOAT, TYPE_SHIFT + 1);
- retp->len = len;
- retp->expo = 0;
- retp->data = (limb_t *)&retp[1];
- retp->vo_full = FLAGS_CONST;
- return (retp);
- }
- object alloc_bigfloat (interpreter *interp, int len)
- {
- auto retp = bigfloat::alloc_raw (len);
- interp->alval = retp->as_obj ();
- gc_register (interp, retp, sizeof (*retp) + len * sizeof (limb_t));
- return (interp->alval);
- }
- static bigfloat*
- get_bigfloat (object obj, int& sign)
- {
- #ifdef KP_ARCH_WIDE
- sign = (obj & SIGN_BIT) / SIGN_BIT;
- #else
- sign = as_bigfloat(obj)->len < 0;
- #endif
- return (as_bigfloat (obj));
- }
- static_assert (sizeof (bigfloat) >= sizeof (fltobj),
- "big floats must be larger or equal in size to floats");
- static object
- ret_F (interpreter *interp, bigfloat *src, int sign)
- {
- if (src->len == 0)
- {
- xfree (src);
- return (FLT_ZERO);
- }
- int b2exp = (src->expo - src->len) * LIMB_BITS;
- if (b2exp >= DBL_MIN_EXP && b2exp <= DBL_MAX_EXP)
- { // May fit in a regular double.
- double dbl = uitodbl (src->data, src->len, src->expo);
- if (finf_p (dbl))
- goto do_bigfloat;
- /* Fits in a regular double - Transmute the type.
- * This wastes some memory, but is preferable to
- * allocating a new object. */
- fltobj tmp;
- tmp.gc_link = src->gc_link;
- tmp.vo_full = src->vo_full;
- tmp.vo_type = typecode::FLOAT;
- tmp.val = sign ? -dbl : dbl;
- memcpy (src, &tmp, sizeof (tmp));
- interp->retval = src->as_obj ();
- }
- else
- do_bigfloat:
- interp->retval = make_bigfloat (src, sign);
- gc_register (interp, src, sizeof (*src) + src->len * sizeof (limb_t));
- return (interp->retval);
- }
- static bigfloat*
- bfloat_add (interpreter *interp, const bigfloat *x, const bigfloat *y)
- {
- if (x->expo < y->expo)
- swap (x, y);
- int ediff = x->expo - y->expo, xl = F_ABS (x->len), yl = F_ABS (y->len);
- int prec = max (xl, yl);
- if (ediff >= prec)
- // Y completely cancelled.
- return ((bigfloat *)x);
- auto ret = bigfloat::alloc_raw (prec + 1);
- ret->expo = x->expo;
- int size;
- limb_t cy;
- if (!(xl > ediff))
- {
- size = yl + ediff - xl;
- memcpy (ret->data, y->data, yl * sizeof (limb_t));
- memset (ret->data + yl, 0, (ediff - xl) * sizeof (limb_t));
- memcpy (ret->data + size, x->data, xl * sizeof (limb_t));
- cy = 0, ret->len = size + xl;
- }
- else if (yl + ediff <= xl)
- {
- size = xl - ediff - yl;
- memcpy (ret->data, x->data, size * sizeof (limb_t));
- cy = uiadd (ret->data + size, x->data + size, xl - size, y->data, yl);
- ret->len = yl;
- }
- else
- {
- size = yl + ediff - xl;
- memcpy (ret->data, y->data, size * sizeof (limb_t));
- cy = uiadd (ret->data + size, x->data, xl, y->data + size, xl - ediff);
- ret->len = yl + ediff;
- }
- ret->data[ret->len] = cy;
- ret->len += cy, ret->expo += cy;
- return (ret);
- }
- static bigfloat*
- bfloat_sub (interpreter *interp, const bigfloat *x,
- const bigfloat *y, int *negp)
- {
- int neg = 0;
- if (x->expo < y->expo)
- {
- swap (x, y);
- neg = 1;
- }
- int expo = x->expo, ediff = x->expo - y->expo;
- int xl = F_ABS (x->len), yl = F_ABS (y->len), len = max (xl, yl);
- limb_t *xp = x->data, *yp = y->data;
- auto ret = bigfloat::alloc_raw (len);
- limb_t cy;
- if (ediff <= 1)
- { /* The numbers are very close. Scan them and
- * ignore the limbs that compare equal. */
- if (ediff == 0)
- {
- if (xp[xl - 1] == yp[yl - 1])
- do
- {
- --yl, --expo;
- if (--xl == 0)
- {
- neg ^= 1;
- cancel:
- for (; yl && !yp[yl - 1]; --yl, --expo) ;
- memcpy (ret->data, y->data, yl * sizeof (limb_t));
- ret->len = yl;
- goto done;
- }
- if (yl == 0)
- {
- yp = xp;
- yl = xl;
- goto cancel;
- }
- }
- while (xp[xl - 1] == yp[yl - 1]);
- if (xp[xl - 1] < yp[yl - 1])
- {
- swap (xp, yp);
- swap (xl, yl);
- neg ^= 1;
- }
- if (xp[xl - 1] != yp[yl - 1] + 1)
- goto main;
-
- --xl, --yl, --expo;
- }
- else // (ediff == 1)
- {
- if (xp[xl - 1] != 1 || yp[yl - 1] != ~(limb_t)0 ||
- (xl >= 2 && xp[xl - 2] != 0))
- goto main;
- --xl, --expo;
- }
- for (; yl && xl && !xp[xl - 1] &&
- yp[yl - 1] == ~(limb_t)0; --xl, --yl, --expo) ;
- if (xl == 0)
- for (; yl && yp[yl - 1] == ~(limb_t)0; --yl, --expo) ;
-
- if (yl == 0)
- {
- memcpy (ret->data, xp, (ret->len = xl) * sizeof (limb_t));
- cy = 0;
- }
- else if (xl == 0)
- {
- for (int i = 0; i < yl; ++i)
- ret->data[i] = ~yp[i];
- cy = 1 - uiadd1 (ret->data, ret->data, yl, 1);
- ret->len = yl;
- }
- else if (xl >= yl)
- {
- int size = xl - yl;
- memcpy (ret->data, xp, size * sizeof (limb_t));
- cy = uisubn (ret->data + size, xp + size, yp, yl);
- ret->len = xl;
- }
- else
- {
- int size = yl - xl;
- for (int i = 0; i < size; ++i)
- ret->data[i] = ~yp[i];
- cy = uisubn (ret->data + size, xp, yp + size, xl);
- cy += uisub1 (ret->data + size, ret->data + size, xl, 1);
- cy -= uiadd1 (ret->data, ret->data, yl, 1);
- ret->len = yl;
- }
- if (cy == 0)
- {
- ret->data[ret->len++] = 1;
- ++expo;
- goto normalize;
- }
- }
-
- main:
- if (ediff >= len)
- {
- memcpy (ret->data, xp, xl);
- ret->len = xl;
- expo = ((const bigfloat *)
- ((char *)xp - sizeof (bigfloat)))->expo;
- goto done;
- }
- // Search for the least significant non-zero limbs in X and Y.
-
- while (true)
- {
- if (yl == 0)
- {
- memcpy (ret->data, xp, xl * sizeof (limb_t));
- ret->len = xl;
- goto done;
- }
-
- if (*yp != 0)
- break;
-
- ++yp, --yl;
- }
- while (true)
- {
- if (xl == 0)
- {
- memcpy (ret->data, yp, yl * sizeof (limb_t));
- ret->len = yl;
- neg ^= 1;
- goto done;
- }
-
- if (*xp != 0)
- break;
-
- ++xp, --xl;
- }
- if (xl > ediff)
- { // Partial overlapping.
- if (!ediff)
- { /* Compare leading limbs to determine whether to
- * do X - Y or Y - X. */
- if (xl >= yl)
- {
- int size = xl - yl;
- memcpy (ret->data, xp, size * sizeof (limb_t));
- uisubn (ret->data + size, xp + size, yp, yl);
- ret->len = xl;
- }
- else
- {
- int size = yl - xl;
- *ret->data = -*yp;
- for (int i = 1; i < size; ++i)
- ret->data[i] = ~yp[i];
-
- uisubn (ret->data + size, xp, yp + size, xl);
- uisub1 (ret->data + size, ret->data + size, xl, 1);
- ret->len = yl;
- }
- }
- else if (yl + ediff <= xl)
- {
- int size = xl - ediff - yl;
- memcpy (ret->data, xp, size * sizeof (limb_t));
- uisub (ret->data + size, xp + size, xl - size, yp, yl);
- ret->len = xl;
- }
- else
- {
- int size = yl + ediff - xl;
- *ret->data = -*yp;
- for (int i = 1; i < size; ++i)
- ret->data[i] = ~yp[i];
- uisub (ret->data + size, xp, xl, yp + size, xl - ediff);
- uisub1 (ret->data + size, ret->data + size, xl, 1);
- ret->len = yl + ediff;
- }
- }
- else
- { // No overlapping - Trivial copy.
- int size = yl + ediff - xl;
- *ret->data = -*yp;
- for (int i = 1; i < yl; ++i)
- ret->data[i] = ~yp[i];
- memset (&ret->data[yl], 0xff,
- (size - yl) * sizeof (limb_t));
- uisub1 (ret->data + size, xp, xl, 1);
- ret->len = size + xl;
- }
- normalize:
- for (; ret->len && !ret->data[ret->len - 1]; --ret->len, --expo) ;
- done:
- ret->expo = !ret->len ? 0 : expo;
- if (negp != nullptr) *negp = neg;
- return (ret);
- }
- static bigfloat*
- bfloat_mul (interpreter *interp, const bigfloat *x, const bigfloat *y)
- {
- int xl = F_ABS (x->len), yl = F_ABS (y->len), tsize = xl + yl;
- auto ret = bigfloat::alloc_raw (tsize);
- limb_t cy = x == y ? uisqr (interp, ret->data, x->data, xl) :
- uimul (interp, ret->data, x->data, xl, y->data, yl);
- cy = !cy;
- ret->len -= cy;
- ret->expo = x->expo + y->expo - cy;
- return (ret);
- }
- static bigfloat*
- bfloat_div (interpreter *interp, const bigfloat *x, const bigfloat *y)
- {
- int xl = F_ABS (x->len), yl = F_ABS (y->len), rsize = xl - yl;
- limb_t *tmp, *xp = x->data;
- tmp_allocator ta { interp };
- if (rsize <= 0)
- { // Copy and extend X.
- uint32_t padding = rsize == 0 ? 1 : (rsize = -rsize) * 2;
- xp = (limb_t *)ta.alloc ((xl * 2 + padding) * sizeof (*xp));
- memset (xp, 0, padding * sizeof (limb_t));
- memcpy (xp + padding, x->data, xl * sizeof (limb_t));
- xl += padding, tmp = xp + xl;
- }
- else
- tmp = (limb_t *)ta.alloc ((xl + 1) * sizeof (*tmp));
- rsize = uidiv (interp, xp, xl, y->data, yl, tmp, false);
- auto ret = bigfloat::alloc_raw (rsize);
- memcpy (ret->data, tmp + yl, ret->len * sizeof (limb_t));
- ret->expo = x->expo - y->expo + 1 -
- (xl - yl + 1 != ret->len); // i.e: msd equals zero.
- return (ret);
- }
- object add_FF (interpreter *interp, object x, object y)
- {
- int s1, s2;
- auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
- bigfloat *ret;
- if (!(s1 ^ s2))
- {
- ret = bfloat_add (interp, v1, v2);
- if (ret == v1)
- kp_return (x);
- else if (ret == v2)
- kp_return (y);
- }
- else
- ret = s1 ? bfloat_sub (interp, v2, v1, &s1) :
- bfloat_sub (interp, v1, v2, &s1);
- return (ret_F (interp, ret, s1));
- }
- object add_fF (interpreter *interp, object x, object y)
- {
- int s1, c1, s2;
- double v1 = get_dbl (x, s1, c1);
- const auto v2 = get_bigfloat (y, s2);
- if (c1 != FP_NORMAL)
- kp_return (x); // NaN or +/- Inf.
- fake_bigfloat tmp (v1);
- bigfloat *ret;
- if (!(s1 ^ s2))
- {
- ret = bfloat_add (interp, &tmp, v2);
- if (ret == (bigfloat *)&tmp)
- kp_return (x);
- else if (ret == v2)
- kp_return (y);
- }
- else
- ret = s1 ? bfloat_sub (interp, v2, &tmp, &s1) :
- bfloat_sub (interp, &tmp, v2, &s1);
- return (ret_F (interp, ret, s1));
- }
- object sub_FF (interpreter *interp, object x, object y)
- {
- #ifdef KP_ARCH_WIDE
- return (add_FF (interp, x, y ^ SIGN_BIT));
- #else
- local_varobj<bigfloat> tmp;
- bigfloat *yp = as_bigfloat (y);
- tmp.local_init (yp->data, -yp->len, yp->expo);
- return (add_FF (interp, x, tmp.as_obj ()));
- #endif
- }
- object sub_fF (interpreter *interp, object x, object y)
- {
- #ifdef KP_ARCH_WIDE
- return (add_fF (interp, x, y ^ SIGN_BIT));
- #else
- local_varobj<bigfloat> tmp;
- bigfloat *yp = as_bigfloat (y);
- tmp.local_init (yp->data, -yp->len, yp->expo);
- return (add_fF (interp, x, tmp.as_obj ()));
- #endif
- }
- object neg_F (interpreter *interp, object obj)
- {
- #ifdef KP_ARCH_WIDE
- kp_return (obj ^ SIGN_BIT);
- #else
- const bigfloat *srcp = as_bigfloat (obj);
- int sl = F_ABS (srcp->len);
- bigfloat *retp = as_bigfloat (alloc_bigfloat (interp, sl));
- memcpy (retp->data, srcp->data, sl * sizeof (limb_t));
- retp->expo = srcp->expo;
- retp->len = -sl;
- kp_return (retp->as_obj ());
- #endif
- }
- object mul_FF (interpreter *interp, object x, object y)
- {
- int s1, s2;
- const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
- return (ret_F (interp, bfloat_mul (interp, v1, v2), s1 ^ s2));
- }
- object mul_fF (interpreter *interp, object x, object y)
- {
- int s1, c1, s2;
- double v1 = get_dbl (x, s1, c1);
- if (c1 != FP_NORMAL)
- kp_return (x); // NaN or +/- Inf.
- fake_bigfloat tmp (v1);
- const auto v2 = get_bigfloat (y, s2);
- return (ret_F (interp, bfloat_mul (interp, &tmp, v2), s1 ^ s2));
- }
- object div_FF (interpreter *interp, object x, object y)
- {
- int s1, s2;
- const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
- return (ret_F (interp, bfloat_div (interp, v1, v2), s1 ^ s2));
- }
- #define DIV_BODY(left, right, op, alt, zero_x) \
- int s1, s2, cl; \
- double v1 = get_dbl (left, s1, cl); \
- const auto v2 = get_bigfloat (right, s2); \
- \
- if (cl != FP_NORMAL) \
- kp_return (alt); \
- else if (v1 == 0) \
- kp_return (zero_x); \
- \
- fake_bigfloat tmp (v1); \
- return (ret_F (interp, op, s1 ^ s2))
- object div_fF (interpreter *interp, object x, object y)
- {
- DIV_BODY (x, y, bfloat_div (interp, &tmp, v2), x, FLT_ZERO);
- }
- object div_Ff (interpreter *interp, object x, object y)
- {
- DIV_BODY (y, x, bfloat_div (interp, v2, &tmp),
- cl == FP_INFINITE ? FLT_ZERO : x, FLT_PINF);
- }
- #undef DIV_BODY
- static bigfloat*
- bfloat_mod (interpreter *interp, const bigfloat *x, const bigfloat *y)
- {
- bigfloat *p, *q, *r;
- p = bfloat_div (interp, x, y); // p := x / y
- // Truncate P in place.
- if (p->expo <= 0)
- p->expo = p->len = 0;
- else
- {
- int len = min (p->len, p->expo);
- limb_t *lp = p->data + p->len - len;
- if (lp != p->data)
- memmove (p->data, lp, (p->len - len) * sizeof (*lp));
-
- p->len = len;
- }
- q = bfloat_mul (interp, p, y); // q := trunc (p) * y
- r = bfloat_sub (interp, x, q, nullptr); // r := x - q
- xfree (p);
- xfree (q);
- return (r);
- }
- object mod_FF (interpreter *interp, object x, object y)
- {
- int s1, s2;
- const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
- return (ret_F (interp, bfloat_mod (interp, v1, v2), s1));
- }
- object mod_fF (interpreter *interp, object x, object y)
- {
- int s1, c1, s2;
- double v1 = get_dbl (x, s1, c1);
- if (c1 != FP_NORMAL)
- kp_return (FLT_QNAN);
- fake_bigfloat tmp (v1);
- const auto v2 = get_bigfloat (y, s2);
- return (ret_F (interp, bfloat_mod (interp, &tmp, v2), s1));
- }
- object mod_Ff (interpreter *interp, object x, object y)
- {
- int s1, c2, s2;
- double v2 = get_dbl (y, s2, c2);
- if (c2 != FP_NORMAL)
- kp_return (FLT_QNAN);
- const auto v1 = get_bigfloat (x, s1);
- fake_bigfloat tmp (v2);
- return (ret_F (interp, bfloat_mod (interp, v1, &tmp), s1));
- }
- static int
- bfloat_cmp (const bigfloat *x, int s1, const bigfloat *y, int s2)
- {
- int ret = 0;
- if (s1 != s2)
- ret = 1;
- else if (!(ret = x->expo - y->expo))
- {
- limb_t *xp = x->data, *yp = y->data;
- int xl, yl;
- for (; *xp == 0; ++xp) ;
- for (; *yp == 0; ++yp) ;
- xl = (int)((x->data + x->len) - xp);
- yl = (int)((y->data + y->len) - yp);
-
- if (xl > yl)
- ret = uicmpn (xp + xl - yl, yp, yl);
- else if (yl > xl)
- ret = -uicmpn (xp, yp + yl - xl, xl);
- else
- ret = uicmpn (xp, yp, yl);
- }
- return (s1 ? -ret : ret);
- }
- int cmp_FF (interpreter *interp, object x, object y)
- {
- int s1, s2;
- const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
- return (bfloat_cmp (v1, s1, v2, s2));
- }
- int cmp_fF (interpreter *interp, object x, object y)
- {
- int s1, c1, s2;
- double v1 = get_dbl (x, s1, c1);
- if (c1 != FP_NORMAL)
- return (s1); // i.e: +/-1 for Inf, +1 for NaN.
- const bigfloat *v2 = get_bigfloat (y, s2);
- fake_bigfloat tmp (v1);
- return (bfloat_cmp (&tmp, s1, v2, s2));
- }
- bool eq_FF (interpreter *interp, object x, object y)
- {
- int s1, s2;
- const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
- // Exponents and signs must match.
- if (s1 != s2 || v1->expo != v2->expo)
- return (false);
- else if (v1->data == v2->data)
- return (true);
- int minlen, maxlen, i, len, diff;
- int xl = F_ABS (v1->len), yl = F_ABS (v2->len);
- const limb_t *xp = v1->data + xl, *yp = v2->data + yl;
- // Test for most significant limb position.
- if ((*(yp - 1) >> (LIMB_BITS - uiclz (*(xp - 1) - 1))) != 1)
- return (0);
- if (xl < yl)
- minlen = xl, maxlen = yl;
- else
- minlen = yl, maxlen = xl;
- // Compare the most significant limbs.
- xp += minlen, yp += minlen;
- for (i = minlen - 1; i > 0; --i)
- if (xp[i] != yp[i])
- return (0);
- if (0 != (len = maxlen - minlen))
- { /* Either X or Y has been completely examined.
- * Test that the other is made of nil limbs. */
- if (*xp != *yp)
- return (0);
-
- const limb_t *p = (xl > yl ? xp : yp) - len;
- for (i = len -1 ; i > 0; --i)
- if (p[i] != 0)
- return (0);
-
- diff = *p;
- }
- else
- diff = *xp ^ *yp;
- return (diff == 0);
- }
- static object
- bfloat_ceilfloor (interpreter *interp, object obj, int dir)
- {
- int osign;
- const auto x = get_bigfloat (obj, osign);
- int sign = osign ? -1 : 1;
- if (x->expo <= 0)
- kp_return ((sign ^ dir) < 0 ? FLT_ZERO : fltobj::make (interp, dir));
- int len = min (F_ABS (x->len), x->expo);
- auto ret = bigfloat::alloc_raw (len);
- limb_t *xp = x->data + x->len - len;
-
- ret->expo = x->expo;
- if ((sign ^ dir) >= 0)
- {
- const limb_t *tp;
- for (tp = x->data; tp != xp; ++tp)
- if (*tp != 0)
- {
- if (uiadd1 (ret->data, xp, len, 1))
- {
- ret->len = *ret->data = 1;
- ++ret->expo;
- }
-
- return (ret_F (interp, ret, osign));
- }
- }
-
- memcpy (ret->data, xp, len * sizeof (limb_t));
- return (ret_F (interp, ret, osign));
- }
- object floor_F (interpreter *interp, object obj)
- {
- return (bfloat_ceilfloor (interp, obj, -1));
- }
- object ceil_F (interpreter *interp, object obj)
- {
- return (bfloat_ceilfloor (interp, obj, +1));
- }
- object trunc_F (interpreter *interp, object obj)
- {
- int sign;
- const auto x = get_bigfloat (obj, sign);
- if (x->expo <= 0)
- kp_return (FLT_ZERO);
- int xl = F_ABS (x->len), len = min (xl, x->expo);
- auto ret = bigfloat::alloc_raw (len);
- memcpy (ret->data, x->data + xl - len, len * sizeof (limb_t));
- ret->expo = x->expo;
- return (ret_F (interp, ret, sign));
- }
- uint32_t hash_F (interpreter *interp, object obj)
- {
- int sign;
- const auto lp = get_bigfloat (obj, sign);
- int fl = F_ABS (lp->len);
- if (fl == lp->expo)
- {
- local_varobj<bigint> bi;
- bi.local_init (lp->data, lp->len);
- #ifdef KP_ARCH_WIDE
- return (hash_I (interp, bi.as_obj () | (obj & SIGN_BIT)));
- #else
- return (hash_I (interp, bi.as_obj ()));
- #endif
- }
- uint32_t ret = mix_hash (2 * fl, lp->expo);
- for (int i = 0; i < fl; ++i)
- ret = mix_hash (ret, lp->data[i]);
- return (sign ? ~ret :ret);
- }
- result<int64_t> write_F (interpreter *interp, stream *strm,
- object obj, io_info& info)
- {
- const auto lp = as_bigfloat (obj);
- int expo = lp->expo, xl = F_ABS (lp->len), ret = 0;
- int size = uibsize (info.radix, lp->data[xl - 1], xl);
- tmp_allocator ta { interp };
- char *ptr = (char *)ta.alloc (size);
- size = lftostr (interp, ptr, &expo, lp->data, xl, size, info.radix);
- #ifdef KP_ARCH_WIDE
- if (obj & SIGN_BIT)
- #else
- if (lp->len < 0)
- #endif
- { ret += KP_TRY (strm->putb (interp, '-')); }
- // Write the decimal separator, followed by the mantissa.
- ret += KP_TRY (strm->write (interp, "0.", 2));
- ret += KP_TRY (strm->write (interp, ptr, size));
- // Write the exponent.
- size = sprintf (ptr, "%c%d", info.radix <= 10 ? 'e' : '@', expo);
- ret += KP_TRY (strm->write (interp, ptr, size));
- return (ret);
- }
- result<int64_t> pack_F (interpreter *interp, stream *strm,
- object obj, pack_info&)
- {
- const bigfloat *lp = as_bigfloat (obj);
- int expo = lp->expo, xl = lp->len;
- int64_t ret = 0;
- #ifdef KP_ARCH_WIDE
- if (obj & SIGN_BIT)
- xl = -xl;
- #endif
- ret += KP_TRY (strm->write (interp, &expo));
- ret += KP_TRY (strm->write (interp, &xl));
- ret += KP_TRY (strm->write (interp, lp->data,
- F_ABS (lp->len) * sizeof (*lp->data)));
- return (ret);
- }
- result<object> unpack_F (interpreter *interp, stream *strm,
- pack_info& info, bool save)
- {
- int vals[2];
- {
- auto tmp = KP_TRY (strm->read (interp, vals, sizeof (vals)));
- if ((size_t)tmp != sizeof (vals))
- return (info.error ("failed to read floating point value"));
- }
- int len = abs (vals[1]);
- bigfloat *lp = as_bigfloat (alloc_bigfloat (interp, len));
- auto tmp = KP_TRY (strm->read (interp, lp->data, len * sizeof (*lp->data)));
- if ((size_t)tmp != (size_t)(len * sizeof (*lp->data)))
- return (info.error ("invalid number of limbs read"));
- else if (save)
- KP_VTRY (info.add_mapping (interp, *info.offset, lp->as_obj ()));
- lp->expo = vals[0];
- kp_return (make_bigfloat (lp, vals[1] < 0));
- }
- // Mixed operations.
- static inline object
- int2flt (object i, local_varobj<fltobj>& out)
- {
- out.vo_full = 0;
- out.vo_type = typecode::FLOAT;
- out.val = as_int (i);
- return (out.as_obj ());
- }
- static inline object
- bigint2flt (object i, local_varobj<bigfloat>& out)
- {
- const auto lp = as_bigint (i);
- out.local_init (lp->data, lp->len, F_ABS (lp->len));
- #ifdef KP_ARCH_WIDE
- return (out.as_obj () | (i & SIGN_BIT));
- #else
- return (out.as_obj ());
- #endif
- }
- static object
- downgrade_float (interpreter *interp, object obj)
- {
- if (float_p (obj))
- {
- double i, val = as_float (obj);
- if (modf (val, &i) == 0)
- {
- int64_t uv = (int64_t)val;
- return (intobj (interp, uv));
- }
- }
- else
- {
- const auto fp = as_bigfloat (obj);
- int fl = F_ABS (fp->len);
- if (fl == fp->expo)
- {
- auto lp = as_bigint (alloc_bigint (interp, fl));
- memcpy (lp->data, fp->data, fl * sizeof (*fp->data));
- lp->len = fp->len;
- #ifdef KP_ARCH_WIDE
- interp->alval |= obj & SIGN_BIT;
- #endif
- return (interp->alval);
- }
- }
- return (obj);
- }
- #define MIXED_FN_NOCHECK() ((void)0)
- #define MIXED_FN_CHECK() \
- if (ret.error_p ()) \
- return (exception ()); \
- else (void)0
- #define MIXED_FN(name, ret_ty, check) \
- ret_ty name##_if (interpreter *interp, object x, object y) \
- { \
- local_varobj<fltobj> tmp; \
- ret_ty ret = name##_ff (interp, int2flt (x, tmp), y); \
- check (); \
- kp_return (downgrade_float (interp, deref (ret))); \
- } \
- \
- object name##_iF (interpreter *interp, object x, object y) \
- { \
- local_varobj<fltobj> tmp; \
- kp_return (downgrade_float (interp, \
- name##_fF (interp, int2flt (x, tmp), y))); \
- }
- MIXED_FN (add, object, MIXED_FN_NOCHECK)
- MIXED_FN (sub, object, MIXED_FN_NOCHECK)
- MIXED_FN (mul, object, MIXED_FN_NOCHECK)
- #define MIXED_FN2(name) \
- MIXED_FN (name, result<object>, MIXED_FN_CHECK) \
- \
- result<object> name##_fi (interpreter *interp, object x, object y) \
- { \
- local_varobj<fltobj> tmp; \
- object ret = KP_TRY (name##_ff (interp, x, int2flt (y, tmp))); \
- kp_return (downgrade_float (interp, ret)); \
- } \
- \
- result<object> name##_Fi (interpreter *interp, object x, object y) \
- { \
- local_varobj<fltobj> tmp; \
- object ret = KP_TRY (name##_Ff (interp, x, int2flt (y, tmp))); \
- kp_return (downgrade_float (interp, ret)); \
- }
- MIXED_FN2 (div)
- MIXED_FN2 (mod)
- #undef MIXED_FN
- #undef MIXED_FN2
- #define MIXED_FN(name, ret_ty, check) \
- ret_ty name##_If (interpreter *interp, object x, object y) \
- { \
- int s2, c2; \
- double v2 = get_dbl (y, s2, c2); \
- if (c2 != FP_NORMAL) \
- kp_return (y); \
- \
- check (v2); \
- local_varobj<bigfloat> tmp; \
- fake_bigfloat ry (v2); \
- \
- kp_return (downgrade_float (interp, \
- name##_FF (interp, bigint2flt (x, tmp), \
- ry.get_obj (s2)))); \
- } \
- \
- object name##_IF (interpreter *interp, object x, object y) \
- { \
- local_varobj<bigfloat> tmp; \
- kp_return (downgrade_float (interp, \
- name##_FF (interp, bigint2flt (x, tmp), y))); \
- } \
- #undef MIXED_FN_NOCHECK
- #define MIXED_FN_NOCHECK(x) ((void)0)
- MIXED_FN (add, object, MIXED_FN_NOCHECK)
- MIXED_FN (sub, object, MIXED_FN_NOCHECK)
- MIXED_FN (mul, object, MIXED_FN_NOCHECK)
- #define MIXED_FN2(name) \
- MIXED_FN (name, result<object>, MIXED_FN_CHECK) \
- \
- object name##_fI (interpreter *interp, object x, object y) \
- { \
- int s1, c1; \
- double v1 = get_dbl (x, s1, c1); \
- if (c1 != FP_NORMAL) \
- kp_return (x); \
- \
- local_varobj<bigfloat> tmp; \
- fake_bigfloat rx (v1); \
- \
- kp_return (downgrade_float (interp, \
- name##_FF (interp, rx.get_obj (s1), \
- bigint2flt (y, tmp)))); \
- } \
- \
- object name##_FI (interpreter *interp, object x, object y) \
- { \
- local_varobj<bigfloat> tmp; \
- kp_return (downgrade_float (interp, \
- name##_FF (interp, x, bigint2flt (y, tmp)))); \
- }
- #undef MIXED_FN_CHECK
- #define MIXED_FN_CHECK(x) \
- if ((x) == 0) \
- return (interp->raise ("arith-error", "division by zero")); \
- else (void)0
- MIXED_FN2 (div)
- MIXED_FN2 (mod)
- // Division of integers that may produce floating point values.
- result<object> div_ii (interpreter *interp, object x, object y)
- {
- int v1 = as_int (x), v2 = as_int (y);
- if (v2 == 0)
- return (interp->raise ("arith-error", "division by zero"));
- kp_return (v1 % v2 == 0 ? fixint (v1 / v2) :
- fltobj::make (interp, (double)v1 / v2));
- }
- object div_iI (interpreter *interp, object x, object y)
- {
- local_varobj<fltobj> flt;
- local_varobj<bigfloat> bigf;
- return (div_fF (interp, int2flt (x, flt), bigint2flt (y, bigf)));
- }
- result<object> div_Ii (interpreter *interp, object x, object y)
- {
- if (as_int (y) == 0)
- return (interp->raise ("arith-error", "division by zero"));
- local_varobj<fltobj> flt;
- local_varobj<bigfloat> bigf;
- return (downgrade_float (interp,
- div_Ff (interp, bigint2flt (x, bigf),
- int2flt (y, flt))));
- }
- object div_II (interpreter *interp, object x, object y)
- {
- local_varobj<bigfloat> lf1, lf2;
- return (downgrade_float (interp,
- div_FF (interp, bigint2flt (x, lf1),
- bigint2flt (y, lf2))));
- }
- #undef MIXED_FN
- #undef MIXED_FN2
- int cmp_if (interpreter *interp, object x, object y)
- {
- local_varobj<fltobj> fv;
- return (cmp_ff (interp, int2flt (x, fv), y));
- }
- int cmp_iF (interpreter *interp, object x, object y)
- {
- local_varobj<fltobj> fv;
- return (cmp_fF (interp, int2flt (x, fv), y));
- }
- int cmp_If (interpreter *interp, object x, object y)
- {
- local_varobj<bigfloat> fv;
- return (-cmp_fF (interp, y, bigint2flt (x, fv)));
- }
- int cmp_IF (interpreter *interp, object x, object y)
- {
- local_varobj<bigfloat> fv;
- return (cmp_FF (interp, bigint2flt (x, fv), y));
- }
- bool eq_if (interpreter *, object x, object y)
- {
- return ((double)as_int (x) == as_float (y));
- }
- bool eq_If (interpreter *, object x, object y)
- {
- const auto lp = as_bigint (x);
- int xl = F_ABS (lp->len);
- if (xl > LNDBL_SIZE)
- return (false);
- double dbl = uitodbl (lp->data, xl, xl);
- #ifdef KP_ARCH_WIDE
- if (x & SIGN_BIT)
- #else
- if (lp->len < 0)
- #endif
- dbl = -dbl;
- return (as_float (y) == dbl);
- }
- bool eq_IF (interpreter *, object x, object y)
- {
- const auto lp = as_bigint (x);
- const auto fp = as_bigfloat (y);
- #ifdef KP_ARCH_WIDE
- if ((x & SIGN_BIT) != (y & SIGN_BIT))
- #else
- if ((lp->len ^ fp->len) < 0)
- #endif
- return (false);
- return (F_ABS (fp->len) == fp->expo && fp->len == lp->len &&
- memcmp (fp->data, lp->data, F_ABS (lp->len) * sizeof (limb_t)) == 0);
- }
- object lsh_fi (interpreter *interp, object x, object y)
- {
- int s1, c1, ival = as_int (y);
- double d1 = get_dbl (x, s1, c1);
- if (kp_unlikely (c1 != FP_NORMAL))
- kp_return (x); // +/- Inf or NaN
- double ret = ldexp (d1, ival);
- if (kp_likely (!finf_p (ret)))
- kp_return (fltobj::make (interp, ret));
- fake_bigfloat tx (d1);
- return (lsh_Fi (interp, tx.as_obj (), y));
- }
- object lsh_Fi (interpreter *interp, object x, object y)
- {
- int sign, ival = as_int (y);
- if (ival < 0)
- return (rsh_Fi (interp, x, -ival));
- const auto flt = get_bigfloat (x, sign);
- int len = F_ABS (flt->len);
- bigfloat *ret;
- if (ival % LIMB_BITS == 0)
- {
- ret = bigfloat::alloc_raw (len);
- memcpy (ret->data, flt->data, len * sizeof (limb_t));
- ret->expo = flt->expo + flt->expo / LIMB_BITS;
- ret->len = len;
- }
- else
- {
- ret = bigfloat::alloc_raw (len + 1);
- limb_t cy = uimul2exp (ret->data, flt->data,
- len, ival % LIMB_BITS);
- ret->data[len] = cy;
- ret->len = len + (cy != 0);
- ret->expo = flt->expo + ival / LIMB_BITS + (cy != 0);
- }
- return (ret_F (interp, ret, sign));
- }
- object rsh_fi (interpreter *interp, object x, object y)
- {
- int s1, c1, ival = as_int (y);
- double d1 = get_dbl (x, s1, c1);
- if (kp_unlikely (c1 != FP_NORMAL))
- kp_return (x); // +/- Inf or NaN
- double ret = ldexp (d1, -ival);
- if (kp_likely (!finf_p (ret)))
- kp_return (fltobj::make (interp, ret));
- fake_bigfloat tx (d1);
- return (rsh_Fi (interp, tx.as_obj (), y));
- }
- object rsh_Fi (interpreter *interp, object x, object y)
- {
- int sign, ival = as_int (y);
- if (ival < 0)
- return (rsh_Fi (interp, x, -ival));
- const auto flt = get_bigfloat (x, sign);
- int len = F_ABS (flt->len);
- bigfloat *ret;
- if (ival % LIMB_BITS == 0)
- {
- ret = bigfloat::alloc_raw (len);
- memcpy (ret->data, flt->data, len * sizeof (limb_t));
- ret->expo = flt->expo - flt->expo / LIMB_BITS;
- ret->len = len;
- }
- else
- {
- ret = bigfloat::alloc_raw (len + 1);
- limb_t cy = uidiv2exp (ret->data, flt->data,
- len, LIMB_BITS - ival % LIMB_BITS);
- ret->data[len] = cy;
- ret->len = len + (cy != 0);
- ret->expo = flt->expo - ival / LIMB_BITS + (cy != 0) - 1;
- }
- return (ret_F (interp, ret, sign));
- }
- // External definitions.
- object FLT_PINF;
- object FLT_NINF;
- object FLT_QNAN;
- object FLT_ZERO;
- static int
- do_init_float (interpreter *)
- {
- static fltobj SF_PINF;
- static fltobj SF_QNAN;
- static fltobj SF_ZERO;
- #ifndef KP_ARCH_WIDE
- static fltobj SF_NINF;
- SF_NINF.vo_type = typecode::FLOAT;
- #endif
- SF_PINF.vo_type = SF_QNAN.vo_type = SF_ZERO.vo_type = typecode::FLOAT;
- SF_PINF.val = std::numeric_limits<double>::infinity ();
- SF_QNAN.val = std::numeric_limits<double>::quiet_NaN ();
- SF_ZERO.val = 0.;
- // Now export the global objects.
- #ifdef KP_ARCH_WIDE
- /* On 64-bit platforms, we use an additional bit to store the sign.
- * This helps us to avoid allocating a full value when negating it. */
- #define DEF_FLOAT(ptr) \
- ptrtype (ensure_mask ((ptr), TYPE_SHIFT - 1), typecode::FLOAT)
- FLT_PINF = DEF_FLOAT (&SF_PINF);
- FLT_NINF = FLT_PINF | SIGN_BIT;
- FLT_QNAN = DEF_FLOAT (&SF_QNAN);
- FLT_ZERO = DEF_FLOAT (&SF_ZERO);
- # undef DEF_FLOAT
- #else
- SF_NINF.val = -SF_PINF.val;
- FLT_PINF = SF_PINF.as_obj ();
- FLT_NINF = SF_NINF.as_obj ();
- FLT_QNAN = SF_QNAN.as_obj ();
- FLT_ZERO = SF_ZERO.as_obj ();
- #endif
- return (init_op::result_ok);
- }
- init_op init_float (do_init_float, "float");
- KP_DECLS_END
|