123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765376637673768376937703771377237733774377537763777377837793780378137823783378437853786378737883789379037913792379337943795379637973798379938003801380238033804380538063807380838093810381138123813381438153816381738183819382038213822382338243825382638273828382938303831383238333834383538363837383838393840384138423843384438453846384738483849385038513852385338543855385638573858385938603861386238633864386538663867386838693870387138723873387438753876387738783879388038813882388338843885388638873888388938903891389238933894389538963897389838993900390139023903390439053906390739083909 |
- /* Common base code for the decNumber C Library.
- Copyright (C) 2007-2018 Free Software Foundation, Inc.
- Contributed by IBM Corporation. Author Mike Cowlishaw.
- This file is part of GCC.
- GCC is free software; you can redistribute it and/or modify it under
- the terms of the GNU General Public License as published by the Free
- Software Foundation; either version 3, or (at your option) any later
- version.
- GCC 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 General Public License
- for more details.
- Under Section 7 of GPL version 3, you are granted additional
- permissions described in the GCC Runtime Library Exception, version
- 3.1, as published by the Free Software Foundation.
- You should have received a copy of the GNU General Public License and
- a copy of the GCC Runtime Library Exception along with this program;
- see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
- <http://www.gnu.org/licenses/>. */
- /* ------------------------------------------------------------------ */
- /* decBasic.c -- common base code for Basic decimal types */
- /* ------------------------------------------------------------------ */
- /* This module comprises code that is shared between decDouble and */
- /* decQuad (but not decSingle). The main arithmetic operations are */
- /* here (Add, Subtract, Multiply, FMA, and Division operators). */
- /* */
- /* Unlike decNumber, parameterization takes place at compile time */
- /* rather than at runtime. The parameters are set in the decDouble.c */
- /* (etc.) files, which then include this one to produce the compiled */
- /* code. The functions here, therefore, are code shared between */
- /* multiple formats. */
- /* */
- /* This must be included after decCommon.c. */
- /* ------------------------------------------------------------------ */
- /* Names here refer to decFloat rather than to decDouble, etc., and */
- /* the functions are in strict alphabetical order. */
- /* The compile-time flags SINGLE, DOUBLE, and QUAD are set up in */
- /* decCommon.c */
- #if !defined(QUAD)
- #error decBasic.c must be included after decCommon.c
- #endif
- #if SINGLE
- #error Routines in decBasic.c are for decDouble and decQuad only
- #endif
- /* Private constants */
- #define DIVIDE 0x80000000 /* Divide operations [as flags] */
- #define REMAINDER 0x40000000 /* .. */
- #define DIVIDEINT 0x20000000 /* .. */
- #define REMNEAR 0x10000000 /* .. */
- /* Private functions (local, used only by routines in this module) */
- static decFloat *decDivide(decFloat *, const decFloat *,
- const decFloat *, decContext *, uInt);
- static decFloat *decCanonical(decFloat *, const decFloat *);
- static void decFiniteMultiply(bcdnum *, uByte *, const decFloat *,
- const decFloat *);
- static decFloat *decInfinity(decFloat *, const decFloat *);
- static decFloat *decInvalid(decFloat *, decContext *);
- static decFloat *decNaNs(decFloat *, const decFloat *, const decFloat *,
- decContext *);
- static Int decNumCompare(const decFloat *, const decFloat *, Flag);
- static decFloat *decToIntegral(decFloat *, const decFloat *, decContext *,
- enum rounding, Flag);
- static uInt decToInt32(const decFloat *, decContext *, enum rounding,
- Flag, Flag);
- /* ------------------------------------------------------------------ */
- /* decCanonical -- copy a decFloat, making canonical */
- /* */
- /* result gets the canonicalized df */
- /* df is the decFloat to copy and make canonical */
- /* returns result */
- /* */
- /* This is exposed via decFloatCanonical for Double and Quad only. */
- /* This works on specials, too; no error or exception is possible. */
- /* ------------------------------------------------------------------ */
- static decFloat * decCanonical(decFloat *result, const decFloat *df) {
- uInt encode, precode, dpd; /* work */
- uInt inword, uoff, canon; /* .. */
- Int n; /* counter (down) */
- if (df!=result) *result=*df; /* effect copy if needed */
- if (DFISSPECIAL(result)) {
- if (DFISINF(result)) return decInfinity(result, df); /* clean Infinity */
- /* is a NaN */
- DFWORD(result, 0)&=~ECONNANMASK; /* clear ECON except selector */
- if (DFISCCZERO(df)) return result; /* coefficient continuation is 0 */
- /* drop through to check payload */
- }
- /* return quickly if the coefficient continuation is canonical */
- { /* declare block */
- #if DOUBLE
- uInt sourhi=DFWORD(df, 0);
- uInt sourlo=DFWORD(df, 1);
- if (CANONDPDOFF(sourhi, 8)
- && CANONDPDTWO(sourhi, sourlo, 30)
- && CANONDPDOFF(sourlo, 20)
- && CANONDPDOFF(sourlo, 10)
- && CANONDPDOFF(sourlo, 0)) return result;
- #elif QUAD
- uInt sourhi=DFWORD(df, 0);
- uInt sourmh=DFWORD(df, 1);
- uInt sourml=DFWORD(df, 2);
- uInt sourlo=DFWORD(df, 3);
- if (CANONDPDOFF(sourhi, 4)
- && CANONDPDTWO(sourhi, sourmh, 26)
- && CANONDPDOFF(sourmh, 16)
- && CANONDPDOFF(sourmh, 6)
- && CANONDPDTWO(sourmh, sourml, 28)
- && CANONDPDOFF(sourml, 18)
- && CANONDPDOFF(sourml, 8)
- && CANONDPDTWO(sourml, sourlo, 30)
- && CANONDPDOFF(sourlo, 20)
- && CANONDPDOFF(sourlo, 10)
- && CANONDPDOFF(sourlo, 0)) return result;
- #endif
- } /* block */
- /* Loop to repair a non-canonical coefficent, as needed */
- inword=DECWORDS-1; /* current input word */
- uoff=0; /* bit offset of declet */
- encode=DFWORD(result, inword);
- for (n=DECLETS-1; n>=0; n--) { /* count down declets of 10 bits */
- dpd=encode>>uoff;
- uoff+=10;
- if (uoff>32) { /* crossed uInt boundary */
- inword--;
- encode=DFWORD(result, inword);
- uoff-=32;
- dpd|=encode<<(10-uoff); /* get pending bits */
- }
- dpd&=0x3ff; /* clear uninteresting bits */
- if (dpd<0x16e) continue; /* must be canonical */
- canon=BIN2DPD[DPD2BIN[dpd]]; /* determine canonical declet */
- if (canon==dpd) continue; /* have canonical declet */
- /* need to replace declet */
- if (uoff>=10) { /* all within current word */
- encode&=~(0x3ff<<(uoff-10)); /* clear the 10 bits ready for replace */
- encode|=canon<<(uoff-10); /* insert the canonical form */
- DFWORD(result, inword)=encode; /* .. and save */
- continue;
- }
- /* straddled words */
- precode=DFWORD(result, inword+1); /* get previous */
- precode&=0xffffffff>>(10-uoff); /* clear top bits */
- DFWORD(result, inword+1)=precode|(canon<<(32-(10-uoff)));
- encode&=0xffffffff<<uoff; /* clear bottom bits */
- encode|=canon>>(10-uoff); /* insert canonical */
- DFWORD(result, inword)=encode; /* .. and save */
- } /* n */
- return result;
- } /* decCanonical */
- /* ------------------------------------------------------------------ */
- /* decDivide -- divide operations */
- /* */
- /* result gets the result of dividing dfl by dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* op is the operation selector */
- /* returns result */
- /* */
- /* op is one of DIVIDE, REMAINDER, DIVIDEINT, or REMNEAR. */
- /* ------------------------------------------------------------------ */
- #define DIVCOUNT 0 /* 1 to instrument subtractions counter */
- #define DIVBASE ((uInt)BILLION) /* the base used for divide */
- #define DIVOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
- #define DIVACCLEN (DIVOPLEN*3) /* accumulator length (ditto) */
- static decFloat * decDivide(decFloat *result, const decFloat *dfl,
- const decFloat *dfr, decContext *set, uInt op) {
- decFloat quotient; /* for remainders */
- bcdnum num; /* for final conversion */
- uInt acc[DIVACCLEN]; /* coefficent in base-billion .. */
- uInt div[DIVOPLEN]; /* divisor in base-billion .. */
- uInt quo[DIVOPLEN+1]; /* quotient in base-billion .. */
- uByte bcdacc[(DIVOPLEN+1)*9+2]; /* for quotient in BCD, +1, +1 */
- uInt *msua, *msud, *msuq; /* -> msu of acc, div, and quo */
- Int divunits, accunits; /* lengths */
- Int quodigits; /* digits in quotient */
- uInt *lsua, *lsuq; /* -> current acc and quo lsus */
- Int length, multiplier; /* work */
- uInt carry, sign; /* .. */
- uInt *ua, *ud, *uq; /* .. */
- uByte *ub; /* .. */
- uInt uiwork; /* for macros */
- uInt divtop; /* top unit of div adjusted for estimating */
- #if DIVCOUNT
- static uInt maxcount=0; /* worst-seen subtractions count */
- uInt divcount=0; /* subtractions count [this divide] */
- #endif
- /* calculate sign */
- num.sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
- if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
- /* NaNs are handled as usual */
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- /* one or two infinities */
- if (DFISINF(dfl)) {
- if (DFISINF(dfr)) return decInvalid(result, set); /* Two infinities bad */
- if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* as is rem */
- /* Infinity/x is infinite and quiet, even if x=0 */
- DFWORD(result, 0)=num.sign;
- return decInfinity(result, result);
- }
- /* must be x/Infinity -- remainders are lhs */
- if (op&(REMAINDER|REMNEAR)) return decCanonical(result, dfl);
- /* divides: return zero with correct sign and exponent depending */
- /* on op (Etiny for divide, 0 for divideInt) */
- decFloatZero(result);
- if (op==DIVIDEINT) DFWORD(result, 0)|=num.sign; /* add sign */
- else DFWORD(result, 0)=num.sign; /* zeros the exponent, too */
- return result;
- }
- /* next, handle zero operands (x/0 and 0/x) */
- if (DFISZERO(dfr)) { /* x/0 */
- if (DFISZERO(dfl)) { /* 0/0 is undefined */
- decFloatZero(result);
- DFWORD(result, 0)=DECFLOAT_qNaN;
- set->status|=DEC_Division_undefined;
- return result;
- }
- if (op&(REMAINDER|REMNEAR)) return decInvalid(result, set); /* bad rem */
- set->status|=DEC_Division_by_zero;
- DFWORD(result, 0)=num.sign;
- return decInfinity(result, result); /* x/0 -> signed Infinity */
- }
- num.exponent=GETEXPUN(dfl)-GETEXPUN(dfr); /* ideal exponent */
- if (DFISZERO(dfl)) { /* 0/x (x!=0) */
- /* if divide, result is 0 with ideal exponent; divideInt has */
- /* exponent=0, remainders give zero with lower exponent */
- if (op&DIVIDEINT) {
- decFloatZero(result);
- DFWORD(result, 0)|=num.sign; /* add sign */
- return result;
- }
- if (!(op&DIVIDE)) { /* a remainder */
- /* exponent is the minimum of the operands */
- num.exponent=MINI(GETEXPUN(dfl), GETEXPUN(dfr));
- /* if the result is zero the sign shall be sign of dfl */
- num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
- }
- bcdacc[0]=0;
- num.msd=bcdacc; /* -> 0 */
- num.lsd=bcdacc; /* .. */
- return decFinalize(result, &num, set); /* [divide may clamp exponent] */
- } /* 0/x */
- /* [here, both operands are known to be finite and non-zero] */
- /* extract the operand coefficents into 'units' which are */
- /* base-billion; the lhs is high-aligned in acc and the msu of both */
- /* acc and div is at the right-hand end of array (offset length-1); */
- /* the quotient can need one more unit than the operands as digits */
- /* in it are not necessarily aligned neatly; further, the quotient */
- /* may not start accumulating until after the end of the initial */
- /* operand in acc if that is small (e.g., 1) so the accumulator */
- /* must have at least that number of units extra (at the ls end) */
- GETCOEFFBILL(dfl, acc+DIVACCLEN-DIVOPLEN);
- GETCOEFFBILL(dfr, div);
- /* zero the low uInts of acc */
- acc[0]=0;
- acc[1]=0;
- acc[2]=0;
- acc[3]=0;
- #if DOUBLE
- #if DIVOPLEN!=2
- #error Unexpected Double DIVOPLEN
- #endif
- #elif QUAD
- acc[4]=0;
- acc[5]=0;
- acc[6]=0;
- acc[7]=0;
- #if DIVOPLEN!=4
- #error Unexpected Quad DIVOPLEN
- #endif
- #endif
- /* set msu and lsu pointers */
- msua=acc+DIVACCLEN-1; /* [leading zeros removed below] */
- msuq=quo+DIVOPLEN;
- /*[loop for div will terminate because operands are non-zero] */
- for (msud=div+DIVOPLEN-1; *msud==0;) msud--;
- /* the initial least-significant unit of acc is set so acc appears */
- /* to have the same length as div. */
- /* This moves one position towards the least possible for each */
- /* iteration */
- divunits=(Int)(msud-div+1); /* precalculate */
- lsua=msua-divunits+1; /* initial working lsu of acc */
- lsuq=msuq; /* and of quo */
- /* set up the estimator for the multiplier; this is the msu of div, */
- /* plus two bits from the unit below (if any) rounded up by one if */
- /* there are any non-zero bits or units below that [the extra two */
- /* bits makes for a much better estimate when the top unit is small] */
- divtop=*msud<<2;
- if (divunits>1) {
- uInt *um=msud-1;
- uInt d=*um;
- if (d>=750000000) {divtop+=3; d-=750000000;}
- else if (d>=500000000) {divtop+=2; d-=500000000;}
- else if (d>=250000000) {divtop++; d-=250000000;}
- if (d) divtop++;
- else for (um--; um>=div; um--) if (*um) {
- divtop++;
- break;
- }
- } /* >1 unit */
- #if DECTRACE
- {Int i;
- printf("----- div=");
- for (i=divunits-1; i>=0; i--) printf("%09ld ", (LI)div[i]);
- printf("\n");}
- #endif
- /* now collect up to DECPMAX+1 digits in the quotient (this may */
- /* need OPLEN+1 uInts if unaligned) */
- quodigits=0; /* no digits yet */
- for (;; lsua--) { /* outer loop -- each input position */
- #if DECCHECK
- if (lsua<acc) {
- printf("Acc underrun...\n");
- break;
- }
- #endif
- #if DECTRACE
- printf("Outer: quodigits=%ld acc=", (LI)quodigits);
- for (ua=msua; ua>=lsua; ua--) printf("%09ld ", (LI)*ua);
- printf("\n");
- #endif
- *lsuq=0; /* default unit result is 0 */
- for (;;) { /* inner loop -- calculate quotient unit */
- /* strip leading zero units from acc (either there initially or */
- /* from subtraction below); this may strip all if exactly 0 */
- for (; *msua==0 && msua>=lsua;) msua--;
- accunits=(Int)(msua-lsua+1); /* [maybe 0] */
- /* subtraction is only necessary and possible if there are as */
- /* least as many units remaining in acc for this iteration as */
- /* there are in div */
- if (accunits<divunits) {
- if (accunits==0) msua++; /* restore */
- break;
- }
- /* If acc is longer than div then subtraction is definitely */
- /* possible (as msu of both is non-zero), but if they are the */
- /* same length a comparison is needed. */
- /* If a subtraction is needed then a good estimate of the */
- /* multiplier for the subtraction is also needed in order to */
- /* minimise the iterations of this inner loop because the */
- /* subtractions needed dominate division performance. */
- if (accunits==divunits) {
- /* compare the high divunits of acc and div: */
- /* acc<div: this quotient unit is unchanged; subtraction */
- /* will be possible on the next iteration */
- /* acc==div: quotient gains 1, set acc=0 */
- /* acc>div: subtraction necessary at this position */
- for (ud=msud, ua=msua; ud>div; ud--, ua--) if (*ud!=*ua) break;
- /* [now at first mismatch or lsu] */
- if (*ud>*ua) break; /* next time... */
- if (*ud==*ua) { /* all compared equal */
- *lsuq+=1; /* increment result */
- msua=lsua; /* collapse acc units */
- *msua=0; /* .. to a zero */
- break;
- }
- /* subtraction necessary; estimate multiplier [see above] */
- /* if both *msud and *msua are small it is cost-effective to */
- /* bring in part of the following units (if any) to get a */
- /* better estimate (assume some other non-zero in div) */
- #define DIVLO 1000000U
- #define DIVHI (DIVBASE/DIVLO)
- #if DECUSE64
- if (divunits>1) {
- /* there cannot be a *(msud-2) for DECDOUBLE so next is */
- /* an exact calculation unless DECQUAD (which needs to */
- /* assume bits out there if divunits>2) */
- uLong mul=(uLong)*msua * DIVBASE + *(msua-1);
- uLong div=(uLong)*msud * DIVBASE + *(msud-1);
- #if QUAD
- if (divunits>2) div++;
- #endif
- mul/=div;
- multiplier=(Int)mul;
- }
- else multiplier=*msua/(*msud);
- #else
- if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
- multiplier=(*msua*DIVHI + *(msua-1)/DIVLO)
- /(*msud*DIVHI + *(msud-1)/DIVLO +1);
- }
- else multiplier=(*msua<<2)/divtop;
- #endif
- }
- else { /* accunits>divunits */
- /* msud is one unit 'lower' than msua, so estimate differently */
- #if DECUSE64
- uLong mul;
- /* as before, bring in extra digits if possible */
- if (divunits>1 && *msua<DIVLO && *msud<DIVLO) {
- mul=((uLong)*msua * DIVHI * DIVBASE) + *(msua-1) * DIVHI
- + *(msua-2)/DIVLO;
- mul/=(*msud*DIVHI + *(msud-1)/DIVLO +1);
- }
- else if (divunits==1) {
- mul=(uLong)*msua * DIVBASE + *(msua-1);
- mul/=*msud; /* no more to the right */
- }
- else {
- mul=(uLong)(*msua) * (uInt)(DIVBASE<<2)
- + (*(msua-1)<<2);
- mul/=divtop; /* [divtop already allows for sticky bits] */
- }
- multiplier=(Int)mul;
- #else
- multiplier=*msua * ((DIVBASE<<2)/divtop);
- #endif
- }
- if (multiplier==0) multiplier=1; /* marginal case */
- *lsuq+=multiplier;
- #if DIVCOUNT
- /* printf("Multiplier: %ld\n", (LI)multiplier); */
- divcount++;
- #endif
- /* Carry out the subtraction acc-(div*multiplier); for each */
- /* unit in div, do the multiply, split to units (see */
- /* decFloatMultiply for the algorithm), and subtract from acc */
- #define DIVMAGIC 2305843009U /* 2**61/10**9 */
- #define DIVSHIFTA 29
- #define DIVSHIFTB 32
- carry=0;
- for (ud=div, ua=lsua; ud<=msud; ud++, ua++) {
- uInt lo, hop;
- #if DECUSE64
- uLong sub=(uLong)multiplier*(*ud)+carry;
- if (sub<DIVBASE) {
- carry=0;
- lo=(uInt)sub;
- }
- else {
- hop=(uInt)(sub>>DIVSHIFTA);
- carry=(uInt)(((uLong)hop*DIVMAGIC)>>DIVSHIFTB);
- /* the estimate is now in hi; now calculate sub-hi*10**9 */
- /* to get the remainder (which will be <DIVBASE)) */
- lo=(uInt)sub;
- lo-=carry*DIVBASE; /* low word of result */
- if (lo>=DIVBASE) {
- lo-=DIVBASE; /* correct by +1 */
- carry++;
- }
- }
- #else /* 32-bit */
- uInt hi;
- /* calculate multiplier*(*ud) into hi and lo */
- LONGMUL32HI(hi, *ud, multiplier); /* get the high word */
- lo=multiplier*(*ud); /* .. and the low */
- lo+=carry; /* add the old hi */
- carry=hi+(lo<carry); /* .. with any carry */
- if (carry || lo>=DIVBASE) { /* split is needed */
- hop=(carry<<3)+(lo>>DIVSHIFTA); /* hi:lo/2**29 */
- LONGMUL32HI(carry, hop, DIVMAGIC); /* only need the high word */
- /* [DIVSHIFTB is 32, so carry can be used directly] */
- /* the estimate is now in carry; now calculate hi:lo-est*10**9; */
- /* happily the top word of the result is irrelevant because it */
- /* will always be zero so this needs only one multiplication */
- lo-=(carry*DIVBASE);
- /* the correction here will be at most +1; do it */
- if (lo>=DIVBASE) {
- lo-=DIVBASE;
- carry++;
- }
- }
- #endif
- if (lo>*ua) { /* borrow needed */
- *ua+=DIVBASE;
- carry++;
- }
- *ua-=lo;
- } /* ud loop */
- if (carry) *ua-=carry; /* accdigits>divdigits [cannot borrow] */
- } /* inner loop */
- /* the outer loop terminates when there is either an exact result */
- /* or enough digits; first update the quotient digit count and */
- /* pointer (if any significant digits) */
- #if DECTRACE
- if (*lsuq || quodigits) printf("*lsuq=%09ld\n", (LI)*lsuq);
- #endif
- if (quodigits) {
- quodigits+=9; /* had leading unit earlier */
- lsuq--;
- if (quodigits>DECPMAX+1) break; /* have enough */
- }
- else if (*lsuq) { /* first quotient digits */
- const uInt *pow;
- for (pow=DECPOWERS; *lsuq>=*pow; pow++) quodigits++;
- lsuq--;
- /* [cannot have >DECPMAX+1 on first unit] */
- }
- if (*msua!=0) continue; /* not an exact result */
- /* acc is zero iff used all of original units and zero down to lsua */
- /* (must also continue to original lsu for correct quotient length) */
- if (lsua>acc+DIVACCLEN-DIVOPLEN) continue;
- for (; msua>lsua && *msua==0;) msua--;
- if (*msua==0 && msua==lsua) break;
- } /* outer loop */
- /* all of the original operand in acc has been covered at this point */
- /* quotient now has at least DECPMAX+2 digits */
- /* *msua is now non-0 if inexact and sticky bits */
- /* lsuq is one below the last uint of the quotient */
- lsuq++; /* set -> true lsu of quo */
- if (*msua) *lsuq|=1; /* apply sticky bit */
- /* quo now holds the (unrounded) quotient in base-billion; one */
- /* base-billion 'digit' per uInt. */
- #if DECTRACE
- printf("DivQuo:");
- for (uq=msuq; uq>=lsuq; uq--) printf(" %09ld", (LI)*uq);
- printf("\n");
- #endif
- /* Now convert to BCD for rounding and cleanup, starting from the */
- /* most significant end [offset by one into bcdacc to leave room */
- /* for a possible carry digit if rounding for REMNEAR is needed] */
- for (uq=msuq, ub=bcdacc+1; uq>=lsuq; uq--, ub+=9) {
- uInt top, mid, rem; /* work */
- if (*uq==0) { /* no split needed */
- UBFROMUI(ub, 0); /* clear 9 BCD8s */
- UBFROMUI(ub+4, 0); /* .. */
- *(ub+8)=0; /* .. */
- continue;
- }
- /* *uq is non-zero -- split the base-billion digit into */
- /* hi, mid, and low three-digits */
- #define divsplit9 1000000 /* divisor */
- #define divsplit6 1000 /* divisor */
- /* The splitting is done by simple divides and remainders, */
- /* assuming the compiler will optimize these [GCC does] */
- top=*uq/divsplit9;
- rem=*uq%divsplit9;
- mid=rem/divsplit6;
- rem=rem%divsplit6;
- /* lay out the nine BCD digits (plus one unwanted byte) */
- UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
- UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
- UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
- } /* BCD conversion loop */
- ub--; /* -> lsu */
- /* complete the bcdnum; quodigits is correct, so the position of */
- /* the first non-zero is known */
- num.msd=bcdacc+1+(msuq-lsuq+1)*9-quodigits;
- num.lsd=ub;
- /* make exponent adjustments, etc */
- if (lsua<acc+DIVACCLEN-DIVOPLEN) { /* used extra digits */
- num.exponent-=(Int)((acc+DIVACCLEN-DIVOPLEN-lsua)*9);
- /* if the result was exact then there may be up to 8 extra */
- /* trailing zeros in the overflowed quotient final unit */
- if (*msua==0) {
- for (; *ub==0;) ub--; /* drop zeros */
- num.exponent+=(Int)(num.lsd-ub); /* and adjust exponent */
- num.lsd=ub;
- }
- } /* adjustment needed */
- #if DIVCOUNT
- if (divcount>maxcount) { /* new high-water nark */
- maxcount=divcount;
- printf("DivNewMaxCount: %ld\n", (LI)maxcount);
- }
- #endif
- if (op&DIVIDE) return decFinalize(result, &num, set); /* all done */
- /* Is DIVIDEINT or a remainder; there is more to do -- first form */
- /* the integer (this is done 'after the fact', unlike as in */
- /* decNumber, so as not to tax DIVIDE) */
- /* The first non-zero digit will be in the first 9 digits, known */
- /* from quodigits and num.msd, so there is always space for DECPMAX */
- /* digits */
- length=(Int)(num.lsd-num.msd+1);
- /*printf("Length exp: %ld %ld\n", (LI)length, (LI)num.exponent); */
- if (length+num.exponent>DECPMAX) { /* cannot fit */
- decFloatZero(result);
- DFWORD(result, 0)=DECFLOAT_qNaN;
- set->status|=DEC_Division_impossible;
- return result;
- }
- if (num.exponent>=0) { /* already an int, or need pad zeros */
- for (ub=num.lsd+1; ub<=num.lsd+num.exponent; ub++) *ub=0;
- num.lsd+=num.exponent;
- }
- else { /* too long: round or truncate needed */
- Int drop=-num.exponent;
- if (!(op&REMNEAR)) { /* simple truncate */
- num.lsd-=drop;
- if (num.lsd<num.msd) { /* truncated all */
- num.lsd=num.msd; /* make 0 */
- *num.lsd=0; /* .. [sign still relevant] */
- }
- }
- else { /* round to nearest even [sigh] */
- /* round-to-nearest, in-place; msd is at or to right of bcdacc+1 */
- /* (this is a special case of Quantize -- q.v. for commentary) */
- uByte *roundat; /* -> re-round digit */
- uByte reround; /* reround value */
- *(num.msd-1)=0; /* in case of left carry, or make 0 */
- if (drop<length) roundat=num.lsd-drop+1;
- else if (drop==length) roundat=num.msd;
- else roundat=num.msd-1; /* [-> 0] */
- reround=*roundat;
- for (ub=roundat+1; ub<=num.lsd; ub++) {
- if (*ub!=0) {
- reround=DECSTICKYTAB[reround];
- break;
- }
- } /* check stickies */
- if (roundat>num.msd) num.lsd=roundat-1;
- else {
- num.msd--; /* use the 0 .. */
- num.lsd=num.msd; /* .. at the new MSD place */
- }
- if (reround!=0) { /* discarding non-zero */
- uInt bump=0;
- /* rounding is DEC_ROUND_HALF_EVEN always */
- if (reround>5) bump=1; /* >0.5 goes up */
- else if (reround==5) /* exactly 0.5000 .. */
- bump=*(num.lsd) & 0x01; /* .. up iff [new] lsd is odd */
- if (bump!=0) { /* need increment */
- /* increment the coefficient; this might end up with 1000... */
- ub=num.lsd;
- for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
- for (; *ub==9; ub--) *ub=0; /* at most 3 more */
- *ub+=1;
- if (ub<num.msd) num.msd--; /* carried */
- } /* bump needed */
- } /* reround!=0 */
- } /* remnear */
- } /* round or truncate needed */
- num.exponent=0; /* all paths */
- /*decShowNum(&num, "int"); */
- if (op&DIVIDEINT) return decFinalize(result, &num, set); /* all done */
- /* Have a remainder to calculate */
- decFinalize("ient, &num, set); /* lay out the integer so far */
- DFWORD("ient, 0)^=DECFLOAT_Sign; /* negate it */
- sign=DFWORD(dfl, 0); /* save sign of dfl */
- decFloatFMA(result, "ient, dfr, dfl, set);
- if (!DFISZERO(result)) return result;
- /* if the result is zero the sign shall be sign of dfl */
- DFWORD("ient, 0)=sign; /* construct decFloat of sign */
- return decFloatCopySign(result, result, "ient);
- } /* decDivide */
- /* ------------------------------------------------------------------ */
- /* decFiniteMultiply -- multiply two finite decFloats */
- /* */
- /* num gets the result of multiplying dfl and dfr */
- /* bcdacc .. with the coefficient in this array */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* */
- /* This effects the multiplication of two decFloats, both known to be */
- /* finite, leaving the result in a bcdnum ready for decFinalize (for */
- /* use in Multiply) or in a following addition (FMA). */
- /* */
- /* bcdacc must have space for at least DECPMAX9*18+1 bytes. */
- /* No error is possible and no status is set. */
- /* ------------------------------------------------------------------ */
- /* This routine has two separate implementations of the core */
- /* multiplication; both using base-billion. One uses only 32-bit */
- /* variables (Ints and uInts) or smaller; the other uses uLongs (for */
- /* multiplication and addition only). Both implementations cover */
- /* both arithmetic sizes (DOUBLE and QUAD) in order to allow timing */
- /* comparisons. In any one compilation only one implementation for */
- /* each size can be used, and if DECUSE64 is 0 then use of the 32-bit */
- /* version is forced. */
- /* */
- /* Historical note: an earlier version of this code also supported the */
- /* 256-bit format and has been preserved. That is somewhat trickier */
- /* during lazy carry splitting because the initial quotient estimate */
- /* (est) can exceed 32 bits. */
- #define MULTBASE ((uInt)BILLION) /* the base used for multiply */
- #define MULOPLEN DECPMAX9 /* operand length ('digits' base 10**9) */
- #define MULACCLEN (MULOPLEN*2) /* accumulator length (ditto) */
- #define LEADZEROS (MULACCLEN*9 - DECPMAX*2) /* leading zeros always */
- /* Assertions: exponent not too large and MULACCLEN is a multiple of 4 */
- #if DECEMAXD>9
- #error Exponent may overflow when doubled for Multiply
- #endif
- #if MULACCLEN!=(MULACCLEN/4)*4
- /* This assumption is used below only for initialization */
- #error MULACCLEN is not a multiple of 4
- #endif
- static void decFiniteMultiply(bcdnum *num, uByte *bcdacc,
- const decFloat *dfl, const decFloat *dfr) {
- uInt bufl[MULOPLEN]; /* left coefficient (base-billion) */
- uInt bufr[MULOPLEN]; /* right coefficient (base-billion) */
- uInt *ui, *uj; /* work */
- uByte *ub; /* .. */
- uInt uiwork; /* for macros */
- #if DECUSE64
- uLong accl[MULACCLEN]; /* lazy accumulator (base-billion+) */
- uLong *pl; /* work -> lazy accumulator */
- uInt acc[MULACCLEN]; /* coefficent in base-billion .. */
- #else
- uInt acc[MULACCLEN*2]; /* accumulator in base-billion .. */
- #endif
- uInt *pa; /* work -> accumulator */
- /*printf("Base10**9: OpLen=%d MulAcclen=%d\n", OPLEN, MULACCLEN); */
- /* Calculate sign and exponent */
- num->sign=(DFWORD(dfl, 0)^DFWORD(dfr, 0)) & DECFLOAT_Sign;
- num->exponent=GETEXPUN(dfl)+GETEXPUN(dfr); /* [see assertion above] */
- /* Extract the coefficients and prepare the accumulator */
- /* the coefficients of the operands are decoded into base-billion */
- /* numbers in uInt arrays (bufl and bufr, LSD at offset 0) of the */
- /* appropriate size. */
- GETCOEFFBILL(dfl, bufl);
- GETCOEFFBILL(dfr, bufr);
- #if DECTRACE && 0
- printf("CoeffbL:");
- for (ui=bufl+MULOPLEN-1; ui>=bufl; ui--) printf(" %08lx", (LI)*ui);
- printf("\n");
- printf("CoeffbR:");
- for (uj=bufr+MULOPLEN-1; uj>=bufr; uj--) printf(" %08lx", (LI)*uj);
- printf("\n");
- #endif
- /* start the 64-bit/32-bit differing paths... */
- #if DECUSE64
- /* zero the accumulator */
- #if MULACCLEN==4
- accl[0]=0; accl[1]=0; accl[2]=0; accl[3]=0;
- #else /* use a loop */
- /* MULACCLEN is a multiple of four, asserted above */
- for (pl=accl; pl<accl+MULACCLEN; pl+=4) {
- *pl=0; *(pl+1)=0; *(pl+2)=0; *(pl+3)=0;/* [reduce overhead] */
- } /* pl */
- #endif
- /* Effect the multiplication */
- /* The multiplcation proceeds using MFC's lazy-carry resolution */
- /* algorithm from decNumber. First, the multiplication is */
- /* effected, allowing accumulation of the partial products (which */
- /* are in base-billion at each column position) into 64 bits */
- /* without resolving back to base=billion after each addition. */
- /* These 64-bit numbers (which may contain up to 19 decimal digits) */
- /* are then split using the Clark & Cowlishaw algorithm (see below). */
- /* [Testing for 0 in the inner loop is not really a 'win'] */
- for (ui=bufr; ui<bufr+MULOPLEN; ui++) { /* over each item in rhs */
- if (*ui==0) continue; /* product cannot affect result */
- pl=accl+(ui-bufr); /* where to add the lhs */
- for (uj=bufl; uj<bufl+MULOPLEN; uj++, pl++) { /* over each item in lhs */
- /* if (*uj==0) continue; // product cannot affect result */
- *pl+=((uLong)*ui)*(*uj);
- } /* uj */
- } /* ui */
- /* The 64-bit carries must now be resolved; this means that a */
- /* quotient/remainder has to be calculated for base-billion (1E+9). */
- /* For this, Clark & Cowlishaw's quotient estimation approach (also */
- /* used in decNumber) is needed, because 64-bit divide is generally */
- /* extremely slow on 32-bit machines, and may be slower than this */
- /* approach even on 64-bit machines. This algorithm splits X */
- /* using: */
- /* */
- /* magic=2**(A+B)/1E+9; // 'magic number' */
- /* hop=X/2**A; // high order part of X (by shift) */
- /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
- /* */
- /* A and B are quite constrained; hop and magic must fit in 32 bits, */
- /* and 2**(A+B) must be as large as possible (which is 2**61 if */
- /* magic is to fit). Further, maxX increases with the length of */
- /* the operands (and hence the number of partial products */
- /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
- /* */
- /* It can be shown that when OPLEN is 2 then the maximum error in */
- /* the estimated quotient is <1, but for larger maximum x the */
- /* maximum error is above 1 so a correction that is >1 may be */
- /* needed. Values of A and B are chosen to satisfy the constraints */
- /* just mentioned while minimizing the maximum error (and hence the */
- /* maximum correction), as shown in the following table: */
- /* */
- /* Type OPLEN A B maxX maxError maxCorrection */
- /* --------------------------------------------------------- */
- /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
- /* QUAD 4 30 31 <4*10**18 1.17 2 */
- /* */
- /* In the OPLEN==2 case there is most choice, but the value for B */
- /* of 32 has a big advantage as then the calculation of the */
- /* estimate requires no shifting; the compiler can extract the high */
- /* word directly after multiplying magic*hop. */
- #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
- #if DOUBLE
- #define MULSHIFTA 29
- #define MULSHIFTB 32
- #elif QUAD
- #define MULSHIFTA 30
- #define MULSHIFTB 31
- #else
- #error Unexpected type
- #endif
- #if DECTRACE
- printf("MulAccl:");
- for (pl=accl+MULACCLEN-1; pl>=accl; pl--)
- printf(" %08lx:%08lx", (LI)(*pl>>32), (LI)(*pl&0xffffffff));
- printf("\n");
- #endif
- for (pl=accl, pa=acc; pl<accl+MULACCLEN; pl++, pa++) { /* each column position */
- uInt lo, hop; /* work */
- uInt est; /* cannot exceed 4E+9 */
- if (*pl>=MULTBASE) {
- /* *pl holds a binary number which needs to be split */
- hop=(uInt)(*pl>>MULSHIFTA);
- est=(uInt)(((uLong)hop*MULMAGIC)>>MULSHIFTB);
- /* the estimate is now in est; now calculate hi:lo-est*10**9; */
- /* happily the top word of the result is irrelevant because it */
- /* will always be zero so this needs only one multiplication */
- lo=(uInt)(*pl-((uLong)est*MULTBASE)); /* low word of result */
- /* If QUAD, the correction here could be +2 */
- if (lo>=MULTBASE) {
- lo-=MULTBASE; /* correct by +1 */
- est++;
- #if QUAD
- /* may need to correct by +2 */
- if (lo>=MULTBASE) {
- lo-=MULTBASE;
- est++;
- }
- #endif
- }
- /* finally place lo as the new coefficient 'digit' and add est to */
- /* the next place up [this is safe because this path is never */
- /* taken on the final iteration as *pl will fit] */
- *pa=lo;
- *(pl+1)+=est;
- } /* *pl needed split */
- else { /* *pl<MULTBASE */
- *pa=(uInt)*pl; /* just copy across */
- }
- } /* pl loop */
- #else /* 32-bit */
- for (pa=acc;; pa+=4) { /* zero the accumulator */
- *pa=0; *(pa+1)=0; *(pa+2)=0; *(pa+3)=0; /* [reduce overhead] */
- if (pa==acc+MULACCLEN*2-4) break; /* multiple of 4 asserted */
- } /* pa */
- /* Effect the multiplication */
- /* uLongs are not available (and in particular, there is no uLong */
- /* divide) but it is still possible to use MFC's lazy-carry */
- /* resolution algorithm from decNumber. First, the multiplication */
- /* is effected, allowing accumulation of the partial products */
- /* (which are in base-billion at each column position) into 64 bits */
- /* [with the high-order 32 bits in each position being held at */
- /* offset +ACCLEN from the low-order 32 bits in the accumulator]. */
- /* These 64-bit numbers (which may contain up to 19 decimal digits) */
- /* are then split using the Clark & Cowlishaw algorithm (see */
- /* below). */
- for (ui=bufr;; ui++) { /* over each item in rhs */
- uInt hi, lo; /* words of exact multiply result */
- pa=acc+(ui-bufr); /* where to add the lhs */
- for (uj=bufl;; uj++, pa++) { /* over each item in lhs */
- LONGMUL32HI(hi, *ui, *uj); /* calculate product of digits */
- lo=(*ui)*(*uj); /* .. */
- *pa+=lo; /* accumulate low bits and .. */
- *(pa+MULACCLEN)+=hi+(*pa<lo); /* .. high bits with any carry */
- if (uj==bufl+MULOPLEN-1) break;
- }
- if (ui==bufr+MULOPLEN-1) break;
- }
- /* The 64-bit carries must now be resolved; this means that a */
- /* quotient/remainder has to be calculated for base-billion (1E+9). */
- /* For this, Clark & Cowlishaw's quotient estimation approach (also */
- /* used in decNumber) is needed, because 64-bit divide is generally */
- /* extremely slow on 32-bit machines. This algorithm splits X */
- /* using: */
- /* */
- /* magic=2**(A+B)/1E+9; // 'magic number' */
- /* hop=X/2**A; // high order part of X (by shift) */
- /* est=magic*hop/2**B // quotient estimate (may be low by 1) */
- /* */
- /* A and B are quite constrained; hop and magic must fit in 32 bits, */
- /* and 2**(A+B) must be as large as possible (which is 2**61 if */
- /* magic is to fit). Further, maxX increases with the length of */
- /* the operands (and hence the number of partial products */
- /* accumulated); maxX is OPLEN*(10**18), which is up to 19 digits. */
- /* */
- /* It can be shown that when OPLEN is 2 then the maximum error in */
- /* the estimated quotient is <1, but for larger maximum x the */
- /* maximum error is above 1 so a correction that is >1 may be */
- /* needed. Values of A and B are chosen to satisfy the constraints */
- /* just mentioned while minimizing the maximum error (and hence the */
- /* maximum correction), as shown in the following table: */
- /* */
- /* Type OPLEN A B maxX maxError maxCorrection */
- /* --------------------------------------------------------- */
- /* DOUBLE 2 29 32 <2*10**18 0.63 1 */
- /* QUAD 4 30 31 <4*10**18 1.17 2 */
- /* */
- /* In the OPLEN==2 case there is most choice, but the value for B */
- /* of 32 has a big advantage as then the calculation of the */
- /* estimate requires no shifting; the high word is simply */
- /* calculated from multiplying magic*hop. */
- #define MULMAGIC 2305843009U /* 2**61/10**9 [both cases] */
- #if DOUBLE
- #define MULSHIFTA 29
- #define MULSHIFTB 32
- #elif QUAD
- #define MULSHIFTA 30
- #define MULSHIFTB 31
- #else
- #error Unexpected type
- #endif
- #if DECTRACE
- printf("MulHiLo:");
- for (pa=acc+MULACCLEN-1; pa>=acc; pa--)
- printf(" %08lx:%08lx", (LI)*(pa+MULACCLEN), (LI)*pa);
- printf("\n");
- #endif
- for (pa=acc;; pa++) { /* each low uInt */
- uInt hi, lo; /* words of exact multiply result */
- uInt hop, estlo; /* work */
- #if QUAD
- uInt esthi; /* .. */
- #endif
- lo=*pa;
- hi=*(pa+MULACCLEN); /* top 32 bits */
- /* hi and lo now hold a binary number which needs to be split */
- #if DOUBLE
- hop=(hi<<3)+(lo>>MULSHIFTA); /* hi:lo/2**29 */
- LONGMUL32HI(estlo, hop, MULMAGIC);/* only need the high word */
- /* [MULSHIFTB is 32, so estlo can be used directly] */
- /* the estimate is now in estlo; now calculate hi:lo-est*10**9; */
- /* happily the top word of the result is irrelevant because it */
- /* will always be zero so this needs only one multiplication */
- lo-=(estlo*MULTBASE);
- /* esthi=0; // high word is ignored below */
- /* the correction here will be at most +1; do it */
- if (lo>=MULTBASE) {
- lo-=MULTBASE;
- estlo++;
- }
- #elif QUAD
- hop=(hi<<2)+(lo>>MULSHIFTA); /* hi:lo/2**30 */
- LONGMUL32HI(esthi, hop, MULMAGIC);/* shift will be 31 .. */
- estlo=hop*MULMAGIC; /* .. so low word needed */
- estlo=(esthi<<1)+(estlo>>MULSHIFTB); /* [just the top bit] */
- /* esthi=0; // high word is ignored below */
- lo-=(estlo*MULTBASE); /* as above */
- /* the correction here could be +1 or +2 */
- if (lo>=MULTBASE) {
- lo-=MULTBASE;
- estlo++;
- }
- if (lo>=MULTBASE) {
- lo-=MULTBASE;
- estlo++;
- }
- #else
- #error Unexpected type
- #endif
- /* finally place lo as the new accumulator digit and add est to */
- /* the next place up; this latter add could cause a carry of 1 */
- /* to the high word of the next place */
- *pa=lo;
- *(pa+1)+=estlo;
- /* esthi is always 0 for DOUBLE and QUAD so this is skipped */
- /* *(pa+1+MULACCLEN)+=esthi; */
- if (*(pa+1)<estlo) *(pa+1+MULACCLEN)+=1; /* carry */
- if (pa==acc+MULACCLEN-2) break; /* [MULACCLEN-1 will never need split] */
- } /* pa loop */
- #endif
- /* At this point, whether using the 64-bit or the 32-bit paths, the */
- /* accumulator now holds the (unrounded) result in base-billion; */
- /* one base-billion 'digit' per uInt. */
- #if DECTRACE
- printf("MultAcc:");
- for (pa=acc+MULACCLEN-1; pa>=acc; pa--) printf(" %09ld", (LI)*pa);
- printf("\n");
- #endif
- /* Now convert to BCD for rounding and cleanup, starting from the */
- /* most significant end */
- pa=acc+MULACCLEN-1;
- if (*pa!=0) num->msd=bcdacc+LEADZEROS;/* drop known lead zeros */
- else { /* >=1 word of leading zeros */
- num->msd=bcdacc; /* known leading zeros are gone */
- pa--; /* skip first word .. */
- for (; *pa==0; pa--) if (pa==acc) break; /* .. and any more leading 0s */
- }
- for (ub=bcdacc;; pa--, ub+=9) {
- if (*pa!=0) { /* split(s) needed */
- uInt top, mid, rem; /* work */
- /* *pa is non-zero -- split the base-billion acc digit into */
- /* hi, mid, and low three-digits */
- #define mulsplit9 1000000 /* divisor */
- #define mulsplit6 1000 /* divisor */
- /* The splitting is done by simple divides and remainders, */
- /* assuming the compiler will optimize these where useful */
- /* [GCC does] */
- top=*pa/mulsplit9;
- rem=*pa%mulsplit9;
- mid=rem/mulsplit6;
- rem=rem%mulsplit6;
- /* lay out the nine BCD digits (plus one unwanted byte) */
- UBFROMUI(ub, UBTOUI(&BIN2BCD8[top*4]));
- UBFROMUI(ub+3, UBTOUI(&BIN2BCD8[mid*4]));
- UBFROMUI(ub+6, UBTOUI(&BIN2BCD8[rem*4]));
- }
- else { /* *pa==0 */
- UBFROMUI(ub, 0); /* clear 9 BCD8s */
- UBFROMUI(ub+4, 0); /* .. */
- *(ub+8)=0; /* .. */
- }
- if (pa==acc) break;
- } /* BCD conversion loop */
- num->lsd=ub+8; /* complete the bcdnum .. */
- #if DECTRACE
- decShowNum(num, "postmult");
- decFloatShow(dfl, "dfl");
- decFloatShow(dfr, "dfr");
- #endif
- return;
- } /* decFiniteMultiply */
- /* ------------------------------------------------------------------ */
- /* decFloatAbs -- absolute value, heeding NaNs, etc. */
- /* */
- /* result gets the canonicalized df with sign 0 */
- /* df is the decFloat to abs */
- /* set is the context */
- /* returns result */
- /* */
- /* This has the same effect as decFloatPlus unless df is negative, */
- /* in which case it has the same effect as decFloatMinus. The */
- /* effect is also the same as decFloatCopyAbs except that NaNs are */
- /* handled normally (the sign of a NaN is not affected, and an sNaN */
- /* will signal) and the result will be canonical. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatAbs(decFloat *result, const decFloat *df,
- decContext *set) {
- if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
- decCanonical(result, df); /* copy and check */
- DFBYTE(result, 0)&=~0x80; /* zero sign bit */
- return result;
- } /* decFloatAbs */
- /* ------------------------------------------------------------------ */
- /* decFloatAdd -- add two decFloats */
- /* */
- /* result gets the result of adding dfl and dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- #if QUAD
- /* Table for testing MSDs for fastpath elimination; returns the MSD of */
- /* a decDouble or decQuad (top 6 bits tested) ignoring the sign. */
- /* Infinities return -32 and NaNs return -128 so that summing the two */
- /* MSDs also allows rapid tests for the Specials (see code below). */
- const Int DECTESTMSD[64]={
- 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
- 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128,
- 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7,
- 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 8, 9, 8, 9, -32, -128};
- #else
- /* The table for testing MSDs is shared between the modules */
- extern const Int DECTESTMSD[64];
- #endif
- decFloat * decFloatAdd(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- bcdnum num; /* for final conversion */
- Int bexpl, bexpr; /* left and right biased exponents */
- uByte *ub, *us, *ut; /* work */
- uInt uiwork; /* for macros */
- #if QUAD
- uShort uswork; /* .. */
- #endif
- uInt sourhil, sourhir; /* top words from source decFloats */
- /* [valid only through end of */
- /* fastpath code -- before swap] */
- uInt diffsign; /* non-zero if signs differ */
- uInt carry; /* carry: 0 or 1 before add loop */
- Int overlap; /* coefficient overlap (if full) */
- Int summ; /* sum of the MSDs */
- /* the following buffers hold coefficients with various alignments */
- /* (see commentary and diagrams below) */
- uByte acc[4+2+DECPMAX*3+8];
- uByte buf[4+2+DECPMAX*2];
- uByte *umsd, *ulsd; /* local MSD and LSD pointers */
- #if DECLITEND
- #define CARRYPAT 0x01000000 /* carry=1 pattern */
- #else
- #define CARRYPAT 0x00000001 /* carry=1 pattern */
- #endif
- /* Start decoding the arguments */
- /* The initial exponents are placed into the opposite Ints to */
- /* that which might be expected; there are two sets of data to */
- /* keep track of (each decFloat and the corresponding exponent), */
- /* and this scheme means that at the swap point (after comparing */
- /* exponents) only one pair of words needs to be swapped */
- /* whichever path is taken (thereby minimising worst-case path). */
- /* The calculated exponents will be nonsense when the arguments are */
- /* Special, but are not used in that path */
- sourhil=DFWORD(dfl, 0); /* LHS top word */
- summ=DECTESTMSD[sourhil>>26]; /* get first MSD for testing */
- bexpr=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
- bexpr+=GETECON(dfl); /* .. + continuation */
- sourhir=DFWORD(dfr, 0); /* RHS top word */
- summ+=DECTESTMSD[sourhir>>26]; /* sum MSDs for testing */
- bexpl=DECCOMBEXP[sourhir>>26];
- bexpl+=GETECON(dfr);
- /* here bexpr has biased exponent from lhs, and vice versa */
- diffsign=(sourhil^sourhir)&DECFLOAT_Sign;
- /* now determine whether to take a fast path or the full-function */
- /* slow path. The slow path must be taken when: */
- /* -- both numbers are finite, and: */
- /* the exponents are different, or */
- /* the signs are different, or */
- /* the sum of the MSDs is >8 (hence might overflow) */
- /* specialness and the sum of the MSDs can be tested at once using */
- /* the summ value just calculated, so the test for specials is no */
- /* longer on the worst-case path (as of 3.60) */
- if (summ<=8) { /* MSD+MSD is good, or there is a special */
- if (summ<0) { /* there is a special */
- /* Inf+Inf would give -64; Inf+finite is -32 or higher */
- if (summ<-64) return decNaNs(result, dfl, dfr, set); /* one or two NaNs */
- /* two infinities with different signs is invalid */
- if (summ==-64 && diffsign) return decInvalid(result, set);
- if (DFISINF(dfl)) return decInfinity(result, dfl); /* LHS is infinite */
- return decInfinity(result, dfr); /* RHS must be Inf */
- }
- /* Here when both arguments are finite; fast path is possible */
- /* (currently only for aligned and same-sign) */
- if (bexpr==bexpl && !diffsign) {
- uInt tac[DECLETS+1]; /* base-1000 coefficient */
- uInt encode; /* work */
- /* Get one coefficient as base-1000 and add the other */
- GETCOEFFTHOU(dfl, tac); /* least-significant goes to [0] */
- ADDCOEFFTHOU(dfr, tac);
- /* here the sum of the MSDs (plus any carry) will be <10 due to */
- /* the fastpath test earlier */
- /* construct the result; low word is the same for both formats */
- encode =BIN2DPD[tac[0]];
- encode|=BIN2DPD[tac[1]]<<10;
- encode|=BIN2DPD[tac[2]]<<20;
- encode|=BIN2DPD[tac[3]]<<30;
- DFWORD(result, (DECBYTES/4)-1)=encode;
- /* collect next two declets (all that remains, for Double) */
- encode =BIN2DPD[tac[3]]>>2;
- encode|=BIN2DPD[tac[4]]<<8;
- #if QUAD
- /* complete and lay out middling words */
- encode|=BIN2DPD[tac[5]]<<18;
- encode|=BIN2DPD[tac[6]]<<28;
- DFWORD(result, 2)=encode;
- encode =BIN2DPD[tac[6]]>>4;
- encode|=BIN2DPD[tac[7]]<<6;
- encode|=BIN2DPD[tac[8]]<<16;
- encode|=BIN2DPD[tac[9]]<<26;
- DFWORD(result, 1)=encode;
- /* and final two declets */
- encode =BIN2DPD[tac[9]]>>6;
- encode|=BIN2DPD[tac[10]]<<4;
- #endif
- /* add exponent continuation and sign (from either argument) */
- encode|=sourhil & (ECONMASK | DECFLOAT_Sign);
- /* create lookup index = MSD + top two bits of biased exponent <<4 */
- tac[DECLETS]|=(bexpl>>DECECONL)<<4;
- encode|=DECCOMBFROM[tac[DECLETS]]; /* add constructed combination field */
- DFWORD(result, 0)=encode; /* complete */
- /* decFloatShow(result, ">"); */
- return result;
- } /* fast path OK */
- /* drop through to slow path */
- } /* low sum or Special(s) */
- /* Slow path required -- arguments are finite and might overflow, */
- /* or require alignment, or might have different signs */
- /* now swap either exponents or argument pointers */
- if (bexpl<=bexpr) {
- /* original left is bigger */
- Int bexpswap=bexpl;
- bexpl=bexpr;
- bexpr=bexpswap;
- /* printf("left bigger\n"); */
- }
- else {
- const decFloat *dfswap=dfl;
- dfl=dfr;
- dfr=dfswap;
- /* printf("right bigger\n"); */
- }
- /* [here dfl and bexpl refer to the datum with the larger exponent, */
- /* of if the exponents are equal then the original LHS argument] */
- /* if lhs is zero then result will be the rhs (now known to have */
- /* the smaller exponent), which also may need to be tested for zero */
- /* for the weird IEEE 754 sign rules */
- if (DFISZERO(dfl)) {
- decCanonical(result, dfr); /* clean copy */
- /* "When the sum of two operands with opposite signs is */
- /* exactly zero, the sign of that sum shall be '+' in all */
- /* rounding modes except round toward -Infinity, in which */
- /* mode that sign shall be '-'." */
- if (diffsign && DFISZERO(result)) {
- DFWORD(result, 0)&=~DECFLOAT_Sign; /* assume sign 0 */
- if (set->round==DEC_ROUND_FLOOR) DFWORD(result, 0)|=DECFLOAT_Sign;
- }
- return result;
- } /* numfl is zero */
- /* [here, LHS is non-zero; code below assumes that] */
- /* Coefficients layout during the calculations to follow: */
- /* */
- /* Overlap case: */
- /* +------------------------------------------------+ */
- /* acc: |0000| coeffa | tail B | | */
- /* +------------------------------------------------+ */
- /* buf: |0000| pad0s | coeffb | | */
- /* +------------------------------------------------+ */
- /* */
- /* Touching coefficients or gap: */
- /* +------------------------------------------------+ */
- /* acc: |0000| coeffa | gap | coeffb | */
- /* +------------------------------------------------+ */
- /* [buf not used or needed; gap clamped to Pmax] */
- /* lay out lhs coefficient into accumulator; this starts at acc+4 */
- /* for decDouble or acc+6 for decQuad so the LSD is word- */
- /* aligned; the top word gap is there only in case a carry digit */
- /* is prefixed after the add -- it does not need to be zeroed */
- #if DOUBLE
- #define COFF 4 /* offset into acc */
- #elif QUAD
- UBFROMUS(acc+4, 0); /* prefix 00 */
- #define COFF 6 /* offset into acc */
- #endif
- GETCOEFF(dfl, acc+COFF); /* decode from decFloat */
- ulsd=acc+COFF+DECPMAX-1;
- umsd=acc+4; /* [having this here avoids */
- #if DECTRACE
- {bcdnum tum;
- tum.msd=umsd;
- tum.lsd=ulsd;
- tum.exponent=bexpl-DECBIAS;
- tum.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
- decShowNum(&tum, "dflx");}
- #endif
- /* if signs differ, take ten's complement of lhs (here the */
- /* coefficient is subtracted from all-nines; the 1 is added during */
- /* the later add cycle -- zeros to the right do not matter because */
- /* the complement of zero is zero); these are fixed-length inverts */
- /* where the lsd is known to be at a 4-byte boundary (so no borrow */
- /* possible) */
- carry=0; /* assume no carry */
- if (diffsign) {
- carry=CARRYPAT; /* for +1 during add */
- UBFROMUI(acc+ 4, 0x09090909-UBTOUI(acc+ 4));
- UBFROMUI(acc+ 8, 0x09090909-UBTOUI(acc+ 8));
- UBFROMUI(acc+12, 0x09090909-UBTOUI(acc+12));
- UBFROMUI(acc+16, 0x09090909-UBTOUI(acc+16));
- #if QUAD
- UBFROMUI(acc+20, 0x09090909-UBTOUI(acc+20));
- UBFROMUI(acc+24, 0x09090909-UBTOUI(acc+24));
- UBFROMUI(acc+28, 0x09090909-UBTOUI(acc+28));
- UBFROMUI(acc+32, 0x09090909-UBTOUI(acc+32));
- UBFROMUI(acc+36, 0x09090909-UBTOUI(acc+36));
- #endif
- } /* diffsign */
- /* now process the rhs coefficient; if it cannot overlap lhs then */
- /* it can be put straight into acc (with an appropriate gap, if */
- /* needed) because no actual addition will be needed (except */
- /* possibly to complete ten's complement) */
- overlap=DECPMAX-(bexpl-bexpr);
- #if DECTRACE
- printf("exps: %ld %ld\n", (LI)(bexpl-DECBIAS), (LI)(bexpr-DECBIAS));
- printf("Overlap=%ld carry=%08lx\n", (LI)overlap, (LI)carry);
- #endif
- if (overlap<=0) { /* no overlap possible */
- uInt gap; /* local work */
- /* since a full addition is not needed, a ten's complement */
- /* calculation started above may need to be completed */
- if (carry) {
- for (ub=ulsd; *ub==9; ub--) *ub=0;
- *ub+=1;
- carry=0; /* taken care of */
- }
- /* up to DECPMAX-1 digits of the final result can extend down */
- /* below the LSD of the lhs, so if the gap is >DECPMAX then the */
- /* rhs will be simply sticky bits. In this case the gap is */
- /* clamped to DECPMAX and the exponent adjusted to suit [this is */
- /* safe because the lhs is non-zero]. */
- gap=-overlap;
- if (gap>DECPMAX) {
- bexpr+=gap-1;
- gap=DECPMAX;
- }
- ub=ulsd+gap+1; /* where MSD will go */
- /* Fill the gap with 0s; note that there is no addition to do */
- ut=acc+COFF+DECPMAX; /* start of gap */
- for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* mind the gap */
- if (overlap<-DECPMAX) { /* gap was > DECPMAX */
- *ub=(uByte)(!DFISZERO(dfr)); /* make sticky digit */
- }
- else { /* need full coefficient */
- GETCOEFF(dfr, ub); /* decode from decFloat */
- ub+=DECPMAX-1; /* new LSD... */
- }
- ulsd=ub; /* save new LSD */
- } /* no overlap possible */
- else { /* overlap>0 */
- /* coefficients overlap (perhaps completely, although also */
- /* perhaps only where zeros) */
- if (overlap==DECPMAX) { /* aligned */
- ub=buf+COFF; /* where msd will go */
- #if QUAD
- UBFROMUS(buf+4, 0); /* clear quad's 00 */
- #endif
- GETCOEFF(dfr, ub); /* decode from decFloat */
- }
- else { /* unaligned */
- ub=buf+COFF+DECPMAX-overlap; /* where MSD will go */
- /* Fill the prefix gap with 0s; 8 will cover most common */
- /* unalignments, so start with direct assignments (a loop is */
- /* then used for any remaining -- the loop (and the one in a */
- /* moment) is not then on the critical path because the number */
- /* of additions is reduced by (at least) two in this case) */
- UBFROMUI(buf+4, 0); /* [clears decQuad 00 too] */
- UBFROMUI(buf+8, 0);
- if (ub>buf+12) {
- ut=buf+12; /* start any remaining */
- for (; ut<ub; ut+=4) UBFROMUI(ut, 0); /* fill them */
- }
- GETCOEFF(dfr, ub); /* decode from decFloat */
- /* now move tail of rhs across to main acc; again use direct */
- /* copies for 8 digits-worth */
- UBFROMUI(acc+COFF+DECPMAX, UBTOUI(buf+COFF+DECPMAX));
- UBFROMUI(acc+COFF+DECPMAX+4, UBTOUI(buf+COFF+DECPMAX+4));
- if (buf+COFF+DECPMAX+8<ub+DECPMAX) {
- us=buf+COFF+DECPMAX+8; /* source */
- ut=acc+COFF+DECPMAX+8; /* target */
- for (; us<ub+DECPMAX; us+=4, ut+=4) UBFROMUI(ut, UBTOUI(us));
- }
- } /* unaligned */
- ulsd=acc+(ub-buf+DECPMAX-1); /* update LSD pointer */
- /* Now do the add of the non-tail; this is all nicely aligned, */
- /* and is over a multiple of four digits (because for Quad two */
- /* zero digits were added on the left); words in both acc and */
- /* buf (buf especially) will often be zero */
- /* [byte-by-byte add, here, is about 15% slower total effect than */
- /* the by-fours] */
- /* Now effect the add; this is harder on a little-endian */
- /* machine as the inter-digit carry cannot use the usual BCD */
- /* addition trick because the bytes are loaded in the wrong order */
- /* [this loop could be unrolled, but probably scarcely worth it] */
- ut=acc+COFF+DECPMAX-4; /* target LSW (acc) */
- us=buf+COFF+DECPMAX-4; /* source LSW (buf, to add to acc) */
- #if !DECLITEND
- for (; ut>=acc+4; ut-=4, us-=4) { /* big-endian add loop */
- /* bcd8 add */
- carry+=UBTOUI(us); /* rhs + carry */
- if (carry==0) continue; /* no-op */
- carry+=UBTOUI(ut); /* lhs */
- /* Big-endian BCD adjust (uses internal carry) */
- carry+=0x76f6f6f6; /* note top nibble not all bits */
- /* apply BCD adjust and save */
- UBFROMUI(ut, (carry & 0x0f0f0f0f) - ((carry & 0x60606060)>>4));
- carry>>=31; /* true carry was at far left */
- } /* add loop */
- #else
- for (; ut>=acc+4; ut-=4, us-=4) { /* little-endian add loop */
- /* bcd8 add */
- carry+=UBTOUI(us); /* rhs + carry */
- if (carry==0) continue; /* no-op [common if unaligned] */
- carry+=UBTOUI(ut); /* lhs */
- /* Little-endian BCD adjust; inter-digit carry must be manual */
- /* because the lsb from the array will be in the most-significant */
- /* byte of carry */
- carry+=0x76767676; /* note no inter-byte carries */
- carry+=(carry & 0x80000000)>>15;
- carry+=(carry & 0x00800000)>>15;
- carry+=(carry & 0x00008000)>>15;
- carry-=(carry & 0x60606060)>>4; /* BCD adjust back */
- UBFROMUI(ut, carry & 0x0f0f0f0f); /* clear debris and save */
- /* here, final carry-out bit is at 0x00000080; move it ready */
- /* for next word-add (i.e., to 0x01000000) */
- carry=(carry & 0x00000080)<<17;
- } /* add loop */
- #endif
- #if DECTRACE
- {bcdnum tum;
- printf("Add done, carry=%08lx, diffsign=%ld\n", (LI)carry, (LI)diffsign);
- tum.msd=umsd; /* acc+4; */
- tum.lsd=ulsd;
- tum.exponent=0;
- tum.sign=0;
- decShowNum(&tum, "dfadd");}
- #endif
- } /* overlap possible */
- /* ordering here is a little strange in order to have slowest path */
- /* first in GCC asm listing */
- if (diffsign) { /* subtraction */
- if (!carry) { /* no carry out means RHS<LHS */
- /* borrowed -- take ten's complement */
- /* sign is lhs sign */
- num.sign=DFWORD(dfl, 0) & DECFLOAT_Sign;
- /* invert the coefficient first by fours, then add one; space */
- /* at the end of the buffer ensures the by-fours is always */
- /* safe, but lsd+1 must be cleared to prevent a borrow */
- /* if big-endian */
- #if !DECLITEND
- *(ulsd+1)=0;
- #endif
- /* there are always at least four coefficient words */
- UBFROMUI(umsd, 0x09090909-UBTOUI(umsd));
- UBFROMUI(umsd+4, 0x09090909-UBTOUI(umsd+4));
- UBFROMUI(umsd+8, 0x09090909-UBTOUI(umsd+8));
- UBFROMUI(umsd+12, 0x09090909-UBTOUI(umsd+12));
- #if DOUBLE
- #define BNEXT 16
- #elif QUAD
- UBFROMUI(umsd+16, 0x09090909-UBTOUI(umsd+16));
- UBFROMUI(umsd+20, 0x09090909-UBTOUI(umsd+20));
- UBFROMUI(umsd+24, 0x09090909-UBTOUI(umsd+24));
- UBFROMUI(umsd+28, 0x09090909-UBTOUI(umsd+28));
- UBFROMUI(umsd+32, 0x09090909-UBTOUI(umsd+32));
- #define BNEXT 36
- #endif
- if (ulsd>=umsd+BNEXT) { /* unaligned */
- /* eight will handle most unaligments for Double; 16 for Quad */
- UBFROMUI(umsd+BNEXT, 0x09090909-UBTOUI(umsd+BNEXT));
- UBFROMUI(umsd+BNEXT+4, 0x09090909-UBTOUI(umsd+BNEXT+4));
- #if DOUBLE
- #define BNEXTY (BNEXT+8)
- #elif QUAD
- UBFROMUI(umsd+BNEXT+8, 0x09090909-UBTOUI(umsd+BNEXT+8));
- UBFROMUI(umsd+BNEXT+12, 0x09090909-UBTOUI(umsd+BNEXT+12));
- #define BNEXTY (BNEXT+16)
- #endif
- if (ulsd>=umsd+BNEXTY) { /* very unaligned */
- ut=umsd+BNEXTY; /* -> continue */
- for (;;ut+=4) {
- UBFROMUI(ut, 0x09090909-UBTOUI(ut)); /* invert four digits */
- if (ut>=ulsd-3) break; /* all done */
- }
- }
- }
- /* complete the ten's complement by adding 1 */
- for (ub=ulsd; *ub==9; ub--) *ub=0;
- *ub+=1;
- } /* borrowed */
- else { /* carry out means RHS>=LHS */
- num.sign=DFWORD(dfr, 0) & DECFLOAT_Sign;
- /* all done except for the special IEEE 754 exact-zero-result */
- /* rule (see above); while testing for zero, strip leading */
- /* zeros (which will save decFinalize doing it) (this is in */
- /* diffsign path, so carry impossible and true umsd is */
- /* acc+COFF) */
- /* Check the initial coefficient area using the fast macro; */
- /* this will often be all that needs to be done (as on the */
- /* worst-case path when the subtraction was aligned and */
- /* full-length) */
- if (ISCOEFFZERO(acc+COFF)) {
- umsd=acc+COFF+DECPMAX-1; /* so far, so zero */
- if (ulsd>umsd) { /* more to check */
- umsd++; /* to align after checked area */
- for (; UBTOUI(umsd)==0 && umsd+3<ulsd;) umsd+=4;
- for (; *umsd==0 && umsd<ulsd;) umsd++;
- }
- if (*umsd==0) { /* must be true zero (and diffsign) */
- num.sign=0; /* assume + */
- if (set->round==DEC_ROUND_FLOOR) num.sign=DECFLOAT_Sign;
- }
- }
- /* [else was not zero, might still have leading zeros] */
- } /* subtraction gave positive result */
- } /* diffsign */
- else { /* same-sign addition */
- num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
- #if DOUBLE
- if (carry) { /* only possible with decDouble */
- *(acc+3)=1; /* [Quad has leading 00] */
- umsd=acc+3;
- }
- #endif
- } /* same sign */
- num.msd=umsd; /* set MSD .. */
- num.lsd=ulsd; /* .. and LSD */
- num.exponent=bexpr-DECBIAS; /* set exponent to smaller, unbiassed */
- #if DECTRACE
- decFloatShow(dfl, "dfl");
- decFloatShow(dfr, "dfr");
- decShowNum(&num, "postadd");
- #endif
- return decFinalize(result, &num, set); /* round, check, and lay out */
- } /* decFloatAdd */
- /* ------------------------------------------------------------------ */
- /* decFloatAnd -- logical digitwise AND of two decFloats */
- /* */
- /* result gets the result of ANDing dfl and dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result, which will be canonical with sign=0 */
- /* */
- /* The operands must be positive, finite with exponent q=0, and */
- /* comprise just zeros and ones; if not, Invalid operation results. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatAnd(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
- || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
- /* the operands are positive finite integers (q=0) with just 0s and 1s */
- #if DOUBLE
- DFWORD(result, 0)=ZEROWORD
- |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04009124);
- DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x49124491;
- #elif QUAD
- DFWORD(result, 0)=ZEROWORD
- |((DFWORD(dfl, 0) & DFWORD(dfr, 0))&0x04000912);
- DFWORD(result, 1)=(DFWORD(dfl, 1) & DFWORD(dfr, 1))&0x44912449;
- DFWORD(result, 2)=(DFWORD(dfl, 2) & DFWORD(dfr, 2))&0x12449124;
- DFWORD(result, 3)=(DFWORD(dfl, 3) & DFWORD(dfr, 3))&0x49124491;
- #endif
- return result;
- } /* decFloatAnd */
- /* ------------------------------------------------------------------ */
- /* decFloatCanonical -- copy a decFloat, making canonical */
- /* */
- /* result gets the canonicalized df */
- /* df is the decFloat to copy and make canonical */
- /* returns result */
- /* */
- /* This works on specials, too; no error or exception is possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCanonical(decFloat *result, const decFloat *df) {
- return decCanonical(result, df);
- } /* decFloatCanonical */
- /* ------------------------------------------------------------------ */
- /* decFloatClass -- return the class of a decFloat */
- /* */
- /* df is the decFloat to test */
- /* returns the decClass that df falls into */
- /* ------------------------------------------------------------------ */
- enum decClass decFloatClass(const decFloat *df) {
- Int exp; /* exponent */
- if (DFISSPECIAL(df)) {
- if (DFISQNAN(df)) return DEC_CLASS_QNAN;
- if (DFISSNAN(df)) return DEC_CLASS_SNAN;
- /* must be an infinity */
- if (DFISSIGNED(df)) return DEC_CLASS_NEG_INF;
- return DEC_CLASS_POS_INF;
- }
- if (DFISZERO(df)) { /* quite common */
- if (DFISSIGNED(df)) return DEC_CLASS_NEG_ZERO;
- return DEC_CLASS_POS_ZERO;
- }
- /* is finite and non-zero; similar code to decFloatIsNormal, here */
- /* [this could be speeded up slightly by in-lining decFloatDigits] */
- exp=GETEXPUN(df) /* get unbiased exponent .. */
- +decFloatDigits(df)-1; /* .. and make adjusted exponent */
- if (exp>=DECEMIN) { /* is normal */
- if (DFISSIGNED(df)) return DEC_CLASS_NEG_NORMAL;
- return DEC_CLASS_POS_NORMAL;
- }
- /* is subnormal */
- if (DFISSIGNED(df)) return DEC_CLASS_NEG_SUBNORMAL;
- return DEC_CLASS_POS_SUBNORMAL;
- } /* decFloatClass */
- /* ------------------------------------------------------------------ */
- /* decFloatClassString -- return the class of a decFloat as a string */
- /* */
- /* df is the decFloat to test */
- /* returns a constant string describing the class df falls into */
- /* ------------------------------------------------------------------ */
- const char *decFloatClassString(const decFloat *df) {
- enum decClass eclass=decFloatClass(df);
- if (eclass==DEC_CLASS_POS_NORMAL) return DEC_ClassString_PN;
- if (eclass==DEC_CLASS_NEG_NORMAL) return DEC_ClassString_NN;
- if (eclass==DEC_CLASS_POS_ZERO) return DEC_ClassString_PZ;
- if (eclass==DEC_CLASS_NEG_ZERO) return DEC_ClassString_NZ;
- if (eclass==DEC_CLASS_POS_SUBNORMAL) return DEC_ClassString_PS;
- if (eclass==DEC_CLASS_NEG_SUBNORMAL) return DEC_ClassString_NS;
- if (eclass==DEC_CLASS_POS_INF) return DEC_ClassString_PI;
- if (eclass==DEC_CLASS_NEG_INF) return DEC_ClassString_NI;
- if (eclass==DEC_CLASS_QNAN) return DEC_ClassString_QN;
- if (eclass==DEC_CLASS_SNAN) return DEC_ClassString_SN;
- return DEC_ClassString_UN; /* Unknown */
- } /* decFloatClassString */
- /* ------------------------------------------------------------------ */
- /* decFloatCompare -- compare two decFloats; quiet NaNs allowed */
- /* */
- /* result gets the result of comparing dfl and dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCompare(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int comp; /* work */
- /* NaNs are handled as usual */
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- /* numeric comparison needed */
- comp=decNumCompare(dfl, dfr, 0);
- decFloatZero(result);
- if (comp==0) return result;
- DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
- if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
- return result;
- } /* decFloatCompare */
- /* ------------------------------------------------------------------ */
- /* decFloatCompareSignal -- compare two decFloats; all NaNs signal */
- /* */
- /* result gets the result of comparing dfl and dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result, which may be -1, 0, 1, or NaN (Unordered) */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCompareSignal(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int comp; /* work */
- /* NaNs are handled as usual, except that all NaNs signal */
- if (DFISNAN(dfl) || DFISNAN(dfr)) {
- set->status|=DEC_Invalid_operation;
- return decNaNs(result, dfl, dfr, set);
- }
- /* numeric comparison needed */
- comp=decNumCompare(dfl, dfr, 0);
- decFloatZero(result);
- if (comp==0) return result;
- DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
- if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
- return result;
- } /* decFloatCompareSignal */
- /* ------------------------------------------------------------------ */
- /* decFloatCompareTotal -- compare two decFloats with total ordering */
- /* */
- /* result gets the result of comparing dfl and dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* returns result, which may be -1, 0, or 1 */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCompareTotal(decFloat *result,
- const decFloat *dfl, const decFloat *dfr) {
- Int comp; /* work */
- uInt uiwork; /* for macros */
- #if QUAD
- uShort uswork; /* .. */
- #endif
- if (DFISNAN(dfl) || DFISNAN(dfr)) {
- Int nanl, nanr; /* work */
- /* morph NaNs to +/- 1 or 2, leave numbers as 0 */
- nanl=DFISSNAN(dfl)+DFISQNAN(dfl)*2; /* quiet > signalling */
- if (DFISSIGNED(dfl)) nanl=-nanl;
- nanr=DFISSNAN(dfr)+DFISQNAN(dfr)*2;
- if (DFISSIGNED(dfr)) nanr=-nanr;
- if (nanl>nanr) comp=+1;
- else if (nanl<nanr) comp=-1;
- else { /* NaNs are the same type and sign .. must compare payload */
- /* buffers need +2 for QUAD */
- uByte bufl[DECPMAX+4]; /* for LHS coefficient + foot */
- uByte bufr[DECPMAX+4]; /* for RHS coefficient + foot */
- uByte *ub, *uc; /* work */
- Int sigl; /* signum of LHS */
- sigl=(DFISSIGNED(dfl) ? -1 : +1);
- /* decode the coefficients */
- /* (shift both right two if Quad to make a multiple of four) */
- #if QUAD
- UBFROMUS(bufl, 0);
- UBFROMUS(bufr, 0);
- #endif
- GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
- GETCOEFF(dfr, bufr+QUAD*2); /* .. */
- /* all multiples of four, here */
- comp=0; /* assume equal */
- for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
- uInt ui=UBTOUI(ub);
- if (ui==UBTOUI(uc)) continue; /* so far so same */
- /* about to find a winner; go by bytes in case little-endian */
- for (;; ub++, uc++) {
- if (*ub==*uc) continue;
- if (*ub>*uc) comp=sigl; /* difference found */
- else comp=-sigl; /* .. */
- break;
- }
- }
- } /* same NaN type and sign */
- }
- else {
- /* numeric comparison needed */
- comp=decNumCompare(dfl, dfr, 1); /* total ordering */
- }
- decFloatZero(result);
- if (comp==0) return result;
- DFBYTE(result, DECBYTES-1)=0x01; /* LSD=1 */
- if (comp<0) DFBYTE(result, 0)|=0x80; /* set sign bit */
- return result;
- } /* decFloatCompareTotal */
- /* ------------------------------------------------------------------ */
- /* decFloatCompareTotalMag -- compare magnitudes with total ordering */
- /* */
- /* result gets the result of comparing abs(dfl) and abs(dfr) */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* returns result, which may be -1, 0, or 1 */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCompareTotalMag(decFloat *result,
- const decFloat *dfl, const decFloat *dfr) {
- decFloat a, b; /* for copy if needed */
- /* copy and redirect signed operand(s) */
- if (DFISSIGNED(dfl)) {
- decFloatCopyAbs(&a, dfl);
- dfl=&a;
- }
- if (DFISSIGNED(dfr)) {
- decFloatCopyAbs(&b, dfr);
- dfr=&b;
- }
- return decFloatCompareTotal(result, dfl, dfr);
- } /* decFloatCompareTotalMag */
- /* ------------------------------------------------------------------ */
- /* decFloatCopy -- copy a decFloat as-is */
- /* */
- /* result gets the copy of dfl */
- /* dfl is the decFloat to copy */
- /* returns result */
- /* */
- /* This is a bitwise operation; no errors or exceptions are possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCopy(decFloat *result, const decFloat *dfl) {
- if (dfl!=result) *result=*dfl; /* copy needed */
- return result;
- } /* decFloatCopy */
- /* ------------------------------------------------------------------ */
- /* decFloatCopyAbs -- copy a decFloat as-is and set sign bit to 0 */
- /* */
- /* result gets the copy of dfl with sign bit 0 */
- /* dfl is the decFloat to copy */
- /* returns result */
- /* */
- /* This is a bitwise operation; no errors or exceptions are possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCopyAbs(decFloat *result, const decFloat *dfl) {
- if (dfl!=result) *result=*dfl; /* copy needed */
- DFBYTE(result, 0)&=~0x80; /* zero sign bit */
- return result;
- } /* decFloatCopyAbs */
- /* ------------------------------------------------------------------ */
- /* decFloatCopyNegate -- copy a decFloat as-is with inverted sign bit */
- /* */
- /* result gets the copy of dfl with sign bit inverted */
- /* dfl is the decFloat to copy */
- /* returns result */
- /* */
- /* This is a bitwise operation; no errors or exceptions are possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCopyNegate(decFloat *result, const decFloat *dfl) {
- if (dfl!=result) *result=*dfl; /* copy needed */
- DFBYTE(result, 0)^=0x80; /* invert sign bit */
- return result;
- } /* decFloatCopyNegate */
- /* ------------------------------------------------------------------ */
- /* decFloatCopySign -- copy a decFloat with the sign of another */
- /* */
- /* result gets the result of copying dfl with the sign of dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* returns result */
- /* */
- /* This is a bitwise operation; no errors or exceptions are possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatCopySign(decFloat *result,
- const decFloat *dfl, const decFloat *dfr) {
- uByte sign=(uByte)(DFBYTE(dfr, 0)&0x80); /* save sign bit */
- if (dfl!=result) *result=*dfl; /* copy needed */
- DFBYTE(result, 0)&=~0x80; /* clear sign .. */
- DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* .. and set saved */
- return result;
- } /* decFloatCopySign */
- /* ------------------------------------------------------------------ */
- /* decFloatDigits -- return the number of digits in a decFloat */
- /* */
- /* df is the decFloat to investigate */
- /* returns the number of significant digits in the decFloat; a */
- /* zero coefficient returns 1 as does an infinity (a NaN returns */
- /* the number of digits in the payload) */
- /* ------------------------------------------------------------------ */
- /* private macro to extract a declet according to provided formula */
- /* (form), and if it is non-zero then return the calculated digits */
- /* depending on the declet number (n), where n=0 for the most */
- /* significant declet; uses uInt dpd for work */
- #define dpdlenchk(n, form) {dpd=(form)&0x3ff; \
- if (dpd) return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
- /* next one is used when it is known that the declet must be */
- /* non-zero, or is the final zero declet */
- #define dpdlendun(n, form) {dpd=(form)&0x3ff; \
- if (dpd==0) return 1; \
- return (DECPMAX-1-3*(n))-(3-DPD2BCD8[dpd*4+3]);}
- uInt decFloatDigits(const decFloat *df) {
- uInt dpd; /* work */
- uInt sourhi=DFWORD(df, 0); /* top word from source decFloat */
- #if QUAD
- uInt sourmh, sourml;
- #endif
- uInt sourlo;
- if (DFISINF(df)) return 1;
- /* A NaN effectively has an MSD of 0; otherwise if non-zero MSD */
- /* then the coefficient is full-length */
- if (!DFISNAN(df) && DECCOMBMSD[sourhi>>26]) return DECPMAX;
- #if DOUBLE
- if (sourhi&0x0003ffff) { /* ends in first */
- dpdlenchk(0, sourhi>>8);
- sourlo=DFWORD(df, 1);
- dpdlendun(1, (sourhi<<2) | (sourlo>>30));
- } /* [cannot drop through] */
- sourlo=DFWORD(df, 1); /* sourhi not involved now */
- if (sourlo&0xfff00000) { /* in one of first two */
- dpdlenchk(1, sourlo>>30); /* very rare */
- dpdlendun(2, sourlo>>20);
- } /* [cannot drop through] */
- dpdlenchk(3, sourlo>>10);
- dpdlendun(4, sourlo);
- /* [cannot drop through] */
- #elif QUAD
- if (sourhi&0x00003fff) { /* ends in first */
- dpdlenchk(0, sourhi>>4);
- sourmh=DFWORD(df, 1);
- dpdlendun(1, ((sourhi)<<6) | (sourmh>>26));
- } /* [cannot drop through] */
- sourmh=DFWORD(df, 1);
- if (sourmh) {
- dpdlenchk(1, sourmh>>26);
- dpdlenchk(2, sourmh>>16);
- dpdlenchk(3, sourmh>>6);
- sourml=DFWORD(df, 2);
- dpdlendun(4, ((sourmh)<<4) | (sourml>>28));
- } /* [cannot drop through] */
- sourml=DFWORD(df, 2);
- if (sourml) {
- dpdlenchk(4, sourml>>28);
- dpdlenchk(5, sourml>>18);
- dpdlenchk(6, sourml>>8);
- sourlo=DFWORD(df, 3);
- dpdlendun(7, ((sourml)<<2) | (sourlo>>30));
- } /* [cannot drop through] */
- sourlo=DFWORD(df, 3);
- if (sourlo&0xfff00000) { /* in one of first two */
- dpdlenchk(7, sourlo>>30); /* very rare */
- dpdlendun(8, sourlo>>20);
- } /* [cannot drop through] */
- dpdlenchk(9, sourlo>>10);
- dpdlendun(10, sourlo);
- /* [cannot drop through] */
- #endif
- } /* decFloatDigits */
- /* ------------------------------------------------------------------ */
- /* decFloatDivide -- divide a decFloat by another */
- /* */
- /* result gets the result of dividing dfl by dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- /* This is just a wrapper. */
- decFloat * decFloatDivide(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- return decDivide(result, dfl, dfr, set, DIVIDE);
- } /* decFloatDivide */
- /* ------------------------------------------------------------------ */
- /* decFloatDivideInteger -- integer divide a decFloat by another */
- /* */
- /* result gets the result of dividing dfl by dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatDivideInteger(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- return decDivide(result, dfl, dfr, set, DIVIDEINT);
- } /* decFloatDivideInteger */
- /* ------------------------------------------------------------------ */
- /* decFloatFMA -- multiply and add three decFloats, fused */
- /* */
- /* result gets the result of (dfl*dfr)+dff with a single rounding */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* dff is the final decFloat (fhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatFMA(decFloat *result, const decFloat *dfl,
- const decFloat *dfr, const decFloat *dff,
- decContext *set) {
- /* The accumulator has the bytes needed for FiniteMultiply, plus */
- /* one byte to the left in case of carry, plus DECPMAX+2 to the */
- /* right for the final addition (up to full fhs + round & sticky) */
- #define FMALEN (ROUNDUP4(1+ (DECPMAX9*18+1) +DECPMAX+2))
- uByte acc[FMALEN]; /* for multiplied coefficient in BCD */
- /* .. and for final result */
- bcdnum mul; /* for multiplication result */
- bcdnum fin; /* for final operand, expanded */
- uByte coe[ROUNDUP4(DECPMAX)]; /* dff coefficient in BCD */
- bcdnum *hi, *lo; /* bcdnum with higher/lower exponent */
- uInt diffsign; /* non-zero if signs differ */
- uInt hipad; /* pad digit for hi if needed */
- Int padding; /* excess exponent */
- uInt carry; /* +1 for ten's complement and during add */
- uByte *ub, *uh, *ul; /* work */
- uInt uiwork; /* for macros */
- /* handle all the special values [any special operand leads to a */
- /* special result] */
- if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr) || DFISSPECIAL(dff)) {
- decFloat proxy; /* multiplication result proxy */
- /* NaNs are handled as usual, giving priority to sNaNs */
- if (DFISSNAN(dfl) || DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- if (DFISSNAN(dff)) return decNaNs(result, dff, NULL, set);
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- if (DFISNAN(dff)) return decNaNs(result, dff, NULL, set);
- /* One or more of the three is infinite */
- /* infinity times zero is bad */
- decFloatZero(&proxy);
- if (DFISINF(dfl)) {
- if (DFISZERO(dfr)) return decInvalid(result, set);
- decInfinity(&proxy, &proxy);
- }
- else if (DFISINF(dfr)) {
- if (DFISZERO(dfl)) return decInvalid(result, set);
- decInfinity(&proxy, &proxy);
- }
- /* compute sign of multiplication and place in proxy */
- DFWORD(&proxy, 0)|=(DFWORD(dfl, 0)^DFWORD(dfr, 0))&DECFLOAT_Sign;
- if (!DFISINF(dff)) return decFloatCopy(result, &proxy);
- /* dff is Infinite */
- if (!DFISINF(&proxy)) return decInfinity(result, dff);
- /* both sides of addition are infinite; different sign is bad */
- if ((DFWORD(dff, 0)&DECFLOAT_Sign)!=(DFWORD(&proxy, 0)&DECFLOAT_Sign))
- return decInvalid(result, set);
- return decFloatCopy(result, &proxy);
- }
- /* Here when all operands are finite */
- /* First multiply dfl*dfr */
- decFiniteMultiply(&mul, acc+1, dfl, dfr);
- /* The multiply is complete, exact and unbounded, and described in */
- /* mul with the coefficient held in acc[1...] */
- /* now add in dff; the algorithm is essentially the same as */
- /* decFloatAdd, but the code is different because the code there */
- /* is highly optimized for adding two numbers of the same size */
- fin.exponent=GETEXPUN(dff); /* get dff exponent and sign */
- fin.sign=DFWORD(dff, 0)&DECFLOAT_Sign;
- diffsign=mul.sign^fin.sign; /* note if signs differ */
- fin.msd=coe;
- fin.lsd=coe+DECPMAX-1;
- GETCOEFF(dff, coe); /* extract the coefficient */
- /* now set hi and lo so that hi points to whichever of mul and fin */
- /* has the higher exponent and lo points to the other [don't care, */
- /* if the same]. One coefficient will be in acc, the other in coe. */
- if (mul.exponent>=fin.exponent) {
- hi=&mul;
- lo=&fin;
- }
- else {
- hi=&fin;
- lo=&mul;
- }
- /* remove leading zeros on both operands; this will save time later */
- /* and make testing for zero trivial (tests are safe because acc */
- /* and coe are rounded up to uInts) */
- for (; UBTOUI(hi->msd)==0 && hi->msd+3<hi->lsd;) hi->msd+=4;
- for (; *hi->msd==0 && hi->msd<hi->lsd;) hi->msd++;
- for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
- for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
- /* if hi is zero then result will be lo (which has the smaller */
- /* exponent), which also may need to be tested for zero for the */
- /* weird IEEE 754 sign rules */
- if (*hi->msd==0) { /* hi is zero */
- /* "When the sum of two operands with opposite signs is */
- /* exactly zero, the sign of that sum shall be '+' in all */
- /* rounding modes except round toward -Infinity, in which */
- /* mode that sign shall be '-'." */
- if (diffsign) {
- if (*lo->msd==0) { /* lo is zero */
- lo->sign=0;
- if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
- } /* diffsign && lo=0 */
- } /* diffsign */
- return decFinalize(result, lo, set); /* may need clamping */
- } /* numfl is zero */
- /* [here, both are minimal length and hi is non-zero] */
- /* (if lo is zero then padding with zeros may be needed, below) */
- /* if signs differ, take the ten's complement of hi (zeros to the */
- /* right do not matter because the complement of zero is zero); the */
- /* +1 is done later, as part of the addition, inserted at the */
- /* correct digit */
- hipad=0;
- carry=0;
- if (diffsign) {
- hipad=9;
- carry=1;
- /* exactly the correct number of digits must be inverted */
- for (uh=hi->msd; uh<hi->lsd-3; uh+=4) UBFROMUI(uh, 0x09090909-UBTOUI(uh));
- for (; uh<=hi->lsd; uh++) *uh=(uByte)(0x09-*uh);
- }
- /* ready to add; note that hi has no leading zeros so gap */
- /* calculation does not have to be as pessimistic as in decFloatAdd */
- /* (this is much more like the arbitrary-precision algorithm in */
- /* Rexx and decNumber) */
- /* padding is the number of zeros that would need to be added to hi */
- /* for its lsd to be aligned with the lsd of lo */
- padding=hi->exponent-lo->exponent;
- /* printf("FMA pad %ld\n", (LI)padding); */
- /* the result of the addition will be built into the accumulator, */
- /* starting from the far right; this could be either hi or lo, and */
- /* will be aligned */
- ub=acc+FMALEN-1; /* where lsd of result will go */
- ul=lo->lsd; /* lsd of rhs */
- if (padding!=0) { /* unaligned */
- /* if the msd of lo is more than DECPMAX+2 digits to the right of */
- /* the original msd of hi then it can be reduced to a single */
- /* digit at the right place, as it stays clear of hi digits */
- /* [it must be DECPMAX+2 because during a subtraction the msd */
- /* could become 0 after a borrow from 1.000 to 0.9999...] */
- Int hilen=(Int)(hi->lsd-hi->msd+1); /* length of hi */
- Int lolen=(Int)(lo->lsd-lo->msd+1); /* and of lo */
- if (hilen+padding-lolen > DECPMAX+2) { /* can reduce lo to single */
- /* make sure it is virtually at least DECPMAX from hi->msd, at */
- /* least to right of hi->lsd (in case of destructive subtract), */
- /* and separated by at least two digits from either of those */
- /* (the tricky DOUBLE case is when hi is a 1 that will become a */
- /* 0.9999... by subtraction: */
- /* hi: 1 E+16 */
- /* lo: .................1000000000000000 E-16 */
- /* which for the addition pads to: */
- /* hi: 1000000000000000000 E-16 */
- /* lo: .................1000000000000000 E-16 */
- Int newexp=MINI(hi->exponent, hi->exponent+hilen-DECPMAX)-3;
- /* printf("FMA reduce: %ld\n", (LI)reduce); */
- lo->lsd=lo->msd; /* to single digit [maybe 0] */
- lo->exponent=newexp; /* new lowest exponent */
- padding=hi->exponent-lo->exponent; /* recalculate */
- ul=lo->lsd; /* .. and repoint */
- }
- /* padding is still > 0, but will fit in acc (less leading carry slot) */
- #if DECCHECK
- if (padding<=0) printf("FMA low padding: %ld\n", (LI)padding);
- if (hilen+padding+1>FMALEN)
- printf("FMA excess hilen+padding: %ld+%ld \n", (LI)hilen, (LI)padding);
- /* printf("FMA padding: %ld\n", (LI)padding); */
- #endif
- /* padding digits can now be set in the result; one or more of */
- /* these will come from lo; others will be zeros in the gap */
- for (; ul-3>=lo->msd && padding>3; padding-=4, ul-=4, ub-=4) {
- UBFROMUI(ub-3, UBTOUI(ul-3)); /* [cannot overlap] */
- }
- for (; ul>=lo->msd && padding>0; padding--, ul--, ub--) *ub=*ul;
- for (;padding>0; padding--, ub--) *ub=0; /* mind the gap */
- }
- /* addition now complete to the right of the rightmost digit of hi */
- uh=hi->lsd;
- /* dow do the add from hi->lsd to the left */
- /* [bytewise, because either operand can run out at any time] */
- /* carry was set up depending on ten's complement above */
- /* first assume both operands have some digits */
- for (;; ub--) {
- if (uh<hi->msd || ul<lo->msd) break;
- *ub=(uByte)(carry+(*uh--)+(*ul--));
- carry=0;
- if (*ub<10) continue;
- *ub-=10;
- carry=1;
- } /* both loop */
- if (ul<lo->msd) { /* to left of lo */
- for (;; ub--) {
- if (uh<hi->msd) break;
- *ub=(uByte)(carry+(*uh--)); /* [+0] */
- carry=0;
- if (*ub<10) continue;
- *ub-=10;
- carry=1;
- } /* hi loop */
- }
- else { /* to left of hi */
- for (;; ub--) {
- if (ul<lo->msd) break;
- *ub=(uByte)(carry+hipad+(*ul--));
- carry=0;
- if (*ub<10) continue;
- *ub-=10;
- carry=1;
- } /* lo loop */
- }
- /* addition complete -- now handle carry, borrow, etc. */
- /* use lo to set up the num (its exponent is already correct, and */
- /* sign usually is) */
- lo->msd=ub+1;
- lo->lsd=acc+FMALEN-1;
- /* decShowNum(lo, "lo"); */
- if (!diffsign) { /* same-sign addition */
- if (carry) { /* carry out */
- *ub=1; /* place the 1 .. */
- lo->msd--; /* .. and update */
- }
- } /* same sign */
- else { /* signs differed (subtraction) */
- if (!carry) { /* no carry out means hi<lo */
- /* borrowed -- take ten's complement of the right digits */
- lo->sign=hi->sign; /* sign is lhs sign */
- for (ul=lo->msd; ul<lo->lsd-3; ul+=4) UBFROMUI(ul, 0x09090909-UBTOUI(ul));
- for (; ul<=lo->lsd; ul++) *ul=(uByte)(0x09-*ul); /* [leaves ul at lsd+1] */
- /* complete the ten's complement by adding 1 [cannot overrun] */
- for (ul--; *ul==9; ul--) *ul=0;
- *ul+=1;
- } /* borrowed */
- else { /* carry out means hi>=lo */
- /* sign to use is lo->sign */
- /* all done except for the special IEEE 754 exact-zero-result */
- /* rule (see above); while testing for zero, strip leading */
- /* zeros (which will save decFinalize doing it) */
- for (; UBTOUI(lo->msd)==0 && lo->msd+3<lo->lsd;) lo->msd+=4;
- for (; *lo->msd==0 && lo->msd<lo->lsd;) lo->msd++;
- if (*lo->msd==0) { /* must be true zero (and diffsign) */
- lo->sign=0; /* assume + */
- if (set->round==DEC_ROUND_FLOOR) lo->sign=DECFLOAT_Sign;
- }
- /* [else was not zero, might still have leading zeros] */
- } /* subtraction gave positive result */
- } /* diffsign */
- #if DECCHECK
- /* assert no left underrun */
- if (lo->msd<acc) {
- printf("FMA underrun by %ld \n", (LI)(acc-lo->msd));
- }
- #endif
- return decFinalize(result, lo, set); /* round, check, and lay out */
- } /* decFloatFMA */
- /* ------------------------------------------------------------------ */
- /* decFloatFromInt -- initialise a decFloat from an Int */
- /* */
- /* result gets the converted Int */
- /* n is the Int to convert */
- /* returns result */
- /* */
- /* The result is Exact; no errors or exceptions are possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatFromInt32(decFloat *result, Int n) {
- uInt u=(uInt)n; /* copy as bits */
- uInt encode; /* work */
- DFWORD(result, 0)=ZEROWORD; /* always */
- #if QUAD
- DFWORD(result, 1)=0;
- DFWORD(result, 2)=0;
- #endif
- if (n<0) { /* handle -n with care */
- /* [This can be done without the test, but is then slightly slower] */
- u=(~u)+1;
- DFWORD(result, 0)|=DECFLOAT_Sign;
- }
- /* Since the maximum value of u now is 2**31, only the low word of */
- /* result is affected */
- encode=BIN2DPD[u%1000];
- u/=1000;
- encode|=BIN2DPD[u%1000]<<10;
- u/=1000;
- encode|=BIN2DPD[u%1000]<<20;
- u/=1000; /* now 0, 1, or 2 */
- encode|=u<<30;
- DFWORD(result, DECWORDS-1)=encode;
- return result;
- } /* decFloatFromInt32 */
- /* ------------------------------------------------------------------ */
- /* decFloatFromUInt -- initialise a decFloat from a uInt */
- /* */
- /* result gets the converted uInt */
- /* n is the uInt to convert */
- /* returns result */
- /* */
- /* The result is Exact; no errors or exceptions are possible. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatFromUInt32(decFloat *result, uInt u) {
- uInt encode; /* work */
- DFWORD(result, 0)=ZEROWORD; /* always */
- #if QUAD
- DFWORD(result, 1)=0;
- DFWORD(result, 2)=0;
- #endif
- encode=BIN2DPD[u%1000];
- u/=1000;
- encode|=BIN2DPD[u%1000]<<10;
- u/=1000;
- encode|=BIN2DPD[u%1000]<<20;
- u/=1000; /* now 0 -> 4 */
- encode|=u<<30;
- DFWORD(result, DECWORDS-1)=encode;
- DFWORD(result, DECWORDS-2)|=u>>2; /* rarely non-zero */
- return result;
- } /* decFloatFromUInt32 */
- /* ------------------------------------------------------------------ */
- /* decFloatInvert -- logical digitwise INVERT of a decFloat */
- /* */
- /* result gets the result of INVERTing df */
- /* df is the decFloat to invert */
- /* set is the context */
- /* returns result, which will be canonical with sign=0 */
- /* */
- /* The operand must be positive, finite with exponent q=0, and */
- /* comprise just zeros and ones; if not, Invalid operation results. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatInvert(decFloat *result, const decFloat *df,
- decContext *set) {
- uInt sourhi=DFWORD(df, 0); /* top word of dfs */
- if (!DFISUINT01(df) || !DFISCC01(df)) return decInvalid(result, set);
- /* the operand is a finite integer (q=0) */
- #if DOUBLE
- DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04009124);
- DFWORD(result, 1)=(~DFWORD(df, 1)) &0x49124491;
- #elif QUAD
- DFWORD(result, 0)=ZEROWORD|((~sourhi)&0x04000912);
- DFWORD(result, 1)=(~DFWORD(df, 1)) &0x44912449;
- DFWORD(result, 2)=(~DFWORD(df, 2)) &0x12449124;
- DFWORD(result, 3)=(~DFWORD(df, 3)) &0x49124491;
- #endif
- return result;
- } /* decFloatInvert */
- /* ------------------------------------------------------------------ */
- /* decFloatIs -- decFloat tests (IsSigned, etc.) */
- /* */
- /* df is the decFloat to test */
- /* returns 0 or 1 in a uInt */
- /* */
- /* Many of these could be macros, but having them as real functions */
- /* is a little cleaner (and they can be referred to here by the */
- /* generic names) */
- /* ------------------------------------------------------------------ */
- uInt decFloatIsCanonical(const decFloat *df) {
- if (DFISSPECIAL(df)) {
- if (DFISINF(df)) {
- if (DFWORD(df, 0)&ECONMASK) return 0; /* exponent continuation */
- if (!DFISCCZERO(df)) return 0; /* coefficient continuation */
- return 1;
- }
- /* is a NaN */
- if (DFWORD(df, 0)&ECONNANMASK) return 0; /* exponent continuation */
- if (DFISCCZERO(df)) return 1; /* coefficient continuation */
- /* drop through to check payload */
- }
- { /* declare block */
- #if DOUBLE
- uInt sourhi=DFWORD(df, 0);
- uInt sourlo=DFWORD(df, 1);
- if (CANONDPDOFF(sourhi, 8)
- && CANONDPDTWO(sourhi, sourlo, 30)
- && CANONDPDOFF(sourlo, 20)
- && CANONDPDOFF(sourlo, 10)
- && CANONDPDOFF(sourlo, 0)) return 1;
- #elif QUAD
- uInt sourhi=DFWORD(df, 0);
- uInt sourmh=DFWORD(df, 1);
- uInt sourml=DFWORD(df, 2);
- uInt sourlo=DFWORD(df, 3);
- if (CANONDPDOFF(sourhi, 4)
- && CANONDPDTWO(sourhi, sourmh, 26)
- && CANONDPDOFF(sourmh, 16)
- && CANONDPDOFF(sourmh, 6)
- && CANONDPDTWO(sourmh, sourml, 28)
- && CANONDPDOFF(sourml, 18)
- && CANONDPDOFF(sourml, 8)
- && CANONDPDTWO(sourml, sourlo, 30)
- && CANONDPDOFF(sourlo, 20)
- && CANONDPDOFF(sourlo, 10)
- && CANONDPDOFF(sourlo, 0)) return 1;
- #endif
- } /* block */
- return 0; /* a declet is non-canonical */
- }
- uInt decFloatIsFinite(const decFloat *df) {
- return !DFISSPECIAL(df);
- }
- uInt decFloatIsInfinite(const decFloat *df) {
- return DFISINF(df);
- }
- uInt decFloatIsInteger(const decFloat *df) {
- return DFISINT(df);
- }
- uInt decFloatIsNaN(const decFloat *df) {
- return DFISNAN(df);
- }
- uInt decFloatIsNormal(const decFloat *df) {
- Int exp; /* exponent */
- if (DFISSPECIAL(df)) return 0;
- if (DFISZERO(df)) return 0;
- /* is finite and non-zero */
- exp=GETEXPUN(df) /* get unbiased exponent .. */
- +decFloatDigits(df)-1; /* .. and make adjusted exponent */
- return (exp>=DECEMIN); /* < DECEMIN is subnormal */
- }
- uInt decFloatIsSignaling(const decFloat *df) {
- return DFISSNAN(df);
- }
- uInt decFloatIsSignalling(const decFloat *df) {
- return DFISSNAN(df);
- }
- uInt decFloatIsSigned(const decFloat *df) {
- return DFISSIGNED(df);
- }
- uInt decFloatIsSubnormal(const decFloat *df) {
- if (DFISSPECIAL(df)) return 0;
- /* is finite */
- if (decFloatIsNormal(df)) return 0;
- /* it is <Nmin, but could be zero */
- if (DFISZERO(df)) return 0;
- return 1; /* is subnormal */
- }
- uInt decFloatIsZero(const decFloat *df) {
- return DFISZERO(df);
- } /* decFloatIs... */
- /* ------------------------------------------------------------------ */
- /* decFloatLogB -- return adjusted exponent, by 754 rules */
- /* */
- /* result gets the adjusted exponent as an integer, or a NaN etc. */
- /* df is the decFloat to be examined */
- /* set is the context */
- /* returns result */
- /* */
- /* Notable cases: */
- /* A<0 -> Use |A| */
- /* A=0 -> -Infinity (Division by zero) */
- /* A=Infinite -> +Infinity (Exact) */
- /* A=1 exactly -> 0 (Exact) */
- /* NaNs are propagated as usual */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatLogB(decFloat *result, const decFloat *df,
- decContext *set) {
- Int ae; /* adjusted exponent */
- if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
- if (DFISINF(df)) {
- DFWORD(result, 0)=0; /* need +ve */
- return decInfinity(result, result); /* canonical +Infinity */
- }
- if (DFISZERO(df)) {
- set->status|=DEC_Division_by_zero; /* as per 754 */
- DFWORD(result, 0)=DECFLOAT_Sign; /* make negative */
- return decInfinity(result, result); /* canonical -Infinity */
- }
- ae=GETEXPUN(df) /* get unbiased exponent .. */
- +decFloatDigits(df)-1; /* .. and make adjusted exponent */
- /* ae has limited range (3 digits for DOUBLE and 4 for QUAD), so */
- /* it is worth using a special case of decFloatFromInt32 */
- DFWORD(result, 0)=ZEROWORD; /* always */
- if (ae<0) {
- DFWORD(result, 0)|=DECFLOAT_Sign; /* -0 so far */
- ae=-ae;
- }
- #if DOUBLE
- DFWORD(result, 1)=BIN2DPD[ae]; /* a single declet */
- #elif QUAD
- DFWORD(result, 1)=0;
- DFWORD(result, 2)=0;
- DFWORD(result, 3)=(ae/1000)<<10; /* is <10, so need no DPD encode */
- DFWORD(result, 3)|=BIN2DPD[ae%1000];
- #endif
- return result;
- } /* decFloatLogB */
- /* ------------------------------------------------------------------ */
- /* decFloatMax -- return maxnum of two operands */
- /* */
- /* result gets the chosen decFloat */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* If just one operand is a quiet NaN it is ignored. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatMax(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int comp;
- if (DFISNAN(dfl)) {
- /* sNaN or both NaNs leads to normal NaN processing */
- if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
- return decCanonical(result, dfr); /* RHS is numeric */
- }
- if (DFISNAN(dfr)) {
- /* sNaN leads to normal NaN processing (both NaN handled above) */
- if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- return decCanonical(result, dfl); /* LHS is numeric */
- }
- /* Both operands are numeric; numeric comparison needed -- use */
- /* total order for a well-defined choice (and +0 > -0) */
- comp=decNumCompare(dfl, dfr, 1);
- if (comp>=0) return decCanonical(result, dfl);
- return decCanonical(result, dfr);
- } /* decFloatMax */
- /* ------------------------------------------------------------------ */
- /* decFloatMaxMag -- return maxnummag of two operands */
- /* */
- /* result gets the chosen decFloat */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* Returns according to the magnitude comparisons if both numeric and */
- /* unequal, otherwise returns maxnum */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatMaxMag(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int comp;
- decFloat absl, absr;
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMax(result, dfl, dfr, set);
- decFloatCopyAbs(&absl, dfl);
- decFloatCopyAbs(&absr, dfr);
- comp=decNumCompare(&absl, &absr, 0);
- if (comp>0) return decCanonical(result, dfl);
- if (comp<0) return decCanonical(result, dfr);
- return decFloatMax(result, dfl, dfr, set);
- } /* decFloatMaxMag */
- /* ------------------------------------------------------------------ */
- /* decFloatMin -- return minnum of two operands */
- /* */
- /* result gets the chosen decFloat */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* If just one operand is a quiet NaN it is ignored. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatMin(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int comp;
- if (DFISNAN(dfl)) {
- /* sNaN or both NaNs leads to normal NaN processing */
- if (DFISNAN(dfr) || DFISSNAN(dfl)) return decNaNs(result, dfl, dfr, set);
- return decCanonical(result, dfr); /* RHS is numeric */
- }
- if (DFISNAN(dfr)) {
- /* sNaN leads to normal NaN processing (both NaN handled above) */
- if (DFISSNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- return decCanonical(result, dfl); /* LHS is numeric */
- }
- /* Both operands are numeric; numeric comparison needed -- use */
- /* total order for a well-defined choice (and +0 > -0) */
- comp=decNumCompare(dfl, dfr, 1);
- if (comp<=0) return decCanonical(result, dfl);
- return decCanonical(result, dfr);
- } /* decFloatMin */
- /* ------------------------------------------------------------------ */
- /* decFloatMinMag -- return minnummag of two operands */
- /* */
- /* result gets the chosen decFloat */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* Returns according to the magnitude comparisons if both numeric and */
- /* unequal, otherwise returns minnum */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatMinMag(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int comp;
- decFloat absl, absr;
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decFloatMin(result, dfl, dfr, set);
- decFloatCopyAbs(&absl, dfl);
- decFloatCopyAbs(&absr, dfr);
- comp=decNumCompare(&absl, &absr, 0);
- if (comp<0) return decCanonical(result, dfl);
- if (comp>0) return decCanonical(result, dfr);
- return decFloatMin(result, dfl, dfr, set);
- } /* decFloatMinMag */
- /* ------------------------------------------------------------------ */
- /* decFloatMinus -- negate value, heeding NaNs, etc. */
- /* */
- /* result gets the canonicalized 0-df */
- /* df is the decFloat to minus */
- /* set is the context */
- /* returns result */
- /* */
- /* This has the same effect as 0-df where the exponent of the zero is */
- /* the same as that of df (if df is finite). */
- /* The effect is also the same as decFloatCopyNegate except that NaNs */
- /* are handled normally (the sign of a NaN is not affected, and an */
- /* sNaN will signal), the result is canonical, and zero gets sign 0. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatMinus(decFloat *result, const decFloat *df,
- decContext *set) {
- if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
- decCanonical(result, df); /* copy and check */
- if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
- else DFBYTE(result, 0)^=0x80; /* flip sign bit */
- return result;
- } /* decFloatMinus */
- /* ------------------------------------------------------------------ */
- /* decFloatMultiply -- multiply two decFloats */
- /* */
- /* result gets the result of multiplying dfl and dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatMultiply(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- bcdnum num; /* for final conversion */
- uByte bcdacc[DECPMAX9*18+1]; /* for coefficent in BCD */
- if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) { /* either is special? */
- /* NaNs are handled as usual */
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- /* infinity times zero is bad */
- if (DFISINF(dfl) && DFISZERO(dfr)) return decInvalid(result, set);
- if (DFISINF(dfr) && DFISZERO(dfl)) return decInvalid(result, set);
- /* both infinite; return canonical infinity with computed sign */
- DFWORD(result, 0)=DFWORD(dfl, 0)^DFWORD(dfr, 0); /* compute sign */
- return decInfinity(result, result);
- }
- /* Here when both operands are finite */
- decFiniteMultiply(&num, bcdacc, dfl, dfr);
- return decFinalize(result, &num, set); /* round, check, and lay out */
- } /* decFloatMultiply */
- /* ------------------------------------------------------------------ */
- /* decFloatNextMinus -- next towards -Infinity */
- /* */
- /* result gets the next lesser decFloat */
- /* dfl is the decFloat to start with */
- /* set is the context */
- /* returns result */
- /* */
- /* This is 754 nextdown; Invalid is the only status possible (from */
- /* an sNaN). */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatNextMinus(decFloat *result, const decFloat *dfl,
- decContext *set) {
- decFloat delta; /* tiny increment */
- uInt savestat; /* saves status */
- enum rounding saveround; /* .. and mode */
- /* +Infinity is the special case */
- if (DFISINF(dfl) && !DFISSIGNED(dfl)) {
- DFSETNMAX(result);
- return result; /* [no status to set] */
- }
- /* other cases are effected by sutracting a tiny delta -- this */
- /* should be done in a wider format as the delta is unrepresentable */
- /* here (but can be done with normal add if the sign of zero is */
- /* treated carefully, because no Inexactitude is interesting); */
- /* rounding to -Infinity then pushes the result to next below */
- decFloatZero(&delta); /* set up tiny delta */
- DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
- DFWORD(&delta, 0)=DECFLOAT_Sign; /* Sign=1 + biased exponent=0 */
- /* set up for the directional round */
- saveround=set->round; /* save mode */
- set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
- savestat=set->status; /* save status */
- decFloatAdd(result, dfl, &delta, set);
- /* Add rules mess up the sign when going from +Ntiny to 0 */
- if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
- set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
- set->status|=savestat; /* restore pending flags */
- set->round=saveround; /* .. and mode */
- return result;
- } /* decFloatNextMinus */
- /* ------------------------------------------------------------------ */
- /* decFloatNextPlus -- next towards +Infinity */
- /* */
- /* result gets the next larger decFloat */
- /* dfl is the decFloat to start with */
- /* set is the context */
- /* returns result */
- /* */
- /* This is 754 nextup; Invalid is the only status possible (from */
- /* an sNaN). */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatNextPlus(decFloat *result, const decFloat *dfl,
- decContext *set) {
- uInt savestat; /* saves status */
- enum rounding saveround; /* .. and mode */
- decFloat delta; /* tiny increment */
- /* -Infinity is the special case */
- if (DFISINF(dfl) && DFISSIGNED(dfl)) {
- DFSETNMAX(result);
- DFWORD(result, 0)|=DECFLOAT_Sign; /* make negative */
- return result; /* [no status to set] */
- }
- /* other cases are effected by sutracting a tiny delta -- this */
- /* should be done in a wider format as the delta is unrepresentable */
- /* here (but can be done with normal add if the sign of zero is */
- /* treated carefully, because no Inexactitude is interesting); */
- /* rounding to +Infinity then pushes the result to next above */
- decFloatZero(&delta); /* set up tiny delta */
- DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
- DFWORD(&delta, 0)=0; /* Sign=0 + biased exponent=0 */
- /* set up for the directional round */
- saveround=set->round; /* save mode */
- set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
- savestat=set->status; /* save status */
- decFloatAdd(result, dfl, &delta, set);
- /* Add rules mess up the sign when going from -Ntiny to -0 */
- if (DFISZERO(result)) DFWORD(result, 0)^=DECFLOAT_Sign; /* correct */
- set->status&=DEC_Invalid_operation; /* preserve only sNaN status */
- set->status|=savestat; /* restore pending flags */
- set->round=saveround; /* .. and mode */
- return result;
- } /* decFloatNextPlus */
- /* ------------------------------------------------------------------ */
- /* decFloatNextToward -- next towards a decFloat */
- /* */
- /* result gets the next decFloat */
- /* dfl is the decFloat to start with */
- /* dfr is the decFloat to move toward */
- /* set is the context */
- /* returns result */
- /* */
- /* This is 754-1985 nextafter, as modified during revision (dropped */
- /* from 754-2008); status may be set unless the result is a normal */
- /* number. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatNextToward(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- decFloat delta; /* tiny increment or decrement */
- decFloat pointone; /* 1e-1 */
- uInt savestat; /* saves status */
- enum rounding saveround; /* .. and mode */
- uInt deltatop; /* top word for delta */
- Int comp; /* work */
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- /* Both are numeric, so Invalid no longer a possibility */
- comp=decNumCompare(dfl, dfr, 0);
- if (comp==0) return decFloatCopySign(result, dfl, dfr); /* equal */
- /* unequal; do NextPlus or NextMinus but with different status rules */
- if (comp<0) { /* lhs<rhs, do NextPlus, see above for commentary */
- if (DFISINF(dfl) && DFISSIGNED(dfl)) { /* -Infinity special case */
- DFSETNMAX(result);
- DFWORD(result, 0)|=DECFLOAT_Sign;
- return result;
- }
- saveround=set->round; /* save mode */
- set->round=DEC_ROUND_CEILING; /* .. round towards +Infinity */
- deltatop=0; /* positive delta */
- }
- else { /* lhs>rhs, do NextMinus, see above for commentary */
- if (DFISINF(dfl) && !DFISSIGNED(dfl)) { /* +Infinity special case */
- DFSETNMAX(result);
- return result;
- }
- saveround=set->round; /* save mode */
- set->round=DEC_ROUND_FLOOR; /* .. round towards -Infinity */
- deltatop=DECFLOAT_Sign; /* negative delta */
- }
- savestat=set->status; /* save status */
- /* Here, Inexact is needed where appropriate (and hence Underflow, */
- /* etc.). Therefore the tiny delta which is otherwise */
- /* unrepresentable (see NextPlus and NextMinus) is constructed */
- /* using the multiplication of FMA. */
- decFloatZero(&delta); /* set up tiny delta */
- DFWORD(&delta, DECWORDS-1)=1; /* coefficient=1 */
- DFWORD(&delta, 0)=deltatop; /* Sign + biased exponent=0 */
- decFloatFromString(&pointone, "1E-1", set); /* set up multiplier */
- decFloatFMA(result, &delta, &pointone, dfl, set);
- /* [Delta is truly tiny, so no need to correct sign of zero] */
- /* use new status unless the result is normal */
- if (decFloatIsNormal(result)) set->status=savestat; /* else goes forward */
- set->round=saveround; /* restore mode */
- return result;
- } /* decFloatNextToward */
- /* ------------------------------------------------------------------ */
- /* decFloatOr -- logical digitwise OR of two decFloats */
- /* */
- /* result gets the result of ORing dfl and dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result, which will be canonical with sign=0 */
- /* */
- /* The operands must be positive, finite with exponent q=0, and */
- /* comprise just zeros and ones; if not, Invalid operation results. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatOr(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
- || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
- /* the operands are positive finite integers (q=0) with just 0s and 1s */
- #if DOUBLE
- DFWORD(result, 0)=ZEROWORD
- |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04009124);
- DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x49124491;
- #elif QUAD
- DFWORD(result, 0)=ZEROWORD
- |((DFWORD(dfl, 0) | DFWORD(dfr, 0))&0x04000912);
- DFWORD(result, 1)=(DFWORD(dfl, 1) | DFWORD(dfr, 1))&0x44912449;
- DFWORD(result, 2)=(DFWORD(dfl, 2) | DFWORD(dfr, 2))&0x12449124;
- DFWORD(result, 3)=(DFWORD(dfl, 3) | DFWORD(dfr, 3))&0x49124491;
- #endif
- return result;
- } /* decFloatOr */
- /* ------------------------------------------------------------------ */
- /* decFloatPlus -- add value to 0, heeding NaNs, etc. */
- /* */
- /* result gets the canonicalized 0+df */
- /* df is the decFloat to plus */
- /* set is the context */
- /* returns result */
- /* */
- /* This has the same effect as 0+df where the exponent of the zero is */
- /* the same as that of df (if df is finite). */
- /* The effect is also the same as decFloatCopy except that NaNs */
- /* are handled normally (the sign of a NaN is not affected, and an */
- /* sNaN will signal), the result is canonical, and zero gets sign 0. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatPlus(decFloat *result, const decFloat *df,
- decContext *set) {
- if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
- decCanonical(result, df); /* copy and check */
- if (DFISZERO(df)) DFBYTE(result, 0)&=~0x80; /* turn off sign bit */
- return result;
- } /* decFloatPlus */
- /* ------------------------------------------------------------------ */
- /* decFloatQuantize -- quantize a decFloat */
- /* */
- /* result gets the result of quantizing dfl to match dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs), which sets the exponent */
- /* set is the context */
- /* returns result */
- /* */
- /* Unless there is an error or the result is infinite, the exponent */
- /* of result is guaranteed to be the same as that of dfr. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatQuantize(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int explb, exprb; /* left and right biased exponents */
- uByte *ulsd; /* local LSD pointer */
- uByte *ub, *uc; /* work */
- Int drop; /* .. */
- uInt dpd; /* .. */
- uInt encode; /* encoding accumulator */
- uInt sourhil, sourhir; /* top words from source decFloats */
- uInt uiwork; /* for macros */
- #if QUAD
- uShort uswork; /* .. */
- #endif
- /* the following buffer holds the coefficient for manipulation */
- uByte buf[4+DECPMAX*3+2*QUAD]; /* + space for zeros to left or right */
- #if DECTRACE
- bcdnum num; /* for trace displays */
- #endif
- /* Start decoding the arguments */
- sourhil=DFWORD(dfl, 0); /* LHS top word */
- explb=DECCOMBEXP[sourhil>>26]; /* get exponent high bits (in place) */
- sourhir=DFWORD(dfr, 0); /* RHS top word */
- exprb=DECCOMBEXP[sourhir>>26];
- if (EXPISSPECIAL(explb | exprb)) { /* either is special? */
- /* NaNs are handled as usual */
- if (DFISNAN(dfl) || DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- /* one infinity but not both is bad */
- if (DFISINF(dfl)!=DFISINF(dfr)) return decInvalid(result, set);
- /* both infinite; return canonical infinity with sign of LHS */
- return decInfinity(result, dfl);
- }
- /* Here when both arguments are finite */
- /* complete extraction of the exponents [no need to unbias] */
- explb+=GETECON(dfl); /* + continuation */
- exprb+=GETECON(dfr); /* .. */
- /* calculate the number of digits to drop from the coefficient */
- drop=exprb-explb; /* 0 if nothing to do */
- if (drop==0) return decCanonical(result, dfl); /* return canonical */
- /* the coefficient is needed; lay it out into buf, offset so zeros */
- /* can be added before or after as needed -- an extra heading is */
- /* added so can safely pad Quad DECPMAX-1 zeros to the left by */
- /* fours */
- #define BUFOFF (buf+4+DECPMAX)
- GETCOEFF(dfl, BUFOFF); /* decode from decFloat */
- /* [now the msd is at BUFOFF and the lsd is at BUFOFF+DECPMAX-1] */
- #if DECTRACE
- num.msd=BUFOFF;
- num.lsd=BUFOFF+DECPMAX-1;
- num.exponent=explb-DECBIAS;
- num.sign=sourhil & DECFLOAT_Sign;
- decShowNum(&num, "dfl");
- #endif
- if (drop>0) { /* [most common case] */
- /* (this code is very similar to that in decFloatFinalize, but */
- /* has many differences so is duplicated here -- so any changes */
- /* may need to be made there, too) */
- uByte *roundat; /* -> re-round digit */
- uByte reround; /* reround value */
- /* printf("Rounding; drop=%ld\n", (LI)drop); */
- /* there is at least one zero needed to the left, in all but one */
- /* exceptional (all-nines) case, so place four zeros now; this is */
- /* needed almost always and makes rounding all-nines by fours safe */
- UBFROMUI(BUFOFF-4, 0);
- /* Three cases here: */
- /* 1. new LSD is in coefficient (almost always) */
- /* 2. new LSD is digit to left of coefficient (so MSD is */
- /* round-for-reround digit) */
- /* 3. new LSD is to left of case 2 (whole coefficient is sticky) */
- /* Note that leading zeros can safely be treated as useful digits */
- /* [duplicate check-stickies code to save a test] */
- /* [by-digit check for stickies as runs of zeros are rare] */
- if (drop<DECPMAX) { /* NB lengths not addresses */
- roundat=BUFOFF+DECPMAX-drop;
- reround=*roundat;
- for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
- if (*ub!=0) { /* non-zero to be discarded */
- reround=DECSTICKYTAB[reround]; /* apply sticky bit */
- break; /* [remainder don't-care] */
- }
- } /* check stickies */
- ulsd=roundat-1; /* set LSD */
- }
- else { /* edge case */
- if (drop==DECPMAX) {
- roundat=BUFOFF;
- reround=*roundat;
- }
- else {
- roundat=BUFOFF-1;
- reround=0;
- }
- for (ub=roundat+1; ub<BUFOFF+DECPMAX; ub++) {
- if (*ub!=0) { /* non-zero to be discarded */
- reround=DECSTICKYTAB[reround]; /* apply sticky bit */
- break; /* [remainder don't-care] */
- }
- } /* check stickies */
- *BUFOFF=0; /* make a coefficient of 0 */
- ulsd=BUFOFF; /* .. at the MSD place */
- }
- if (reround!=0) { /* discarding non-zero */
- uInt bump=0;
- set->status|=DEC_Inexact;
- /* next decide whether to increment the coefficient */
- if (set->round==DEC_ROUND_HALF_EVEN) { /* fastpath slowest case */
- if (reround>5) bump=1; /* >0.5 goes up */
- else if (reround==5) /* exactly 0.5000 .. */
- bump=*ulsd & 0x01; /* .. up iff [new] lsd is odd */
- } /* r-h-e */
- else switch (set->round) {
- case DEC_ROUND_DOWN: {
- /* no change */
- break;} /* r-d */
- case DEC_ROUND_HALF_DOWN: {
- if (reround>5) bump=1;
- break;} /* r-h-d */
- case DEC_ROUND_HALF_UP: {
- if (reround>=5) bump=1;
- break;} /* r-h-u */
- case DEC_ROUND_UP: {
- if (reround>0) bump=1;
- break;} /* r-u */
- case DEC_ROUND_CEILING: {
- /* same as _UP for positive numbers, and as _DOWN for negatives */
- if (!(sourhil&DECFLOAT_Sign) && reround>0) bump=1;
- break;} /* r-c */
- case DEC_ROUND_FLOOR: {
- /* same as _UP for negative numbers, and as _DOWN for positive */
- /* [negative reround cannot occur on 0] */
- if (sourhil&DECFLOAT_Sign && reround>0) bump=1;
- break;} /* r-f */
- case DEC_ROUND_05UP: {
- if (reround>0) { /* anything out there is 'sticky' */
- /* bump iff lsd=0 or 5; this cannot carry so it could be */
- /* effected immediately with no bump -- but the code */
- /* is clearer if this is done the same way as the others */
- if (*ulsd==0 || *ulsd==5) bump=1;
- }
- break;} /* r-r */
- default: { /* e.g., DEC_ROUND_MAX */
- set->status|=DEC_Invalid_context;
- #if DECCHECK
- printf("Unknown rounding mode: %ld\n", (LI)set->round);
- #endif
- break;}
- } /* switch (not r-h-e) */
- /* printf("ReRound: %ld bump: %ld\n", (LI)reround, (LI)bump); */
- if (bump!=0) { /* need increment */
- /* increment the coefficient; this could give 1000... (after */
- /* the all nines case) */
- ub=ulsd;
- for (; UBTOUI(ub-3)==0x09090909; ub-=4) UBFROMUI(ub-3, 0);
- /* now at most 3 digits left to non-9 (usually just the one) */
- for (; *ub==9; ub--) *ub=0;
- *ub+=1;
- /* [the all-nines case will have carried one digit to the */
- /* left of the original MSD -- just where it is needed] */
- } /* bump needed */
- } /* inexact rounding */
- /* now clear zeros to the left so exactly DECPMAX digits will be */
- /* available in the coefficent -- the first word to the left was */
- /* cleared earlier for safe carry; now add any more needed */
- if (drop>4) {
- UBFROMUI(BUFOFF-8, 0); /* must be at least 5 */
- for (uc=BUFOFF-12; uc>ulsd-DECPMAX-3; uc-=4) UBFROMUI(uc, 0);
- }
- } /* need round (drop>0) */
- else { /* drop<0; padding with -drop digits is needed */
- /* This is the case where an error can occur if the padded */
- /* coefficient will not fit; checking for this can be done in the */
- /* same loop as padding for zeros if the no-hope and zero cases */
- /* are checked first */
- if (-drop>DECPMAX-1) { /* cannot fit unless 0 */
- if (!ISCOEFFZERO(BUFOFF)) return decInvalid(result, set);
- /* a zero can have any exponent; just drop through and use it */
- ulsd=BUFOFF+DECPMAX-1;
- }
- else { /* padding will fit (but may still be too long) */
- /* final-word mask depends on endianess */
- #if DECLITEND
- static const uInt dmask[]={0, 0x000000ff, 0x0000ffff, 0x00ffffff};
- #else
- static const uInt dmask[]={0, 0xff000000, 0xffff0000, 0xffffff00};
- #endif
- /* note that here zeros to the right are added by fours, so in */
- /* the Quad case this could write 36 zeros if the coefficient has */
- /* fewer than three significant digits (hence the +2*QUAD for buf) */
- for (uc=BUFOFF+DECPMAX;; uc+=4) {
- UBFROMUI(uc, 0);
- if (UBTOUI(uc-DECPMAX)!=0) { /* could be bad */
- /* if all four digits should be zero, definitely bad */
- if (uc<=BUFOFF+DECPMAX+(-drop)-4)
- return decInvalid(result, set);
- /* must be a 1- to 3-digit sequence; check more carefully */
- if ((UBTOUI(uc-DECPMAX)&dmask[(-drop)%4])!=0)
- return decInvalid(result, set);
- break; /* no need for loop end test */
- }
- if (uc>=BUFOFF+DECPMAX+(-drop)-4) break; /* done */
- }
- ulsd=BUFOFF+DECPMAX+(-drop)-1;
- } /* pad and check leading zeros */
- } /* drop<0 */
- #if DECTRACE
- num.msd=ulsd-DECPMAX+1;
- num.lsd=ulsd;
- num.exponent=explb-DECBIAS;
- num.sign=sourhil & DECFLOAT_Sign;
- decShowNum(&num, "res");
- #endif
- /*------------------------------------------------------------------*/
- /* At this point the result is DECPMAX digits, ending at ulsd, so */
- /* fits the encoding exactly; there is no possibility of error */
- /*------------------------------------------------------------------*/
- encode=((exprb>>DECECONL)<<4) + *(ulsd-DECPMAX+1); /* make index */
- encode=DECCOMBFROM[encode]; /* indexed by (0-2)*16+msd */
- /* the exponent continuation can be extracted from the original RHS */
- encode|=sourhir & ECONMASK;
- encode|=sourhil&DECFLOAT_Sign; /* add the sign from LHS */
- /* finally encode the coefficient */
- /* private macro to encode a declet; this version can be used */
- /* because all coefficient digits exist */
- #define getDPD3q(dpd, n) ub=ulsd-(3*(n))-2; \
- dpd=BCD2DPD[(*ub*256)+(*(ub+1)*16)+*(ub+2)];
- #if DOUBLE
- getDPD3q(dpd, 4); encode|=dpd<<8;
- getDPD3q(dpd, 3); encode|=dpd>>2;
- DFWORD(result, 0)=encode;
- encode=dpd<<30;
- getDPD3q(dpd, 2); encode|=dpd<<20;
- getDPD3q(dpd, 1); encode|=dpd<<10;
- getDPD3q(dpd, 0); encode|=dpd;
- DFWORD(result, 1)=encode;
- #elif QUAD
- getDPD3q(dpd,10); encode|=dpd<<4;
- getDPD3q(dpd, 9); encode|=dpd>>6;
- DFWORD(result, 0)=encode;
- encode=dpd<<26;
- getDPD3q(dpd, 8); encode|=dpd<<16;
- getDPD3q(dpd, 7); encode|=dpd<<6;
- getDPD3q(dpd, 6); encode|=dpd>>4;
- DFWORD(result, 1)=encode;
- encode=dpd<<28;
- getDPD3q(dpd, 5); encode|=dpd<<18;
- getDPD3q(dpd, 4); encode|=dpd<<8;
- getDPD3q(dpd, 3); encode|=dpd>>2;
- DFWORD(result, 2)=encode;
- encode=dpd<<30;
- getDPD3q(dpd, 2); encode|=dpd<<20;
- getDPD3q(dpd, 1); encode|=dpd<<10;
- getDPD3q(dpd, 0); encode|=dpd;
- DFWORD(result, 3)=encode;
- #endif
- return result;
- } /* decFloatQuantize */
- /* ------------------------------------------------------------------ */
- /* decFloatReduce -- reduce finite coefficient to minimum length */
- /* */
- /* result gets the reduced decFloat */
- /* df is the source decFloat */
- /* set is the context */
- /* returns result, which will be canonical */
- /* */
- /* This removes all possible trailing zeros from the coefficient; */
- /* some may remain when the number is very close to Nmax. */
- /* Special values are unchanged and no status is set unless df=sNaN. */
- /* Reduced zero has an exponent q=0. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatReduce(decFloat *result, const decFloat *df,
- decContext *set) {
- bcdnum num; /* work */
- uByte buf[DECPMAX], *ub; /* coefficient and pointer */
- if (df!=result) *result=*df; /* copy, if needed */
- if (DFISNAN(df)) return decNaNs(result, df, NULL, set); /* sNaN */
- /* zeros and infinites propagate too */
- if (DFISINF(df)) return decInfinity(result, df); /* canonical */
- if (DFISZERO(df)) {
- uInt sign=DFWORD(df, 0)&DECFLOAT_Sign;
- decFloatZero(result);
- DFWORD(result, 0)|=sign;
- return result; /* exponent dropped, sign OK */
- }
- /* non-zero finite */
- GETCOEFF(df, buf);
- ub=buf+DECPMAX-1; /* -> lsd */
- if (*ub) return result; /* no trailing zeros */
- for (ub--; *ub==0;) ub--; /* terminates because non-zero */
- /* *ub is the first non-zero from the right */
- num.sign=DFWORD(df, 0)&DECFLOAT_Sign; /* set up number... */
- num.exponent=GETEXPUN(df)+(Int)(buf+DECPMAX-1-ub); /* adjusted exponent */
- num.msd=buf;
- num.lsd=ub;
- return decFinalize(result, &num, set);
- } /* decFloatReduce */
- /* ------------------------------------------------------------------ */
- /* decFloatRemainder -- integer divide and return remainder */
- /* */
- /* result gets the remainder of dividing dfl by dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatRemainder(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- return decDivide(result, dfl, dfr, set, REMAINDER);
- } /* decFloatRemainder */
- /* ------------------------------------------------------------------ */
- /* decFloatRemainderNear -- integer divide to nearest and remainder */
- /* */
- /* result gets the remainder of dividing dfl by dfr: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* This is the IEEE remainder, where the nearest integer is used. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatRemainderNear(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- return decDivide(result, dfl, dfr, set, REMNEAR);
- } /* decFloatRemainderNear */
- /* ------------------------------------------------------------------ */
- /* decFloatRotate -- rotate the coefficient of a decFloat left/right */
- /* */
- /* result gets the result of rotating dfl */
- /* dfl is the source decFloat to rotate */
- /* dfr is the count of digits to rotate, an integer (with q=0) */
- /* set is the context */
- /* returns result */
- /* */
- /* The digits of the coefficient of dfl are rotated to the left (if */
- /* dfr is positive) or to the right (if dfr is negative) without */
- /* adjusting the exponent or the sign of dfl. */
- /* */
- /* dfr must be in the range -DECPMAX through +DECPMAX. */
- /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
- /* dfr must be valid). No status is set unless dfr is invalid or an */
- /* operand is an sNaN. The result is canonical. */
- /* ------------------------------------------------------------------ */
- #define PHALF (ROUNDUP(DECPMAX/2, 4)) /* half length, rounded up */
- decFloat * decFloatRotate(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int rotate; /* dfr as an Int */
- uByte buf[DECPMAX+PHALF]; /* coefficient + half */
- uInt digits, savestat; /* work */
- bcdnum num; /* .. */
- uByte *ub; /* .. */
- if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- if (!DFISINT(dfr)) return decInvalid(result, set);
- digits=decFloatDigits(dfr); /* calculate digits */
- if (digits>2) return decInvalid(result, set); /* definitely out of range */
- rotate=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
- if (rotate>DECPMAX) return decInvalid(result, set); /* too big */
- /* [from here on no error or status change is possible] */
- if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
- /* handle no-rotate cases */
- if (rotate==0 || rotate==DECPMAX) return decCanonical(result, dfl);
- /* a real rotate is needed: 0 < rotate < DECPMAX */
- /* reduce the rotation to no more than half to reduce copying later */
- /* (for QUAD in fact half + 2 digits) */
- if (DFISSIGNED(dfr)) rotate=-rotate;
- if (abs(rotate)>PHALF) {
- if (rotate<0) rotate=DECPMAX+rotate;
- else rotate=rotate-DECPMAX;
- }
- /* now lay out the coefficient, leaving room to the right or the */
- /* left depending on the direction of rotation */
- ub=buf;
- if (rotate<0) ub+=PHALF; /* rotate right, so space to left */
- GETCOEFF(dfl, ub);
- /* copy half the digits to left or right, and set num.msd */
- if (rotate<0) {
- memcpy(buf, buf+DECPMAX, PHALF);
- num.msd=buf+PHALF+rotate;
- }
- else {
- memcpy(buf+DECPMAX, buf, PHALF);
- num.msd=buf+rotate;
- }
- /* fill in rest of num */
- num.lsd=num.msd+DECPMAX-1;
- num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
- num.exponent=GETEXPUN(dfl);
- savestat=set->status; /* record */
- decFinalize(result, &num, set);
- set->status=savestat; /* restore */
- return result;
- } /* decFloatRotate */
- /* ------------------------------------------------------------------ */
- /* decFloatSameQuantum -- test decFloats for same quantum */
- /* */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* returns 1 if the operands have the same quantum, 0 otherwise */
- /* */
- /* No error is possible and no status results. */
- /* ------------------------------------------------------------------ */
- uInt decFloatSameQuantum(const decFloat *dfl, const decFloat *dfr) {
- if (DFISSPECIAL(dfl) || DFISSPECIAL(dfr)) {
- if (DFISNAN(dfl) && DFISNAN(dfr)) return 1;
- if (DFISINF(dfl) && DFISINF(dfr)) return 1;
- return 0; /* any other special mixture gives false */
- }
- if (GETEXP(dfl)==GETEXP(dfr)) return 1; /* biased exponents match */
- return 0;
- } /* decFloatSameQuantum */
- /* ------------------------------------------------------------------ */
- /* decFloatScaleB -- multiply by a power of 10, as per 754 */
- /* */
- /* result gets the result of the operation */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs), am integer (with q=0) */
- /* set is the context */
- /* returns result */
- /* */
- /* This computes result=dfl x 10**dfr where dfr is an integer in the */
- /* range +/-2*(emax+pmax), typically resulting from LogB. */
- /* Underflow and Overflow (with Inexact) may occur. NaNs propagate */
- /* as usual. */
- /* ------------------------------------------------------------------ */
- #define SCALEBMAX 2*(DECEMAX+DECPMAX) /* D=800, Q=12356 */
- decFloat * decFloatScaleB(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- uInt digits; /* work */
- Int expr; /* dfr as an Int */
- if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- if (!DFISINT(dfr)) return decInvalid(result, set);
- digits=decFloatDigits(dfr); /* calculate digits */
- #if DOUBLE
- if (digits>3) return decInvalid(result, set); /* definitely out of range */
- expr=DPD2BIN[DFWORD(dfr, 1)&0x3ff]; /* must be in bottom declet */
- #elif QUAD
- if (digits>5) return decInvalid(result, set); /* definitely out of range */
- expr=DPD2BIN[DFWORD(dfr, 3)&0x3ff] /* in bottom 2 declets .. */
- +DPD2BIN[(DFWORD(dfr, 3)>>10)&0x3ff]*1000; /* .. */
- #endif
- if (expr>SCALEBMAX) return decInvalid(result, set); /* oops */
- /* [from now on no error possible] */
- if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
- if (DFISSIGNED(dfr)) expr=-expr;
- /* dfl is finite and expr is valid */
- *result=*dfl; /* copy to target */
- return decFloatSetExponent(result, set, GETEXPUN(result)+expr);
- } /* decFloatScaleB */
- /* ------------------------------------------------------------------ */
- /* decFloatShift -- shift the coefficient of a decFloat left or right */
- /* */
- /* result gets the result of shifting dfl */
- /* dfl is the source decFloat to shift */
- /* dfr is the count of digits to shift, an integer (with q=0) */
- /* set is the context */
- /* returns result */
- /* */
- /* The digits of the coefficient of dfl are shifted to the left (if */
- /* dfr is positive) or to the right (if dfr is negative) without */
- /* adjusting the exponent or the sign of dfl. */
- /* */
- /* dfr must be in the range -DECPMAX through +DECPMAX. */
- /* NaNs are propagated as usual. An infinite dfl is unaffected (but */
- /* dfr must be valid). No status is set unless dfr is invalid or an */
- /* operand is an sNaN. The result is canonical. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatShift(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- Int shift; /* dfr as an Int */
- uByte buf[DECPMAX*2]; /* coefficient + padding */
- uInt digits, savestat; /* work */
- bcdnum num; /* .. */
- uInt uiwork; /* for macros */
- if (DFISNAN(dfl)||DFISNAN(dfr)) return decNaNs(result, dfl, dfr, set);
- if (!DFISINT(dfr)) return decInvalid(result, set);
- digits=decFloatDigits(dfr); /* calculate digits */
- if (digits>2) return decInvalid(result, set); /* definitely out of range */
- shift=DPD2BIN[DFWORD(dfr, DECWORDS-1)&0x3ff]; /* is in bottom declet */
- if (shift>DECPMAX) return decInvalid(result, set); /* too big */
- /* [from here on no error or status change is possible] */
- if (DFISINF(dfl)) return decInfinity(result, dfl); /* canonical */
- /* handle no-shift and all-shift (clear to zero) cases */
- if (shift==0) return decCanonical(result, dfl);
- if (shift==DECPMAX) { /* zero with sign */
- uByte sign=(uByte)(DFBYTE(dfl, 0)&0x80); /* save sign bit */
- decFloatZero(result); /* make +0 */
- DFBYTE(result, 0)=(uByte)(DFBYTE(result, 0)|sign); /* and set sign */
- /* [cannot safely use CopySign] */
- return result;
- }
- /* a real shift is needed: 0 < shift < DECPMAX */
- num.sign=DFWORD(dfl, 0)&DECFLOAT_Sign;
- num.exponent=GETEXPUN(dfl);
- num.msd=buf;
- GETCOEFF(dfl, buf);
- if (DFISSIGNED(dfr)) { /* shift right */
- /* edge cases are taken care of, so this is easy */
- num.lsd=buf+DECPMAX-shift-1;
- }
- else { /* shift left -- zero padding needed to right */
- UBFROMUI(buf+DECPMAX, 0); /* 8 will handle most cases */
- UBFROMUI(buf+DECPMAX+4, 0); /* .. */
- if (shift>8) memset(buf+DECPMAX+8, 0, 8+QUAD*18); /* all other cases */
- num.msd+=shift;
- num.lsd=num.msd+DECPMAX-1;
- }
- savestat=set->status; /* record */
- decFinalize(result, &num, set);
- set->status=savestat; /* restore */
- return result;
- } /* decFloatShift */
- /* ------------------------------------------------------------------ */
- /* decFloatSubtract -- subtract a decFloat from another */
- /* */
- /* result gets the result of subtracting dfr from dfl: */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result */
- /* */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatSubtract(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- decFloat temp;
- /* NaNs must propagate without sign change */
- if (DFISNAN(dfr)) return decFloatAdd(result, dfl, dfr, set);
- temp=*dfr; /* make a copy */
- DFBYTE(&temp, 0)^=0x80; /* flip sign */
- return decFloatAdd(result, dfl, &temp, set); /* and add to the lhs */
- } /* decFloatSubtract */
- /* ------------------------------------------------------------------ */
- /* decFloatToInt -- round to 32-bit binary integer (4 flavours) */
- /* */
- /* df is the decFloat to round */
- /* set is the context */
- /* round is the rounding mode to use */
- /* returns a uInt or an Int, rounded according to the name */
- /* */
- /* Invalid will always be signaled if df is a NaN, is Infinite, or is */
- /* outside the range of the target; Inexact will not be signaled for */
- /* simple rounding unless 'Exact' appears in the name. */
- /* ------------------------------------------------------------------ */
- uInt decFloatToUInt32(const decFloat *df, decContext *set,
- enum rounding round) {
- return decToInt32(df, set, round, 0, 1);}
- uInt decFloatToUInt32Exact(const decFloat *df, decContext *set,
- enum rounding round) {
- return decToInt32(df, set, round, 1, 1);}
- Int decFloatToInt32(const decFloat *df, decContext *set,
- enum rounding round) {
- return (Int)decToInt32(df, set, round, 0, 0);}
- Int decFloatToInt32Exact(const decFloat *df, decContext *set,
- enum rounding round) {
- return (Int)decToInt32(df, set, round, 1, 0);}
- /* ------------------------------------------------------------------ */
- /* decFloatToIntegral -- round to integral value (two flavours) */
- /* */
- /* result gets the result */
- /* df is the decFloat to round */
- /* set is the context */
- /* round is the rounding mode to use */
- /* returns result */
- /* */
- /* No exceptions, even Inexact, are raised except for sNaN input, or */
- /* if 'Exact' appears in the name. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatToIntegralValue(decFloat *result, const decFloat *df,
- decContext *set, enum rounding round) {
- return decToIntegral(result, df, set, round, 0);}
- decFloat * decFloatToIntegralExact(decFloat *result, const decFloat *df,
- decContext *set) {
- return decToIntegral(result, df, set, set->round, 1);}
- /* ------------------------------------------------------------------ */
- /* decFloatXor -- logical digitwise XOR of two decFloats */
- /* */
- /* result gets the result of XORing dfl and dfr */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) */
- /* set is the context */
- /* returns result, which will be canonical with sign=0 */
- /* */
- /* The operands must be positive, finite with exponent q=0, and */
- /* comprise just zeros and ones; if not, Invalid operation results. */
- /* ------------------------------------------------------------------ */
- decFloat * decFloatXor(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- if (!DFISUINT01(dfl) || !DFISUINT01(dfr)
- || !DFISCC01(dfl) || !DFISCC01(dfr)) return decInvalid(result, set);
- /* the operands are positive finite integers (q=0) with just 0s and 1s */
- #if DOUBLE
- DFWORD(result, 0)=ZEROWORD
- |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04009124);
- DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x49124491;
- #elif QUAD
- DFWORD(result, 0)=ZEROWORD
- |((DFWORD(dfl, 0) ^ DFWORD(dfr, 0))&0x04000912);
- DFWORD(result, 1)=(DFWORD(dfl, 1) ^ DFWORD(dfr, 1))&0x44912449;
- DFWORD(result, 2)=(DFWORD(dfl, 2) ^ DFWORD(dfr, 2))&0x12449124;
- DFWORD(result, 3)=(DFWORD(dfl, 3) ^ DFWORD(dfr, 3))&0x49124491;
- #endif
- return result;
- } /* decFloatXor */
- /* ------------------------------------------------------------------ */
- /* decInvalid -- set Invalid_operation result */
- /* */
- /* result gets a canonical NaN */
- /* set is the context */
- /* returns result */
- /* */
- /* status has Invalid_operation added */
- /* ------------------------------------------------------------------ */
- static decFloat *decInvalid(decFloat *result, decContext *set) {
- decFloatZero(result);
- DFWORD(result, 0)=DECFLOAT_qNaN;
- set->status|=DEC_Invalid_operation;
- return result;
- } /* decInvalid */
- /* ------------------------------------------------------------------ */
- /* decInfinity -- set canonical Infinity with sign from a decFloat */
- /* */
- /* result gets a canonical Infinity */
- /* df is source decFloat (only the sign is used) */
- /* returns result */
- /* */
- /* df may be the same as result */
- /* ------------------------------------------------------------------ */
- static decFloat *decInfinity(decFloat *result, const decFloat *df) {
- uInt sign=DFWORD(df, 0); /* save source signword */
- decFloatZero(result); /* clear everything */
- DFWORD(result, 0)=DECFLOAT_Inf | (sign & DECFLOAT_Sign);
- return result;
- } /* decInfinity */
- /* ------------------------------------------------------------------ */
- /* decNaNs -- handle NaN argument(s) */
- /* */
- /* result gets the result of handling dfl and dfr, one or both of */
- /* which is a NaN */
- /* dfl is the first decFloat (lhs) */
- /* dfr is the second decFloat (rhs) -- may be NULL for a single- */
- /* operand operation */
- /* set is the context */
- /* returns result */
- /* */
- /* Called when one or both operands is a NaN, and propagates the */
- /* appropriate result to res. When an sNaN is found, it is changed */
- /* to a qNaN and Invalid operation is set. */
- /* ------------------------------------------------------------------ */
- static decFloat *decNaNs(decFloat *result,
- const decFloat *dfl, const decFloat *dfr,
- decContext *set) {
- /* handle sNaNs first */
- if (dfr!=NULL && DFISSNAN(dfr) && !DFISSNAN(dfl)) dfl=dfr; /* use RHS */
- if (DFISSNAN(dfl)) {
- decCanonical(result, dfl); /* propagate canonical sNaN */
- DFWORD(result, 0)&=~(DECFLOAT_qNaN ^ DECFLOAT_sNaN); /* quiet */
- set->status|=DEC_Invalid_operation;
- return result;
- }
- /* one or both is a quiet NaN */
- if (!DFISNAN(dfl)) dfl=dfr; /* RHS must be NaN, use it */
- return decCanonical(result, dfl); /* propagate canonical qNaN */
- } /* decNaNs */
- /* ------------------------------------------------------------------ */
- /* decNumCompare -- numeric comparison of two decFloats */
- /* */
- /* dfl is the left-hand decFloat, which is not a NaN */
- /* dfr is the right-hand decFloat, which is not a NaN */
- /* tot is 1 for total order compare, 0 for simple numeric */
- /* returns -1, 0, or +1 for dfl<dfr, dfl=dfr, dfl>dfr */
- /* */
- /* No error is possible; status and mode are unchanged. */
- /* ------------------------------------------------------------------ */
- static Int decNumCompare(const decFloat *dfl, const decFloat *dfr, Flag tot) {
- Int sigl, sigr; /* LHS and RHS non-0 signums */
- Int shift; /* shift needed to align operands */
- uByte *ub, *uc; /* work */
- uInt uiwork; /* for macros */
- /* buffers +2 if Quad (36 digits), need double plus 4 for safe padding */
- uByte bufl[DECPMAX*2+QUAD*2+4]; /* for LHS coefficient + padding */
- uByte bufr[DECPMAX*2+QUAD*2+4]; /* for RHS coefficient + padding */
- sigl=1;
- if (DFISSIGNED(dfl)) {
- if (!DFISSIGNED(dfr)) { /* -LHS +RHS */
- if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
- return -1; /* RHS wins */
- }
- sigl=-1;
- }
- if (DFISSIGNED(dfr)) {
- if (!DFISSIGNED(dfl)) { /* +LHS -RHS */
- if (DFISZERO(dfl) && DFISZERO(dfr) && !tot) return 0;
- return +1; /* LHS wins */
- }
- }
- /* signs are the same; operand(s) could be zero */
- sigr=-sigl; /* sign to return if abs(RHS) wins */
- if (DFISINF(dfl)) {
- if (DFISINF(dfr)) return 0; /* both infinite & same sign */
- return sigl; /* inf > n */
- }
- if (DFISINF(dfr)) return sigr; /* n < inf [dfl is finite] */
- /* here, both are same sign and finite; calculate their offset */
- shift=GETEXP(dfl)-GETEXP(dfr); /* [0 means aligned] */
- /* [bias can be ignored -- the absolute exponent is not relevant] */
- if (DFISZERO(dfl)) {
- if (!DFISZERO(dfr)) return sigr; /* LHS=0, RHS!=0 */
- /* both are zero, return 0 if both same exponent or numeric compare */
- if (shift==0 || !tot) return 0;
- if (shift>0) return sigl;
- return sigr; /* [shift<0] */
- }
- else { /* LHS!=0 */
- if (DFISZERO(dfr)) return sigl; /* LHS!=0, RHS=0 */
- }
- /* both are known to be non-zero at this point */
- /* if the exponents are so different that the coefficients do not */
- /* overlap (by even one digit) then a full comparison is not needed */
- if (abs(shift)>=DECPMAX) { /* no overlap */
- /* coefficients are known to be non-zero */
- if (shift>0) return sigl;
- return sigr; /* [shift<0] */
- }
- /* decode the coefficients */
- /* (shift both right two if Quad to make a multiple of four) */
- #if QUAD
- UBFROMUI(bufl, 0);
- UBFROMUI(bufr, 0);
- #endif
- GETCOEFF(dfl, bufl+QUAD*2); /* decode from decFloat */
- GETCOEFF(dfr, bufr+QUAD*2); /* .. */
- if (shift==0) { /* aligned; common and easy */
- /* all multiples of four, here */
- for (ub=bufl, uc=bufr; ub<bufl+DECPMAX+QUAD*2; ub+=4, uc+=4) {
- uInt ui=UBTOUI(ub);
- if (ui==UBTOUI(uc)) continue; /* so far so same */
- /* about to find a winner; go by bytes in case little-endian */
- for (;; ub++, uc++) {
- if (*ub>*uc) return sigl; /* difference found */
- if (*ub<*uc) return sigr; /* .. */
- }
- }
- } /* aligned */
- else if (shift>0) { /* lhs to left */
- ub=bufl; /* RHS pointer */
- /* pad bufl so right-aligned; most shifts will fit in 8 */
- UBFROMUI(bufl+DECPMAX+QUAD*2, 0); /* add eight zeros */
- UBFROMUI(bufl+DECPMAX+QUAD*2+4, 0); /* .. */
- if (shift>8) {
- /* more than eight; fill the rest, and also worth doing the */
- /* lead-in by fours */
- uByte *up; /* work */
- uByte *upend=bufl+DECPMAX+QUAD*2+shift;
- for (up=bufl+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
- /* [pads up to 36 in all for Quad] */
- for (;; ub+=4) {
- if (UBTOUI(ub)!=0) return sigl;
- if (ub+4>bufl+shift-4) break;
- }
- }
- /* check remaining leading digits */
- for (; ub<bufl+shift; ub++) if (*ub!=0) return sigl;
- /* now start the overlapped part; bufl has been padded, so the */
- /* comparison can go for the full length of bufr, which is a */
- /* multiple of 4 bytes */
- for (uc=bufr; ; uc+=4, ub+=4) {
- uInt ui=UBTOUI(ub);
- if (ui!=UBTOUI(uc)) { /* mismatch found */
- for (;; uc++, ub++) { /* check from left [little-endian?] */
- if (*ub>*uc) return sigl; /* difference found */
- if (*ub<*uc) return sigr; /* .. */
- }
- } /* mismatch */
- if (uc==bufr+QUAD*2+DECPMAX-4) break; /* all checked */
- }
- } /* shift>0 */
- else { /* shift<0) .. RHS is to left of LHS; mirror shift>0 */
- uc=bufr; /* RHS pointer */
- /* pad bufr so right-aligned; most shifts will fit in 8 */
- UBFROMUI(bufr+DECPMAX+QUAD*2, 0); /* add eight zeros */
- UBFROMUI(bufr+DECPMAX+QUAD*2+4, 0); /* .. */
- if (shift<-8) {
- /* more than eight; fill the rest, and also worth doing the */
- /* lead-in by fours */
- uByte *up; /* work */
- uByte *upend=bufr+DECPMAX+QUAD*2-shift;
- for (up=bufr+DECPMAX+QUAD*2+8; up<upend; up+=4) UBFROMUI(up, 0);
- /* [pads up to 36 in all for Quad] */
- for (;; uc+=4) {
- if (UBTOUI(uc)!=0) return sigr;
- if (uc+4>bufr-shift-4) break;
- }
- }
- /* check remaining leading digits */
- for (; uc<bufr-shift; uc++) if (*uc!=0) return sigr;
- /* now start the overlapped part; bufr has been padded, so the */
- /* comparison can go for the full length of bufl, which is a */
- /* multiple of 4 bytes */
- for (ub=bufl; ; ub+=4, uc+=4) {
- uInt ui=UBTOUI(ub);
- if (ui!=UBTOUI(uc)) { /* mismatch found */
- for (;; ub++, uc++) { /* check from left [little-endian?] */
- if (*ub>*uc) return sigl; /* difference found */
- if (*ub<*uc) return sigr; /* .. */
- }
- } /* mismatch */
- if (ub==bufl+QUAD*2+DECPMAX-4) break; /* all checked */
- }
- } /* shift<0 */
- /* Here when compare equal */
- if (!tot) return 0; /* numerically equal */
- /* total ordering .. exponent matters */
- if (shift>0) return sigl; /* total order by exponent */
- if (shift<0) return sigr; /* .. */
- return 0;
- } /* decNumCompare */
- /* ------------------------------------------------------------------ */
- /* decToInt32 -- local routine to effect ToInteger conversions */
- /* */
- /* df is the decFloat to convert */
- /* set is the context */
- /* rmode is the rounding mode to use */
- /* exact is 1 if Inexact should be signalled */
- /* unsign is 1 if the result a uInt, 0 if an Int (cast to uInt) */
- /* returns 32-bit result as a uInt */
- /* */
- /* Invalid is set is df is a NaN, is infinite, or is out-of-range; in */
- /* these cases 0 is returned. */
- /* ------------------------------------------------------------------ */
- static uInt decToInt32(const decFloat *df, decContext *set,
- enum rounding rmode, Flag exact, Flag unsign) {
- Int exp; /* exponent */
- uInt sourhi, sourpen, sourlo; /* top word from source decFloat .. */
- uInt hi, lo; /* .. penultimate, least, etc. */
- decFloat zero, result; /* work */
- Int i; /* .. */
- /* Start decoding the argument */
- sourhi=DFWORD(df, 0); /* top word */
- exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
- if (EXPISSPECIAL(exp)) { /* is special? */
- set->status|=DEC_Invalid_operation; /* signal */
- return 0;
- }
- /* Here when the argument is finite */
- if (GETEXPUN(df)==0) result=*df; /* already a true integer */
- else { /* need to round to integer */
- enum rounding saveround; /* saver */
- uInt savestatus; /* .. */
- saveround=set->round; /* save rounding mode .. */
- savestatus=set->status; /* .. and status */
- set->round=rmode; /* set mode */
- decFloatZero(&zero); /* make 0E+0 */
- set->status=0; /* clear */
- decFloatQuantize(&result, df, &zero, set); /* [this may fail] */
- set->round=saveround; /* restore rounding mode .. */
- if (exact) set->status|=savestatus; /* include Inexact */
- else set->status=savestatus; /* .. or just original status */
- }
- /* only the last four declets of the coefficient can contain */
- /* non-zero; check for others (and also NaN or Infinity from the */
- /* Quantize) first (see DFISZERO for explanation): */
- /* decFloatShow(&result, "sofar"); */
- #if DOUBLE
- if ((DFWORD(&result, 0)&0x1c03ff00)!=0
- || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
- #elif QUAD
- if ((DFWORD(&result, 2)&0xffffff00)!=0
- || DFWORD(&result, 1)!=0
- || (DFWORD(&result, 0)&0x1c003fff)!=0
- || (DFWORD(&result, 0)&0x60000000)==0x60000000) {
- #endif
- set->status|=DEC_Invalid_operation; /* Invalid or out of range */
- return 0;
- }
- /* get last twelve digits of the coefficent into hi & ho, base */
- /* 10**9 (see GETCOEFFBILL): */
- sourlo=DFWORD(&result, DECWORDS-1);
- lo=DPD2BIN0[sourlo&0x3ff]
- +DPD2BINK[(sourlo>>10)&0x3ff]
- +DPD2BINM[(sourlo>>20)&0x3ff];
- sourpen=DFWORD(&result, DECWORDS-2);
- hi=DPD2BIN0[((sourpen<<2) | (sourlo>>30))&0x3ff];
- /* according to request, check range carefully */
- if (unsign) {
- if (hi>4 || (hi==4 && lo>294967295) || (hi+lo!=0 && DFISSIGNED(&result))) {
- set->status|=DEC_Invalid_operation; /* out of range */
- return 0;
- }
- return hi*BILLION+lo;
- }
- /* signed */
- if (hi>2 || (hi==2 && lo>147483647)) {
- /* handle the usual edge case */
- if (lo==147483648 && hi==2 && DFISSIGNED(&result)) return 0x80000000;
- set->status|=DEC_Invalid_operation; /* truly out of range */
- return 0;
- }
- i=hi*BILLION+lo;
- if (DFISSIGNED(&result)) i=-i;
- return (uInt)i;
- } /* decToInt32 */
- /* ------------------------------------------------------------------ */
- /* decToIntegral -- local routine to effect ToIntegral value */
- /* */
- /* result gets the result */
- /* df is the decFloat to round */
- /* set is the context */
- /* rmode is the rounding mode to use */
- /* exact is 1 if Inexact should be signalled */
- /* returns result */
- /* ------------------------------------------------------------------ */
- static decFloat * decToIntegral(decFloat *result, const decFloat *df,
- decContext *set, enum rounding rmode,
- Flag exact) {
- Int exp; /* exponent */
- uInt sourhi; /* top word from source decFloat */
- enum rounding saveround; /* saver */
- uInt savestatus; /* .. */
- decFloat zero; /* work */
- /* Start decoding the argument */
- sourhi=DFWORD(df, 0); /* top word */
- exp=DECCOMBEXP[sourhi>>26]; /* get exponent high bits (in place) */
- if (EXPISSPECIAL(exp)) { /* is special? */
- /* NaNs are handled as usual */
- if (DFISNAN(df)) return decNaNs(result, df, NULL, set);
- /* must be infinite; return canonical infinity with sign of df */
- return decInfinity(result, df);
- }
- /* Here when the argument is finite */
- /* complete extraction of the exponent */
- exp+=GETECON(df)-DECBIAS; /* .. + continuation and unbias */
- if (exp>=0) return decCanonical(result, df); /* already integral */
- saveround=set->round; /* save rounding mode .. */
- savestatus=set->status; /* .. and status */
- set->round=rmode; /* set mode */
- decFloatZero(&zero); /* make 0E+0 */
- decFloatQuantize(result, df, &zero, set); /* 'integrate'; cannot fail */
- set->round=saveround; /* restore rounding mode .. */
- if (!exact) set->status=savestatus; /* .. and status, unless exact */
- return result;
- } /* decToIntegral */
|