dtoa.c 86 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453345434553456345734583459346034613462346334643465346634673468346934703471347234733474347534763477347834793480348134823483348434853486348734883489349034913492349334943495349634973498349935003501350235033504350535063507350835093510351135123513351435153516351735183519352035213522352335243525352635273528352935303531353235333534353535363537353835393540354135423543354435453546354735483549355035513552355335543555355635573558355935603561356235633564356535663567356835693570357135723573357435753576357735783579358035813582358335843585358635873588358935903591359235933594359535963597359835993600360136023603360436053606360736083609361036113612361336143615361636173618361936203621362236233624362536263627362836293630363136323633363436353636363736383639364036413642364336443645364636473648364936503651365236533654365536563657365836593660366136623663366436653666366736683669367036713672367336743675367636773678367936803681368236833684368536863687368836893690369136923693369436953696369736983699370037013702370337043705370637073708370937103711371237133714371537163717371837193720372137223723372437253726372737283729373037313732373337343735373637373738373937403741374237433744374537463747374837493750375137523753375437553756375737583759376037613762376337643765376637673768376937703771377237733774377537763777377837793780378137823783378437853786378737883789379037913792379337943795379637973798379938003801380238033804380538063807380838093810381138123813381438153816381738183819382038213822382338243825382638273828382938303831383238333834383538363837383838393840384138423843384438453846384738483849385038513852385338543855385638573858385938603861386238633864386538663867386838693870387138723873387438753876387738783879388038813882388338843885388638873888388938903891389238933894389538963897389838993900390139023903390439053906390739083909391039113912391339143915391639173918391939203921392239233924392539263927392839293930393139323933393439353936393739383939394039413942394339443945394639473948394939503951395239533954395539563957395839593960396139623963396439653966396739683969397039713972397339743975397639773978397939803981398239833984398539863987398839893990399139923993399439953996399739983999400040014002400340044005400640074008400940104011401240134014401540164017401840194020402140224023402440254026402740284029403040314032403340344035403640374038403940404041404240434044404540464047404840494050405140524053405440554056405740584059406040614062406340644065406640674068406940704071407240734074407540764077407840794080408140824083408440854086408740884089409040914092409340944095409640974098409941004101410241034104410541064107410841094110411141124113411441154116411741184119412041214122412341244125412641274128412941304131413241334134413541364137413841394140414141424143414441454146414741484149415041514152415341544155415641574158415941604161416241634164416541664167416841694170417141724173417441754176417741784179418041814182418341844185418641874188418941904191419241934194419541964197419841994200420142024203420442054206420742084209421042114212421342144215421642174218421942204221422242234224422542264227422842294230423142324233423442354236423742384239424042414242424342444245424642474248424942504251425242534254425542564257425842594260426142624263426442654266426742684269427042714272427342744275427642774278427942804281428242834284428542864287428842894290429142924293429442954296429742984299430043014302430343044305430643074308430943104311431243134314431543164317431843194320432143224323432443254326432743284329433043314332433343344335433643374338433943404341434243434344434543464347434843494350435143524353435443554356435743584359
  1. /****************************************************************
  2. *
  3. * The author of this software is David M. Gay.
  4. *
  5. * Copyright (c) 1991, 2000, 2001 by Lucent Technologies.
  6. *
  7. * Permission to use, copy, modify, and distribute this software for any
  8. * purpose without fee is hereby granted, provided that this entire notice
  9. * is included in all copies of any software which is or includes a copy
  10. * or modification of this software and in all copies of the supporting
  11. * documentation for such software.
  12. *
  13. * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  14. * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR LUCENT MAKES ANY
  15. * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  16. * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  17. *
  18. ***************************************************************/
  19. /* Please send bug reports to David M. Gay (dmg at acm dot org,
  20. * with " at " changed at "@" and " dot " changed to "."). */
  21. /* On a machine with IEEE extended-precision registers, it is
  22. * necessary to specify double-precision (53-bit) rounding precision
  23. * before invoking strtod or dtoa. If the machine uses (the equivalent
  24. * of) Intel 80x87 arithmetic, the call
  25. * _control87(PC_53, MCW_PC);
  26. * does this with many compilers. Whether this or another call is
  27. * appropriate depends on the compiler; for this to work, it may be
  28. * necessary to #include "float.h" or another system-dependent header
  29. * file.
  30. */
  31. /* strtod for IEEE-, VAX-, and IBM-arithmetic machines.
  32. *
  33. * This strtod returns a nearest machine number to the input decimal
  34. * string (or sets errno to ERANGE). With IEEE arithmetic, ties are
  35. * broken by the IEEE round-even rule. Otherwise ties are broken by
  36. * biased rounding (add half and chop).
  37. *
  38. * Inspired loosely by William D. Clinger's paper "How to Read Floating
  39. * Point Numbers Accurately" [Proc. ACM SIGPLAN '90, pp. 92-101].
  40. *
  41. * Modifications:
  42. *
  43. * 1. We only require IEEE, IBM, or VAX double-precision
  44. * arithmetic (not IEEE double-extended).
  45. * 2. We get by with floating-point arithmetic in a case that
  46. * Clinger missed -- when we're computing d * 10^n
  47. * for a small integer d and the integer n is not too
  48. * much larger than 22 (the maximum integer k for which
  49. * we can represent 10^k exactly), we may be able to
  50. * compute (d*10^k) * 10^(e-k) with just one roundoff.
  51. * 3. Rather than a bit-at-a-time adjustment of the binary
  52. * result in the hard case, we use floating-point
  53. * arithmetic to determine the adjustment to within
  54. * one bit; only in really hard cases do we need to
  55. * compute a second residual.
  56. * 4. Because of 3., we don't need a large table of powers of 10
  57. * for ten-to-e (just some small tables, e.g. of 10^k
  58. * for 0 <= k <= 22).
  59. */
  60. /*
  61. * #define IEEE_8087 for IEEE-arithmetic machines where the least
  62. * significant byte has the lowest address.
  63. * #define IEEE_MC68k for IEEE-arithmetic machines where the most
  64. * significant byte has the lowest address.
  65. * #define Long int on machines with 32-bit ints and 64-bit longs.
  66. * #define IBM for IBM mainframe-style floating-point arithmetic.
  67. * #define VAX for VAX-style floating-point arithmetic (D_floating).
  68. * #define No_leftright to omit left-right logic in fast floating-point
  69. * computation of dtoa. This will cause dtoa modes 4 and 5 to be
  70. * treated the same as modes 2 and 3 for some inputs.
  71. * #define Honor_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  72. * and strtod and dtoa should round accordingly. Unless Trust_FLT_ROUNDS
  73. * is also #defined, fegetround() will be queried for the rounding mode.
  74. * Note that both FLT_ROUNDS and fegetround() are specified by the C99
  75. * standard (and are specified to be consistent, with fesetround()
  76. * affecting the value of FLT_ROUNDS), but that some (Linux) systems
  77. * do not work correctly in this regard, so using fegetround() is more
  78. * portable than using FLT_ROUNDS directly.
  79. * #define Check_FLT_ROUNDS if FLT_ROUNDS can assume the values 2 or 3
  80. * and Honor_FLT_ROUNDS is not #defined.
  81. * #define RND_PRODQUOT to use rnd_prod and rnd_quot (assembly routines
  82. * that use extended-precision instructions to compute rounded
  83. * products and quotients) with IBM.
  84. * #define ROUND_BIASED for IEEE-format with biased rounding and arithmetic
  85. * that rounds toward +Infinity.
  86. * #define ROUND_BIASED_without_Round_Up for IEEE-format with biased
  87. * rounding when the underlying floating-point arithmetic uses
  88. * unbiased rounding. This prevent using ordinary floating-point
  89. * arithmetic when the result could be computed with one rounding error.
  90. * #define Inaccurate_Divide for IEEE-format with correctly rounded
  91. * products but inaccurate quotients, e.g., for Intel i860.
  92. * #define NO_LONG_LONG on machines that do not have a "long long"
  93. * integer type (of >= 64 bits). On such machines, you can
  94. * #define Just_16 to store 16 bits per 32-bit Long when doing
  95. * high-precision integer arithmetic. Whether this speeds things
  96. * up or slows things down depends on the machine and the number
  97. * being converted. If long long is available and the name is
  98. * something other than "long long", #define Llong to be the name,
  99. * and if "unsigned Llong" does not work as an unsigned version of
  100. * Llong, #define #ULLong to be the corresponding unsigned type.
  101. * #define KR_headers for old-style C function headers.
  102. * #define Bad_float_h if your system lacks a float.h or if it does not
  103. * define some or all of DBL_DIG, DBL_MAX_10_EXP, DBL_MAX_EXP,
  104. * FLT_RADIX, FLT_ROUNDS, and DBL_MAX.
  105. * #define MALLOC your_malloc, where your_malloc(n) acts like malloc(n)
  106. * if memory is available and otherwise does something you deem
  107. * appropriate. If MALLOC is undefined, malloc will be invoked
  108. * directly -- and assumed always to succeed. Similarly, if you
  109. * want something other than the system's free() to be called to
  110. * recycle memory acquired from MALLOC, #define FREE to be the
  111. * name of the alternate routine. (FREE or free is only called in
  112. * pathological cases, e.g., in a dtoa call after a dtoa return in
  113. * mode 3 with thousands of digits requested.)
  114. * #define Omit_Private_Memory to omit logic (added Jan. 1998) for making
  115. * memory allocations from a private pool of memory when possible.
  116. * When used, the private pool is PRIVATE_MEM bytes long: 2304 bytes,
  117. * unless #defined to be a different length. This default length
  118. * suffices to get rid of MALLOC calls except for unusual cases,
  119. * such as decimal-to-binary conversion of a very long string of
  120. * digits. The longest string dtoa can return is about 751 bytes
  121. * long. For conversions by strtod of strings of 800 digits and
  122. * all dtoa conversions in single-threaded executions with 8-byte
  123. * pointers, PRIVATE_MEM >= 7400 appears to suffice; with 4-byte
  124. * pointers, PRIVATE_MEM >= 7112 appears adequate.
  125. * #define NO_INFNAN_CHECK if you do not wish to have INFNAN_CHECK
  126. * #defined automatically on IEEE systems. On such systems,
  127. * when INFNAN_CHECK is #defined, strtod checks
  128. * for Infinity and NaN (case insensitively). On some systems
  129. * (e.g., some HP systems), it may be necessary to #define NAN_WORD0
  130. * appropriately -- to the most significant word of a quiet NaN.
  131. * (On HP Series 700/800 machines, -DNAN_WORD0=0x7ff40000 works.)
  132. * When INFNAN_CHECK is #defined and No_Hex_NaN is not #defined,
  133. * strtod also accepts (case insensitively) strings of the form
  134. * NaN(x), where x is a string of hexadecimal digits and spaces;
  135. * if there is only one string of hexadecimal digits, it is taken
  136. * for the 52 fraction bits of the resulting NaN; if there are two
  137. * or more strings of hex digits, the first is for the high 20 bits,
  138. * the second and subsequent for the low 32 bits, with intervening
  139. * white space ignored; but if this results in none of the 52
  140. * fraction bits being on (an IEEE Infinity symbol), then NAN_WORD0
  141. * and NAN_WORD1 are used instead.
  142. * #define MULTIPLE_THREADS if the system offers preemptively scheduled
  143. * multiple threads. In this case, you must provide (or suitably
  144. * #define) two locks, acquired by ACQUIRE_DTOA_LOCK(n) and freed
  145. * by FREE_DTOA_LOCK(n) for n = 0 or 1. (The second lock, accessed
  146. * in pow5mult, ensures lazy evaluation of only one copy of high
  147. * powers of 5; omitting this lock would introduce a small
  148. * probability of wasting memory, but would otherwise be harmless.)
  149. * You must also invoke freedtoa(s) to free the value s returned by
  150. * dtoa. You may do so whether or not MULTIPLE_THREADS is #defined.
  151. * #define NO_IEEE_Scale to disable new (Feb. 1997) logic in strtod that
  152. * avoids underflows on inputs whose result does not underflow.
  153. * If you #define NO_IEEE_Scale on a machine that uses IEEE-format
  154. * floating-point numbers and flushes underflows to zero rather
  155. * than implementing gradual underflow, then you must also #define
  156. * Sudden_Underflow.
  157. * #define USE_LOCALE to use the current locale's decimal_point value.
  158. * #define SET_INEXACT if IEEE arithmetic is being used and extra
  159. * computation should be done to set the inexact flag when the
  160. * result is inexact and avoid setting inexact when the result
  161. * is exact. In this case, dtoa.c must be compiled in
  162. * an environment, perhaps provided by #include "dtoa.c" in a
  163. * suitable wrapper, that defines two functions,
  164. * int get_inexact(void);
  165. * void clear_inexact(void);
  166. * such that get_inexact() returns a nonzero value if the
  167. * inexact bit is already set, and clear_inexact() sets the
  168. * inexact bit to 0. When SET_INEXACT is #defined, strtod
  169. * also does extra computations to set the underflow and overflow
  170. * flags when appropriate (i.e., when the result is tiny and
  171. * inexact or when it is a numeric value rounded to +-infinity).
  172. * #define NO_ERRNO if strtod should not assign errno = ERANGE when
  173. * the result overflows to +-Infinity or underflows to 0.
  174. * #define NO_HEX_FP to omit recognition of hexadecimal floating-point
  175. * values by strtod.
  176. * #define NO_STRTOD_BIGCOMP (on IEEE-arithmetic systems only for now)
  177. * to disable logic for "fast" testing of very long input strings
  178. * to strtod. This testing proceeds by initially truncating the
  179. * input string, then if necessary comparing the whole string with
  180. * a decimal expansion to decide close cases. This logic is only
  181. * used for input more than STRTOD_DIGLIM digits long (default 40).
  182. */
  183. #if 0
  184. #include "dtoa_config.h"
  185. #ifndef Long
  186. #define Long long
  187. #endif
  188. #ifndef ULong
  189. typedef unsigned Long ULong;
  190. #endif
  191. #ifdef DEBUG
  192. #include "stdio.h"
  193. #define Bug(x) {fprintf(stderr, "%s\n", x); exit(1);}
  194. #endif
  195. #include "stdlib.h"
  196. #include "string.h"
  197. #ifdef USE_LOCALE
  198. #include "locale.h"
  199. #endif
  200. #ifdef Honor_FLT_ROUNDS
  201. #ifndef Trust_FLT_ROUNDS
  202. #include <fenv.h>
  203. #endif
  204. #endif
  205. #ifdef MALLOC
  206. #ifdef KR_headers
  207. extern char *MALLOC();
  208. #else
  209. extern void *MALLOC(size_t);
  210. #endif
  211. #else
  212. #define MALLOC malloc
  213. #endif
  214. #ifndef Omit_Private_Memory
  215. #ifndef PRIVATE_MEM
  216. #define PRIVATE_MEM 2304
  217. #endif
  218. #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
  219. static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
  220. #endif
  221. #undef IEEE_Arith
  222. #undef Avoid_Underflow
  223. #ifdef IEEE_MC68k
  224. #define IEEE_Arith
  225. #endif
  226. #ifdef IEEE_8087
  227. #define IEEE_Arith
  228. #endif
  229. #ifdef IEEE_Arith
  230. #ifndef NO_INFNAN_CHECK
  231. #undef INFNAN_CHECK
  232. #define INFNAN_CHECK
  233. #endif
  234. #else
  235. #undef INFNAN_CHECK
  236. #define NO_STRTOD_BIGCOMP
  237. #endif
  238. #include "errno.h"
  239. #ifdef Bad_float_h
  240. #ifdef IEEE_Arith
  241. #define DBL_DIG 15
  242. #define DBL_MAX_10_EXP 308
  243. #define DBL_MAX_EXP 1024
  244. #define FLT_RADIX 2
  245. #endif /*IEEE_Arith*/
  246. #ifdef IBM
  247. #define DBL_DIG 16
  248. #define DBL_MAX_10_EXP 75
  249. #define DBL_MAX_EXP 63
  250. #define FLT_RADIX 16
  251. #define DBL_MAX 7.2370055773322621e+75
  252. #endif
  253. #ifdef VAX
  254. #define DBL_DIG 16
  255. #define DBL_MAX_10_EXP 38
  256. #define DBL_MAX_EXP 127
  257. #define FLT_RADIX 2
  258. #define DBL_MAX 1.7014118346046923e+38
  259. #endif
  260. #ifndef LONG_MAX
  261. #define LONG_MAX 2147483647
  262. #endif
  263. #else /* ifndef Bad_float_h */
  264. #include "float.h"
  265. #endif /* Bad_float_h */
  266. #ifndef __MATH_H__
  267. #include "math.h"
  268. #endif
  269. #ifdef __cplusplus
  270. extern "C" {
  271. #endif
  272. #ifndef CONST
  273. #ifdef KR_headers
  274. #define CONST /* blank */
  275. #else
  276. #define CONST const
  277. #endif
  278. #endif
  279. #if defined(IEEE_8087) + defined(IEEE_MC68k) + defined(VAX) + defined(IBM) != 1
  280. Exactly one of IEEE_8087, IEEE_MC68k, VAX, or IBM should be defined.
  281. #endif
  282. typedef union { double d; ULong L[2]; } U;
  283. #ifdef IEEE_8087
  284. #define word0(x) (x)->L[1]
  285. #define word1(x) (x)->L[0]
  286. #else
  287. #define word0(x) (x)->L[0]
  288. #define word1(x) (x)->L[1]
  289. #endif
  290. #define dval(x) (x)->d
  291. #ifndef STRTOD_DIGLIM
  292. #define STRTOD_DIGLIM 40
  293. #endif
  294. #ifdef DIGLIM_DEBUG
  295. extern int strtod_diglim;
  296. #else
  297. #define strtod_diglim STRTOD_DIGLIM
  298. #endif
  299. /* The following definition of Storeinc is appropriate for MIPS processors.
  300. * An alternative that might be better on some machines is
  301. * #define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)
  302. */
  303. #if defined(IEEE_8087) + defined(VAX)
  304. #define Storeinc(a,b,c) (((unsigned short *)a)[1] = (unsigned short)b, \
  305. ((unsigned short *)a)[0] = (unsigned short)c, a++)
  306. #else
  307. #define Storeinc(a,b,c) (((unsigned short *)a)[0] = (unsigned short)b, \
  308. ((unsigned short *)a)[1] = (unsigned short)c, a++)
  309. #endif
  310. /* #define P DBL_MANT_DIG */
  311. /* Ten_pmax = floor(P*log(2)/log(5)) */
  312. /* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */
  313. /* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */
  314. /* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */
  315. #ifdef IEEE_Arith
  316. #define Exp_shift 20
  317. #define Exp_shift1 20
  318. #define Exp_msk1 0x100000
  319. #define Exp_msk11 0x100000
  320. #define Exp_mask 0x7ff00000
  321. #define P 53
  322. #define Nbits 53
  323. #define Bias 1023
  324. #define Emax 1023
  325. #define Emin (-1022)
  326. #define Exp_1 0x3ff00000
  327. #define Exp_11 0x3ff00000
  328. #define Ebits 11
  329. #define Frac_mask 0xfffff
  330. #define Frac_mask1 0xfffff
  331. #define Ten_pmax 22
  332. #define Bletch 0x10
  333. #define Bndry_mask 0xfffff
  334. #define Bndry_mask1 0xfffff
  335. #define LSB 1
  336. #define Sign_bit 0x80000000
  337. #define Log2P 1
  338. #define Tiny0 0
  339. #define Tiny1 1
  340. #define Quick_max 14
  341. #define Int_max 14
  342. #ifndef NO_IEEE_Scale
  343. #define Avoid_Underflow
  344. #ifdef Flush_Denorm /* debugging option */
  345. #undef Sudden_Underflow
  346. #endif
  347. #endif
  348. #ifndef Flt_Rounds
  349. #ifdef FLT_ROUNDS
  350. #define Flt_Rounds FLT_ROUNDS
  351. #else
  352. #define Flt_Rounds 1
  353. #endif
  354. #endif /*Flt_Rounds*/
  355. #ifdef Honor_FLT_ROUNDS
  356. #undef Check_FLT_ROUNDS
  357. #define Check_FLT_ROUNDS
  358. #else
  359. #define Rounding Flt_Rounds
  360. #endif
  361. #else /* ifndef IEEE_Arith */
  362. #undef Check_FLT_ROUNDS
  363. #undef Honor_FLT_ROUNDS
  364. #undef SET_INEXACT
  365. #undef Sudden_Underflow
  366. #define Sudden_Underflow
  367. #ifdef IBM
  368. #undef Flt_Rounds
  369. #define Flt_Rounds 0
  370. #define Exp_shift 24
  371. #define Exp_shift1 24
  372. #define Exp_msk1 0x1000000
  373. #define Exp_msk11 0x1000000
  374. #define Exp_mask 0x7f000000
  375. #define P 14
  376. #define Nbits 56
  377. #define Bias 65
  378. #define Emax 248
  379. #define Emin (-260)
  380. #define Exp_1 0x41000000
  381. #define Exp_11 0x41000000
  382. #define Ebits 8 /* exponent has 7 bits, but 8 is the right value in b2d */
  383. #define Frac_mask 0xffffff
  384. #define Frac_mask1 0xffffff
  385. #define Bletch 4
  386. #define Ten_pmax 22
  387. #define Bndry_mask 0xefffff
  388. #define Bndry_mask1 0xffffff
  389. #define LSB 1
  390. #define Sign_bit 0x80000000
  391. #define Log2P 4
  392. #define Tiny0 0x100000
  393. #define Tiny1 0
  394. #define Quick_max 14
  395. #define Int_max 15
  396. #else /* VAX */
  397. #undef Flt_Rounds
  398. #define Flt_Rounds 1
  399. #define Exp_shift 23
  400. #define Exp_shift1 7
  401. #define Exp_msk1 0x80
  402. #define Exp_msk11 0x800000
  403. #define Exp_mask 0x7f80
  404. #define P 56
  405. #define Nbits 56
  406. #define Bias 129
  407. #define Emax 126
  408. #define Emin (-129)
  409. #define Exp_1 0x40800000
  410. #define Exp_11 0x4080
  411. #define Ebits 8
  412. #define Frac_mask 0x7fffff
  413. #define Frac_mask1 0xffff007f
  414. #define Ten_pmax 24
  415. #define Bletch 2
  416. #define Bndry_mask 0xffff007f
  417. #define Bndry_mask1 0xffff007f
  418. #define LSB 0x10000
  419. #define Sign_bit 0x8000
  420. #define Log2P 1
  421. #define Tiny0 0x80
  422. #define Tiny1 0
  423. #define Quick_max 15
  424. #define Int_max 15
  425. #endif /* IBM, VAX */
  426. #endif /* IEEE_Arith */
  427. #ifndef IEEE_Arith
  428. #define ROUND_BIASED
  429. #else
  430. #ifdef ROUND_BIASED_without_Round_Up
  431. #undef ROUND_BIASED
  432. #define ROUND_BIASED
  433. #endif
  434. #endif
  435. #ifdef RND_PRODQUOT
  436. #define rounded_product(a,b) a = rnd_prod(a, b)
  437. #define rounded_quotient(a,b) a = rnd_quot(a, b)
  438. #ifdef KR_headers
  439. extern double rnd_prod(), rnd_quot();
  440. #else
  441. extern double rnd_prod(double, double), rnd_quot(double, double);
  442. #endif
  443. #else
  444. #define rounded_product(a,b) a *= b
  445. #define rounded_quotient(a,b) a /= b
  446. #endif
  447. #define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))
  448. #define Big1 0xffffffff
  449. #ifndef Pack_32
  450. #define Pack_32
  451. #endif
  452. typedef struct BCinfo BCinfo;
  453. struct
  454. BCinfo { int dp0, dp1, dplen, dsign, e0, inexact, nd, nd0, rounding, scale, uflchk; };
  455. #ifdef KR_headers
  456. #define FFFFFFFF ((((unsigned long)0xffff)<<16)|(unsigned long)0xffff)
  457. #else
  458. #define FFFFFFFF 0xffffffffUL
  459. #endif
  460. #ifdef NO_LONG_LONG
  461. #undef ULLong
  462. #ifdef Just_16
  463. #undef Pack_32
  464. /* When Pack_32 is not defined, we store 16 bits per 32-bit Long.
  465. * This makes some inner loops simpler and sometimes saves work
  466. * during multiplications, but it often seems to make things slightly
  467. * slower. Hence the default is now to store 32 bits per Long.
  468. */
  469. #endif
  470. #else /* long long available */
  471. #ifndef Llong
  472. #define Llong long long
  473. #endif
  474. #ifndef ULLong
  475. #define ULLong unsigned Llong
  476. #endif
  477. #endif /* NO_LONG_LONG */
  478. #ifndef MULTIPLE_THREADS
  479. #define ACQUIRE_DTOA_LOCK(n) /*nothing*/
  480. #define FREE_DTOA_LOCK(n) /*nothing*/
  481. #endif
  482. #define Kmax 7
  483. #ifdef __cplusplus
  484. extern "C" double fpconv_strtod(const char *s00, char **se);
  485. extern "C" char *dtoa(double d, int mode, int ndigits,
  486. int *decpt, int *sign, char **rve);
  487. #endif
  488. struct
  489. Bigint {
  490. struct Bigint *next;
  491. int k, maxwds, sign, wds;
  492. ULong x[1];
  493. };
  494. typedef struct Bigint Bigint;
  495. static Bigint *freelist[Kmax+1];
  496. static Bigint *
  497. Balloc
  498. #ifdef KR_headers
  499. (k) int k;
  500. #else
  501. (int k)
  502. #endif
  503. {
  504. int x;
  505. Bigint *rv;
  506. #ifndef Omit_Private_Memory
  507. unsigned int len;
  508. #endif
  509. ACQUIRE_DTOA_LOCK(0);
  510. /* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
  511. /* but this case seems very unlikely. */
  512. if (k <= Kmax && (rv = freelist[k]))
  513. freelist[k] = rv->next;
  514. else {
  515. x = 1 << k;
  516. #ifdef Omit_Private_Memory
  517. rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
  518. #else
  519. len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
  520. /sizeof(double);
  521. if (k <= Kmax && pmem_next - private_mem + len <= PRIVATE_mem) {
  522. rv = (Bigint*)pmem_next;
  523. pmem_next += len;
  524. }
  525. else
  526. rv = (Bigint*)MALLOC(len*sizeof(double));
  527. #endif
  528. rv->k = k;
  529. rv->maxwds = x;
  530. }
  531. FREE_DTOA_LOCK(0);
  532. rv->sign = rv->wds = 0;
  533. return rv;
  534. }
  535. static void
  536. Bfree
  537. #ifdef KR_headers
  538. (v) Bigint *v;
  539. #else
  540. (Bigint *v)
  541. #endif
  542. {
  543. if (v) {
  544. if (v->k > Kmax)
  545. #ifdef FREE
  546. FREE((void*)v);
  547. #else
  548. free((void*)v);
  549. #endif
  550. else {
  551. ACQUIRE_DTOA_LOCK(0);
  552. v->next = freelist[v->k];
  553. freelist[v->k] = v;
  554. FREE_DTOA_LOCK(0);
  555. }
  556. }
  557. }
  558. #define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \
  559. y->wds*sizeof(Long) + 2*sizeof(int))
  560. static Bigint *
  561. multadd
  562. #ifdef KR_headers
  563. (b, m, a) Bigint *b; int m, a;
  564. #else
  565. (Bigint *b, int m, int a) /* multiply by m and add a */
  566. #endif
  567. {
  568. int i, wds;
  569. #ifdef ULLong
  570. ULong *x;
  571. ULLong carry, y;
  572. #else
  573. ULong carry, *x, y;
  574. #ifdef Pack_32
  575. ULong xi, z;
  576. #endif
  577. #endif
  578. Bigint *b1;
  579. wds = b->wds;
  580. x = b->x;
  581. i = 0;
  582. carry = a;
  583. do {
  584. #ifdef ULLong
  585. y = *x * (ULLong)m + carry;
  586. carry = y >> 32;
  587. *x++ = y & FFFFFFFF;
  588. #else
  589. #ifdef Pack_32
  590. xi = *x;
  591. y = (xi & 0xffff) * m + carry;
  592. z = (xi >> 16) * m + (y >> 16);
  593. carry = z >> 16;
  594. *x++ = (z << 16) + (y & 0xffff);
  595. #else
  596. y = *x * m + carry;
  597. carry = y >> 16;
  598. *x++ = y & 0xffff;
  599. #endif
  600. #endif
  601. }
  602. while(++i < wds);
  603. if (carry) {
  604. if (wds >= b->maxwds) {
  605. b1 = Balloc(b->k+1);
  606. Bcopy(b1, b);
  607. Bfree(b);
  608. b = b1;
  609. }
  610. b->x[wds++] = carry;
  611. b->wds = wds;
  612. }
  613. return b;
  614. }
  615. static Bigint *
  616. s2b
  617. #ifdef KR_headers
  618. (s, nd0, nd, y9, dplen) CONST char *s; int nd0, nd, dplen; ULong y9;
  619. #else
  620. (const char *s, int nd0, int nd, ULong y9, int dplen)
  621. #endif
  622. {
  623. Bigint *b;
  624. int i, k;
  625. Long x, y;
  626. x = (nd + 8) / 9;
  627. for(k = 0, y = 1; x > y; y <<= 1, k++) ;
  628. #ifdef Pack_32
  629. b = Balloc(k);
  630. b->x[0] = y9;
  631. b->wds = 1;
  632. #else
  633. b = Balloc(k+1);
  634. b->x[0] = y9 & 0xffff;
  635. b->wds = (b->x[1] = y9 >> 16) ? 2 : 1;
  636. #endif
  637. i = 9;
  638. if (9 < nd0) {
  639. s += 9;
  640. do b = multadd(b, 10, *s++ - '0');
  641. while(++i < nd0);
  642. s += dplen;
  643. }
  644. else
  645. s += dplen + 9;
  646. for(; i < nd; i++)
  647. b = multadd(b, 10, *s++ - '0');
  648. return b;
  649. }
  650. static int
  651. hi0bits
  652. #ifdef KR_headers
  653. (x) ULong x;
  654. #else
  655. (ULong x)
  656. #endif
  657. {
  658. int k = 0;
  659. if (!(x & 0xffff0000)) {
  660. k = 16;
  661. x <<= 16;
  662. }
  663. if (!(x & 0xff000000)) {
  664. k += 8;
  665. x <<= 8;
  666. }
  667. if (!(x & 0xf0000000)) {
  668. k += 4;
  669. x <<= 4;
  670. }
  671. if (!(x & 0xc0000000)) {
  672. k += 2;
  673. x <<= 2;
  674. }
  675. if (!(x & 0x80000000)) {
  676. k++;
  677. if (!(x & 0x40000000))
  678. return 32;
  679. }
  680. return k;
  681. }
  682. static int
  683. lo0bits
  684. #ifdef KR_headers
  685. (y) ULong *y;
  686. #else
  687. (ULong *y)
  688. #endif
  689. {
  690. int k;
  691. ULong x = *y;
  692. if (x & 7) {
  693. if (x & 1)
  694. return 0;
  695. if (x & 2) {
  696. *y = x >> 1;
  697. return 1;
  698. }
  699. *y = x >> 2;
  700. return 2;
  701. }
  702. k = 0;
  703. if (!(x & 0xffff)) {
  704. k = 16;
  705. x >>= 16;
  706. }
  707. if (!(x & 0xff)) {
  708. k += 8;
  709. x >>= 8;
  710. }
  711. if (!(x & 0xf)) {
  712. k += 4;
  713. x >>= 4;
  714. }
  715. if (!(x & 0x3)) {
  716. k += 2;
  717. x >>= 2;
  718. }
  719. if (!(x & 1)) {
  720. k++;
  721. x >>= 1;
  722. if (!x)
  723. return 32;
  724. }
  725. *y = x;
  726. return k;
  727. }
  728. static Bigint *
  729. i2b
  730. #ifdef KR_headers
  731. (i) int i;
  732. #else
  733. (int i)
  734. #endif
  735. {
  736. Bigint *b;
  737. b = Balloc(1);
  738. b->x[0] = i;
  739. b->wds = 1;
  740. return b;
  741. }
  742. static Bigint *
  743. mult
  744. #ifdef KR_headers
  745. (a, b) Bigint *a, *b;
  746. #else
  747. (Bigint *a, Bigint *b)
  748. #endif
  749. {
  750. Bigint *c;
  751. int k, wa, wb, wc;
  752. ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
  753. ULong y;
  754. #ifdef ULLong
  755. ULLong carry, z;
  756. #else
  757. ULong carry, z;
  758. #ifdef Pack_32
  759. ULong z2;
  760. #endif
  761. #endif
  762. if (a->wds < b->wds) {
  763. c = a;
  764. a = b;
  765. b = c;
  766. }
  767. k = a->k;
  768. wa = a->wds;
  769. wb = b->wds;
  770. wc = wa + wb;
  771. if (wc > a->maxwds)
  772. k++;
  773. c = Balloc(k);
  774. for(x = c->x, xa = x + wc; x < xa; x++)
  775. *x = 0;
  776. xa = a->x;
  777. xae = xa + wa;
  778. xb = b->x;
  779. xbe = xb + wb;
  780. xc0 = c->x;
  781. #ifdef ULLong
  782. for(; xb < xbe; xc0++) {
  783. if ((y = *xb++)) {
  784. x = xa;
  785. xc = xc0;
  786. carry = 0;
  787. do {
  788. z = *x++ * (ULLong)y + *xc + carry;
  789. carry = z >> 32;
  790. *xc++ = z & FFFFFFFF;
  791. }
  792. while(x < xae);
  793. *xc = carry;
  794. }
  795. }
  796. #else
  797. #ifdef Pack_32
  798. for(; xb < xbe; xb++, xc0++) {
  799. if (y = *xb & 0xffff) {
  800. x = xa;
  801. xc = xc0;
  802. carry = 0;
  803. do {
  804. z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;
  805. carry = z >> 16;
  806. z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;
  807. carry = z2 >> 16;
  808. Storeinc(xc, z2, z);
  809. }
  810. while(x < xae);
  811. *xc = carry;
  812. }
  813. if (y = *xb >> 16) {
  814. x = xa;
  815. xc = xc0;
  816. carry = 0;
  817. z2 = *xc;
  818. do {
  819. z = (*x & 0xffff) * y + (*xc >> 16) + carry;
  820. carry = z >> 16;
  821. Storeinc(xc, z, z2);
  822. z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;
  823. carry = z2 >> 16;
  824. }
  825. while(x < xae);
  826. *xc = z2;
  827. }
  828. }
  829. #else
  830. for(; xb < xbe; xc0++) {
  831. if (y = *xb++) {
  832. x = xa;
  833. xc = xc0;
  834. carry = 0;
  835. do {
  836. z = *x++ * y + *xc + carry;
  837. carry = z >> 16;
  838. *xc++ = z & 0xffff;
  839. }
  840. while(x < xae);
  841. *xc = carry;
  842. }
  843. }
  844. #endif
  845. #endif
  846. for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
  847. c->wds = wc;
  848. return c;
  849. }
  850. static Bigint *p5s;
  851. static Bigint *
  852. pow5mult
  853. #ifdef KR_headers
  854. (b, k) Bigint *b; int k;
  855. #else
  856. (Bigint *b, int k)
  857. #endif
  858. {
  859. Bigint *b1, *p5, *p51;
  860. int i;
  861. static int p05[3] = { 5, 25, 125 };
  862. if ((i = k & 3))
  863. b = multadd(b, p05[i-1], 0);
  864. if (!(k >>= 2))
  865. return b;
  866. if (!(p5 = p5s)) {
  867. /* first time */
  868. #ifdef MULTIPLE_THREADS
  869. ACQUIRE_DTOA_LOCK(1);
  870. if (!(p5 = p5s)) {
  871. p5 = p5s = i2b(625);
  872. p5->next = 0;
  873. }
  874. FREE_DTOA_LOCK(1);
  875. #else
  876. p5 = p5s = i2b(625);
  877. p5->next = 0;
  878. #endif
  879. }
  880. for(;;) {
  881. if (k & 1) {
  882. b1 = mult(b, p5);
  883. Bfree(b);
  884. b = b1;
  885. }
  886. if (!(k >>= 1))
  887. break;
  888. if (!(p51 = p5->next)) {
  889. #ifdef MULTIPLE_THREADS
  890. ACQUIRE_DTOA_LOCK(1);
  891. if (!(p51 = p5->next)) {
  892. p51 = p5->next = mult(p5,p5);
  893. p51->next = 0;
  894. }
  895. FREE_DTOA_LOCK(1);
  896. #else
  897. p51 = p5->next = mult(p5,p5);
  898. p51->next = 0;
  899. #endif
  900. }
  901. p5 = p51;
  902. }
  903. return b;
  904. }
  905. static Bigint *
  906. lshift
  907. #ifdef KR_headers
  908. (b, k) Bigint *b; int k;
  909. #else
  910. (Bigint *b, int k)
  911. #endif
  912. {
  913. int i, k1, n, n1;
  914. Bigint *b1;
  915. ULong *x, *x1, *xe, z;
  916. #ifdef Pack_32
  917. n = k >> 5;
  918. #else
  919. n = k >> 4;
  920. #endif
  921. k1 = b->k;
  922. n1 = n + b->wds + 1;
  923. for(i = b->maxwds; n1 > i; i <<= 1)
  924. k1++;
  925. b1 = Balloc(k1);
  926. x1 = b1->x;
  927. for(i = 0; i < n; i++)
  928. *x1++ = 0;
  929. x = b->x;
  930. xe = x + b->wds;
  931. #ifdef Pack_32
  932. if (k &= 0x1f) {
  933. k1 = 32 - k;
  934. z = 0;
  935. do {
  936. *x1++ = *x << k | z;
  937. z = *x++ >> k1;
  938. }
  939. while(x < xe);
  940. if ((*x1 = z))
  941. ++n1;
  942. }
  943. #else
  944. if (k &= 0xf) {
  945. k1 = 16 - k;
  946. z = 0;
  947. do {
  948. *x1++ = *x << k & 0xffff | z;
  949. z = *x++ >> k1;
  950. }
  951. while(x < xe);
  952. if (*x1 = z)
  953. ++n1;
  954. }
  955. #endif
  956. else do
  957. *x1++ = *x++;
  958. while(x < xe);
  959. b1->wds = n1 - 1;
  960. Bfree(b);
  961. return b1;
  962. }
  963. static int
  964. cmp
  965. #ifdef KR_headers
  966. (a, b) Bigint *a, *b;
  967. #else
  968. (Bigint *a, Bigint *b)
  969. #endif
  970. {
  971. ULong *xa, *xa0, *xb, *xb0;
  972. int i, j;
  973. i = a->wds;
  974. j = b->wds;
  975. #ifdef DEBUG
  976. if (i > 1 && !a->x[i-1])
  977. Bug("cmp called with a->x[a->wds-1] == 0");
  978. if (j > 1 && !b->x[j-1])
  979. Bug("cmp called with b->x[b->wds-1] == 0");
  980. #endif
  981. if (i -= j)
  982. return i;
  983. xa0 = a->x;
  984. xa = xa0 + j;
  985. xb0 = b->x;
  986. xb = xb0 + j;
  987. for(;;) {
  988. if (*--xa != *--xb)
  989. return *xa < *xb ? -1 : 1;
  990. if (xa <= xa0)
  991. break;
  992. }
  993. return 0;
  994. }
  995. static Bigint *
  996. diff
  997. #ifdef KR_headers
  998. (a, b) Bigint *a, *b;
  999. #else
  1000. (Bigint *a, Bigint *b)
  1001. #endif
  1002. {
  1003. Bigint *c;
  1004. int i, wa, wb;
  1005. ULong *xa, *xae, *xb, *xbe, *xc;
  1006. #ifdef ULLong
  1007. ULLong borrow, y;
  1008. #else
  1009. ULong borrow, y;
  1010. #ifdef Pack_32
  1011. ULong z;
  1012. #endif
  1013. #endif
  1014. i = cmp(a,b);
  1015. if (!i) {
  1016. c = Balloc(0);
  1017. c->wds = 1;
  1018. c->x[0] = 0;
  1019. return c;
  1020. }
  1021. if (i < 0) {
  1022. c = a;
  1023. a = b;
  1024. b = c;
  1025. i = 1;
  1026. }
  1027. else
  1028. i = 0;
  1029. c = Balloc(a->k);
  1030. c->sign = i;
  1031. wa = a->wds;
  1032. xa = a->x;
  1033. xae = xa + wa;
  1034. wb = b->wds;
  1035. xb = b->x;
  1036. xbe = xb + wb;
  1037. xc = c->x;
  1038. borrow = 0;
  1039. #ifdef ULLong
  1040. do {
  1041. y = (ULLong)*xa++ - *xb++ - borrow;
  1042. borrow = y >> 32 & (ULong)1;
  1043. *xc++ = y & FFFFFFFF;
  1044. }
  1045. while(xb < xbe);
  1046. while(xa < xae) {
  1047. y = *xa++ - borrow;
  1048. borrow = y >> 32 & (ULong)1;
  1049. *xc++ = y & FFFFFFFF;
  1050. }
  1051. #else
  1052. #ifdef Pack_32
  1053. do {
  1054. y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;
  1055. borrow = (y & 0x10000) >> 16;
  1056. z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;
  1057. borrow = (z & 0x10000) >> 16;
  1058. Storeinc(xc, z, y);
  1059. }
  1060. while(xb < xbe);
  1061. while(xa < xae) {
  1062. y = (*xa & 0xffff) - borrow;
  1063. borrow = (y & 0x10000) >> 16;
  1064. z = (*xa++ >> 16) - borrow;
  1065. borrow = (z & 0x10000) >> 16;
  1066. Storeinc(xc, z, y);
  1067. }
  1068. #else
  1069. do {
  1070. y = *xa++ - *xb++ - borrow;
  1071. borrow = (y & 0x10000) >> 16;
  1072. *xc++ = y & 0xffff;
  1073. }
  1074. while(xb < xbe);
  1075. while(xa < xae) {
  1076. y = *xa++ - borrow;
  1077. borrow = (y & 0x10000) >> 16;
  1078. *xc++ = y & 0xffff;
  1079. }
  1080. #endif
  1081. #endif
  1082. while(!*--xc)
  1083. wa--;
  1084. c->wds = wa;
  1085. return c;
  1086. }
  1087. static double
  1088. ulp
  1089. #ifdef KR_headers
  1090. (x) U *x;
  1091. #else
  1092. (U *x)
  1093. #endif
  1094. {
  1095. Long L;
  1096. U u;
  1097. L = (word0(x) & Exp_mask) - (P-1)*Exp_msk1;
  1098. #ifndef Avoid_Underflow
  1099. #ifndef Sudden_Underflow
  1100. if (L > 0) {
  1101. #endif
  1102. #endif
  1103. #ifdef IBM
  1104. L |= Exp_msk1 >> 4;
  1105. #endif
  1106. word0(&u) = L;
  1107. word1(&u) = 0;
  1108. #ifndef Avoid_Underflow
  1109. #ifndef Sudden_Underflow
  1110. }
  1111. else {
  1112. L = -L >> Exp_shift;
  1113. if (L < Exp_shift) {
  1114. word0(&u) = 0x80000 >> L;
  1115. word1(&u) = 0;
  1116. }
  1117. else {
  1118. word0(&u) = 0;
  1119. L -= Exp_shift;
  1120. word1(&u) = L >= 31 ? 1 : 1 << 31 - L;
  1121. }
  1122. }
  1123. #endif
  1124. #endif
  1125. return dval(&u);
  1126. }
  1127. static double
  1128. b2d
  1129. #ifdef KR_headers
  1130. (a, e) Bigint *a; int *e;
  1131. #else
  1132. (Bigint *a, int *e)
  1133. #endif
  1134. {
  1135. ULong *xa, *xa0, w, y, z;
  1136. int k;
  1137. U d;
  1138. #ifdef VAX
  1139. ULong d0, d1;
  1140. #else
  1141. #define d0 word0(&d)
  1142. #define d1 word1(&d)
  1143. #endif
  1144. xa0 = a->x;
  1145. xa = xa0 + a->wds;
  1146. y = *--xa;
  1147. #ifdef DEBUG
  1148. if (!y) Bug("zero y in b2d");
  1149. #endif
  1150. k = hi0bits(y);
  1151. *e = 32 - k;
  1152. #ifdef Pack_32
  1153. if (k < Ebits) {
  1154. d0 = Exp_1 | y >> (Ebits - k);
  1155. w = xa > xa0 ? *--xa : 0;
  1156. d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
  1157. goto ret_d;
  1158. }
  1159. z = xa > xa0 ? *--xa : 0;
  1160. if (k -= Ebits) {
  1161. d0 = Exp_1 | y << k | z >> (32 - k);
  1162. y = xa > xa0 ? *--xa : 0;
  1163. d1 = z << k | y >> (32 - k);
  1164. }
  1165. else {
  1166. d0 = Exp_1 | y;
  1167. d1 = z;
  1168. }
  1169. #else
  1170. if (k < Ebits + 16) {
  1171. z = xa > xa0 ? *--xa : 0;
  1172. d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
  1173. w = xa > xa0 ? *--xa : 0;
  1174. y = xa > xa0 ? *--xa : 0;
  1175. d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
  1176. goto ret_d;
  1177. }
  1178. z = xa > xa0 ? *--xa : 0;
  1179. w = xa > xa0 ? *--xa : 0;
  1180. k -= Ebits + 16;
  1181. d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
  1182. y = xa > xa0 ? *--xa : 0;
  1183. d1 = w << k + 16 | y << k;
  1184. #endif
  1185. ret_d:
  1186. #ifdef VAX
  1187. word0(&d) = d0 >> 16 | d0 << 16;
  1188. word1(&d) = d1 >> 16 | d1 << 16;
  1189. #else
  1190. #undef d0
  1191. #undef d1
  1192. #endif
  1193. return dval(&d);
  1194. }
  1195. static Bigint *
  1196. d2b
  1197. #ifdef KR_headers
  1198. (d, e, bits) U *d; int *e, *bits;
  1199. #else
  1200. (U *d, int *e, int *bits)
  1201. #endif
  1202. {
  1203. Bigint *b;
  1204. int de, k;
  1205. ULong *x, y, z;
  1206. #ifndef Sudden_Underflow
  1207. int i;
  1208. #endif
  1209. #ifdef VAX
  1210. ULong d0, d1;
  1211. d0 = word0(d) >> 16 | word0(d) << 16;
  1212. d1 = word1(d) >> 16 | word1(d) << 16;
  1213. #else
  1214. #define d0 word0(d)
  1215. #define d1 word1(d)
  1216. #endif
  1217. #ifdef Pack_32
  1218. b = Balloc(1);
  1219. #else
  1220. b = Balloc(2);
  1221. #endif
  1222. x = b->x;
  1223. z = d0 & Frac_mask;
  1224. d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
  1225. #ifdef Sudden_Underflow
  1226. de = (int)(d0 >> Exp_shift);
  1227. #ifndef IBM
  1228. z |= Exp_msk11;
  1229. #endif
  1230. #else
  1231. if ((de = (int)(d0 >> Exp_shift)))
  1232. z |= Exp_msk1;
  1233. #endif
  1234. #ifdef Pack_32
  1235. if ((y = d1)) {
  1236. if ((k = lo0bits(&y))) {
  1237. x[0] = y | z << (32 - k);
  1238. z >>= k;
  1239. }
  1240. else
  1241. x[0] = y;
  1242. #ifndef Sudden_Underflow
  1243. i =
  1244. #endif
  1245. b->wds = (x[1] = z) ? 2 : 1;
  1246. }
  1247. else {
  1248. k = lo0bits(&z);
  1249. x[0] = z;
  1250. #ifndef Sudden_Underflow
  1251. i =
  1252. #endif
  1253. b->wds = 1;
  1254. k += 32;
  1255. }
  1256. #else
  1257. if (y = d1) {
  1258. if (k = lo0bits(&y))
  1259. if (k >= 16) {
  1260. x[0] = y | z << 32 - k & 0xffff;
  1261. x[1] = z >> k - 16 & 0xffff;
  1262. x[2] = z >> k;
  1263. i = 2;
  1264. }
  1265. else {
  1266. x[0] = y & 0xffff;
  1267. x[1] = y >> 16 | z << 16 - k & 0xffff;
  1268. x[2] = z >> k & 0xffff;
  1269. x[3] = z >> k+16;
  1270. i = 3;
  1271. }
  1272. else {
  1273. x[0] = y & 0xffff;
  1274. x[1] = y >> 16;
  1275. x[2] = z & 0xffff;
  1276. x[3] = z >> 16;
  1277. i = 3;
  1278. }
  1279. }
  1280. else {
  1281. #ifdef DEBUG
  1282. if (!z)
  1283. Bug("Zero passed to d2b");
  1284. #endif
  1285. k = lo0bits(&z);
  1286. if (k >= 16) {
  1287. x[0] = z;
  1288. i = 0;
  1289. }
  1290. else {
  1291. x[0] = z & 0xffff;
  1292. x[1] = z >> 16;
  1293. i = 1;
  1294. }
  1295. k += 32;
  1296. }
  1297. while(!x[i])
  1298. --i;
  1299. b->wds = i + 1;
  1300. #endif
  1301. #ifndef Sudden_Underflow
  1302. if (de) {
  1303. #endif
  1304. #ifdef IBM
  1305. *e = (de - Bias - (P-1) << 2) + k;
  1306. *bits = 4*P + 8 - k - hi0bits(word0(d) & Frac_mask);
  1307. #else
  1308. *e = de - Bias - (P-1) + k;
  1309. *bits = P - k;
  1310. #endif
  1311. #ifndef Sudden_Underflow
  1312. }
  1313. else {
  1314. *e = de - Bias - (P-1) + 1 + k;
  1315. #ifdef Pack_32
  1316. *bits = 32*i - hi0bits(x[i-1]);
  1317. #else
  1318. *bits = (i+2)*16 - hi0bits(x[i]);
  1319. #endif
  1320. }
  1321. #endif
  1322. return b;
  1323. }
  1324. #undef d0
  1325. #undef d1
  1326. static double
  1327. ratio
  1328. #ifdef KR_headers
  1329. (a, b) Bigint *a, *b;
  1330. #else
  1331. (Bigint *a, Bigint *b)
  1332. #endif
  1333. {
  1334. U da, db;
  1335. int k, ka, kb;
  1336. dval(&da) = b2d(a, &ka);
  1337. dval(&db) = b2d(b, &kb);
  1338. #ifdef Pack_32
  1339. k = ka - kb + 32*(a->wds - b->wds);
  1340. #else
  1341. k = ka - kb + 16*(a->wds - b->wds);
  1342. #endif
  1343. #ifdef IBM
  1344. if (k > 0) {
  1345. word0(&da) += (k >> 2)*Exp_msk1;
  1346. if (k &= 3)
  1347. dval(&da) *= 1 << k;
  1348. }
  1349. else {
  1350. k = -k;
  1351. word0(&db) += (k >> 2)*Exp_msk1;
  1352. if (k &= 3)
  1353. dval(&db) *= 1 << k;
  1354. }
  1355. #else
  1356. if (k > 0)
  1357. word0(&da) += k*Exp_msk1;
  1358. else {
  1359. k = -k;
  1360. word0(&db) += k*Exp_msk1;
  1361. }
  1362. #endif
  1363. return dval(&da) / dval(&db);
  1364. }
  1365. static CONST double
  1366. tens[] = {
  1367. 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
  1368. 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
  1369. 1e20, 1e21, 1e22
  1370. #ifdef VAX
  1371. , 1e23, 1e24
  1372. #endif
  1373. };
  1374. static CONST double
  1375. #ifdef IEEE_Arith
  1376. bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
  1377. static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
  1378. #ifdef Avoid_Underflow
  1379. 9007199254740992.*9007199254740992.e-256
  1380. /* = 2^106 * 1e-256 */
  1381. #else
  1382. 1e-256
  1383. #endif
  1384. };
  1385. /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
  1386. /* flag unnecessarily. It leads to a song and dance at the end of strtod. */
  1387. #define Scale_Bit 0x10
  1388. #define n_bigtens 5
  1389. #else
  1390. #ifdef IBM
  1391. bigtens[] = { 1e16, 1e32, 1e64 };
  1392. static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64 };
  1393. #define n_bigtens 3
  1394. #else
  1395. bigtens[] = { 1e16, 1e32 };
  1396. static CONST double tinytens[] = { 1e-16, 1e-32 };
  1397. #define n_bigtens 2
  1398. #endif
  1399. #endif
  1400. #undef Need_Hexdig
  1401. #ifdef INFNAN_CHECK
  1402. #ifndef No_Hex_NaN
  1403. #define Need_Hexdig
  1404. #endif
  1405. #endif
  1406. #ifndef Need_Hexdig
  1407. #ifndef NO_HEX_FP
  1408. #define Need_Hexdig
  1409. #endif
  1410. #endif
  1411. #ifdef Need_Hexdig /*{*/
  1412. static unsigned char hexdig[256];
  1413. static void
  1414. #ifdef KR_headers
  1415. htinit(h, s, inc) unsigned char *h; unsigned char *s; int inc;
  1416. #else
  1417. htinit(unsigned char *h, unsigned char *s, int inc)
  1418. #endif
  1419. {
  1420. int i, j;
  1421. for(i = 0; (j = s[i]) !=0; i++)
  1422. h[j] = i + inc;
  1423. }
  1424. static void
  1425. #ifdef KR_headers
  1426. hexdig_init()
  1427. #else
  1428. hexdig_init(void)
  1429. #endif
  1430. {
  1431. #define USC (unsigned char *)
  1432. htinit(hexdig, USC "0123456789", 0x10);
  1433. htinit(hexdig, USC "abcdef", 0x10 + 10);
  1434. htinit(hexdig, USC "ABCDEF", 0x10 + 10);
  1435. }
  1436. #endif /* } Need_Hexdig */
  1437. #ifdef INFNAN_CHECK
  1438. #ifndef NAN_WORD0
  1439. #define NAN_WORD0 0x7ff80000
  1440. #endif
  1441. #ifndef NAN_WORD1
  1442. #define NAN_WORD1 0
  1443. #endif
  1444. static int
  1445. match
  1446. #ifdef KR_headers
  1447. (sp, t) char **sp, *t;
  1448. #else
  1449. (const char **sp, const char *t)
  1450. #endif
  1451. {
  1452. int c, d;
  1453. CONST char *s = *sp;
  1454. while((d = *t++)) {
  1455. if ((c = *++s) >= 'A' && c <= 'Z')
  1456. c += 'a' - 'A';
  1457. if (c != d)
  1458. return 0;
  1459. }
  1460. *sp = s + 1;
  1461. return 1;
  1462. }
  1463. #ifndef No_Hex_NaN
  1464. static void
  1465. hexnan
  1466. #ifdef KR_headers
  1467. (rvp, sp) U *rvp; CONST char **sp;
  1468. #else
  1469. (U *rvp, const char **sp)
  1470. #endif
  1471. {
  1472. ULong c, x[2];
  1473. CONST char *s;
  1474. int c1, havedig, udx0, xshift;
  1475. if (!hexdig['0'])
  1476. hexdig_init();
  1477. x[0] = x[1] = 0;
  1478. havedig = xshift = 0;
  1479. udx0 = 1;
  1480. s = *sp;
  1481. /* allow optional initial 0x or 0X */
  1482. while((c = *(CONST unsigned char*)(s+1)) && c <= ' ')
  1483. ++s;
  1484. if (s[1] == '0' && (s[2] == 'x' || s[2] == 'X'))
  1485. s += 2;
  1486. while((c = *(CONST unsigned char*)++s)) {
  1487. if ((c1 = hexdig[c]))
  1488. c = c1 & 0xf;
  1489. else if (c <= ' ') {
  1490. if (udx0 && havedig) {
  1491. udx0 = 0;
  1492. xshift = 1;
  1493. }
  1494. continue;
  1495. }
  1496. #ifdef GDTOA_NON_PEDANTIC_NANCHECK
  1497. else if (/*(*/ c == ')' && havedig) {
  1498. *sp = s + 1;
  1499. break;
  1500. }
  1501. else
  1502. return; /* invalid form: don't change *sp */
  1503. #else
  1504. else {
  1505. do {
  1506. if (/*(*/ c == ')') {
  1507. *sp = s + 1;
  1508. break;
  1509. }
  1510. } while((c = *++s));
  1511. break;
  1512. }
  1513. #endif
  1514. havedig = 1;
  1515. if (xshift) {
  1516. xshift = 0;
  1517. x[0] = x[1];
  1518. x[1] = 0;
  1519. }
  1520. if (udx0)
  1521. x[0] = (x[0] << 4) | (x[1] >> 28);
  1522. x[1] = (x[1] << 4) | c;
  1523. }
  1524. if ((x[0] &= 0xfffff) || x[1]) {
  1525. word0(rvp) = Exp_mask | x[0];
  1526. word1(rvp) = x[1];
  1527. }
  1528. }
  1529. #endif /*No_Hex_NaN*/
  1530. #endif /* INFNAN_CHECK */
  1531. #ifdef Pack_32
  1532. #define ULbits 32
  1533. #define kshift 5
  1534. #define kmask 31
  1535. #else
  1536. #define ULbits 16
  1537. #define kshift 4
  1538. #define kmask 15
  1539. #endif
  1540. #if !defined(NO_HEX_FP) || defined(Honor_FLT_ROUNDS) /*{*/
  1541. static Bigint *
  1542. #ifdef KR_headers
  1543. increment(b) Bigint *b;
  1544. #else
  1545. increment(Bigint *b)
  1546. #endif
  1547. {
  1548. ULong *x, *xe;
  1549. Bigint *b1;
  1550. x = b->x;
  1551. xe = x + b->wds;
  1552. do {
  1553. if (*x < (ULong)0xffffffffL) {
  1554. ++*x;
  1555. return b;
  1556. }
  1557. *x++ = 0;
  1558. } while(x < xe);
  1559. {
  1560. if (b->wds >= b->maxwds) {
  1561. b1 = Balloc(b->k+1);
  1562. Bcopy(b1,b);
  1563. Bfree(b);
  1564. b = b1;
  1565. }
  1566. b->x[b->wds++] = 1;
  1567. }
  1568. return b;
  1569. }
  1570. #endif /*}*/
  1571. #ifndef NO_HEX_FP /*{*/
  1572. static void
  1573. #ifdef KR_headers
  1574. rshift(b, k) Bigint *b; int k;
  1575. #else
  1576. rshift(Bigint *b, int k)
  1577. #endif
  1578. {
  1579. ULong *x, *x1, *xe, y;
  1580. int n;
  1581. x = x1 = b->x;
  1582. n = k >> kshift;
  1583. if (n < b->wds) {
  1584. xe = x + b->wds;
  1585. x += n;
  1586. if (k &= kmask) {
  1587. n = 32 - k;
  1588. y = *x++ >> k;
  1589. while(x < xe) {
  1590. *x1++ = (y | (*x << n)) & 0xffffffff;
  1591. y = *x++ >> k;
  1592. }
  1593. if ((*x1 = y) !=0)
  1594. x1++;
  1595. }
  1596. else
  1597. while(x < xe)
  1598. *x1++ = *x++;
  1599. }
  1600. if ((b->wds = x1 - b->x) == 0)
  1601. b->x[0] = 0;
  1602. }
  1603. static ULong
  1604. #ifdef KR_headers
  1605. any_on(b, k) Bigint *b; int k;
  1606. #else
  1607. any_on(Bigint *b, int k)
  1608. #endif
  1609. {
  1610. int n, nwds;
  1611. ULong *x, *x0, x1, x2;
  1612. x = b->x;
  1613. nwds = b->wds;
  1614. n = k >> kshift;
  1615. if (n > nwds)
  1616. n = nwds;
  1617. else if (n < nwds && (k &= kmask)) {
  1618. x1 = x2 = x[n];
  1619. x1 >>= k;
  1620. x1 <<= k;
  1621. if (x1 != x2)
  1622. return 1;
  1623. }
  1624. x0 = x;
  1625. x += n;
  1626. while(x > x0)
  1627. if (*--x)
  1628. return 1;
  1629. return 0;
  1630. }
  1631. enum { /* rounding values: same as FLT_ROUNDS */
  1632. Round_zero = 0,
  1633. Round_near = 1,
  1634. Round_up = 2,
  1635. Round_down = 3
  1636. };
  1637. void
  1638. #ifdef KR_headers
  1639. gethex(sp, rvp, rounding, sign)
  1640. CONST char **sp; U *rvp; int rounding, sign;
  1641. #else
  1642. gethex( CONST char **sp, U *rvp, int rounding, int sign)
  1643. #endif
  1644. {
  1645. Bigint *b;
  1646. CONST unsigned char *decpt, *s0, *s, *s1;
  1647. Long e, e1;
  1648. ULong L, lostbits, *x;
  1649. int big, denorm, esign, havedig, k, n, nbits, up, zret;
  1650. #ifdef IBM
  1651. int j;
  1652. #endif
  1653. enum {
  1654. #ifdef IEEE_Arith /*{{*/
  1655. emax = 0x7fe - Bias - P + 1,
  1656. emin = Emin - P + 1
  1657. #else /*}{*/
  1658. emin = Emin - P,
  1659. #ifdef VAX
  1660. emax = 0x7ff - Bias - P + 1
  1661. #endif
  1662. #ifdef IBM
  1663. emax = 0x7f - Bias - P
  1664. #endif
  1665. #endif /*}}*/
  1666. };
  1667. #ifdef USE_LOCALE
  1668. int i;
  1669. #ifdef NO_LOCALE_CACHE
  1670. const unsigned char *decimalpoint = (unsigned char*)
  1671. localeconv()->decimal_point;
  1672. #else
  1673. const unsigned char *decimalpoint;
  1674. static unsigned char *decimalpoint_cache;
  1675. if (!(s0 = decimalpoint_cache)) {
  1676. s0 = (unsigned char*)localeconv()->decimal_point;
  1677. if ((decimalpoint_cache = (unsigned char*)
  1678. MALLOC(strlen((CONST char*)s0) + 1))) {
  1679. strcpy((char*)decimalpoint_cache, (CONST char*)s0);
  1680. s0 = decimalpoint_cache;
  1681. }
  1682. }
  1683. decimalpoint = s0;
  1684. #endif
  1685. #endif
  1686. if (!hexdig['0'])
  1687. hexdig_init();
  1688. havedig = 0;
  1689. s0 = *(CONST unsigned char **)sp + 2;
  1690. while(s0[havedig] == '0')
  1691. havedig++;
  1692. s0 += havedig;
  1693. s = s0;
  1694. decpt = 0;
  1695. zret = 0;
  1696. e = 0;
  1697. if (hexdig[*s])
  1698. havedig++;
  1699. else {
  1700. zret = 1;
  1701. #ifdef USE_LOCALE
  1702. for(i = 0; decimalpoint[i]; ++i) {
  1703. if (s[i] != decimalpoint[i])
  1704. goto pcheck;
  1705. }
  1706. decpt = s += i;
  1707. #else
  1708. if (*s != '.')
  1709. goto pcheck;
  1710. decpt = ++s;
  1711. #endif
  1712. if (!hexdig[*s])
  1713. goto pcheck;
  1714. while(*s == '0')
  1715. s++;
  1716. if (hexdig[*s])
  1717. zret = 0;
  1718. havedig = 1;
  1719. s0 = s;
  1720. }
  1721. while(hexdig[*s])
  1722. s++;
  1723. #ifdef USE_LOCALE
  1724. if (*s == *decimalpoint && !decpt) {
  1725. for(i = 1; decimalpoint[i]; ++i) {
  1726. if (s[i] != decimalpoint[i])
  1727. goto pcheck;
  1728. }
  1729. decpt = s += i;
  1730. #else
  1731. if (*s == '.' && !decpt) {
  1732. decpt = ++s;
  1733. #endif
  1734. while(hexdig[*s])
  1735. s++;
  1736. }/*}*/
  1737. if (decpt)
  1738. e = -(((Long)(s-decpt)) << 2);
  1739. pcheck:
  1740. s1 = s;
  1741. big = esign = 0;
  1742. switch(*s) {
  1743. case 'p':
  1744. case 'P':
  1745. switch(*++s) {
  1746. case '-':
  1747. esign = 1;
  1748. /* no break */
  1749. case '+':
  1750. s++;
  1751. }
  1752. if ((n = hexdig[*s]) == 0 || n > 0x19) {
  1753. s = s1;
  1754. break;
  1755. }
  1756. e1 = n - 0x10;
  1757. while((n = hexdig[*++s]) !=0 && n <= 0x19) {
  1758. if (e1 & 0xf8000000)
  1759. big = 1;
  1760. e1 = 10*e1 + n - 0x10;
  1761. }
  1762. if (esign)
  1763. e1 = -e1;
  1764. e += e1;
  1765. }
  1766. *sp = (char*)s;
  1767. if (!havedig)
  1768. *sp = (char*)s0 - 1;
  1769. if (zret)
  1770. goto retz1;
  1771. if (big) {
  1772. if (esign) {
  1773. #ifdef IEEE_Arith
  1774. switch(rounding) {
  1775. case Round_up:
  1776. if (sign)
  1777. break;
  1778. goto ret_tiny;
  1779. case Round_down:
  1780. if (!sign)
  1781. break;
  1782. goto ret_tiny;
  1783. }
  1784. #endif
  1785. goto retz;
  1786. #ifdef IEEE_Arith
  1787. ret_tiny:
  1788. #ifndef NO_ERRNO
  1789. errno = ERANGE;
  1790. #endif
  1791. word0(rvp) = 0;
  1792. word1(rvp) = 1;
  1793. return;
  1794. #endif /* IEEE_Arith */
  1795. }
  1796. switch(rounding) {
  1797. case Round_near:
  1798. goto ovfl1;
  1799. case Round_up:
  1800. if (!sign)
  1801. goto ovfl1;
  1802. goto ret_big;
  1803. case Round_down:
  1804. if (sign)
  1805. goto ovfl1;
  1806. goto ret_big;
  1807. }
  1808. ret_big:
  1809. word0(rvp) = Big0;
  1810. word1(rvp) = Big1;
  1811. return;
  1812. }
  1813. n = s1 - s0 - 1;
  1814. for(k = 0; n > (1 << (kshift-2)) - 1; n >>= 1)
  1815. k++;
  1816. b = Balloc(k);
  1817. x = b->x;
  1818. n = 0;
  1819. L = 0;
  1820. #ifdef USE_LOCALE
  1821. for(i = 0; decimalpoint[i+1]; ++i);
  1822. #endif
  1823. while(s1 > s0) {
  1824. #ifdef USE_LOCALE
  1825. if (*--s1 == decimalpoint[i]) {
  1826. s1 -= i;
  1827. continue;
  1828. }
  1829. #else
  1830. if (*--s1 == '.')
  1831. continue;
  1832. #endif
  1833. if (n == ULbits) {
  1834. *x++ = L;
  1835. L = 0;
  1836. n = 0;
  1837. }
  1838. L |= (hexdig[*s1] & 0x0f) << n;
  1839. n += 4;
  1840. }
  1841. *x++ = L;
  1842. b->wds = n = x - b->x;
  1843. n = ULbits*n - hi0bits(L);
  1844. nbits = Nbits;
  1845. lostbits = 0;
  1846. x = b->x;
  1847. if (n > nbits) {
  1848. n -= nbits;
  1849. if (any_on(b,n)) {
  1850. lostbits = 1;
  1851. k = n - 1;
  1852. if (x[k>>kshift] & 1 << (k & kmask)) {
  1853. lostbits = 2;
  1854. if (k > 0 && any_on(b,k))
  1855. lostbits = 3;
  1856. }
  1857. }
  1858. rshift(b, n);
  1859. e += n;
  1860. }
  1861. else if (n < nbits) {
  1862. n = nbits - n;
  1863. b = lshift(b, n);
  1864. e -= n;
  1865. x = b->x;
  1866. }
  1867. if (e > Emax) {
  1868. ovfl:
  1869. Bfree(b);
  1870. ovfl1:
  1871. #ifndef NO_ERRNO
  1872. errno = ERANGE;
  1873. #endif
  1874. word0(rvp) = Exp_mask;
  1875. word1(rvp) = 0;
  1876. return;
  1877. }
  1878. denorm = 0;
  1879. if (e < emin) {
  1880. denorm = 1;
  1881. n = emin - e;
  1882. if (n >= nbits) {
  1883. #ifdef IEEE_Arith /*{*/
  1884. switch (rounding) {
  1885. case Round_near:
  1886. if (n == nbits && (n < 2 || any_on(b,n-1)))
  1887. goto ret_tiny;
  1888. break;
  1889. case Round_up:
  1890. if (!sign)
  1891. goto ret_tiny;
  1892. break;
  1893. case Round_down:
  1894. if (sign)
  1895. goto ret_tiny;
  1896. }
  1897. #endif /* } IEEE_Arith */
  1898. Bfree(b);
  1899. retz:
  1900. #ifndef NO_ERRNO
  1901. errno = ERANGE;
  1902. #endif
  1903. retz1:
  1904. rvp->d = 0.;
  1905. return;
  1906. }
  1907. k = n - 1;
  1908. if (lostbits)
  1909. lostbits = 1;
  1910. else if (k > 0)
  1911. lostbits = any_on(b,k);
  1912. if (x[k>>kshift] & 1 << (k & kmask))
  1913. lostbits |= 2;
  1914. nbits -= n;
  1915. rshift(b,n);
  1916. e = emin;
  1917. }
  1918. if (lostbits) {
  1919. up = 0;
  1920. switch(rounding) {
  1921. case Round_zero:
  1922. break;
  1923. case Round_near:
  1924. if (lostbits & 2
  1925. && (lostbits & 1) | (x[0] & 1))
  1926. up = 1;
  1927. break;
  1928. case Round_up:
  1929. up = 1 - sign;
  1930. break;
  1931. case Round_down:
  1932. up = sign;
  1933. }
  1934. if (up) {
  1935. k = b->wds;
  1936. b = increment(b);
  1937. x = b->x;
  1938. if (denorm) {
  1939. #if 0
  1940. if (nbits == Nbits - 1
  1941. && x[nbits >> kshift] & 1 << (nbits & kmask))
  1942. denorm = 0; /* not currently used */
  1943. #endif
  1944. }
  1945. else if (b->wds > k
  1946. || ((n = nbits & kmask) !=0
  1947. && hi0bits(x[k-1]) < 32-n)) {
  1948. rshift(b,1);
  1949. if (++e > Emax)
  1950. goto ovfl;
  1951. }
  1952. }
  1953. }
  1954. #ifdef IEEE_Arith
  1955. if (denorm)
  1956. word0(rvp) = b->wds > 1 ? b->x[1] & ~0x100000 : 0;
  1957. else
  1958. word0(rvp) = (b->x[1] & ~0x100000) | ((e + 0x3ff + 52) << 20);
  1959. word1(rvp) = b->x[0];
  1960. #endif
  1961. #ifdef IBM
  1962. if ((j = e & 3)) {
  1963. k = b->x[0] & ((1 << j) - 1);
  1964. rshift(b,j);
  1965. if (k) {
  1966. switch(rounding) {
  1967. case Round_up:
  1968. if (!sign)
  1969. increment(b);
  1970. break;
  1971. case Round_down:
  1972. if (sign)
  1973. increment(b);
  1974. break;
  1975. case Round_near:
  1976. j = 1 << (j-1);
  1977. if (k & j && ((k & (j-1)) | lostbits))
  1978. increment(b);
  1979. }
  1980. }
  1981. }
  1982. e >>= 2;
  1983. word0(rvp) = b->x[1] | ((e + 65 + 13) << 24);
  1984. word1(rvp) = b->x[0];
  1985. #endif
  1986. #ifdef VAX
  1987. /* The next two lines ignore swap of low- and high-order 2 bytes. */
  1988. /* word0(rvp) = (b->x[1] & ~0x800000) | ((e + 129 + 55) << 23); */
  1989. /* word1(rvp) = b->x[0]; */
  1990. word0(rvp) = ((b->x[1] & ~0x800000) >> 16) | ((e + 129 + 55) << 7) | (b->x[1] << 16);
  1991. word1(rvp) = (b->x[0] >> 16) | (b->x[0] << 16);
  1992. #endif
  1993. Bfree(b);
  1994. }
  1995. #endif /*!NO_HEX_FP}*/
  1996. static int
  1997. #ifdef KR_headers
  1998. dshift(b, p2) Bigint *b; int p2;
  1999. #else
  2000. dshift(Bigint *b, int p2)
  2001. #endif
  2002. {
  2003. int rv = hi0bits(b->x[b->wds-1]) - 4;
  2004. if (p2 > 0)
  2005. rv -= p2;
  2006. return rv & kmask;
  2007. }
  2008. static int
  2009. quorem
  2010. #ifdef KR_headers
  2011. (b, S) Bigint *b, *S;
  2012. #else
  2013. (Bigint *b, Bigint *S)
  2014. #endif
  2015. {
  2016. int n;
  2017. ULong *bx, *bxe, q, *sx, *sxe;
  2018. #ifdef ULLong
  2019. ULLong borrow, carry, y, ys;
  2020. #else
  2021. ULong borrow, carry, y, ys;
  2022. #ifdef Pack_32
  2023. ULong si, z, zs;
  2024. #endif
  2025. #endif
  2026. n = S->wds;
  2027. #ifdef DEBUG
  2028. /*debug*/ if (b->wds > n)
  2029. /*debug*/ Bug("oversize b in quorem");
  2030. #endif
  2031. if (b->wds < n)
  2032. return 0;
  2033. sx = S->x;
  2034. sxe = sx + --n;
  2035. bx = b->x;
  2036. bxe = bx + n;
  2037. q = *bxe / (*sxe + 1); /* ensure q <= true quotient */
  2038. #ifdef DEBUG
  2039. #ifdef NO_STRTOD_BIGCOMP
  2040. /*debug*/ if (q > 9)
  2041. #else
  2042. /* An oversized q is possible when quorem is called from bigcomp and */
  2043. /* the input is near, e.g., twice the smallest denormalized number. */
  2044. /*debug*/ if (q > 15)
  2045. #endif
  2046. /*debug*/ Bug("oversized quotient in quorem");
  2047. #endif
  2048. if (q) {
  2049. borrow = 0;
  2050. carry = 0;
  2051. do {
  2052. #ifdef ULLong
  2053. ys = *sx++ * (ULLong)q + carry;
  2054. carry = ys >> 32;
  2055. y = *bx - (ys & FFFFFFFF) - borrow;
  2056. borrow = y >> 32 & (ULong)1;
  2057. *bx++ = y & FFFFFFFF;
  2058. #else
  2059. #ifdef Pack_32
  2060. si = *sx++;
  2061. ys = (si & 0xffff) * q + carry;
  2062. zs = (si >> 16) * q + (ys >> 16);
  2063. carry = zs >> 16;
  2064. y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
  2065. borrow = (y & 0x10000) >> 16;
  2066. z = (*bx >> 16) - (zs & 0xffff) - borrow;
  2067. borrow = (z & 0x10000) >> 16;
  2068. Storeinc(bx, z, y);
  2069. #else
  2070. ys = *sx++ * q + carry;
  2071. carry = ys >> 16;
  2072. y = *bx - (ys & 0xffff) - borrow;
  2073. borrow = (y & 0x10000) >> 16;
  2074. *bx++ = y & 0xffff;
  2075. #endif
  2076. #endif
  2077. }
  2078. while(sx <= sxe);
  2079. if (!*bxe) {
  2080. bx = b->x;
  2081. while(--bxe > bx && !*bxe)
  2082. --n;
  2083. b->wds = n;
  2084. }
  2085. }
  2086. if (cmp(b, S) >= 0) {
  2087. q++;
  2088. borrow = 0;
  2089. carry = 0;
  2090. bx = b->x;
  2091. sx = S->x;
  2092. do {
  2093. #ifdef ULLong
  2094. ys = *sx++ + carry;
  2095. carry = ys >> 32;
  2096. y = *bx - (ys & FFFFFFFF) - borrow;
  2097. borrow = y >> 32 & (ULong)1;
  2098. *bx++ = y & FFFFFFFF;
  2099. #else
  2100. #ifdef Pack_32
  2101. si = *sx++;
  2102. ys = (si & 0xffff) + carry;
  2103. zs = (si >> 16) + (ys >> 16);
  2104. carry = zs >> 16;
  2105. y = (*bx & 0xffff) - (ys & 0xffff) - borrow;
  2106. borrow = (y & 0x10000) >> 16;
  2107. z = (*bx >> 16) - (zs & 0xffff) - borrow;
  2108. borrow = (z & 0x10000) >> 16;
  2109. Storeinc(bx, z, y);
  2110. #else
  2111. ys = *sx++ + carry;
  2112. carry = ys >> 16;
  2113. y = *bx - (ys & 0xffff) - borrow;
  2114. borrow = (y & 0x10000) >> 16;
  2115. *bx++ = y & 0xffff;
  2116. #endif
  2117. #endif
  2118. }
  2119. while(sx <= sxe);
  2120. bx = b->x;
  2121. bxe = bx + n;
  2122. if (!*bxe) {
  2123. while(--bxe > bx && !*bxe)
  2124. --n;
  2125. b->wds = n;
  2126. }
  2127. }
  2128. return q;
  2129. }
  2130. #if defined(Avoid_Underflow) || !defined(NO_STRTOD_BIGCOMP) /*{*/
  2131. static double
  2132. sulp
  2133. #ifdef KR_headers
  2134. (x, bc) U *x; BCinfo *bc;
  2135. #else
  2136. (U *x, BCinfo *bc)
  2137. #endif
  2138. {
  2139. U u;
  2140. double rv;
  2141. int i;
  2142. rv = ulp(x);
  2143. if (!bc->scale || (i = 2*P + 1 - ((word0(x) & Exp_mask) >> Exp_shift)) <= 0)
  2144. return rv; /* Is there an example where i <= 0 ? */
  2145. word0(&u) = Exp_1 + (i << Exp_shift);
  2146. word1(&u) = 0;
  2147. return rv * u.d;
  2148. }
  2149. #endif /*}*/
  2150. #ifndef NO_STRTOD_BIGCOMP
  2151. static void
  2152. bigcomp
  2153. #ifdef KR_headers
  2154. (rv, s0, bc)
  2155. U *rv; CONST char *s0; BCinfo *bc;
  2156. #else
  2157. (U *rv, const char *s0, BCinfo *bc)
  2158. #endif
  2159. {
  2160. Bigint *b, *d;
  2161. int b2, bbits, d2, dd, dig, dsign, i, j, nd, nd0, p2, p5, speccase;
  2162. dsign = bc->dsign;
  2163. nd = bc->nd;
  2164. nd0 = bc->nd0;
  2165. p5 = nd + bc->e0 - 1;
  2166. speccase = 0;
  2167. #ifndef Sudden_Underflow
  2168. if (rv->d == 0.) { /* special case: value near underflow-to-zero */
  2169. /* threshold was rounded to zero */
  2170. b = i2b(1);
  2171. p2 = Emin - P + 1;
  2172. bbits = 1;
  2173. #ifdef Avoid_Underflow
  2174. word0(rv) = (P+2) << Exp_shift;
  2175. #else
  2176. word1(rv) = 1;
  2177. #endif
  2178. i = 0;
  2179. #ifdef Honor_FLT_ROUNDS
  2180. if (bc->rounding == 1)
  2181. #endif
  2182. {
  2183. speccase = 1;
  2184. --p2;
  2185. dsign = 0;
  2186. goto have_i;
  2187. }
  2188. }
  2189. else
  2190. #endif
  2191. b = d2b(rv, &p2, &bbits);
  2192. #ifdef Avoid_Underflow
  2193. p2 -= bc->scale;
  2194. #endif
  2195. /* floor(log2(rv)) == bbits - 1 + p2 */
  2196. /* Check for denormal case. */
  2197. i = P - bbits;
  2198. if (i > (j = P - Emin - 1 + p2)) {
  2199. #ifdef Sudden_Underflow
  2200. Bfree(b);
  2201. b = i2b(1);
  2202. p2 = Emin;
  2203. i = P - 1;
  2204. #ifdef Avoid_Underflow
  2205. word0(rv) = (1 + bc->scale) << Exp_shift;
  2206. #else
  2207. word0(rv) = Exp_msk1;
  2208. #endif
  2209. word1(rv) = 0;
  2210. #else
  2211. i = j;
  2212. #endif
  2213. }
  2214. #ifdef Honor_FLT_ROUNDS
  2215. if (bc->rounding != 1) {
  2216. if (i > 0)
  2217. b = lshift(b, i);
  2218. if (dsign)
  2219. b = increment(b);
  2220. }
  2221. else
  2222. #endif
  2223. {
  2224. b = lshift(b, ++i);
  2225. b->x[0] |= 1;
  2226. }
  2227. #ifndef Sudden_Underflow
  2228. have_i:
  2229. #endif
  2230. p2 -= p5 + i;
  2231. d = i2b(1);
  2232. /* Arrange for convenient computation of quotients:
  2233. * shift left if necessary so divisor has 4 leading 0 bits.
  2234. */
  2235. if (p5 > 0)
  2236. d = pow5mult(d, p5);
  2237. else if (p5 < 0)
  2238. b = pow5mult(b, -p5);
  2239. if (p2 > 0) {
  2240. b2 = p2;
  2241. d2 = 0;
  2242. }
  2243. else {
  2244. b2 = 0;
  2245. d2 = -p2;
  2246. }
  2247. i = dshift(d, d2);
  2248. if ((b2 += i) > 0)
  2249. b = lshift(b, b2);
  2250. if ((d2 += i) > 0)
  2251. d = lshift(d, d2);
  2252. /* Now b/d = exactly half-way between the two floating-point values */
  2253. /* on either side of the input string. Compute first digit of b/d. */
  2254. if (!(dig = quorem(b,d))) {
  2255. b = multadd(b, 10, 0); /* very unlikely */
  2256. dig = quorem(b,d);
  2257. }
  2258. /* Compare b/d with s0 */
  2259. for(i = 0; i < nd0; ) {
  2260. if ((dd = s0[i++] - '0' - dig))
  2261. goto ret;
  2262. if (!b->x[0] && b->wds == 1) {
  2263. if (i < nd)
  2264. dd = 1;
  2265. goto ret;
  2266. }
  2267. b = multadd(b, 10, 0);
  2268. dig = quorem(b,d);
  2269. }
  2270. for(j = bc->dp1; i++ < nd;) {
  2271. if ((dd = s0[j++] - '0' - dig))
  2272. goto ret;
  2273. if (!b->x[0] && b->wds == 1) {
  2274. if (i < nd)
  2275. dd = 1;
  2276. goto ret;
  2277. }
  2278. b = multadd(b, 10, 0);
  2279. dig = quorem(b,d);
  2280. }
  2281. if (b->x[0] || b->wds > 1)
  2282. dd = -1;
  2283. ret:
  2284. Bfree(b);
  2285. Bfree(d);
  2286. #ifdef Honor_FLT_ROUNDS
  2287. if (bc->rounding != 1) {
  2288. if (dd < 0) {
  2289. if (bc->rounding == 0) {
  2290. if (!dsign)
  2291. goto retlow1;
  2292. }
  2293. else if (dsign)
  2294. goto rethi1;
  2295. }
  2296. else if (dd > 0) {
  2297. if (bc->rounding == 0) {
  2298. if (dsign)
  2299. goto rethi1;
  2300. goto ret1;
  2301. }
  2302. if (!dsign)
  2303. goto rethi1;
  2304. dval(rv) += 2.*sulp(rv,bc);
  2305. }
  2306. else {
  2307. bc->inexact = 0;
  2308. if (dsign)
  2309. goto rethi1;
  2310. }
  2311. }
  2312. else
  2313. #endif
  2314. if (speccase) {
  2315. if (dd <= 0)
  2316. rv->d = 0.;
  2317. }
  2318. else if (dd < 0) {
  2319. if (!dsign) /* does not happen for round-near */
  2320. retlow1:
  2321. dval(rv) -= sulp(rv,bc);
  2322. }
  2323. else if (dd > 0) {
  2324. if (dsign) {
  2325. rethi1:
  2326. dval(rv) += sulp(rv,bc);
  2327. }
  2328. }
  2329. else {
  2330. /* Exact half-way case: apply round-even rule. */
  2331. if ((j = ((word0(rv) & Exp_mask) >> Exp_shift) - bc->scale) <= 0) {
  2332. i = 1 - j;
  2333. if (i <= 31) {
  2334. if (word1(rv) & (0x1 << i))
  2335. goto odd;
  2336. }
  2337. else if (word0(rv) & (0x1 << (i-32)))
  2338. goto odd;
  2339. }
  2340. else if (word1(rv) & 1) {
  2341. odd:
  2342. if (dsign)
  2343. goto rethi1;
  2344. goto retlow1;
  2345. }
  2346. }
  2347. #ifdef Honor_FLT_ROUNDS
  2348. ret1:
  2349. #endif
  2350. return;
  2351. }
  2352. #endif /* NO_STRTOD_BIGCOMP */
  2353. double
  2354. fpconv_strtod
  2355. #ifdef KR_headers
  2356. (s00, se) CONST char *s00; char **se;
  2357. #else
  2358. (const char *s00, char **se)
  2359. #endif
  2360. {
  2361. int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, e, e1;
  2362. int esign, i, j, k, nd, nd0, nf, nz, nz0, nz1, sign;
  2363. CONST char *s, *s0, *s1;
  2364. double aadj, aadj1;
  2365. Long L;
  2366. U aadj2, adj, rv, rv0;
  2367. ULong y, z;
  2368. BCinfo bc;
  2369. Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
  2370. #ifdef Avoid_Underflow
  2371. ULong Lsb, Lsb1;
  2372. #endif
  2373. #ifdef SET_INEXACT
  2374. int oldinexact;
  2375. #endif
  2376. #ifndef NO_STRTOD_BIGCOMP
  2377. int req_bigcomp = 0;
  2378. #endif
  2379. #ifdef Honor_FLT_ROUNDS /*{*/
  2380. #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
  2381. bc.rounding = Flt_Rounds;
  2382. #else /*}{*/
  2383. bc.rounding = 1;
  2384. switch(fegetround()) {
  2385. case FE_TOWARDZERO: bc.rounding = 0; break;
  2386. case FE_UPWARD: bc.rounding = 2; break;
  2387. case FE_DOWNWARD: bc.rounding = 3;
  2388. }
  2389. #endif /*}}*/
  2390. #endif /*}*/
  2391. #ifdef USE_LOCALE
  2392. CONST char *s2;
  2393. #endif
  2394. sign = nz0 = nz1 = nz = bc.dplen = bc.uflchk = 0;
  2395. dval(&rv) = 0.;
  2396. for(s = s00;;s++) switch(*s) {
  2397. case '-':
  2398. sign = 1;
  2399. /* no break */
  2400. case '+':
  2401. if (*++s)
  2402. goto break2;
  2403. /* no break */
  2404. case 0:
  2405. goto ret0;
  2406. case '\t':
  2407. case '\n':
  2408. case '\v':
  2409. case '\f':
  2410. case '\r':
  2411. case ' ':
  2412. continue;
  2413. default:
  2414. goto break2;
  2415. }
  2416. break2:
  2417. if (*s == '0') {
  2418. #ifndef NO_HEX_FP /*{*/
  2419. switch(s[1]) {
  2420. case 'x':
  2421. case 'X':
  2422. #ifdef Honor_FLT_ROUNDS
  2423. gethex(&s, &rv, bc.rounding, sign);
  2424. #else
  2425. gethex(&s, &rv, 1, sign);
  2426. #endif
  2427. goto ret;
  2428. }
  2429. #endif /*}*/
  2430. nz0 = 1;
  2431. while(*++s == '0') ;
  2432. if (!*s)
  2433. goto ret;
  2434. }
  2435. s0 = s;
  2436. y = z = 0;
  2437. for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
  2438. if (nd < 9)
  2439. y = 10*y + c - '0';
  2440. else if (nd < 16)
  2441. z = 10*z + c - '0';
  2442. nd0 = nd;
  2443. bc.dp0 = bc.dp1 = s - s0;
  2444. for(s1 = s; s1 > s0 && *--s1 == '0'; )
  2445. ++nz1;
  2446. #ifdef USE_LOCALE
  2447. s1 = localeconv()->decimal_point;
  2448. if (c == *s1) {
  2449. c = '.';
  2450. if (*++s1) {
  2451. s2 = s;
  2452. for(;;) {
  2453. if (*++s2 != *s1) {
  2454. c = 0;
  2455. break;
  2456. }
  2457. if (!*++s1) {
  2458. s = s2;
  2459. break;
  2460. }
  2461. }
  2462. }
  2463. }
  2464. #endif
  2465. if (c == '.') {
  2466. c = *++s;
  2467. bc.dp1 = s - s0;
  2468. bc.dplen = bc.dp1 - bc.dp0;
  2469. if (!nd) {
  2470. for(; c == '0'; c = *++s)
  2471. nz++;
  2472. if (c > '0' && c <= '9') {
  2473. bc.dp0 = s0 - s;
  2474. bc.dp1 = bc.dp0 + bc.dplen;
  2475. s0 = s;
  2476. nf += nz;
  2477. nz = 0;
  2478. goto have_dig;
  2479. }
  2480. goto dig_done;
  2481. }
  2482. for(; c >= '0' && c <= '9'; c = *++s) {
  2483. have_dig:
  2484. nz++;
  2485. if (c -= '0') {
  2486. nf += nz;
  2487. for(i = 1; i < nz; i++)
  2488. if (nd++ < 9)
  2489. y *= 10;
  2490. else if (nd <= DBL_DIG + 1)
  2491. z *= 10;
  2492. if (nd++ < 9)
  2493. y = 10*y + c;
  2494. else if (nd <= DBL_DIG + 1)
  2495. z = 10*z + c;
  2496. nz = nz1 = 0;
  2497. }
  2498. }
  2499. }
  2500. dig_done:
  2501. e = 0;
  2502. if (c == 'e' || c == 'E') {
  2503. if (!nd && !nz && !nz0) {
  2504. goto ret0;
  2505. }
  2506. s00 = s;
  2507. esign = 0;
  2508. switch(c = *++s) {
  2509. case '-':
  2510. esign = 1;
  2511. case '+':
  2512. c = *++s;
  2513. }
  2514. if (c >= '0' && c <= '9') {
  2515. while(c == '0')
  2516. c = *++s;
  2517. if (c > '0' && c <= '9') {
  2518. L = c - '0';
  2519. s1 = s;
  2520. while((c = *++s) >= '0' && c <= '9')
  2521. L = 10*L + c - '0';
  2522. if (s - s1 > 8 || L > 19999)
  2523. /* Avoid confusion from exponents
  2524. * so large that e might overflow.
  2525. */
  2526. e = 19999; /* safe for 16 bit ints */
  2527. else
  2528. e = (int)L;
  2529. if (esign)
  2530. e = -e;
  2531. }
  2532. else
  2533. e = 0;
  2534. }
  2535. else
  2536. s = s00;
  2537. }
  2538. if (!nd) {
  2539. if (!nz && !nz0) {
  2540. #ifdef INFNAN_CHECK
  2541. /* Check for Nan and Infinity */
  2542. if (!bc.dplen)
  2543. switch(c) {
  2544. case 'i':
  2545. case 'I':
  2546. if (match(&s,"nf")) {
  2547. --s;
  2548. if (!match(&s,"inity"))
  2549. ++s;
  2550. word0(&rv) = 0x7ff00000;
  2551. word1(&rv) = 0;
  2552. goto ret;
  2553. }
  2554. break;
  2555. case 'n':
  2556. case 'N':
  2557. if (match(&s, "an")) {
  2558. word0(&rv) = NAN_WORD0;
  2559. word1(&rv) = NAN_WORD1;
  2560. #ifndef No_Hex_NaN
  2561. if (*s == '(') /*)*/
  2562. hexnan(&rv, &s);
  2563. #endif
  2564. goto ret;
  2565. }
  2566. }
  2567. #endif /* INFNAN_CHECK */
  2568. ret0:
  2569. s = s00;
  2570. sign = 0;
  2571. }
  2572. goto ret;
  2573. }
  2574. bc.e0 = e1 = e -= nf;
  2575. /* Now we have nd0 digits, starting at s0, followed by a
  2576. * decimal point, followed by nd-nd0 digits. The number we're
  2577. * after is the integer represented by those digits times
  2578. * 10**e */
  2579. if (!nd0)
  2580. nd0 = nd;
  2581. k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
  2582. dval(&rv) = y;
  2583. if (k > 9) {
  2584. #ifdef SET_INEXACT
  2585. if (k > DBL_DIG)
  2586. oldinexact = get_inexact();
  2587. #endif
  2588. dval(&rv) = tens[k - 9] * dval(&rv) + z;
  2589. }
  2590. bd0 = 0;
  2591. if (nd <= DBL_DIG
  2592. #ifndef RND_PRODQUOT
  2593. #ifndef Honor_FLT_ROUNDS
  2594. && Flt_Rounds == 1
  2595. #endif
  2596. #endif
  2597. ) {
  2598. if (!e)
  2599. goto ret;
  2600. #ifndef ROUND_BIASED_without_Round_Up
  2601. if (e > 0) {
  2602. if (e <= Ten_pmax) {
  2603. #ifdef VAX
  2604. goto vax_ovfl_check;
  2605. #else
  2606. #ifdef Honor_FLT_ROUNDS
  2607. /* round correctly FLT_ROUNDS = 2 or 3 */
  2608. if (sign) {
  2609. rv.d = -rv.d;
  2610. sign = 0;
  2611. }
  2612. #endif
  2613. /* rv = */ rounded_product(dval(&rv), tens[e]);
  2614. goto ret;
  2615. #endif
  2616. }
  2617. i = DBL_DIG - nd;
  2618. if (e <= Ten_pmax + i) {
  2619. /* A fancier test would sometimes let us do
  2620. * this for larger i values.
  2621. */
  2622. #ifdef Honor_FLT_ROUNDS
  2623. /* round correctly FLT_ROUNDS = 2 or 3 */
  2624. if (sign) {
  2625. rv.d = -rv.d;
  2626. sign = 0;
  2627. }
  2628. #endif
  2629. e -= i;
  2630. dval(&rv) *= tens[i];
  2631. #ifdef VAX
  2632. /* VAX exponent range is so narrow we must
  2633. * worry about overflow here...
  2634. */
  2635. vax_ovfl_check:
  2636. word0(&rv) -= P*Exp_msk1;
  2637. /* rv = */ rounded_product(dval(&rv), tens[e]);
  2638. if ((word0(&rv) & Exp_mask)
  2639. > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
  2640. goto ovfl;
  2641. word0(&rv) += P*Exp_msk1;
  2642. #else
  2643. /* rv = */ rounded_product(dval(&rv), tens[e]);
  2644. #endif
  2645. goto ret;
  2646. }
  2647. }
  2648. #ifndef Inaccurate_Divide
  2649. else if (e >= -Ten_pmax) {
  2650. #ifdef Honor_FLT_ROUNDS
  2651. /* round correctly FLT_ROUNDS = 2 or 3 */
  2652. if (sign) {
  2653. rv.d = -rv.d;
  2654. sign = 0;
  2655. }
  2656. #endif
  2657. /* rv = */ rounded_quotient(dval(&rv), tens[-e]);
  2658. goto ret;
  2659. }
  2660. #endif
  2661. #endif /* ROUND_BIASED_without_Round_Up */
  2662. }
  2663. e1 += nd - k;
  2664. #ifdef IEEE_Arith
  2665. #ifdef SET_INEXACT
  2666. bc.inexact = 1;
  2667. if (k <= DBL_DIG)
  2668. oldinexact = get_inexact();
  2669. #endif
  2670. #ifdef Avoid_Underflow
  2671. bc.scale = 0;
  2672. #endif
  2673. #ifdef Honor_FLT_ROUNDS
  2674. if (bc.rounding >= 2) {
  2675. if (sign)
  2676. bc.rounding = bc.rounding == 2 ? 0 : 2;
  2677. else
  2678. if (bc.rounding != 2)
  2679. bc.rounding = 0;
  2680. }
  2681. #endif
  2682. #endif /*IEEE_Arith*/
  2683. /* Get starting approximation = rv * 10**e1 */
  2684. if (e1 > 0) {
  2685. if ((i = e1 & 15))
  2686. dval(&rv) *= tens[i];
  2687. if (e1 &= ~15) {
  2688. if (e1 > DBL_MAX_10_EXP) {
  2689. ovfl:
  2690. /* Can't trust HUGE_VAL */
  2691. #ifdef IEEE_Arith
  2692. #ifdef Honor_FLT_ROUNDS
  2693. switch(bc.rounding) {
  2694. case 0: /* toward 0 */
  2695. case 3: /* toward -infinity */
  2696. word0(&rv) = Big0;
  2697. word1(&rv) = Big1;
  2698. break;
  2699. default:
  2700. word0(&rv) = Exp_mask;
  2701. word1(&rv) = 0;
  2702. }
  2703. #else /*Honor_FLT_ROUNDS*/
  2704. word0(&rv) = Exp_mask;
  2705. word1(&rv) = 0;
  2706. #endif /*Honor_FLT_ROUNDS*/
  2707. #ifdef SET_INEXACT
  2708. /* set overflow bit */
  2709. dval(&rv0) = 1e300;
  2710. dval(&rv0) *= dval(&rv0);
  2711. #endif
  2712. #else /*IEEE_Arith*/
  2713. word0(&rv) = Big0;
  2714. word1(&rv) = Big1;
  2715. #endif /*IEEE_Arith*/
  2716. range_err:
  2717. if (bd0) {
  2718. Bfree(bb);
  2719. Bfree(bd);
  2720. Bfree(bs);
  2721. Bfree(bd0);
  2722. Bfree(delta);
  2723. }
  2724. #ifndef NO_ERRNO
  2725. errno = ERANGE;
  2726. #endif
  2727. goto ret;
  2728. }
  2729. e1 >>= 4;
  2730. for(j = 0; e1 > 1; j++, e1 >>= 1)
  2731. if (e1 & 1)
  2732. dval(&rv) *= bigtens[j];
  2733. /* The last multiplication could overflow. */
  2734. word0(&rv) -= P*Exp_msk1;
  2735. dval(&rv) *= bigtens[j];
  2736. if ((z = word0(&rv) & Exp_mask)
  2737. > Exp_msk1*(DBL_MAX_EXP+Bias-P))
  2738. goto ovfl;
  2739. if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
  2740. /* set to largest number */
  2741. /* (Can't trust DBL_MAX) */
  2742. word0(&rv) = Big0;
  2743. word1(&rv) = Big1;
  2744. }
  2745. else
  2746. word0(&rv) += P*Exp_msk1;
  2747. }
  2748. }
  2749. else if (e1 < 0) {
  2750. e1 = -e1;
  2751. if ((i = e1 & 15))
  2752. dval(&rv) /= tens[i];
  2753. if (e1 >>= 4) {
  2754. if (e1 >= 1 << n_bigtens)
  2755. goto undfl;
  2756. #ifdef Avoid_Underflow
  2757. if (e1 & Scale_Bit)
  2758. bc.scale = 2*P;
  2759. for(j = 0; e1 > 0; j++, e1 >>= 1)
  2760. if (e1 & 1)
  2761. dval(&rv) *= tinytens[j];
  2762. if (bc.scale && (j = 2*P + 1 - ((word0(&rv) & Exp_mask)
  2763. >> Exp_shift)) > 0) {
  2764. /* scaled rv is denormal; clear j low bits */
  2765. if (j >= 32) {
  2766. if (j > 54)
  2767. goto undfl;
  2768. word1(&rv) = 0;
  2769. if (j >= 53)
  2770. word0(&rv) = (P+2)*Exp_msk1;
  2771. else
  2772. word0(&rv) &= 0xffffffff << (j-32);
  2773. }
  2774. else
  2775. word1(&rv) &= 0xffffffff << j;
  2776. }
  2777. #else
  2778. for(j = 0; e1 > 1; j++, e1 >>= 1)
  2779. if (e1 & 1)
  2780. dval(&rv) *= tinytens[j];
  2781. /* The last multiplication could underflow. */
  2782. dval(&rv0) = dval(&rv);
  2783. dval(&rv) *= tinytens[j];
  2784. if (!dval(&rv)) {
  2785. dval(&rv) = 2.*dval(&rv0);
  2786. dval(&rv) *= tinytens[j];
  2787. #endif
  2788. if (!dval(&rv)) {
  2789. undfl:
  2790. dval(&rv) = 0.;
  2791. goto range_err;
  2792. }
  2793. #ifndef Avoid_Underflow
  2794. word0(&rv) = Tiny0;
  2795. word1(&rv) = Tiny1;
  2796. /* The refinement below will clean
  2797. * this approximation up.
  2798. */
  2799. }
  2800. #endif
  2801. }
  2802. }
  2803. /* Now the hard part -- adjusting rv to the correct value.*/
  2804. /* Put digits into bd: true value = bd * 10^e */
  2805. bc.nd = nd - nz1;
  2806. #ifndef NO_STRTOD_BIGCOMP
  2807. bc.nd0 = nd0; /* Only needed if nd > strtod_diglim, but done here */
  2808. /* to silence an erroneous warning about bc.nd0 */
  2809. /* possibly not being initialized. */
  2810. if (nd > strtod_diglim) {
  2811. /* ASSERT(strtod_diglim >= 18); 18 == one more than the */
  2812. /* minimum number of decimal digits to distinguish double values */
  2813. /* in IEEE arithmetic. */
  2814. i = j = 18;
  2815. if (i > nd0)
  2816. j += bc.dplen;
  2817. for(;;) {
  2818. if (--j < bc.dp1 && j >= bc.dp0)
  2819. j = bc.dp0 - 1;
  2820. if (s0[j] != '0')
  2821. break;
  2822. --i;
  2823. }
  2824. e += nd - i;
  2825. nd = i;
  2826. if (nd0 > nd)
  2827. nd0 = nd;
  2828. if (nd < 9) { /* must recompute y */
  2829. y = 0;
  2830. for(i = 0; i < nd0; ++i)
  2831. y = 10*y + s0[i] - '0';
  2832. for(j = bc.dp1; i < nd; ++i)
  2833. y = 10*y + s0[j++] - '0';
  2834. }
  2835. }
  2836. #endif
  2837. bd0 = s2b(s0, nd0, nd, y, bc.dplen);
  2838. for(;;) {
  2839. bd = Balloc(bd0->k);
  2840. Bcopy(bd, bd0);
  2841. bb = d2b(&rv, &bbe, &bbbits); /* rv = bb * 2^bbe */
  2842. bs = i2b(1);
  2843. if (e >= 0) {
  2844. bb2 = bb5 = 0;
  2845. bd2 = bd5 = e;
  2846. }
  2847. else {
  2848. bb2 = bb5 = -e;
  2849. bd2 = bd5 = 0;
  2850. }
  2851. if (bbe >= 0)
  2852. bb2 += bbe;
  2853. else
  2854. bd2 -= bbe;
  2855. bs2 = bb2;
  2856. #ifdef Honor_FLT_ROUNDS
  2857. if (bc.rounding != 1)
  2858. bs2++;
  2859. #endif
  2860. #ifdef Avoid_Underflow
  2861. Lsb = LSB;
  2862. Lsb1 = 0;
  2863. j = bbe - bc.scale;
  2864. i = j + bbbits - 1; /* logb(rv) */
  2865. j = P + 1 - bbbits;
  2866. if (i < Emin) { /* denormal */
  2867. i = Emin - i;
  2868. j -= i;
  2869. if (i < 32)
  2870. Lsb <<= i;
  2871. else if (i < 52)
  2872. Lsb1 = Lsb << (i-32);
  2873. else
  2874. Lsb1 = Exp_mask;
  2875. }
  2876. #else /*Avoid_Underflow*/
  2877. #ifdef Sudden_Underflow
  2878. #ifdef IBM
  2879. j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
  2880. #else
  2881. j = P + 1 - bbbits;
  2882. #endif
  2883. #else /*Sudden_Underflow*/
  2884. j = bbe;
  2885. i = j + bbbits - 1; /* logb(rv) */
  2886. if (i < Emin) /* denormal */
  2887. j += P - Emin;
  2888. else
  2889. j = P + 1 - bbbits;
  2890. #endif /*Sudden_Underflow*/
  2891. #endif /*Avoid_Underflow*/
  2892. bb2 += j;
  2893. bd2 += j;
  2894. #ifdef Avoid_Underflow
  2895. bd2 += bc.scale;
  2896. #endif
  2897. i = bb2 < bd2 ? bb2 : bd2;
  2898. if (i > bs2)
  2899. i = bs2;
  2900. if (i > 0) {
  2901. bb2 -= i;
  2902. bd2 -= i;
  2903. bs2 -= i;
  2904. }
  2905. if (bb5 > 0) {
  2906. bs = pow5mult(bs, bb5);
  2907. bb1 = mult(bs, bb);
  2908. Bfree(bb);
  2909. bb = bb1;
  2910. }
  2911. if (bb2 > 0)
  2912. bb = lshift(bb, bb2);
  2913. if (bd5 > 0)
  2914. bd = pow5mult(bd, bd5);
  2915. if (bd2 > 0)
  2916. bd = lshift(bd, bd2);
  2917. if (bs2 > 0)
  2918. bs = lshift(bs, bs2);
  2919. delta = diff(bb, bd);
  2920. bc.dsign = delta->sign;
  2921. delta->sign = 0;
  2922. i = cmp(delta, bs);
  2923. #ifndef NO_STRTOD_BIGCOMP /*{*/
  2924. if (bc.nd > nd && i <= 0) {
  2925. if (bc.dsign) {
  2926. /* Must use bigcomp(). */
  2927. req_bigcomp = 1;
  2928. break;
  2929. }
  2930. #ifdef Honor_FLT_ROUNDS
  2931. if (bc.rounding != 1) {
  2932. if (i < 0) {
  2933. req_bigcomp = 1;
  2934. break;
  2935. }
  2936. }
  2937. else
  2938. #endif
  2939. i = -1; /* Discarded digits make delta smaller. */
  2940. }
  2941. #endif /*}*/
  2942. #ifdef Honor_FLT_ROUNDS /*{*/
  2943. if (bc.rounding != 1) {
  2944. if (i < 0) {
  2945. /* Error is less than an ulp */
  2946. if (!delta->x[0] && delta->wds <= 1) {
  2947. /* exact */
  2948. #ifdef SET_INEXACT
  2949. bc.inexact = 0;
  2950. #endif
  2951. break;
  2952. }
  2953. if (bc.rounding) {
  2954. if (bc.dsign) {
  2955. adj.d = 1.;
  2956. goto apply_adj;
  2957. }
  2958. }
  2959. else if (!bc.dsign) {
  2960. adj.d = -1.;
  2961. if (!word1(&rv)
  2962. && !(word0(&rv) & Frac_mask)) {
  2963. y = word0(&rv) & Exp_mask;
  2964. #ifdef Avoid_Underflow
  2965. if (!bc.scale || y > 2*P*Exp_msk1)
  2966. #else
  2967. if (y)
  2968. #endif
  2969. {
  2970. delta = lshift(delta,Log2P);
  2971. if (cmp(delta, bs) <= 0)
  2972. adj.d = -0.5;
  2973. }
  2974. }
  2975. apply_adj:
  2976. #ifdef Avoid_Underflow /*{*/
  2977. if (bc.scale && (y = word0(&rv) & Exp_mask)
  2978. <= 2*P*Exp_msk1)
  2979. word0(&adj) += (2*P+1)*Exp_msk1 - y;
  2980. #else
  2981. #ifdef Sudden_Underflow
  2982. if ((word0(&rv) & Exp_mask) <=
  2983. P*Exp_msk1) {
  2984. word0(&rv) += P*Exp_msk1;
  2985. dval(&rv) += adj.d*ulp(dval(&rv));
  2986. word0(&rv) -= P*Exp_msk1;
  2987. }
  2988. else
  2989. #endif /*Sudden_Underflow*/
  2990. #endif /*Avoid_Underflow}*/
  2991. dval(&rv) += adj.d*ulp(&rv);
  2992. }
  2993. break;
  2994. }
  2995. adj.d = ratio(delta, bs);
  2996. if (adj.d < 1.)
  2997. adj.d = 1.;
  2998. if (adj.d <= 0x7ffffffe) {
  2999. /* adj = rounding ? ceil(adj) : floor(adj); */
  3000. y = adj.d;
  3001. if (y != adj.d) {
  3002. if (!((bc.rounding>>1) ^ bc.dsign))
  3003. y++;
  3004. adj.d = y;
  3005. }
  3006. }
  3007. #ifdef Avoid_Underflow /*{*/
  3008. if (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
  3009. word0(&adj) += (2*P+1)*Exp_msk1 - y;
  3010. #else
  3011. #ifdef Sudden_Underflow
  3012. if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
  3013. word0(&rv) += P*Exp_msk1;
  3014. adj.d *= ulp(dval(&rv));
  3015. if (bc.dsign)
  3016. dval(&rv) += adj.d;
  3017. else
  3018. dval(&rv) -= adj.d;
  3019. word0(&rv) -= P*Exp_msk1;
  3020. goto cont;
  3021. }
  3022. #endif /*Sudden_Underflow*/
  3023. #endif /*Avoid_Underflow}*/
  3024. adj.d *= ulp(&rv);
  3025. if (bc.dsign) {
  3026. if (word0(&rv) == Big0 && word1(&rv) == Big1)
  3027. goto ovfl;
  3028. dval(&rv) += adj.d;
  3029. }
  3030. else
  3031. dval(&rv) -= adj.d;
  3032. goto cont;
  3033. }
  3034. #endif /*}Honor_FLT_ROUNDS*/
  3035. if (i < 0) {
  3036. /* Error is less than half an ulp -- check for
  3037. * special case of mantissa a power of two.
  3038. */
  3039. if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask
  3040. #ifdef IEEE_Arith /*{*/
  3041. #ifdef Avoid_Underflow
  3042. || (word0(&rv) & Exp_mask) <= (2*P+1)*Exp_msk1
  3043. #else
  3044. || (word0(&rv) & Exp_mask) <= Exp_msk1
  3045. #endif
  3046. #endif /*}*/
  3047. ) {
  3048. #ifdef SET_INEXACT
  3049. if (!delta->x[0] && delta->wds <= 1)
  3050. bc.inexact = 0;
  3051. #endif
  3052. break;
  3053. }
  3054. if (!delta->x[0] && delta->wds <= 1) {
  3055. /* exact result */
  3056. #ifdef SET_INEXACT
  3057. bc.inexact = 0;
  3058. #endif
  3059. break;
  3060. }
  3061. delta = lshift(delta,Log2P);
  3062. if (cmp(delta, bs) > 0)
  3063. goto drop_down;
  3064. break;
  3065. }
  3066. if (i == 0) {
  3067. /* exactly half-way between */
  3068. if (bc.dsign) {
  3069. if ((word0(&rv) & Bndry_mask1) == Bndry_mask1
  3070. && word1(&rv) == (
  3071. #ifdef Avoid_Underflow
  3072. (bc.scale && (y = word0(&rv) & Exp_mask) <= 2*P*Exp_msk1)
  3073. ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
  3074. #endif
  3075. 0xffffffff)) {
  3076. /*boundary case -- increment exponent*/
  3077. if (word0(&rv) == Big0 && word1(&rv) == Big1)
  3078. goto ovfl;
  3079. word0(&rv) = (word0(&rv) & Exp_mask)
  3080. + Exp_msk1
  3081. #ifdef IBM
  3082. | Exp_msk1 >> 4
  3083. #endif
  3084. ;
  3085. word1(&rv) = 0;
  3086. #ifdef Avoid_Underflow
  3087. bc.dsign = 0;
  3088. #endif
  3089. break;
  3090. }
  3091. }
  3092. else if (!(word0(&rv) & Bndry_mask) && !word1(&rv)) {
  3093. drop_down:
  3094. /* boundary case -- decrement exponent */
  3095. #ifdef Sudden_Underflow /*{{*/
  3096. L = word0(&rv) & Exp_mask;
  3097. #ifdef IBM
  3098. if (L < Exp_msk1)
  3099. #else
  3100. #ifdef Avoid_Underflow
  3101. if (L <= (bc.scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
  3102. #else
  3103. if (L <= Exp_msk1)
  3104. #endif /*Avoid_Underflow*/
  3105. #endif /*IBM*/
  3106. {
  3107. if (bc.nd >nd) {
  3108. bc.uflchk = 1;
  3109. break;
  3110. }
  3111. goto undfl;
  3112. }
  3113. L -= Exp_msk1;
  3114. #else /*Sudden_Underflow}{*/
  3115. #ifdef Avoid_Underflow
  3116. if (bc.scale) {
  3117. L = word0(&rv) & Exp_mask;
  3118. if (L <= (2*P+1)*Exp_msk1) {
  3119. if (L > (P+2)*Exp_msk1)
  3120. /* round even ==> */
  3121. /* accept rv */
  3122. break;
  3123. /* rv = smallest denormal */
  3124. if (bc.nd >nd) {
  3125. bc.uflchk = 1;
  3126. break;
  3127. }
  3128. goto undfl;
  3129. }
  3130. }
  3131. #endif /*Avoid_Underflow*/
  3132. L = (word0(&rv) & Exp_mask) - Exp_msk1;
  3133. #endif /*Sudden_Underflow}}*/
  3134. word0(&rv) = L | Bndry_mask1;
  3135. word1(&rv) = 0xffffffff;
  3136. #ifdef IBM
  3137. goto cont;
  3138. #else
  3139. #ifndef NO_STRTOD_BIGCOMP
  3140. if (bc.nd > nd)
  3141. goto cont;
  3142. #endif
  3143. break;
  3144. #endif
  3145. }
  3146. #ifndef ROUND_BIASED
  3147. #ifdef Avoid_Underflow
  3148. if (Lsb1) {
  3149. if (!(word0(&rv) & Lsb1))
  3150. break;
  3151. }
  3152. else if (!(word1(&rv) & Lsb))
  3153. break;
  3154. #else
  3155. if (!(word1(&rv) & LSB))
  3156. break;
  3157. #endif
  3158. #endif
  3159. if (bc.dsign)
  3160. #ifdef Avoid_Underflow
  3161. dval(&rv) += sulp(&rv, &bc);
  3162. #else
  3163. dval(&rv) += ulp(&rv);
  3164. #endif
  3165. #ifndef ROUND_BIASED
  3166. else {
  3167. #ifdef Avoid_Underflow
  3168. dval(&rv) -= sulp(&rv, &bc);
  3169. #else
  3170. dval(&rv) -= ulp(&rv);
  3171. #endif
  3172. #ifndef Sudden_Underflow
  3173. if (!dval(&rv)) {
  3174. if (bc.nd >nd) {
  3175. bc.uflchk = 1;
  3176. break;
  3177. }
  3178. goto undfl;
  3179. }
  3180. #endif
  3181. }
  3182. #ifdef Avoid_Underflow
  3183. bc.dsign = 1 - bc.dsign;
  3184. #endif
  3185. #endif
  3186. break;
  3187. }
  3188. if ((aadj = ratio(delta, bs)) <= 2.) {
  3189. if (bc.dsign)
  3190. aadj = aadj1 = 1.;
  3191. else if (word1(&rv) || word0(&rv) & Bndry_mask) {
  3192. #ifndef Sudden_Underflow
  3193. if (word1(&rv) == Tiny1 && !word0(&rv)) {
  3194. if (bc.nd >nd) {
  3195. bc.uflchk = 1;
  3196. break;
  3197. }
  3198. goto undfl;
  3199. }
  3200. #endif
  3201. aadj = 1.;
  3202. aadj1 = -1.;
  3203. }
  3204. else {
  3205. /* special case -- power of FLT_RADIX to be */
  3206. /* rounded down... */
  3207. if (aadj < 2./FLT_RADIX)
  3208. aadj = 1./FLT_RADIX;
  3209. else
  3210. aadj *= 0.5;
  3211. aadj1 = -aadj;
  3212. }
  3213. }
  3214. else {
  3215. aadj *= 0.5;
  3216. aadj1 = bc.dsign ? aadj : -aadj;
  3217. #ifdef Check_FLT_ROUNDS
  3218. switch(bc.rounding) {
  3219. case 2: /* towards +infinity */
  3220. aadj1 -= 0.5;
  3221. break;
  3222. case 0: /* towards 0 */
  3223. case 3: /* towards -infinity */
  3224. aadj1 += 0.5;
  3225. }
  3226. #else
  3227. if (Flt_Rounds == 0)
  3228. aadj1 += 0.5;
  3229. #endif /*Check_FLT_ROUNDS*/
  3230. }
  3231. y = word0(&rv) & Exp_mask;
  3232. /* Check for overflow */
  3233. if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
  3234. dval(&rv0) = dval(&rv);
  3235. word0(&rv) -= P*Exp_msk1;
  3236. adj.d = aadj1 * ulp(&rv);
  3237. dval(&rv) += adj.d;
  3238. if ((word0(&rv) & Exp_mask) >=
  3239. Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
  3240. if (word0(&rv0) == Big0 && word1(&rv0) == Big1)
  3241. goto ovfl;
  3242. word0(&rv) = Big0;
  3243. word1(&rv) = Big1;
  3244. goto cont;
  3245. }
  3246. else
  3247. word0(&rv) += P*Exp_msk1;
  3248. }
  3249. else {
  3250. #ifdef Avoid_Underflow
  3251. if (bc.scale && y <= 2*P*Exp_msk1) {
  3252. if (aadj <= 0x7fffffff) {
  3253. if ((z = aadj) <= 0)
  3254. z = 1;
  3255. aadj = z;
  3256. aadj1 = bc.dsign ? aadj : -aadj;
  3257. }
  3258. dval(&aadj2) = aadj1;
  3259. word0(&aadj2) += (2*P+1)*Exp_msk1 - y;
  3260. aadj1 = dval(&aadj2);
  3261. adj.d = aadj1 * ulp(&rv);
  3262. dval(&rv) += adj.d;
  3263. if (rv.d == 0.)
  3264. #ifdef NO_STRTOD_BIGCOMP
  3265. goto undfl;
  3266. #else
  3267. {
  3268. if (bc.nd > nd)
  3269. bc.dsign = 1;
  3270. break;
  3271. }
  3272. #endif
  3273. }
  3274. else {
  3275. adj.d = aadj1 * ulp(&rv);
  3276. dval(&rv) += adj.d;
  3277. }
  3278. #else
  3279. #ifdef Sudden_Underflow
  3280. if ((word0(&rv) & Exp_mask) <= P*Exp_msk1) {
  3281. dval(&rv0) = dval(&rv);
  3282. word0(&rv) += P*Exp_msk1;
  3283. adj.d = aadj1 * ulp(&rv);
  3284. dval(&rv) += adj.d;
  3285. #ifdef IBM
  3286. if ((word0(&rv) & Exp_mask) < P*Exp_msk1)
  3287. #else
  3288. if ((word0(&rv) & Exp_mask) <= P*Exp_msk1)
  3289. #endif
  3290. {
  3291. if (word0(&rv0) == Tiny0
  3292. && word1(&rv0) == Tiny1) {
  3293. if (bc.nd >nd) {
  3294. bc.uflchk = 1;
  3295. break;
  3296. }
  3297. goto undfl;
  3298. }
  3299. word0(&rv) = Tiny0;
  3300. word1(&rv) = Tiny1;
  3301. goto cont;
  3302. }
  3303. else
  3304. word0(&rv) -= P*Exp_msk1;
  3305. }
  3306. else {
  3307. adj.d = aadj1 * ulp(&rv);
  3308. dval(&rv) += adj.d;
  3309. }
  3310. #else /*Sudden_Underflow*/
  3311. /* Compute adj so that the IEEE rounding rules will
  3312. * correctly round rv + adj in some half-way cases.
  3313. * If rv * ulp(rv) is denormalized (i.e.,
  3314. * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
  3315. * trouble from bits lost to denormalization;
  3316. * example: 1.2e-307 .
  3317. */
  3318. if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
  3319. aadj1 = (double)(int)(aadj + 0.5);
  3320. if (!bc.dsign)
  3321. aadj1 = -aadj1;
  3322. }
  3323. adj.d = aadj1 * ulp(&rv);
  3324. dval(&rv) += adj.d;
  3325. #endif /*Sudden_Underflow*/
  3326. #endif /*Avoid_Underflow*/
  3327. }
  3328. z = word0(&rv) & Exp_mask;
  3329. #ifndef SET_INEXACT
  3330. if (bc.nd == nd) {
  3331. #ifdef Avoid_Underflow
  3332. if (!bc.scale)
  3333. #endif
  3334. if (y == z) {
  3335. /* Can we stop now? */
  3336. L = (Long)aadj;
  3337. aadj -= L;
  3338. /* The tolerances below are conservative. */
  3339. if (bc.dsign || word1(&rv) || word0(&rv) & Bndry_mask) {
  3340. if (aadj < .4999999 || aadj > .5000001)
  3341. break;
  3342. }
  3343. else if (aadj < .4999999/FLT_RADIX)
  3344. break;
  3345. }
  3346. }
  3347. #endif
  3348. cont:
  3349. Bfree(bb);
  3350. Bfree(bd);
  3351. Bfree(bs);
  3352. Bfree(delta);
  3353. }
  3354. Bfree(bb);
  3355. Bfree(bd);
  3356. Bfree(bs);
  3357. Bfree(bd0);
  3358. Bfree(delta);
  3359. #ifndef NO_STRTOD_BIGCOMP
  3360. if (req_bigcomp) {
  3361. bd0 = 0;
  3362. bc.e0 += nz1;
  3363. bigcomp(&rv, s0, &bc);
  3364. y = word0(&rv) & Exp_mask;
  3365. if (y == Exp_mask)
  3366. goto ovfl;
  3367. if (y == 0 && rv.d == 0.)
  3368. goto undfl;
  3369. }
  3370. #endif
  3371. #ifdef SET_INEXACT
  3372. if (bc.inexact) {
  3373. if (!oldinexact) {
  3374. word0(&rv0) = Exp_1 + (70 << Exp_shift);
  3375. word1(&rv0) = 0;
  3376. dval(&rv0) += 1.;
  3377. }
  3378. }
  3379. else if (!oldinexact)
  3380. clear_inexact();
  3381. #endif
  3382. #ifdef Avoid_Underflow
  3383. if (bc.scale) {
  3384. word0(&rv0) = Exp_1 - 2*P*Exp_msk1;
  3385. word1(&rv0) = 0;
  3386. dval(&rv) *= dval(&rv0);
  3387. #ifndef NO_ERRNO
  3388. /* try to avoid the bug of testing an 8087 register value */
  3389. #ifdef IEEE_Arith
  3390. if (!(word0(&rv) & Exp_mask))
  3391. #else
  3392. if (word0(&rv) == 0 && word1(&rv) == 0)
  3393. #endif
  3394. errno = ERANGE;
  3395. #endif
  3396. }
  3397. #endif /* Avoid_Underflow */
  3398. #ifdef SET_INEXACT
  3399. if (bc.inexact && !(word0(&rv) & Exp_mask)) {
  3400. /* set underflow bit */
  3401. dval(&rv0) = 1e-300;
  3402. dval(&rv0) *= dval(&rv0);
  3403. }
  3404. #endif
  3405. ret:
  3406. if (se)
  3407. *se = (char *)s;
  3408. return sign ? -dval(&rv) : dval(&rv);
  3409. }
  3410. #ifndef MULTIPLE_THREADS
  3411. static char *dtoa_result;
  3412. #endif
  3413. static char *
  3414. #ifdef KR_headers
  3415. rv_alloc(i) int i;
  3416. #else
  3417. rv_alloc(int i)
  3418. #endif
  3419. {
  3420. int j, k, *r;
  3421. j = sizeof(ULong);
  3422. for(k = 0;
  3423. sizeof(Bigint) - sizeof(ULong) - sizeof(int) + j <= i;
  3424. j <<= 1)
  3425. k++;
  3426. r = (int*)Balloc(k);
  3427. *r = k;
  3428. return
  3429. #ifndef MULTIPLE_THREADS
  3430. dtoa_result =
  3431. #endif
  3432. (char *)(r+1);
  3433. }
  3434. static char *
  3435. #ifdef KR_headers
  3436. nrv_alloc(s, rve, n) char *s, **rve; int n;
  3437. #else
  3438. nrv_alloc(const char *s, char **rve, int n)
  3439. #endif
  3440. {
  3441. char *rv, *t;
  3442. t = rv = rv_alloc(n);
  3443. while((*t = *s++)) t++;
  3444. if (rve)
  3445. *rve = t;
  3446. return rv;
  3447. }
  3448. /* freedtoa(s) must be used to free values s returned by dtoa
  3449. * when MULTIPLE_THREADS is #defined. It should be used in all cases,
  3450. * but for consistency with earlier versions of dtoa, it is optional
  3451. * when MULTIPLE_THREADS is not defined.
  3452. */
  3453. void
  3454. #ifdef KR_headers
  3455. freedtoa(s) char *s;
  3456. #else
  3457. freedtoa(char *s)
  3458. #endif
  3459. {
  3460. Bigint *b = (Bigint *)((int *)s - 1);
  3461. b->maxwds = 1 << (b->k = *(int*)b);
  3462. Bfree(b);
  3463. #ifndef MULTIPLE_THREADS
  3464. if (s == dtoa_result)
  3465. dtoa_result = 0;
  3466. #endif
  3467. }
  3468. /* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.
  3469. *
  3470. * Inspired by "How to Print Floating-Point Numbers Accurately" by
  3471. * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 112-126].
  3472. *
  3473. * Modifications:
  3474. * 1. Rather than iterating, we use a simple numeric overestimate
  3475. * to determine k = floor(log10(d)). We scale relevant
  3476. * quantities using O(log2(k)) rather than O(k) multiplications.
  3477. * 2. For some modes > 2 (corresponding to ecvt and fcvt), we don't
  3478. * try to generate digits strictly left to right. Instead, we
  3479. * compute with fewer bits and propagate the carry if necessary
  3480. * when rounding the final digit up. This is often faster.
  3481. * 3. Under the assumption that input will be rounded nearest,
  3482. * mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.
  3483. * That is, we allow equality in stopping tests when the
  3484. * round-nearest rule will give the same floating-point value
  3485. * as would satisfaction of the stopping test with strict
  3486. * inequality.
  3487. * 4. We remove common factors of powers of 2 from relevant
  3488. * quantities.
  3489. * 5. When converting floating-point integers less than 1e16,
  3490. * we use floating-point arithmetic rather than resorting
  3491. * to multiple-precision integers.
  3492. * 6. When asked to produce fewer than 15 digits, we first try
  3493. * to get by with floating-point arithmetic; we resort to
  3494. * multiple-precision integer arithmetic only if we cannot
  3495. * guarantee that the floating-point calculation has given
  3496. * the correctly rounded result. For k requested digits and
  3497. * "uniformly" distributed input, the probability is
  3498. * something like 10^(k-15) that we must resort to the Long
  3499. * calculation.
  3500. */
  3501. char *
  3502. dtoa
  3503. #ifdef KR_headers
  3504. (dd, mode, ndigits, decpt, sign, rve)
  3505. double dd; int mode, ndigits, *decpt, *sign; char **rve;
  3506. #else
  3507. (double dd, int mode, int ndigits, int *decpt, int *sign, char **rve)
  3508. #endif
  3509. {
  3510. /* Arguments ndigits, decpt, sign are similar to those
  3511. of ecvt and fcvt; trailing zeros are suppressed from
  3512. the returned string. If not null, *rve is set to point
  3513. to the end of the return value. If d is +-Infinity or NaN,
  3514. then *decpt is set to 9999.
  3515. mode:
  3516. 0 ==> shortest string that yields d when read in
  3517. and rounded to nearest.
  3518. 1 ==> like 0, but with Steele & White stopping rule;
  3519. e.g. with IEEE P754 arithmetic , mode 0 gives
  3520. 1e23 whereas mode 1 gives 9.999999999999999e22.
  3521. 2 ==> max(1,ndigits) significant digits. This gives a
  3522. return value similar to that of ecvt, except
  3523. that trailing zeros are suppressed.
  3524. 3 ==> through ndigits past the decimal point. This
  3525. gives a return value similar to that from fcvt,
  3526. except that trailing zeros are suppressed, and
  3527. ndigits can be negative.
  3528. 4,5 ==> similar to 2 and 3, respectively, but (in
  3529. round-nearest mode) with the tests of mode 0 to
  3530. possibly return a shorter string that rounds to d.
  3531. With IEEE arithmetic and compilation with
  3532. -DHonor_FLT_ROUNDS, modes 4 and 5 behave the same
  3533. as modes 2 and 3 when FLT_ROUNDS != 1.
  3534. 6-9 ==> Debugging modes similar to mode - 4: don't try
  3535. fast floating-point estimate (if applicable).
  3536. Values of mode other than 0-9 are treated as mode 0.
  3537. Sufficient space is allocated to the return value
  3538. to hold the suppressed trailing zeros.
  3539. */
  3540. int bbits, b2, b5, be, dig, i, ieps, ilim, ilim0, ilim1,
  3541. j, j1, k, k0, k_check, leftright, m2, m5, s2, s5,
  3542. spec_case, try_quick;
  3543. Long L;
  3544. #ifndef Sudden_Underflow
  3545. int denorm;
  3546. ULong x;
  3547. #endif
  3548. Bigint *b, *b1, *delta, *mlo, *mhi, *S;
  3549. U d2, eps, u;
  3550. double ds;
  3551. char *s, *s0;
  3552. #ifndef No_leftright
  3553. #ifdef IEEE_Arith
  3554. U eps1;
  3555. #endif
  3556. #endif
  3557. #ifdef SET_INEXACT
  3558. int inexact, oldinexact;
  3559. #endif
  3560. #ifdef Honor_FLT_ROUNDS /*{*/
  3561. int Rounding;
  3562. #ifdef Trust_FLT_ROUNDS /*{{ only define this if FLT_ROUNDS really works! */
  3563. Rounding = Flt_Rounds;
  3564. #else /*}{*/
  3565. Rounding = 1;
  3566. switch(fegetround()) {
  3567. case FE_TOWARDZERO: Rounding = 0; break;
  3568. case FE_UPWARD: Rounding = 2; break;
  3569. case FE_DOWNWARD: Rounding = 3;
  3570. }
  3571. #endif /*}}*/
  3572. #endif /*}*/
  3573. #ifndef MULTIPLE_THREADS
  3574. if (dtoa_result) {
  3575. freedtoa(dtoa_result);
  3576. dtoa_result = 0;
  3577. }
  3578. #endif
  3579. u.d = dd;
  3580. if (word0(&u) & Sign_bit) {
  3581. /* set sign for everything, including 0's and NaNs */
  3582. *sign = 1;
  3583. word0(&u) &= ~Sign_bit; /* clear sign bit */
  3584. }
  3585. else
  3586. *sign = 0;
  3587. #if defined(IEEE_Arith) + defined(VAX)
  3588. #ifdef IEEE_Arith
  3589. if ((word0(&u) & Exp_mask) == Exp_mask)
  3590. #else
  3591. if (word0(&u) == 0x8000)
  3592. #endif
  3593. {
  3594. /* Infinity or NaN */
  3595. *decpt = 9999;
  3596. #ifdef IEEE_Arith
  3597. if (!word1(&u) && !(word0(&u) & 0xfffff))
  3598. return nrv_alloc("Infinity", rve, 8);
  3599. #endif
  3600. return nrv_alloc("NaN", rve, 3);
  3601. }
  3602. #endif
  3603. #ifdef IBM
  3604. dval(&u) += 0; /* normalize */
  3605. #endif
  3606. if (!dval(&u)) {
  3607. *decpt = 1;
  3608. return nrv_alloc("0", rve, 1);
  3609. }
  3610. #ifdef SET_INEXACT
  3611. try_quick = oldinexact = get_inexact();
  3612. inexact = 1;
  3613. #endif
  3614. #ifdef Honor_FLT_ROUNDS
  3615. if (Rounding >= 2) {
  3616. if (*sign)
  3617. Rounding = Rounding == 2 ? 0 : 2;
  3618. else
  3619. if (Rounding != 2)
  3620. Rounding = 0;
  3621. }
  3622. #endif
  3623. b = d2b(&u, &be, &bbits);
  3624. #ifdef Sudden_Underflow
  3625. i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1));
  3626. #else
  3627. if ((i = (int)(word0(&u) >> Exp_shift1 & (Exp_mask>>Exp_shift1)))) {
  3628. #endif
  3629. dval(&d2) = dval(&u);
  3630. word0(&d2) &= Frac_mask1;
  3631. word0(&d2) |= Exp_11;
  3632. #ifdef IBM
  3633. if (j = 11 - hi0bits(word0(&d2) & Frac_mask))
  3634. dval(&d2) /= 1 << j;
  3635. #endif
  3636. /* log(x) ~=~ log(1.5) + (x-1.5)/1.5
  3637. * log10(x) = log(x) / log(10)
  3638. * ~=~ log(1.5)/log(10) + (x-1.5)/(1.5*log(10))
  3639. * log10(d) = (i-Bias)*log(2)/log(10) + log10(d2)
  3640. *
  3641. * This suggests computing an approximation k to log10(d) by
  3642. *
  3643. * k = (i - Bias)*0.301029995663981
  3644. * + ( (d2-1.5)*0.289529654602168 + 0.176091259055681 );
  3645. *
  3646. * We want k to be too large rather than too small.
  3647. * The error in the first-order Taylor series approximation
  3648. * is in our favor, so we just round up the constant enough
  3649. * to compensate for any error in the multiplication of
  3650. * (i - Bias) by 0.301029995663981; since |i - Bias| <= 1077,
  3651. * and 1077 * 0.30103 * 2^-52 ~=~ 7.2e-14,
  3652. * adding 1e-13 to the constant term more than suffices.
  3653. * Hence we adjust the constant term to 0.1760912590558.
  3654. * (We could get a more accurate k by invoking log10,
  3655. * but this is probably not worthwhile.)
  3656. */
  3657. i -= Bias;
  3658. #ifdef IBM
  3659. i <<= 2;
  3660. i += j;
  3661. #endif
  3662. #ifndef Sudden_Underflow
  3663. denorm = 0;
  3664. }
  3665. else {
  3666. /* d is denormalized */
  3667. i = bbits + be + (Bias + (P-1) - 1);
  3668. x = i > 32 ? word0(&u) << (64 - i) | word1(&u) >> (i - 32)
  3669. : word1(&u) << (32 - i);
  3670. dval(&d2) = x;
  3671. word0(&d2) -= 31*Exp_msk1; /* adjust exponent */
  3672. i -= (Bias + (P-1) - 1) + 1;
  3673. denorm = 1;
  3674. }
  3675. #endif
  3676. ds = (dval(&d2)-1.5)*0.289529654602168 + 0.1760912590558 + i*0.301029995663981;
  3677. k = (int)ds;
  3678. if (ds < 0. && ds != k)
  3679. k--; /* want k = floor(ds) */
  3680. k_check = 1;
  3681. if (k >= 0 && k <= Ten_pmax) {
  3682. if (dval(&u) < tens[k])
  3683. k--;
  3684. k_check = 0;
  3685. }
  3686. j = bbits - i - 1;
  3687. if (j >= 0) {
  3688. b2 = 0;
  3689. s2 = j;
  3690. }
  3691. else {
  3692. b2 = -j;
  3693. s2 = 0;
  3694. }
  3695. if (k >= 0) {
  3696. b5 = 0;
  3697. s5 = k;
  3698. s2 += k;
  3699. }
  3700. else {
  3701. b2 -= k;
  3702. b5 = -k;
  3703. s5 = 0;
  3704. }
  3705. if (mode < 0 || mode > 9)
  3706. mode = 0;
  3707. #ifndef SET_INEXACT
  3708. #ifdef Check_FLT_ROUNDS
  3709. try_quick = Rounding == 1;
  3710. #else
  3711. try_quick = 1;
  3712. #endif
  3713. #endif /*SET_INEXACT*/
  3714. if (mode > 5) {
  3715. mode -= 4;
  3716. try_quick = 0;
  3717. }
  3718. leftright = 1;
  3719. ilim = ilim1 = -1; /* Values for cases 0 and 1; done here to */
  3720. /* silence erroneous "gcc -Wall" warning. */
  3721. switch(mode) {
  3722. case 0:
  3723. case 1:
  3724. i = 18;
  3725. ndigits = 0;
  3726. break;
  3727. case 2:
  3728. leftright = 0;
  3729. /* no break */
  3730. case 4:
  3731. if (ndigits <= 0)
  3732. ndigits = 1;
  3733. ilim = ilim1 = i = ndigits;
  3734. break;
  3735. case 3:
  3736. leftright = 0;
  3737. /* no break */
  3738. case 5:
  3739. i = ndigits + k + 1;
  3740. ilim = i;
  3741. ilim1 = i - 1;
  3742. if (i <= 0)
  3743. i = 1;
  3744. }
  3745. s = s0 = rv_alloc(i);
  3746. #ifdef Honor_FLT_ROUNDS
  3747. if (mode > 1 && Rounding != 1)
  3748. leftright = 0;
  3749. #endif
  3750. if (ilim >= 0 && ilim <= Quick_max && try_quick) {
  3751. /* Try to get by with floating-point arithmetic. */
  3752. i = 0;
  3753. dval(&d2) = dval(&u);
  3754. k0 = k;
  3755. ilim0 = ilim;
  3756. ieps = 2; /* conservative */
  3757. if (k > 0) {
  3758. ds = tens[k&0xf];
  3759. j = k >> 4;
  3760. if (j & Bletch) {
  3761. /* prevent overflows */
  3762. j &= Bletch - 1;
  3763. dval(&u) /= bigtens[n_bigtens-1];
  3764. ieps++;
  3765. }
  3766. for(; j; j >>= 1, i++)
  3767. if (j & 1) {
  3768. ieps++;
  3769. ds *= bigtens[i];
  3770. }
  3771. dval(&u) /= ds;
  3772. }
  3773. else if ((j1 = -k)) {
  3774. dval(&u) *= tens[j1 & 0xf];
  3775. for(j = j1 >> 4; j; j >>= 1, i++)
  3776. if (j & 1) {
  3777. ieps++;
  3778. dval(&u) *= bigtens[i];
  3779. }
  3780. }
  3781. if (k_check && dval(&u) < 1. && ilim > 0) {
  3782. if (ilim1 <= 0)
  3783. goto fast_failed;
  3784. ilim = ilim1;
  3785. k--;
  3786. dval(&u) *= 10.;
  3787. ieps++;
  3788. }
  3789. dval(&eps) = ieps*dval(&u) + 7.;
  3790. word0(&eps) -= (P-1)*Exp_msk1;
  3791. if (ilim == 0) {
  3792. S = mhi = 0;
  3793. dval(&u) -= 5.;
  3794. if (dval(&u) > dval(&eps))
  3795. goto one_digit;
  3796. if (dval(&u) < -dval(&eps))
  3797. goto no_digits;
  3798. goto fast_failed;
  3799. }
  3800. #ifndef No_leftright
  3801. if (leftright) {
  3802. /* Use Steele & White method of only
  3803. * generating digits needed.
  3804. */
  3805. dval(&eps) = 0.5/tens[ilim-1] - dval(&eps);
  3806. #ifdef IEEE_Arith
  3807. if (k0 < 0 && j1 >= 307) {
  3808. eps1.d = 1.01e256; /* 1.01 allows roundoff in the next few lines */
  3809. word0(&eps1) -= Exp_msk1 * (Bias+P-1);
  3810. dval(&eps1) *= tens[j1 & 0xf];
  3811. for(i = 0, j = (j1-256) >> 4; j; j >>= 1, i++)
  3812. if (j & 1)
  3813. dval(&eps1) *= bigtens[i];
  3814. if (eps.d < eps1.d)
  3815. eps.d = eps1.d;
  3816. }
  3817. #endif
  3818. for(i = 0;;) {
  3819. L = dval(&u);
  3820. dval(&u) -= L;
  3821. *s++ = '0' + (int)L;
  3822. if (1. - dval(&u) < dval(&eps))
  3823. goto bump_up;
  3824. if (dval(&u) < dval(&eps))
  3825. goto ret1;
  3826. if (++i >= ilim)
  3827. break;
  3828. dval(&eps) *= 10.;
  3829. dval(&u) *= 10.;
  3830. }
  3831. }
  3832. else {
  3833. #endif
  3834. /* Generate ilim digits, then fix them up. */
  3835. dval(&eps) *= tens[ilim-1];
  3836. for(i = 1;; i++, dval(&u) *= 10.) {
  3837. L = (Long)(dval(&u));
  3838. if (!(dval(&u) -= L))
  3839. ilim = i;
  3840. *s++ = '0' + (int)L;
  3841. if (i == ilim) {
  3842. if (dval(&u) > 0.5 + dval(&eps))
  3843. goto bump_up;
  3844. else if (dval(&u) < 0.5 - dval(&eps)) {
  3845. while(*--s == '0');
  3846. s++;
  3847. goto ret1;
  3848. }
  3849. break;
  3850. }
  3851. }
  3852. #ifndef No_leftright
  3853. }
  3854. #endif
  3855. fast_failed:
  3856. s = s0;
  3857. dval(&u) = dval(&d2);
  3858. k = k0;
  3859. ilim = ilim0;
  3860. }
  3861. /* Do we have a "small" integer? */
  3862. if (be >= 0 && k <= Int_max) {
  3863. /* Yes. */
  3864. ds = tens[k];
  3865. if (ndigits < 0 && ilim <= 0) {
  3866. S = mhi = 0;
  3867. if (ilim < 0 || dval(&u) <= 5*ds)
  3868. goto no_digits;
  3869. goto one_digit;
  3870. }
  3871. for(i = 1;; i++, dval(&u) *= 10.) {
  3872. L = (Long)(dval(&u) / ds);
  3873. dval(&u) -= L*ds;
  3874. #ifdef Check_FLT_ROUNDS
  3875. /* If FLT_ROUNDS == 2, L will usually be high by 1 */
  3876. if (dval(&u) < 0) {
  3877. L--;
  3878. dval(&u) += ds;
  3879. }
  3880. #endif
  3881. *s++ = '0' + (int)L;
  3882. if (!dval(&u)) {
  3883. #ifdef SET_INEXACT
  3884. inexact = 0;
  3885. #endif
  3886. break;
  3887. }
  3888. if (i == ilim) {
  3889. #ifdef Honor_FLT_ROUNDS
  3890. if (mode > 1)
  3891. switch(Rounding) {
  3892. case 0: goto ret1;
  3893. case 2: goto bump_up;
  3894. }
  3895. #endif
  3896. dval(&u) += dval(&u);
  3897. #ifdef ROUND_BIASED
  3898. if (dval(&u) >= ds)
  3899. #else
  3900. if (dval(&u) > ds || (dval(&u) == ds && L & 1))
  3901. #endif
  3902. {
  3903. bump_up:
  3904. while(*--s == '9')
  3905. if (s == s0) {
  3906. k++;
  3907. *s = '0';
  3908. break;
  3909. }
  3910. ++*s++;
  3911. }
  3912. break;
  3913. }
  3914. }
  3915. goto ret1;
  3916. }
  3917. m2 = b2;
  3918. m5 = b5;
  3919. mhi = mlo = 0;
  3920. if (leftright) {
  3921. i =
  3922. #ifndef Sudden_Underflow
  3923. denorm ? be + (Bias + (P-1) - 1 + 1) :
  3924. #endif
  3925. #ifdef IBM
  3926. 1 + 4*P - 3 - bbits + ((bbits + be - 1) & 3);
  3927. #else
  3928. 1 + P - bbits;
  3929. #endif
  3930. b2 += i;
  3931. s2 += i;
  3932. mhi = i2b(1);
  3933. }
  3934. if (m2 > 0 && s2 > 0) {
  3935. i = m2 < s2 ? m2 : s2;
  3936. b2 -= i;
  3937. m2 -= i;
  3938. s2 -= i;
  3939. }
  3940. if (b5 > 0) {
  3941. if (leftright) {
  3942. if (m5 > 0) {
  3943. mhi = pow5mult(mhi, m5);
  3944. b1 = mult(mhi, b);
  3945. Bfree(b);
  3946. b = b1;
  3947. }
  3948. if ((j = b5 - m5))
  3949. b = pow5mult(b, j);
  3950. }
  3951. else
  3952. b = pow5mult(b, b5);
  3953. }
  3954. S = i2b(1);
  3955. if (s5 > 0)
  3956. S = pow5mult(S, s5);
  3957. /* Check for special case that d is a normalized power of 2. */
  3958. spec_case = 0;
  3959. if ((mode < 2 || leftright)
  3960. #ifdef Honor_FLT_ROUNDS
  3961. && Rounding == 1
  3962. #endif
  3963. ) {
  3964. if (!word1(&u) && !(word0(&u) & Bndry_mask)
  3965. #ifndef Sudden_Underflow
  3966. && word0(&u) & (Exp_mask & ~Exp_msk1)
  3967. #endif
  3968. ) {
  3969. /* The special case */
  3970. b2 += Log2P;
  3971. s2 += Log2P;
  3972. spec_case = 1;
  3973. }
  3974. }
  3975. /* Arrange for convenient computation of quotients:
  3976. * shift left if necessary so divisor has 4 leading 0 bits.
  3977. *
  3978. * Perhaps we should just compute leading 28 bits of S once
  3979. * and for all and pass them and a shift to quorem, so it
  3980. * can do shifts and ors to compute the numerator for q.
  3981. */
  3982. i = dshift(S, s2);
  3983. b2 += i;
  3984. m2 += i;
  3985. s2 += i;
  3986. if (b2 > 0)
  3987. b = lshift(b, b2);
  3988. if (s2 > 0)
  3989. S = lshift(S, s2);
  3990. if (k_check) {
  3991. if (cmp(b,S) < 0) {
  3992. k--;
  3993. b = multadd(b, 10, 0); /* we botched the k estimate */
  3994. if (leftright)
  3995. mhi = multadd(mhi, 10, 0);
  3996. ilim = ilim1;
  3997. }
  3998. }
  3999. if (ilim <= 0 && (mode == 3 || mode == 5)) {
  4000. if (ilim < 0 || cmp(b,S = multadd(S,5,0)) <= 0) {
  4001. /* no digits, fcvt style */
  4002. no_digits:
  4003. k = -1 - ndigits;
  4004. goto ret;
  4005. }
  4006. one_digit:
  4007. *s++ = '1';
  4008. k++;
  4009. goto ret;
  4010. }
  4011. if (leftright) {
  4012. if (m2 > 0)
  4013. mhi = lshift(mhi, m2);
  4014. /* Compute mlo -- check for special case
  4015. * that d is a normalized power of 2.
  4016. */
  4017. mlo = mhi;
  4018. if (spec_case) {
  4019. mhi = Balloc(mhi->k);
  4020. Bcopy(mhi, mlo);
  4021. mhi = lshift(mhi, Log2P);
  4022. }
  4023. for(i = 1;;i++) {
  4024. dig = quorem(b,S) + '0';
  4025. /* Do we yet have the shortest decimal string
  4026. * that will round to d?
  4027. */
  4028. j = cmp(b, mlo);
  4029. delta = diff(S, mhi);
  4030. j1 = delta->sign ? 1 : cmp(b, delta);
  4031. Bfree(delta);
  4032. #ifndef ROUND_BIASED
  4033. if (j1 == 0 && mode != 1 && !(word1(&u) & 1)
  4034. #ifdef Honor_FLT_ROUNDS
  4035. && Rounding >= 1
  4036. #endif
  4037. ) {
  4038. if (dig == '9')
  4039. goto round_9_up;
  4040. if (j > 0)
  4041. dig++;
  4042. #ifdef SET_INEXACT
  4043. else if (!b->x[0] && b->wds <= 1)
  4044. inexact = 0;
  4045. #endif
  4046. *s++ = dig;
  4047. goto ret;
  4048. }
  4049. #endif
  4050. if (j < 0 || (j == 0 && mode != 1
  4051. #ifndef ROUND_BIASED
  4052. && !(word1(&u) & 1)
  4053. #endif
  4054. )) {
  4055. if (!b->x[0] && b->wds <= 1) {
  4056. #ifdef SET_INEXACT
  4057. inexact = 0;
  4058. #endif
  4059. goto accept_dig;
  4060. }
  4061. #ifdef Honor_FLT_ROUNDS
  4062. if (mode > 1)
  4063. switch(Rounding) {
  4064. case 0: goto accept_dig;
  4065. case 2: goto keep_dig;
  4066. }
  4067. #endif /*Honor_FLT_ROUNDS*/
  4068. if (j1 > 0) {
  4069. b = lshift(b, 1);
  4070. j1 = cmp(b, S);
  4071. #ifdef ROUND_BIASED
  4072. if (j1 >= 0 /*)*/
  4073. #else
  4074. if ((j1 > 0 || (j1 == 0 && dig & 1))
  4075. #endif
  4076. && dig++ == '9')
  4077. goto round_9_up;
  4078. }
  4079. accept_dig:
  4080. *s++ = dig;
  4081. goto ret;
  4082. }
  4083. if (j1 > 0) {
  4084. #ifdef Honor_FLT_ROUNDS
  4085. if (!Rounding)
  4086. goto accept_dig;
  4087. #endif
  4088. if (dig == '9') { /* possible if i == 1 */
  4089. round_9_up:
  4090. *s++ = '9';
  4091. goto roundoff;
  4092. }
  4093. *s++ = dig + 1;
  4094. goto ret;
  4095. }
  4096. #ifdef Honor_FLT_ROUNDS
  4097. keep_dig:
  4098. #endif
  4099. *s++ = dig;
  4100. if (i == ilim)
  4101. break;
  4102. b = multadd(b, 10, 0);
  4103. if (mlo == mhi)
  4104. mlo = mhi = multadd(mhi, 10, 0);
  4105. else {
  4106. mlo = multadd(mlo, 10, 0);
  4107. mhi = multadd(mhi, 10, 0);
  4108. }
  4109. }
  4110. }
  4111. else
  4112. for(i = 1;; i++) {
  4113. *s++ = dig = quorem(b,S) + '0';
  4114. if (!b->x[0] && b->wds <= 1) {
  4115. #ifdef SET_INEXACT
  4116. inexact = 0;
  4117. #endif
  4118. goto ret;
  4119. }
  4120. if (i >= ilim)
  4121. break;
  4122. b = multadd(b, 10, 0);
  4123. }
  4124. /* Round off last digit */
  4125. #ifdef Honor_FLT_ROUNDS
  4126. switch(Rounding) {
  4127. case 0: goto trimzeros;
  4128. case 2: goto roundoff;
  4129. }
  4130. #endif
  4131. b = lshift(b, 1);
  4132. j = cmp(b, S);
  4133. #ifdef ROUND_BIASED
  4134. if (j >= 0)
  4135. #else
  4136. if (j > 0 || (j == 0 && dig & 1))
  4137. #endif
  4138. {
  4139. roundoff:
  4140. while(*--s == '9')
  4141. if (s == s0) {
  4142. k++;
  4143. *s++ = '1';
  4144. goto ret;
  4145. }
  4146. ++*s++;
  4147. }
  4148. else {
  4149. #ifdef Honor_FLT_ROUNDS
  4150. trimzeros:
  4151. #endif
  4152. while(*--s == '0');
  4153. s++;
  4154. }
  4155. ret:
  4156. Bfree(S);
  4157. if (mhi) {
  4158. if (mlo && mlo != mhi)
  4159. Bfree(mlo);
  4160. Bfree(mhi);
  4161. }
  4162. ret1:
  4163. #ifdef SET_INEXACT
  4164. if (inexact) {
  4165. if (!oldinexact) {
  4166. word0(&u) = Exp_1 + (70 << Exp_shift);
  4167. word1(&u) = 0;
  4168. dval(&u) += 1.;
  4169. }
  4170. }
  4171. else if (!oldinexact)
  4172. clear_inexact();
  4173. #endif
  4174. Bfree(b);
  4175. *s = 0;
  4176. *decpt = k + 1;
  4177. if (rve)
  4178. *rve = s;
  4179. return s0;
  4180. }
  4181. #ifdef __cplusplus
  4182. }
  4183. #endif
  4184. #endif