l coNTRIβU↑IoN TQ THE THEoRY oF NUMERIcAレ
INTEGRATION OF NON−LINEAR 、
、DIFFERENTIAL EQUATIONS（vl）ぺr l：
’ 『 ” BY
TATsuJIRo SHIMIZU t ．． t’』 一 Theory of士ound］㎞g errors iS not yet discuSsed f“11y． In 22 we shall show that the rounding errors． can be determined uniquely under somC conditions．． ． ． ： il123 we sha皿show some叫ethods available to obtain the exact rounding errors fbr additipn and for multiplication under described conditions． In 24 we discuss the true local roundi皿g errors and the so−ca皿ed local rounding errors・ ． ． ．・ In 25 we shall show that we obtain the rounding errors exactly on colnputer by single preciSion aritlmetic in numerical integration of differential equatioris． 22． Determination of ro皿ding errors ． ． ． Let the differential equation y…＝∫（x，γ）， and Euler’s methodγカ＋1＝Yn十hf（xn，γヵ） withアo rγ（xo）and the step size乃fbr the equation be given． We consider round㎞g errors fbr such a method in且oating point arithmetic． Denoti lg by・アπthe values obtahled by rounding yヵ99）（the true values fbエγ蕗＋1＝γ外一＋hf（x。，γ。））we．have ・ ’
ア・＋・一ア・＋hl（X。，γ。）， ア・−y・， where∫（Xn， Yn）denotes the value of∫（xn，ア。）computed with xn andア竺， consideri皿g rounding errors． Hence ・ ” アヵ＋1＝：V＃十hゾてXn， yヵ）十ρ蕗（1｝ 一 ． ＝アヵ十乃プてxヵ， ア鈴）十ρヵ（1｝十ρ（2） ＝ア外十乃∫（X”， ア鈴）十ρ蕗（1｝十ρヵ（2｝十ρヵ（3｝． 1皿these expressions・ノてxn， y外）can not be determined uniquely． It depends㎝t／he processes， by which r皿ndi皿g is perfbmed． ． Sin㏄round辻1g errors．have not yet been discussed，桓general， exactly， we shall discuss under what conditions roundmg errors shall be determined uniquely． ＊ 3eceived September 4，1974．， 88）ラ鈴are not exactly the values obtained by’ ro皿ding yカ． With respect toア． and jn see p．92． ［83］84
T．SHIMIZU
First we restrict to the arit−metic operation（addition， subtraction， mUltiplication or diviSion）performed only once． ［・41ユ RePτes◎ntation Qf耳umbers w地ba…埠（典dbζ）β・ Evi（1旬唾．耳o騨diPg errors ’depend on壌se一βarithmetiq． In thN fd叫qWing w皇郷u媒埠β手20rβ＝10（binary or d㏄㎞al n靱ber§yste叫 ［A2］F蝋po血t aritimetic or fioating po桓t arithmetic． In the肋唖g we assume 且oat祖9 Point number syst㎝， that is， x＝㎜βex， mx mantissa， ex exponent． ［A3］Let∫be the set of姐real n卿ers． ap．d R， be the set of machine nmubers： R−←1…ff・・舌くmx＜1・−P＜ex＜・｝・ In the fb皿owing we assulnC．that all given numhers belong to R．（Operated nUm− bers， fbr instapoes， x十ン， xy do not nQCe． ssarily belong to R）1一 ［A，］ The mantiSsa has’b輌nary（d㏄imal）digits． L45コ The roundi 19 methQds are not u旦iq嬰e・thqt is・tru皿cation（dlopPing or roUnd− ing down）， round㎞g up or rOun（ling Proper1ン 麺 ordinary §ense and others， and there 1nay be the面xed type of these methqds． In the fb耳owing we con− sider truncation or proper roundh19． L46］ There are two cases，垣one of which roUnding or truncation is perfbrmed a丘er cadl oP顕tio興nd then nom｛晦ation is per㊤mle4． In the otheτ邸e normaliZation is perfb血ed after each operation， and then roundmg is perfbrmed． In the fb皿owing we assume the second case only． ［．47］，Arittmetic operation is per．fbrmed in differellt ways by different computers fbr instance， different methods of multiplication to shorten the eXecuting t㎞e， methods of division w仙respect to the remamdeL For． subtraction two．㎞1ds of complem㎝lt are used． These methods do not give．exacdy the same round一 桓gerrors，（esp㏄血11y fbr subtraction see Knuth：Art of computation program− ming， Vol．2， p．197）． In the fbllowing multiplication is perfbrmed hl ordinary way’by the repetition of additiOn，．and． divisiOn retaining the remahlder w仙the salne扇gn with the diyis．Qr． In all calSes we qss．ume tha． ．t the Cρ其1p嘩ation is perfbrmed With numbers of the signed magnitude notation（w鯉hQ誕． Us桓g n硬一 bers of complement notation）．［A8］Some computers have hard warest pQrfqmiing double preciSion system and
etc． The round㎞g e宜ors卿be d征er賦．． In．the following we assume that s桓gle pr㏄ision syStem only．おPerfb頂ed、 ［Ag］ Different manipulations are perfbrmed by．◎omp嘩《？rs when exponent over丑ow or underflow has o㏄ured． The absolute zerQ and relat・ive zero arg． al「so Mallil）− ulated桓different ways． In the fb110wing． we assume that QXponent eve描ow， ご underfidw and the zeroes do not occ町in the◎ourse of◎omputatkm． Mantissa over且ow， under且ow or that caused by roUnding is manipurated almostiP the．．s姻e W by㎝y C・mPUt… ln．．th・fo11・Wing噸袈S皿鋤・t・・鋤n
’ ove｝且ow is manipulated．as on the ordinary computer． ．．．．∴ ’ 、t
CONTRIBUTION TO THE THE．ORYセ）F．NUMERICAL INTEGRA’1［EON
額 ［Alb］ There are twO eases． The lf迂st o血e is to peffofin the 樋thnletjC ep’i r，atiofi； considering the mantissa of’digits or manipulation of computer． The． Sε￠ond ’das6 姪 tO perfb】励 aritlmetie OperatiOn （血athe血at’iCa／皿y） WithoUt cOn8ider㎞9 computer， and lastly to give the血oati且g point expression with SOni6 humber of digits． In the fbllow垣g we perfbml arith血etic頭壱fation after．◎omputer manipulation． ［Ai1コ As：to the singl6 preCision syStem Some eomputerS hiWe single preciSion accumUlaters， and others have doUble preCision aeCU血ulaters（the roun（ling errors are diffefent aCCotd・i血g to the co血1）uterS）． Some COrbputerS， havi血g si血gle pr㏄isiOh aCCu血UlatOr，、血iay haVe ’a．dertain． number of guard d垣ts． Here we ass迦e that the computer has nO gUard digit．89） For the ratiQnal expression（arit㎞etic operations being perfbmed successively in af桓ite llumber of times）we must assume the fb皿owing． ［B，］ The order of arit㎞etic operotions perfbnned must be determined， si lce the distributive law and the associative one do llot hold fbr、 the computation， con− sidering rounding errors， and the roinding errors depend on the order of operations． ［瑚After ea（；h arittmetic operation the number is stored as a single precision floating−point number（the mantissa being of’digits）． Some computer put the intermediate resUlts in some registers to shorten the executing time． We do・not consider such a computer・［B3］ Some computers have sp㏄ial methods fbr the product−sum or repeated
multiplication， in order to be more accurqte or to shorten the executing t㎞e． We do not◎onsider such a special method， but ordinary multiplication and addition only； By［．41｝［An］，［B，｝［B3］we can detetmine’the rounding errOrs u血iquely． But in order tO obta血the roullding errrors 6n computer we must《）hoo…re some numer− ical methodl unde才the．cOnditiOns［A1］一［An］，［β1］一［β3］ aboVe血entioned．90） On the stand point to obtain the founding errors we muSt consider the fb11◎wing．［C1］We must choose FORTRAN， ALGOL or assembler etc． in order to輿fbml
the computation．’FORTRAN， ALGOL use two types integer−type and rea1−type． If we use both types in the course of computation， the rounding errors may depend on them（using integer−type to obtain some result other than computa− tion does not aff㏄t the round血g errors of course）． in the fbllowing we use FORTRAN or ALGOL and real type only． 89）W皿kinson：Roundi皿g Errors ln Algebraic hoceSses， pp．7−i1． 90）Without◎onsidering a new type computer in future， the roundi皿g errors are deter’ mhled uhi袈ely uエlder the￠OhditionS 241−、411， Bl−Bi R）r tlie cOmpUtef at preSeht．h 86
T．・SHIMIZU
゜・， ［C，］ Prope’．method must be chosen a㏄ording to the Cρ．n（Utiqns l［ム｝［A，i］，［B．，］ii［β31三’．1’． 『『 ．． 1』’ tt．一 ’・
We sh姐．・9ive here sgme examples◎onCerning to the ConditionS［A，］一［ムコ，1［B，］一 ［B3］，［C1］一［C2］． ・、It is well known that for uy v∈R 馳． ノ1（〃±〃）＝（u：ヒv）（1＋ε） 1ε1＝θβi−t，0≦θ＜1．91）． The above relation姪true under the conditions ［、41］一［An］，［Bl］一［β3］，［！15］．truncation，［Al1］dOuble precision accumUlator． ．．1．、
On the other hand， the above relation姪not true桓general uPder the conditions ［1【1］一［An］，一［Bl］二［β3］，［A5］truncation，［An］ single pre9姪ign aeCumulator，「．but it 誌true that ． ・ ． ． 「・・ ∵’ 一・ ノ7（u：ヒv）＝u（1＋e、）±〃（1＋ε2） 1ε11，1ε21＝θiβi←t，0≦θ＜1， 、・ ’・ e1，ε2 are differe皿t m general under the conditions given here． ExAMPLE．’ti4，［A5］truncation，［ノ111］single precision a《x㎜垣ator ’ ．『’u＝・1㎜×100 「 ’ 一’
’ v＝．9999×10−1 u−v＝・rooo×10−5， fl（u−v）．＝・1000×10−41馳 If fl（u−v）＝（u−v）（1十ε）， thenε＝9＞10−3， which is a contradiction． We have ／7（u−v）＝（〃−v）（1十ε），ε＝0，under ［A5］ truncation， ［ノln］ double pr㏄ision accu頃 mulator． 「 』 」． The method＞which gives the exact rounding errorsδr〃（u−v）一（u十v） f（）r addition is discussed by many writers． The method given by Dekker is the fbllowing：Let u， v∈R and ti≧lvl w＝fl（〃＋v） 、・ 〆−fl（w−u） ・ 、
δ＝fl（v−vt） which gives w十δ＝〃十v， andδis the rounding error． This is true皿der the condi− tions ［ノ41］一［Zl11］， ［β1］一［β3］， ［C1］一［C2］， ［A5］ truncation，’［An］ dbuble I）reciSion accumUlator． Computation of w−〃must be done with at least t十1digit mantissa． The method of Dekker fails under the c6nditi6ns［Al］一［An］，［B1］一［B3］，［C1］一 ［C2］，［／15］truncation，［∠111］ single pr㏄姪ion accumulator・・ 、． EXAMPLE・ t＝4，β＝2・ ” u＝＝・1111×20， v＝＝・1111×2−3， w＝・1000×2−1 nt−u＝・001×20 and w−u＞v． ． Hen㏄the ro皿d鉦1g error can not be obtained． The method given by Mφ皿er is the fbllowing：Let〃，≡R．and u≧lvl， w ・fJ（〃＋の． If ew≦eu， ew， eu being the exponents Of w and u reSp㏄tively， thpn 91）Wilkinson．1㏄． cit．‘‘fl”denotes且oatmg・point n㎎1快r with’digit mantissa．CONTRIBUTION TO THE T耳耳QRY O耳NUMERICAL INTEGRATION
87∴ 、、’．．，．・〆−fl（Wrの． ．・ ・・．：δr卿一の・1「1・．□ apdδ・おthe ro㎜ding．eロOr・ If ew＞切， then we must consider fUrther ・ ・． 一・． ． ．u’ ・fl（w一の vn＝fl（〃一の u”＝fl（u一め． ．
δrノソ（1〆’十t〆’）十．ノ’1（ρ一ノ7（2〆十ノ7（〃−2〆））），． andδis the roundhlg error．92） ． Mφll。， 9。ve th・m・th・d und・・thd。・nditi6・・（diffe・e・t−ff・m the c・nditi…gi・eP ・b・vの，［A、コbi・ary・血・ber W・t㎝，［∠・コt・unp・ti…［Aiコ「c・mp1㎝・nt「i・t・ti・n’ ［んコ・血91・p・ed・i・n accum・1・t…B・t we can・h・w th・t M輌’・m・th・d蜘eS th6 exact round口1g errors f（）r additi6n under． the conditions［ンll］一［ノ1．11］，［β1］一［β3］乏 ［C1］一［C2］，［A5］truncation，［ノlil］ single precision accumUlator・ W。，h姐即。 he，e・a m。・e・㎞P1・m・th・d u・d・・th・6・n砒i・i・［A・コー［A・・コ，［3・コー ［β3］，［（な］一［C』］，［∠45］truncation，［∠tll］s垣91e precision accumulator・ Let〃， v∈」R and u≧lvl w・−fl（u＋v）． If ew≦elちthen v’＝fl（w一の 』 層 δ＝プ7（〃−t〆）． ． 、 If e｝〃＞e〃， eμ一ev＝K＞O then 『 一 ・ ． ・ 』． 〆’＝ノ7（〆−z）， where z＝．〆7（〃−fl（ノ7（4＋β・・）一βeu））， eu． be血g obtained by◎onsiderhlg logarithm of づ u ＝rnuβeu， and δ＝∫fl（v−t〆’）・ 23．Ro皿ding errors for addition． a血d for伽晦1icati皿 In this 23 we assume th。t［A、コ且・ating P・int aritluneti6 sy・t・m［A・コthe set・f machin・・umb・・s ・nd ・7 v∈R・［A・コt・digit m・nti・…．［A・コ・i・gl・p・ecisi・n、system・ qnd u≧1び1・ w．Kahati gave a method to obtain the exact rounding e町orsδfbr addition w十δ＝u十び⑩be桓g the computed value）， which is equivalent to the method given． by D・k㎞血22． H・did n・t・give any・9・diti・n（・xcept［A・コ）皿d・・whi・h．hi＆ method is avaiable． O．M￠ll・・ 9・v・｛h・．f611・w血9・m・th・d t・・bt・in・th…xact・。unding err・・s a f・・ addition w十δ＝u十v， under the cond五ions ［A，］base−2 arithr etic［A5］truncatiQn and［、47］oomplement notation三fbr s曲traction． He gave w−fl（叶の， vLfl（w一の・． u’efl（w一のヒ・‘ −92）Dekker， M姐1et． See．23， p．92． i 、・ 、00 、8
T．S㎜U
1〆’＝＝ノel（〃−t〆）， tit「’三ンel（a一二友り， δ’＝ノ7（1〆，十t〆つ δ＝fl｛ノ7（〃’十1〆ノ）十ノ7（ρ一ノ7cρ’十2〆ノ））｝ and applied his method to reduce the errors血numerical solutions of d遜erential equations． ． ： ’ He did nOt remark that hiS methOd iS aVa丑able←r．U’ ndef the COnditiOn［All］Singleprecision accumUlator． ’1馳’
D．E． Knuth gave the following method w・・fl（u＋の，’〆⇒ノ’1（w 一一 u）；ti一fl（w．−v’） Z〆’＝＝ノ7（V−t〆）， 〃”＝ノ7（〃一〃ノ）， δ＝ノ1（1〆’十t〆’）， ！ under the Cond冠io血s［AI］．base−20r base−10 aritlmetic ［A5］ proper rouqd並ig and ［ノln］double pr㏄ision accumulator． ・ ・ His method is．not『uMcient when the◎omputer has the single precision accい tnulatof， but is su伍cient when it has double precision accumulator． LBabu§ka gave the method equivalent to the血ethod given by Dekker， without ．stating any conditions under which the method is available． He gave a formU la which sum up the roundhlg errors when addition is perfbrmedヵtimes． ． ．’ E．Vitisek applied the method of Bψu§ka to redu㏄the round㎞g errors血numer− ical solutions of diffierential equations．‘ M．Pichat gave the method equivalent to She method given by Dekker under the conditions［An］single precision accym． Ulator with one guard digit（base−16 arithlnetic） and［．45］tmncation［、47］complement notatiOn． Yasuhiko Hirano also gave the method equivalent to the method given by Dekker， without stating any◎onditions皿der which the method is qVailable． Dekker gave two theorems： ． For base−20r basek3 aritimetic and［A5］proper rounding the method given h1 229ives the exact roundi皿g errors fbr addition． For base一βarltlmetic（βbeロ1g arbitrary）and［A・］．truncation the method given in 22 gives the−exact roun（1ing errors fbr addition． He did not give e文plicitly the conditionS under．which the method is ava五able， but in his．proof we can．see that the method is available when the co症1puter has dOtible preciSion． accumUlator or at least single precision aocumulator with one or two guard digitsl ・『． ： As shown血22 M￠ller，s meth6d is the on］y one， which is available when the com・ puter has single preciSion accumUlator Without guard digits． It seems to us that all the above nlethods， co血sidering the conditions under which they are ava丑able are not fu皿y discussed and not exactly used．In the following we shall give here two theorenis： 一
THEoREM． By the血ethod given by Kahan
1〃＝ノ7（u十〃）， t〆＝ノ7（w−〃）， δ＝ノ7（〃一〃’） we can obtain the exaCt roundilg errors 6 for addition w十δ・＝u十v， under the condトt
CONrRIBUTION TO THE T宜EORY OF NUMERICM INTEGRA皿ON
tions［A、］base一βarithmetic（βbeing arbitra ry）， tA，コttrencatiofi or proper rO’and一 血9 （eμ一ev ＝K＞O if proper roundi血g and ew＞eの， and ［An］ s垣gle precision accumUlator with one guard digit， PRooF． Let ul，〃2，…， ut、and vl， v2ピ・r， vt be the mantissas of u and v respec− tively， andρ一q＝K， P and q be the exponents of u and〃．respectively・ First we assume． that、［．45］trUncation， and no mantissa over且ow．does not−occur・， and consider addition（subtraction）． UIU2 ・・．〃石【＋1〃1【＋2．・・Uξ ± Vl V2白・◆i〔・Vt−k Vt．＿亙＋1・ふ◆抄t ・ WlW2°・°Wk十1°°、°°°b°WtThen
mv’−mw−mu−±0∫を、0〃、v、…v，−K．Hen㏄
．±m“一（±血が）一±o蜘仇．冗．、…vi． When mantissa overdow oecurs we have UIU2 ・・°… UK＋1”X十2・・◆Ut →一 砂1 V2◆・・… ρ卜互 Vt．−K＋ゴ・・Vt 1伊oWI｝レ2°・°ハ伊K十1°°・も、°°“Wt wbei皿g of’digits， we have ．WoWIW2ぺ．・‘・・．6W｝＿1 − UIU2 ・◆・．・．◆Ut一IUt m〆−0｛8惑0“、v、…Vt．K．、“、−K−03…10ぬ’Hence
励一m〆−0・烈・Ovトκ．i…“汁O・｛fZOw， δ＝0・（f！（ WtVt“K＋1…仇 mw＝WoWIψ2’”｝レま一1． Under the condition［．45］proper rounding the proof is the same as aboVe， When rounded down． When rounded up we are to put Wlw2…｝Vt−lw，十 1 m〆一±ol写∼Ov、〃、…vトK．、VI．π＋1 ±m〃一（±mめ一．．±0∫ξ10〃，．κ。、…σ声0恩）・01、 when mantissa over且ow does not occur． “Wh㎝血antiSsa o寸er且ow oCcurs we have 『； Wow1｝v2…wト2協一1十1 卿’−oJ烈0“、“、…Vt一冗＋1−O・（P・ow， ， 鋤一m〆＝0・（氾）・OVt．K。、…Vt−0∫ξLOwま一〇・｛ξ二’？．IO 1． Since K＞0， m〃−m〆6an be calculated completely， and δ＝0・（ξ二1）・01十〇∫ξ≧0｝晦仇＿fi＋1…Vt． ‘．： 1「mUt ・＝｝WoWlw2…Wt−2Wt−1十1． 、．90
T．S．HIMIZU
・． sHEoRBMゴ．By the method gvCn in 22 ・ ．「 』』・ 1，11：．∵1・．・ ’ ・ ・、1．・！・−fl（u＋の・．v’ ・＝fl（w−〃）・δrfl（v−〆）．一一． if ew≦e〃． If ew＞eμand eμ一ev＝k＞0， then 1〆’＝．fl（〃’−z） ．z−＝∫fl∼u−∫1（ノソ（u．＋βeり一’： Beり） and δ＝ノ7（V L t〆’）・ 9㌧ ． ．．． we can obta血the exact roundmg errorSδf（）r additiori w十δ＝ti十v’under the condi・ tions［IAi］baSe」βarith血etic（βbei皿g arbitrary）［ノ15］tr皿cation［An］ single pr㏄ト sion aoeumUlator． ， PRooF． If e〃≦e〃、 the proof is the same as above． If・ew＞eμas shown ill the proof of the theorem given above the last digit ut Of m〃fbr／7（u十のis dropped by s垣gle pr㏄ision accumulator． ． ・ Since K＞0， wo＝1， wl＝O and mw−m4 is ofτ一1 digits． Hence fi（〆−z）is nothing加t vlv2… ρt．＿κ．−OSP・owt ．．． if［A，］ trullcation Thus the proof is the same as above． ’ Dekker gaVe a methOd tO Obtain the．eXaCt rOUndi皿g errOrSδf（）r mUltipliCatiOn w＋δ一・×v（wb・ing th・valu・ ・gmpPte．d）・ But his method depends on the sp㏄ial properties of base−2、 arithmetic and of proper rounding． For any base一βarithmetic．and fbr tmncation we shall give here amethod to obta血δwhere｝v十δ＝u；×v。、 First let’be even（’being the numbeエ． gf digits of man砧sa）， then we take the upper’／2 digits of the mantissa of u aロd of． v as the head of u， denoted by hu．，．and the head of v， denoted by h〃rCsp㏄tively・ u−h酩and v−hv are denoted by t〃and t〃resp㏄tively and are called tails of u and of v． Now h〃×hv， hu×tv，加×h〃，、tu×佃all 1）elong to’digit mantissa．We put ．． ’ ‘
、ρ＝ノ1（hu×hの q＝ノ7（ノ7（hr×tの十ノ7（tu×hの） ’＝ノ7似×tの ． then hμ×tv十伍×hv may．be g￡’十1 digits． If the computer has single pr㏄ision a㏄u皿ulator with one guard digit（or double preciSion aocumUlator）， then we obtain the last digit．γ10f h〃×．−tv十tu×hv by cpm− putation to obta血the rounding erτors fbr additiOn meptioned above， andγ1≠O orγ、＝0． ．．・．．・．．・
Computhlg z＝．ノ’1（P十∂by single pr㏄ision arith血etic we have．δ’r〆7（刀（P−z）十の ＠being of’digit mantissa）． 1． ・ ． 、・・ Sinceδ’十r may be of’十1 digits， we can obta血．the・・1ast digitγ20fδ’十r by computation to obtahl the rounding errorS fbr addition mentioned above． iCONTRIBUTION TO THE THEORY OF NUME】紅CAL INTEGRATION
91 Thus we haveδ＝ノ7（fl（δ’＋り＋：γ1×βP−r（峠÷−1））＋・γ2×βρ一ぽ一1）．
…’・‥d，・h・n・We・・ke・h・upP㏄τ
P1 d垣・…m紐・iss・…一・…ぼ
’the head of．㍑and the head of v resp㏄tively． In this case hμ×h〃may be of’十1 digit・・H・nce， w・m・・t’bb煤Em th61・・t・digit・r・、 Of．
Dby・hg垣the卿記、・f頑座：
plication，−that is，．we must compute hu×〃i， h〃×〃2，、．…，加×．び坦and h〃×ρ1十h躍．×倣 カ
†…＋』×〃乎・Th・醐・nm・・tb・繭m・⑪・Wt』．・0輌9娠S：劔
addition， thel1γ0≠O orγ〇−0， where v＝’・ vlv2・づ担．×βeワ，’．bi，、 v2’．↓．・r．，’ρ担1）e垣gl工hε ’ロ ド て ttニ man白sa．． and． e〃・being the Cxponent．of v・．．． ．．∠．．㌃．．．：．・・’ 元．．∴‘：．・ ・』∴r：N・WW・p噸・曲・v・． 、、．．，、1’i．1．…：≡．、．．
ク＝ノ7（hμ×hの， q＝．fl（hu×tv．十t況×hの，、 r≒fl（tu tt tv）；、” ’． ． 1∵’ then hu×ω十伍×h〃may be of’十1 digits． If the computer has single pr㏄ision accumUlator with on『guard・．digit， then we obtam the last d埴tγ10f hu×tv十加×h〃by computation to Obtain the rbt血dihg errors fbr addition mentioned above， andγ1≠O orγ1＝0． Computi皿g z＝ノ1（）十∂by Single pregision．aritimetic we haveδ’rfl（fl（P−z）十・蚕）・ （亘b・ing．・gf・．t． d寧m・ntlssa）・・’ ・ Sin㏄δ’十r may be’of’十1‘digits we can obtain the last digiレγ20fδ，十’−btt Cqm・・…i・n．t・gb・・血th・’・…S f・・additi°n ment’°ned ab°ve・・，． ll Thus we have馳 』 ．． ・一双刀（・’＋・）＋・γ・×β・−t†…×β・三（t＋号一’））＋…×β・一⑭・．．・・．∵’ The heads of u and of〃and the tais of u and〃can be obtahled by the fbl一 オ 10w血g ways． t エ ロ オ tail of u・一 fl（mu一ノ7（刀（mu＋・1×β丁）一・1xβ丁）） when’odd， ．、 t t tail of u＝．〆7（mu−fl（刀（mu＋・1×声）一・1×β司） when’ev㎝， head of u＝〃−ta丑 of u， so similarly fbr ta皿びand hea（I gf〃． v、，v、，…，坦can be Obtained by the fbH・wing ways． ’ ニ v1＝ノ7（fl（m〃＋・1×βま一1）一・1×β’−1）・ 巧＝ノ7（fl（m4＋・1．×βt−2）一・1×βt−2） ・v・rfl（陸一の V3＝ノ7（刀（mu＋・1×βt−3）一・1 ＞f igt“’3） 〃、rfl（Pt3−〃、）．、 ’ t ■ ● ● ● ● ● ． ・ ． ． ． ● ● ． ・ … ■ ■ ● ● ● ● ● ． ． ・ ・ ・ ・ ・ ・ … 廻イ1（防＋1−v卜1）． 2 ． 2 2 、 v ．t r ’92 T．s’H’㎜蕎
：R直FERENCES
D．E． Knuth：W，Kahan：
O．Mφ皿er： 1．Bab．nsk江： E、V諏s6k： Futther remarks ’on reduc血9 trlmeation errors， Co㎜． ACM，8，1（19651t 40・ Quasi double pτ㏄ision i血且oa血g pohlt addition， BIT（1965），37−50． Art of computhlg pτogram血ing， V61．2（1968），20i−204．皿mp CongtesS，68（1968）．． ． ：
The nU血eriCal s毎b皿ity垣solution of different｛al equations． Con㊤rellce《）f the numeriCal solutio皿of differential equations（1969），101． ． ・ T」．Dekker：A且・ating−P・垣t tec㎞iqu・fb・extending th・qvail・b1・p・㏄i・i・n， Numer・ Math．，18（1971），224−242． M．Pichat：Correction d’皿e somme en arithmetique a virgule fiottante， N迦er． Math、う 19（1972），400−406． Yasuhik・肚an・・R・und・ff・err・・s m h・t・9rati・n・f・・dinary differpnti・1 equati・n・（in Japanese）， Jyooh60 ShOri，14，2（1．973），95−97． 24．On local ro皿ding errors We consider the rounding errors負）r EUler’s method． Since ro1皿d血g errors are treated roughly血ordinary books， we sha皿consider as fb皿ows． The truncation e汀ors fbr the equatio11γ’−f（x， y）by Euler’s method万＋1＝み 十hf（xn，γ悟）may be expresses byξらwhere y（X＃＋1）＝y（x．）十抑（x，）），ア（xヵ））十ξ竺，γ（xn） being the true solution of the equationア’＝ノてx，ア）． ξπis called the loCal trunca− tion error． It is the difference ofア（x＃＋1）andγ％＋1＊＝γ（xn）十hf（xn，ア（xヵ））starting at the pointツ（xヵ）． Now we shall consider the roundhlg errors． Letγカbe the true solutions of the d征eren㏄equation−that is−（Euler’s method）． The true local rounding㎝三〇rsηπcorrespoding to the local trumcatiOn errorsξヵ are the fbl10wing．93） γ。．・＊＊一ア。＋碩κ。，ア。） ＝ア。＋hf（x。，γ。）＋η。， where yn denote the rounded values of the true solutionsγヵ94）（ア． be血1g Of i lfinite number of digits， yカis the伽es tOtmded to t digits nU血ber）and元形・・denote the rounded values of f， hf… ．『．う7ヵiS the diffefence ofア外＋1 and Yn＋1＊＊． S祖㏄we can not obtain exactly y．， we￠an．．not obtam exactlyア。． 且en㏄the Iocal rounding errors are consklered to be the fbllbwhlgρヵmordinaly books． ア。＋・一ア。＋乃∫てx。，ア。）一ア・＋hPt（Xs， Yn）＋ρ・・ whereア． denote the values obtained by the fbllowhig recurtent relatiolls． ア1＝γo十砺（xo， Yo） ifアo＝アo otherwiseア1＝アo．十乃ノてxo，アo）・Thus
ア・一ア・＋h7てX、，ア、），ア・，…，ア・，…・ 93）We assume thatヌis uniquely defined． Yn andアπare not equal kl generaL 94）ノてx”，yヵ）can no†be exactly computed， so we takeアπfbrア外．CONTRIBU皿ON TO THE TH聰躍C暉NUMERICAL INTEGRAT ［ON
93． ア露＋1can be obtained on computer・ Yo＝Yo 烏菰
ri ： ￨’」’ 一yぬ珪 二苫．．・ ＼ h ． 噤f ’ @ ’ @ ’ @ ’ﾏ1／ ％∠’ ＿・ ｨ ，’／ ’ 潤^ ツ〔。。）・象イ㌻ 乞 ﾑ． @ 〉 η暖．先＋1簗πy︵㌔n ツ y（エ2） ∀ 頚工1） 〉 ’ へ ’1 r ラ 25．Computation of tOunding． errors We sh訓1 give here exact evaluation of round垣g errors． In order to dQ敏m祖e the rounding elTors uniquely we qsSume．Φ邸，［ A．・，］一［、411］，［B，］一［」B3］，［C1］一［C2］， ［A5ゴtn皿cation，［、411］sillgle preqision a￠田皿Ulator’with one guard digits， and the n㎜ber of digits of m㎝tissaオis⑰（of Q。Ur’唐f we can assume that［∠11］singlep⑩lon㈱㎜司atpr and Qdd t， but thc餌Q曝of cg⑩tation b㏄o鵬s sqme−
what complicated）。 As shown m 23 the round辻lg er葺orS fol頭d垣on胆4 f6：multipHcatiOn’Can be obta麺ed胆der thQ cond．itiop．s above・95） Fqr・m辿iゆ1おation ノ『Z（μ1μ9＝μ142十δ1 fl（fl（u・μ・）ぬr∫7（μユμ2）U・．＋δ・＝」脳2μ・＋a・U3＋δ・ fl（fl（fl（UIU・）翼淑・）弓U・U・U．3．U・＋δ・U・U・＋δ・U・＋δ・ ’◆°°‥’．’．’一．”°’’’’”◆°’”°”＾” ・ 〉’’”‥°°°°¶°’°°”°’’’”．’”°”． For ad戯iq亘 fl（U、＋μD戸μ、＋n2＋δ・’ fl（fl（U、＋U、）＋U、）＝−fl（U、＋μ・）＋U・＋δ・＝μ・＋U・＋〃・＋＆’＋ゼ fl（ノ7（fl（〃1十〃2）十〃3）十U4＞一＝＝〃1十U2十U3十U4．十δ1’十δ2’十δ3’ ■ ■ ◆◆ ● ● ● ● ●● ■ ・ ● ． ・ ・ ・ … ． ・ ？ ． ． ・ … ◆ ■ ■ ● ● ● ● ● ● ● ● ● ● ● ■ ● ● ● ● ・ ・ 白 ・ ・ ・ ・ ・ ・ … wl r・δ・，δ・，…，δ・’． a・’，・…⑪dmg箪grs・e・p㎞申・96i N・ww・・h・11卵・am・th・d t・・btaiP．・g皿φ・騨r㏄・Qxactly・・d・・蝕・d wh・n we computeア％byア．．1＝ア弘十hf（x．，ア。）十ρ外，ア頒the、 C◎頁Ψ砿ed values andρ。 the 95）By・ingl・pP㏄i・i・n・y・t・m th・…nding C・…S can・be 6btai・・d・㏄砥・tely・・far・・⑩i麟． 』 ． ． ．． ．
96）These．eπ・rs are・equivalent・t・th・se given by・ Wil・kins⑨頁（1・。・cit・）but・we．・ean c・m・ P幟．cxactbe t．he8e e廻r◎彰S． eln◎Q頓PU鯨． by tb金above、輿e隻hρd竺 ．：194
T．SHIM］IZU
so・caned local round垣g errors． ．・ ． 、 ・ 1 ．． ． In order to be able to define the rounding errors uniquely we assume thatノてx，γ）i・ap・1y・・mi・1・（・，y）．97・∵ ’ i
We take an exampleρ（x；y）−ifx‘十aツx2十bx，os）the entors bei 19 able to be obtained si荘血arly fbr the general p（x，γ）． ρ（為夕。）＝＝ア。3Xn＋妙み2＋bx． x。，ア。∈R 4（ア・3・・）一ア・3・・＋ρ・ア…†ρ…＋ρ・’ ノT（砂み2）＝妙郷X”2十ρ4Xヵ2十ρ5X餌十ρ6 fl（bx。）⇒㎞。＋ρ7 whereρ1…ρ7 are ro皿ding errors， and乃， xカ， a，6，γo∈R． Repladngρ1アpt鵬，ρ2xカ， P4Xn，ρ5xヵby／7（ρ1ア＃xn），ノ7（ρ病），…from fl（fl（ρ・ア。）X。｝一ρ・アみ＋ρ、’X。＋ρ、’， fl（ρ、’X。）一ρ、’X。＋ρ、” fl（ノ「」（ρ4Xn）Xヵ）＝・ρ4Xπ2→一ρ4’X”十ρ5’， ノ7（ρ4’Xヵ）＝ρ4「Xn十ρ2” ’ プ7（ρ2Xヵ）＝ρ2Xn十ρ3’， ノ7（ρsXπ）＝＝ρ5X外十ρ6’，we obtain
：→1 ’、、 ．』 『』 一’ ノel（アn3Xn）一アn3x”＋ξノ ． ノz（aアhXn’2）一のヲnX．2＋ξ外” ・ ノ7（bxヵ）＝bXn十ρ7， whereξノ，ξπ”are linear combinations ofρ1，ρ2，…，ρ1’，ρ2’…，ρ1”，…and fl（ρ， アバ蕗），ゾ7（ρ1’xn）， fl（ρ4x蕗2），…． Adding these three temls we have l‘ ’ ノ7（ノ7（ア。3x。）＋ノ7（のアPtヵ2）一ア蕗3Xπ＋ξヵ’＋の戸ぱn2＋ξ．”＋ηn’ ． 1 ’・ ノ7（fl（ノ7（Y”3Xn）＋ゾ1（dy−rpc．2））＋fl（bx認））一ア外3X．＋ajnX”2＋bnX＋ξヵ’＋ξ．”＋ρ7＋ηカ’＋ηカ” whe「e OP・t・η・”ar・th…un・血9・rr・・s bロ“diti・n・Th・t・㎝・ρ・…ρ・’…P・”… ノ1（ρIYnXヵ），…ξ霧’…ξ外”… can be classified in the groups of tenns of’ digit mantiSsa 孤d・fexp・n・・t・fTpM P．一（ノ’＋1）t・ρ一（j＋1）’《ノー0・1・・2・…）・e・pecti・・1y・P−・be− ing the highest exponent among the temls considered． We must sum the telms of the same exponents−group reSpectively by the methdd above mentioned， that is， the method to obta血the rounding errors fbr addition．』 qepeating the addition several times the tota1 rounding errors wi皿be oblained．、Denotmg p（x”， Yn）＝＝pn十ω路， where P・・ a「ethe v垣｝es．？bt・in・d・…Pap・t…㎝dωr th…mPding・π0・S．・ Now fbrπ」0， IS 2，…we Proceed as fb110ws γFy。＋hp（X。，γ。） y。＝ア。 夕・＝ア・＋hp。＋乃ω。一δ。 ア2＝ツ1＋妙（x。，ッ、） ＝・n・十δe十hρ（Xl，ア1十δ0）＝ア1十δ0十hp1十hω1』 ・ ’・・’・f”．『
@“”ア、一ア、＋δ。＋hpi＋乃ω、一δ、 γ3＝γ2＋伽（X2，γ2） 、 一ア・＋δ・＋hp（X・・ア・＋δ・）一ア・＋δ・＋hp・＋乃ω・ ．一．．’．一・・− 97） If∫（x， y）iS a transcendendal function， the a㏄uracy of the table used or of th6 trun− ・・D．’Dcated．tems：and．etc． must be ・ considered． 、 ； ．∴』、 ．．：．．’tt 98）The order to be◎omputed of lnultiplicatio皿and・addition must be giVen五τst．’．CONTRIBUTION TO THE THEORY OF NUMERICAL INTEGR．ATION
95
ア3一ア、＋δ、＋吻・＋乃ω2一δ2 ■●●●●●●●●●■d●●●●●●●●●‘●■●■．・・． ボ●・．・・．．・．．‘．・．・・『．‘… ．‥・・．・● ’ ・ ：A6 n血creases the terms bf’digit ’mantissa and of exponents from p−（nt十1）−to lρ一（・＋1）t・in・・ea・e・， ahd the c・m輿t・ti・n b㏄・m・・v・ry・・頑・qt・9・』 But by the above method all rounding errors can be obtaihed exactly as’desired at・1east・theoretica11y．ぷ 』 『 ・ Since：the abbve niethod「9ives the com餌ted values aiid：the roundhlg errors， the valu・・ア。， P・，δ・， ce・・ar・th・v・1…and th…und血9・r・g・s・bゆ・dl99） The a㏄umulated round祖g errors can be obta桓ed as fb皿ows． ・By recurrencerelations we have・ ・ ．＼ ・’．、
ア、一γ。＋hp（x。，ア。）一δ。・ifア・一ア・』． ・ 』・・㌘ ＝yo十吻（xo，γo）十アo一γo十h｛P（xo，アo）一」ρ（κe，γo）｝一δo if yo≒kア6…1 ＝γ0十ε0十εeqo一δOwhere
ア。一γ。一ε。，q。−h｛P（X。，ア・）−P（X・，ア・一ε・）｝・Thus
ア2一ア、＋（ε。＋ε。4rδ・）＋（ε・＋ε・q・一δ・）q・一δ・ ア3＝γ3十｛（ε0十εOqO一δ0）十（ε0十eOqO一δ0）q1一δ1｝ ＋｛（・・＋・・q・一δ・）＋（・・＋・・q・7δ・）q・一δ・｝σ・一δ・・ Putth19 ア。＝み＋ρ。 夕。＋、一ア。＋、一ρ。一ρ。q。一δ。≡（2。＋・ （∼”＝ρ”＿1十（］，三＿lqn一δ露＿1 q。−h｛P（X。．、，ア。．、）−P（X。．・，ア。．・−0。．・）｝ ＝h｛−3ρ。．、Vn．、2x。一、＋3ρ。．・ア。一・x。．一ρ。一・3Xn−・−aOn−・x・一・2｝・ Since the terms（2。アr−−12x亮一1，3ρ．−1り％．Ix．−1，…are equivalent hl fbrms to the terms ア。・x．，卿。・，…・b・v・m・nti・ned， th…皿ding・rr・rs can b・・bt・in・d・imilarly・・ befbre． Now we consider a computer with eight decimal di it mantissa and double preciSion accumulator． We assume that 10−8＜固， lbl，…，国，1ッ1，…，＜102． and me十〃o＜10，． K K十1＜10，whereピ・ys・is the teml of the highest degree of」ρ（x， y）＝Σα，xゴ均乃”」． フエ The rounding errors do not exceed 100∼1000 t㎞es the・quadruple precision temls， that is， the terms of exponentクー4’． The total rounding errors can not reach thedouble pr㏄ision tems， though we computeア”up to〃＝10000． If we compute
exactly with triple pr㏄誌ion by．the above s五1gle pregiSion arithmetic（that is， the rounding errors are less than 10ρ一24，ρbehlg the highest exponent of the terms 99）Si且ce double pl㏄ision（triple precision）system is eas皿y perfbmed on recent◎om・ puters， the single precision technique seems to be of litde use． But if we apply the same technique for double pr㏄ision system， then we Can obta垣quadrup1C‘pr㏄iSio耳 system and so fburth． ．・：．：％ T．SI迎MI2日」・ considered）then we have’ ．、 1ア。−y。1＜10♪−1C．．n≦10000．． Asomewhat practical method is as fb皿gws．．．F畑t． wQ cQmpute the rounding・errors up tQ the 9【oups of the㊤rms of’dig t ma邸三ssa and of expollent from P−（’十1） to P−2’ fbr 姐addition・and、 mU ltiplication・ Then，ω弧in the expressio皿 ク（x，ハ， y“） ≡ρ垣十ω蕗aEe of double pr㏄姪ion． The computation ofアn is perfbrmed such that／7（み十乃幻ヵ）＝ぷ外，ノソ（s．−y．）＝・勿汚， ノ〕（h㌧←＝みσヵ）＝δ潟，and？∼ヵ十δ牧iS the r◎qU鉱ed SOIution〔 It誌of double precision but the rounding㎝rC曝s can nOt be．ぱactly est加ated亀 As皿ple method is as fb110ws． First we compute the rounding errors by simple pr㏄ision system andω鰺in the expression p（xカ，γ，．）−p．十ωカare omitted． Consider− ingび。＝Σδ， we considerア。十〃． in single pr㏄おion system． This method was given ごヨロ by E． Vi縞sek．100） ﹀ 100） Vitas6k． See 23， p．92．
