123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817 |
- COMMENT Module for Arbitrary Precision Real Arithmetic;
- SYMBOLIC;
- COMMENT *** Tables for Bigfloats ***;
- GLOBAL '(DOMAINLIST!*);
- DOMAINLIST!* := UNION('(!:BF!:),DOMAINLIST!*);
- PUT('BIGFLOAT,'TAG,'!:BF!:);
- PUT('!:BF!:,'DNAME,'BIGFLOAT);
- FLAG('(!:BF!:),'FIELD);
- PUT('!:BF!:,'I2D,'I2BF!:);
- PUT('!:FT!:,'!:BF!:,'!*FT2BF);
- PUT('!:RN!:,'!:BF!:,'!*RN2BF);
- PUT('!:BF!:,'MINUSP,'MINUSP!:);
- PUT('!:BF!:,'PLUS,'BFPLUS!:);
- PUT('!:BF!:,'TIMES,'TTIMES!:);
- PUT('!:BF!:,'DIFFERENCE,'TDIFFERENCE!:);
- PUT('!:BF!:,'QUOTIENT,'BFQUOTIENT!:);
- PUT('!:BF!:,'ZEROP,'ZEROP!:);
- PUT('!:BF!:,'PREPFN,'BFPREP!:);
- PUT('!:BF!:,'SPECPRN,'BFPRIN);
- COMMENT SMACROS needed;
- SYMBOLIC SMACRO PROCEDURE MT!: U; CADR U;
- SYMBOLIC SMACRO PROCEDURE EP!: U; CDDR U;
- SYMBOLIC PROCEDURE I2BF!: U; '!:BF!: . U . 0;
- SYMBOLIC PROCEDURE !*RN2BF U;
- BEGIN SCALAR X;
- X := GET('!:BF!:,'I2D);
- RETURN APPLY(GET('!:BF!:,'QUOTIENT),
- LIST(APPLY(X,LIST CADR U),APPLY(X,LIST CDDR U)))
- END;
- SYMBOLIC PROCEDURE !*FT2BF U; CONV!:A2BF CDR U;
- GLOBAL '(!:PREC!:);
- SYMBOLIC PROCEDURE BFPLUS!:(U,V);
- %value is sum of U and V, or zero (NIL) if outside precision;
- BEGIN SCALAR X,Y;
- X := TPLUS!:(U,V);
- Y := '!:BF!: . ABS MT!: X . (EP!: X+!:PREC!:-1);
- RETURN IF LESSP!:(Y,ABS!: U) AND LESSP!:(Y,ABS!: V) THEN NIL
- ELSE X
- END;
- SYMBOLIC PROCEDURE BFQUOTIENT!:(U,V);
- DIVIDE!:(U,V,!:PREC!:);
- SYMBOLIC PROCEDURE BFPREP!: U; U;
- SYMBOLIC PROCEDURE BFPRIN NMBR;
- %prints a big-float in a variety of formats. Still needs work
- %for fortran output;
- BEGIN INTEGER J,K; SCALAR U,V,W;
- NMBR := ROUND!:MT('!:BF!: . NMBR,!:PREC!:-2);
- IF ZEROP!:(NMBR) THEN RETURN PRIN2!* '!0;
- U := EXPLODE ABS(J := MT!: NMBR);
- K := EP!: NMBR;
- IF K>=0 THEN IF K>5 THEN GO TO ETYPE
- ELSE <<V := LIST('!.,'!0);
- WHILE (K := K-1)>=0 DO V := '!0 . V;
- U := NCONC(U,V)>>
- ELSE IF (K := ORDER!:(NMBR)+1)>0
- THEN <<V := U;
- WHILE (K := K-1)>0 DO V := CDR V;
- RPLACD(V,'!. . CDR V)>>
- ELSE IF K<-10 THEN GO TO ETYPE
- ELSE <<WHILE (K := K+1)<=0 DO U := '!0 . U;
- U := '!0 . '!. . U>>;
- BFPRIN1(U,J);
- RETURN NMBR;
- ETYPE:
- IF NULL( CDR(U)) THEN RPLACD(U , LIST('!0));
- U:= CAR U . '!. . CDR U;
- J := BFPRIN1(U,J);
- IF J=0 THEN <<PRIN2!*("E " ); J:=2>> ELSE
- IF J=1 THEN <<PRIN2!*(" E " ); J:=4>> ELSE
- IF J=2 THEN <<PRIN2!*(" E "); J:=0>> ELSE
- IF J=3 THEN <<PRIN2!*(" E " ); J:=0>> ELSE
- IF J=4 THEN <<PRIN2!*(" E "); J:=2>>;
- U:=EXPLODE( K:=ORDER!:(NMBR));
- IF K>=0 THEN U:=CONS('!+,U);
- WHILE U DO <<PRIN2!*( CAR(U)); U:=CDR(U); J:=J+1;
- IF J=5 THEN <<PRIN2!*(" "); J:=0>> >>;
- RETURN NMBR
- END;
- SYMBOLIC PROCEDURE BFPRIN1(U,J);
- BEGIN SCALAR V,W;
- IF J<0 THEN U := '!- . U;
- %suppress trailing zeros;
- V := U;
- WHILE NOT(CAR V EQ '!.) DO V := CDR V;
- V := CDR V;
- L: WHILE CDR V AND NOT(CADR V EQ '!0) DO V := CDR V;
- W := CDR V;
- WHILE W AND CAR W EQ '!0 DO W := CDR W;
- IF NULL W THEN RPLACD(V,NIL) ELSE <<V := W; GO TO L>>;
- %now print the number;
- J := 0;
- FOR EACH CHAR IN U DO <<PRIN2!* CHAR; J := J+1;
- IF J=5 THEN <<IF !*NAT THEN PRIN2!* '! ;
- J := 0>>>>;
- RETURN J
- END;
- SYMBOLIC PROCEDURE BFLERRMSG U;
- %Standard error message for BFLOAT module;
- REDERR LIST("Invalid argument to",U);
- COMMENT Simp property for !:BF!: since PREP is identity;
- SYMBOLIC PROCEDURE !:BF!:SIMP U; ('!:BF!: . U) ./ 1;
- PUT('!:BF!:,'SIMPFN,'!:BF!:SIMP);
- !:PREC!: := 12; %default value;
- INITDMODE 'BIGFLOAT;
- SYMBOLIC PROCEDURE PRECISION N;
- IF N=0 THEN !:PREC!:-2 ELSE <<!:PREC!: := N+2; N>>;
- SYMBOLIC OPERATOR PRECISION;
- COMMENT *** Tables for Elementary Function Numerical Values ***;
- DEFLIST('((EXP BIGFLOAT) (LOG BIGFLOAT) (SIN BIGFLOAT) (COS BIGFLOAT)
- (TAN BIGFLOAT) (ASIN BIGFLOAT) (ACOS BIGFLOAT)
- (ATAN BIGFLOAT) (SQRT BIGFLOAT)),
- 'TARGETMODE);
- PUT('EXP,'DOMAINFN,'EXP!*);
- SYMBOLIC PROCEDURE EXP!* U; EXP!:(U,!:PREC!:);
- PUT('LOG,'DOMAINFN,'LOG!*);
- SYMBOLIC PROCEDURE LOG!* U; LOG!:(U,!:PREC!:);
- PUT('SIN,'DOMAINFN,'SIN!*);
- SYMBOLIC PROCEDURE SIN!* U; SIN!:(U,!:PREC!:);
- PUT('COS,'DOMAINFN,'COS!*);
- SYMBOLIC PROCEDURE COS!* U; COS!:(U,!:PREC!:);
- PUT('TAN,'DOMAINFN,'TAN!*);
- SYMBOLIC PROCEDURE TAN!* U; TAN!:(U,!:PREC!:);
- PUT('ASIN,'DOMAINFN,'ASIN!*);
- SYMBOLIC PROCEDURE ASIN!* U; ASIN!:(U,!:PREC!:);
- PUT('ACOS,'DOMAINFN,'ACOS!*);
- SYMBOLIC PROCEDURE ACOS!* U; ACOS!:(U,!:PREC!:);
- PUT('ATAN,'DOMAINFN,'ATAN!*);
- SYMBOLIC PROCEDURE ATAN!* U; ATAN!:(U,!:PREC!:);
- PUT('SQRT,'DOMAINFN,'SQRT!*);
- SYMBOLIC PROCEDURE SQRT!* U; SQRT!:(U,!:PREC!:);
- COMMENT *** Tables for constants with numerical values ***;
- DEFLIST('((E BIGFLOAT) (PI BIGFLOAT)),'TARGETMODE);
- PUT('E,'DOMAINFN,'E!*);
- PUT('PI,'DOMAINFN,'PI!*);
- SYMBOLIC PROCEDURE PI!*;
- IF !:PREC!:>1000 THEN !:BIGPI !:PREC!: ELSE !:PI !:PREC!:;
- SYMBOLIC PROCEDURE E!*; !:E !:PREC!:;
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** ARBITRARY PRECISION REAL ARITHMETIC SYSTEM **$
- %** machine-independent version **$
- %** **$
- %** made by **$
- %** **$
- %** Tateaki Sasaki **$
- %** **$
- %** The University of Utah, March 1979 **$
- %** **$
- %**=========================================================**$
- %** **$
- %** For design philosophy and characteristics of this **$
- %** system, see T. Sasaki, "An Arbitrary Precision **$
- %** Real Arithmetic Package in REDUCE," Proceedings **$
- %** of EUROSAM '79, Marseille (France), June 1979. **$
- %** **$
- %** For implementing and using this system, see T. Sasaki, **$
- %** "Manual for Arbitrary Precision Real Arithmetic **$
- %** System in REDUCE," Operating Report of Utah Sym- **$
- %** bolic Computation Group. **$
- %** **$
- %**=========================================================**$
- %** **$
- %** In order to speed up this system, you have only to **$
- %** rewrite four routines (DECPREC!:, INCPREC!:, **$
- %** PRECI!:, and ROUND!:LAST) machine-dependently. **$
- %** **$
- %**=========================================================**$
- %** **$
- %** Table of Contents **$
- %** **$
- %** 1-1. Initialization. **$
- %** 1-2. Constructor, selectors and basic predicate. **$
- %** 1-3. Temporary routines for rational number arithmetic. **$
- %** 1-4. Counters. **$
- %** 1-5. Routines for converting the numeric type. **$
- %** 1-6. Routines for converting a big-float number. **$
- %** 1-7. Routines for reading/printing numbers. **$
- %** 2-1. Arithmetic manipulation routines. **$
- %** 2-2. Arithmetic predicates. **$
- %** 3-1. Elementary constants. **$
- %** 3-2. Routines for saving constants. **$
- %** 4-1. Elementary functions. **$
- %** 5-1. Appendix: routines for defining infix operators. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-1. Initialization. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC$ % Mode ====> SYMBOLIC mode $
- GLOBAL '(!:PREC!:)$ % For the global precision $
- %!:PREC!: := NIL$ % Default value of !:PREC!:$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-2. CONSTRUCTOR, SELECTORS and basic PREDICATE. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC SMACRO PROCEDURE MAKE!:BF(MT,EP); %****************$
- %========================================================$
- % This function constructs an internal representation of $
- % a number "n" composed of the mantissa MT and the $
- % exponent EP with the base 10. The magnitude of $
- % the number thus constructed is hence MT*10**EP. $
- % **** CAUTION! MT and EP are integers. So, EP denotes $
- % **** the order of the last figure in "n", where $
- % **** ORDER(n)=k if 10**k <= ABS(n) < 10**(k+1), $
- % **** with the exception ORDER(0)=0. $
- % The number "n" is said to be of precision "k" if its $
- % mantissa is a k-figure number. $
- % MT and EP are any integers (positive or negative). So,$
- % you can handle any big or small numbers. In this $
- % sense, "BF" denotes a BIG-FLOATING-POINT number. $
- % Hereafter, an internal representation of a number $
- % constructed by MAKE!:BF is referred to as a $
- % BIG-FLOAT representation. $
- %========================================================$
- CONS('!:BF!: , CONS(MT,EP))$
- %*************************************************************$
- SYMBOLIC PROCEDURE BFP!:(X); %******************************$
- %==============================================$
- % This function returns T if X is a BIG-FLOAT $
- % representation, else it returns NIL. $
- % X is any LISP entity. $
- %==============================================$
- IF ATOM(X) THEN NIL ELSE
- IF CAR(X) EQ '!:BF!: THEN T ELSE NIL$
- %*************************************************************$
- SYMBOLIC SMACRO PROCEDURE MT!:(NMBR); %*********************$
- %====================================================$
- % This function selects the mantissa of a number "n".$
- % NMBR is a BIG-FLOAT representation of "n". $
- %====================================================$
- CADR(NMBR)$
- %*************************************************************$
- SYMBOLIC SMACRO PROCEDURE EP!:(NMBR); %*********************$
- %====================================================$
- % This function selects the exponent of a number "n".$
- % NMBR is a BIG-FLOAT representation of "n". $
- %====================================================$
- CDDR(NMBR)$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-3. Temporary routines for rational number arithmetic. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE MAKE!:RATNUM(NM,DN); %*******************$
- %=====================================================$
- % This function constructs an internal representation $
- % of a rational number composed of the numerator $
- % NM and the denominator DN. $
- % NM and DN are any integers (positive or negative). $
- % **** Four routines in this section are temporary. $
- % **** That is, if your system has own routines $
- % **** for rational number arithmetic, you can $
- % **** accommodate our system to yours only by $
- % **** redefining these four routines. $
- %=====================================================$
- IF DN=0 THEN REDERR
- ("ZERO DENOMINATOR IN MAKE!:RATNUM") ELSE
- IF DN>0 THEN CONS('!:RATNUM!: , CONS( NM, DN))
- ELSE CONS('!:RATNUM!: , CONS(-NM,-DN))$
- %*************************************************************$
- SYMBOLIC PROCEDURE RATNUMP!:(X); %**************************$
- %===================================================$
- % This function returns T if X is a rational number $
- % representation, else it returns NIL. $
- % X is any LISP entity. $
- %===================================================$
- IF ATOM(X) THEN NIL ELSE
- IF CAR(X) EQ '!:RATNUM!: THEN T ELSE NIL$
- %*************************************************************$
- SYMBOLIC SMACRO PROCEDURE NUMR!:(RNMBR); %******************$
- %===================================================$
- % This function selects the numerator of a rational $
- % number "n". $
- % RNMBR is a rational number representation of "n". $
- %===================================================$
- CADR(RNMBR)$
- %*************************************************************$
- SYMBOLIC SMACRO PROCEDURE DENM!:(RNMBR); %******************$
- %=====================================================$
- % This function selects the denominator of a rational $
- % number "n". $
- % RNMBR is a rational number representation of "n". $
- %=====================================================$
- CDDR(RNMBR)$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-4. COUNTERS. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC SMACRO PROCEDURE PRECI!:(NMBR); %******************$
- %====================================================$
- % This function counts the precision of a number "n".$
- % NMBR is a BIG-FLOAT representation of "n". $
- %====================================================$
- LENGTH( EXPLODE( ABS( MT!:(NMBR))))$
- %*************************************************************$
- SYMBOLIC PROCEDURE ORDER!:(NMBR); %*************************$
- %================================================$
- % This function counts the order of a number "n".$
- % NMBR is a BIG-FLOAT representation of "n". $
- % **** ORDER(n)=k if 10**k <= ABS(n) < 10**(k+1) $
- % **** when n is not 0, and ORDER(0)=0. $
- %================================================$
- IF MT!:(NMBR)=0 THEN 0
- ELSE PRECI!:(NMBR) + EP!:(NMBR) - 1$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-5. Routines for converting the numeric type. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:A2BF(N); %*************************$
- %======================================================$
- % This function converts a number N or a number-like $
- % entity N to a <BIG-FLOAT>, i.e., a BIG-FLOAT $
- % representation of N. $
- % N is either an integer, a floating-point number, $
- % a string representing a number, a rational $
- % number, or a <BIG-FLOAT>. $
- % **** This function is the most general conversion $
- % **** function to get a BIG-FLOAT representation.$
- % **** In this sense, A means an Arbitrary number.$
- % **** A rational number is converted to a <BIG-FLOAT> $
- % **** of precision !:PREC!: if !:PREC!: is not $
- % **** NIL, else the precision is set 50. $
- %======================================================$
- IF BFP!:(N) THEN N ELSE
- IF FIXP(N) THEN MAKE!:BF(N,0) ELSE
- IF FLOATP(N) THEN READ!:NUM(N) ELSE
- IF STRINGP(N) THEN READ!:NUM(N) ELSE
- IF RATNUMP!:(N) THEN CONV!:R2BF(N ,
- (IF !:PREC!: THEN !:PREC!: ELSE 50) )
- ELSE BFLERRMSG 'CONV!:A2BF$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:F2BF(FNMBR); %*********************$
- %================================================$
- % This function converts a floating-point number $
- % FNMBR to a <BIG-FLOAT>, i.e., a BIG-FLOAT $
- % representation. $
- % FNMBR is a floating-point number. $
- % **** CAUSION!. If you input a number, say, 0.1,$
- % **** some systems do not accept it as 0.1 $
- % **** but may accept it as 0.09999999. $
- % **** In such a case, you had better use $
- % **** CONV!:S2BF than to use CONV!:F2BF. $
- %================================================$
- IF FLOATP(FNMBR) THEN READ!:NUM(FNMBR)
- ELSE BFLERRMSG 'CONV!:F2BF$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:I2BF(INTGR); %*********************$
- %====================================================$
- % This function converts an integer INTGR to a <BIG- $
- % FLOAT>, i.e., a BIG-FLOAT representation. $
- % INTGR is an integer. $
- %====================================================$
- IF FIXP(INTGR) THEN MAKE!:BF(INTGR,0)
- ELSE BFLERRMSG 'CONV!:I2BF$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:R2BF(RNMBR,K); %*******************$
- %=====================================================$
- % This function converts a rational number RNMBR to a $
- % <BIG-FLOAT> of precision K, i.e., a BIG-FLOAT $
- % representation with a given precision. $
- % RNMBR is a rational number representation. $
- % K is a positive integer. $
- %=====================================================$
- IF RATNUMP!:(RNMBR) AND FIXP(K) AND K>0 THEN
- DIVIDE!:( MAKE!:BF( NUMR!:(RNMBR),0) ,
- MAKE!:BF( DENM!:(RNMBR),0) , K)
- ELSE BFLERRMSG 'CONV!:R2BF$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:S2BF(STRNG); %*********************$
- %==============================================$
- % This function converts a string representing $
- % a number "n" to a <BIG-FLOAT>, i.e., $
- % a BIG-FLOAT representation. $
- % STRNG is a string representing "n". "n" may $
- % be an integer, a floating-point number $
- % of any precision, or a rational number. $
- % **** CAUTION! Some systems may set the $
- % **** maximum size of string. $
- %==============================================$
- IF STRINGP(STRNG) THEN READ!:NUM(STRNG)
- ELSE BFLERRMSG 'CONV!:S2BF$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:BF2F(NMBR); %**********************$
- %=========================================================$
- % This function converts a <BIG-FLOAT>, i.e., a BIG-FLOAT $
- % representation of "n", to a floating-point number. $
- % NMBR is a BIG-FLOAT representation of the number "n". $
- %=========================================================$
- IF BFP!:(NMBR) THEN
- TIMES( FLOAT( MT!:(NMBR)) ,
- FLOAT( EXPT(10 , EP!:(NMBR))) )
- ELSE BFLERRMSG 'CONV!:BF2F$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:BF2I(NMBR); %**********************$
- %=========================================================$
- % This function converts a <BIG-FLOAT>, i.e., a BIG-FLOAT $
- % representation of "n", to an integer. The result $
- % is the integer part of "n". $
- % **** For getting the nearest integer to "n", please use $
- % **** the combination MT!:( CONV!:EP(NMBR,0)). $
- % NMBR is a BIG-FLOAT representation of the number "n". $
- %=========================================================$
- IF BFP!:(NMBR) THEN
- IF EP!:(NMBR:=CUT!:EP(NMBR,0)) = 0 THEN MT!:(NMBR)
- ELSE MT!:(NMBR)*EXPT(10 , EP!:(NMBR))
- ELSE BFLERRMSG 'CONV!:BF2I$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:BF2R(NMBR); %**********************$
- %=========================================================$
- % This function converts a <BIG-FLOAT>, i.e., a BIG-FLOAT $
- % representation of "n", to a rational number. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % **** The numerator and the denominator of the result $
- % **** have no common divisor. $
- %=========================================================$
- IF BFP!:(NMBR) THEN
- BEGIN INTEGER NN,ND,M,N,Q;
- IF (Q:=EP!:(NMBR)) >= 0 THEN
- <<NN:=MT!:(NMBR)*EXPT(10,Q); ND:=1; M:=1>>
- ELSE <<NN:=MT!:(NMBR); ND:=EXPT(10,-Q);
- IF ABS(NN) > ABS(ND) THEN <<M:=NN; N:=ND>>
- ELSE <<M:=ND; N:=NN>>;
- WHILE NOT(N=0) DO
- <<Q:=REMAINDER(M,N); M:=N; N:=Q>> >>;
- RETURN MAKE!:RATNUM( NN/M , ND/M);
- END
- ELSE BFLERRMSG 'CONV!:BF2R$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-6. Routines for converting a BIG-FLOAT number. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE DECPREC!:(NMBR,K); %*********************$
- %======================================================$
- % This function converts a number "n" to an equivalent $
- % number the precision of which is decreased by K.$
- % **** CAUTION! No rounding is made. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is a positive integer. $
- %======================================================$
- MAKE!:BF( MT!:(NMBR)/EXPT(10,K) , EP!:(NMBR)+K)$
- %*************************************************************$
- SYMBOLIC PROCEDURE INCPREC!:(NMBR,K); %*********************$
- %======================================================$
- % This function converts a number "n" to an equivalent $
- % number the precision of which is increased by K.$
- % **** CAUTION! No rounding is made. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is a positive integer. $
- %======================================================$
- MAKE!:BF( MT!:(NMBR)*EXPT(10,K) , EP!:(NMBR)-K)$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:MT(NMBR,K); %**********************$
- %===========================================$
- % This function converts a number "n" to an $
- % equivalent number of precision K by $
- % rounding "n" or adding "0"s to "n". $
- % NMBR is a BIG-FLOAT representation of "n".$
- % K is a positive integer. $
- %===========================================$
- IF BFP!:(NMBR) AND FIXP(K) AND K>0 THEN
- IF (K:=PRECI!:(NMBR)-K) = 0 THEN NMBR
- ELSE IF K<0 THEN INCPREC!:(NMBR,-K)
- ELSE ROUND!:LAST( DECPREC!:(NMBR,K-1))
- ELSE BFLERRMSG 'CONV!:MT$
- %*************************************************************$
- SYMBOLIC PROCEDURE CONV!:EP(NMBR,K); %**********************$
- %==============================================$
- % This function converts a number "n" to an $
- % equivalent number having the exponent K $
- % by rounding "n" or adding "0"s to "n". $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is an integer (positive or negative). $
- %==============================================$
- IF BFP!:(NMBR) AND FIXP(K) THEN
- IF (K:=K-EP!:(NMBR)) = 0 THEN NMBR
- ELSE IF K<0 THEN INCPREC!:(NMBR,-K)
- ELSE ROUND!:LAST( DECPREC!:(NMBR,K-1))
- ELSE BFLERRMSG 'CONV!:EP$
- %*************************************************************$
- SYMBOLIC PROCEDURE CUT!:MT(NMBR,K); %***********************$
- %======================================================$
- % This function returns a given number "n" unchanged $
- % if its precision is not greater than K, else it $
- % cuts off its mantissa at the (K+1)th place and $
- % returns an equivalent number of precision K. $
- % **** CAUTION! No rounding is made. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is a positive integer. $
- %======================================================$
- IF BFP!:(NMBR) AND FIXP(K) AND K>0 THEN
- IF (K:=PRECI!:(NMBR)-K) <= 0 THEN NMBR
- ELSE DECPREC!:(NMBR,K)
- ELSE BFLERRMSG 'CUT!:MT$
- %*************************************************************$
- SYMBOLIC PROCEDURE CUT!:EP(NMBR,K); %***********************$
- %======================================================$
- % This function returns a given number "n" unchanged $
- % if its exponent is not less than K, else it $
- % cuts off its mantissa and returns an equivalent $
- % number of exponent K. $
- % **** CAUTION! No rounding is made. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is an integer (positive or negative). $
- %======================================================$
- IF BFP!:(NMBR) AND FIXP(K) THEN
- IF (K:=K-EP!:(NMBR)) <= 0 THEN NMBR
- ELSE DECPREC!:(NMBR,K)
- ELSE BFLERRMSG 'CUT!:EP$
- %*************************************************************$
- SYMBOLIC PROCEDURE MATCH!:(N1,N2); %************************$
- %==========================================================$
- % This function converts either "n1" or "n2" so that they $
- % have the same exponent, which is the smaller of $
- % the exponents of "n1" and "n2". $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- % **** CAUTION! Using this function, one of the previous $
- % **** expressions of "n1" and "n2" is lost. $
- %==========================================================$
- IF BFP!:(N1) AND BFP!:(N2) THEN
- BEGIN INTEGER E1,E2; SCALAR N;
- IF (E1:=EP!:(N1)) = (E2:=EP!:(N2)) THEN RETURN T;
- IF E1>E2 THEN <<RPLACA(N1 , CAR(N:=CONV!:EP(N1,E2)));
- RPLACD(N1 , CDR(N)) >>
- ELSE <<RPLACA(N2 , CAR(N:=CONV!:EP(N2,E1)));
- RPLACD(N2 , CDR(N)) >>; RETURN T;
- END
- ELSE BFLERRMSG 'MATCH!:$
- %*************************************************************$
- SYMBOLIC PROCEDURE ROUND!:MT(NMBR,K); %*********************$
- %========================================================$
- % This function rounds a number "n" at the (K+1)th place $
- % and returns an equivalent number of precision K $
- % if the precision of "n" is greater than K, else $
- % it returns the given number unchanged. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is a positive integer. $
- %========================================================$
- IF BFP!:(NMBR) AND FIXP(K) AND K>0 THEN
- IF (K:=PRECI!:(NMBR)-K-1) < 0 THEN NMBR
- ELSE IF K=0 THEN ROUND!:LAST(NMBR)
- ELSE ROUND!:LAST( DECPREC!:(NMBR,K))
- ELSE BFLERRMSG 'ROUND!:MT$
- %*************************************************************$
- SYMBOLIC PROCEDURE ROUND!:EP(NMBR,K); %*********************$
- %==================================================$
- % This function rounds a number "n" and returns an $
- % equivalent number having the exponent K if $
- % the exponent of "n" is less than K, else $
- % it returns the given number unchanged. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is an integer (positive or negative). $
- %==================================================$
- IF BFP!:(NMBR) AND FIXP(K) THEN
- IF (K:=K-1-EP!:(NMBR)) < 0 THEN NMBR
- ELSE IF K=0 THEN ROUND!:LAST(NMBR)
- ELSE ROUND!:LAST( DECPREC!:(NMBR,K))
- ELSE BFLERRMSG 'ROUND!:EP$
- %*************************************************************$
- SYMBOLIC PROCEDURE ROUND!:LAST(NMBR); %*********************$
- %=====================================================$
- % This function rounds a number "n" at its last place.$
- % NMBR is a BIG-FLOAT representation of "n". $
- %=====================================================$
- BEGIN SCALAR N;
- N := DIVIDE(ABS(MT!:(NMBR)),10);
- IF CDR N<5 THEN N := CAR N ELSE N := CAR N+1;
- IF MT!:(NMBR) < 0 THEN N := -N;
- RETURN MAKE!:BF(N , EP!:(NMBR)+1);
- END$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 1-7. Routines for reading/printing numbers. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE READ!:LNUM(L); %*************************$
- %=======================================================$
- % This function reads a long number "n" represented by $
- % a list in a way described below, and constructs $
- % a BIG-FLOAT representation of "n". $
- % **** Using this function, you can input any long $
- % **** floating-point numbers without difficulty. $
- % L is a list of integers, the first element of which $
- % gives the order of "n" and all the next elements $
- % when concatenated give the mantissa of "n". $
- % **** ORDER(n)=k if 10**k <= ABS(n) < 10**(k+1). $
- % **** Except for the first element, all integers in L $
- % **** should not begin with "0" because some $
- % **** systems suppress leading zeros. $
- %=======================================================$
- IF MEMBER(NIL , MAPCAR(L,'FIXP)) THEN BFLERRMSG
- 'READ!:LNUM ELSE
- BEGIN INTEGER MT,EP,K,SIGN; SCALAR U,V;
- MT:=0;
- EP:=CAR( U:=L)+1;
- IF CADR(L)>0 THEN SIGN:=1 ELSE SIGN:=-1;
- WHILE U:=CDR(U) DO
- <<V:=EXPLODE( ABS( CAR(U))); K:=0;
- WHILE V DO <<K:=K+1; V:=CDR(V) >>;
- MT:=MT*EXPT(10,K)+ABS( CAR(U)); EP:=EP-K>>;
- RETURN MAKE!:BF(SIGN*MT,EP);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE READ!:NUM(N); %**************************$
- %========================================================$
- % This function reads a number or a number-like entity N $
- % and constructs a BIG-FLOAT representation of it. $
- % N is an integer, a floating-point number, or a string $
- % representing a number. $
- % **** If the system does not accept or may incorrectly $
- % **** accept the floating-point numbers, you can $
- % **** input them as strings such as "1.234E-56", $
- % **** "-78.90 D+12" , "+3456 B -78", or "901/234". $
- % **** A rational number in a string form is converted $
- % **** to a <BIG-FLOAT> of precision !:PREC!: if $
- % **** !:PREC!: is not NIL, else the precision of $
- % **** the result is set 50. $
- % **** Some systems set the maximum size of strings. If $
- % **** you want to input long numbers exceeding $
- % **** such a maximum size, please use READ!:LNUM. $
- %========================================================$
- IF FIXP(N) THEN MAKE!:BF(N,0) ELSE
- IF NOT( NUMBERP(N) OR STRINGP(N)) THEN BFLERRMSG
- 'READ!:NUM ELSE
- BEGIN INTEGER J,M,SIGN; SCALAR CH,U,V,L,APPEAR!.,APPEAR!/;
- J:=M:=0;
- SIGN:=1;
- U:=V:=APPEAR!.:=APPEAR!/:=NIL;
- L:=EXPLODE(N);
- LOOP: CH:=CAR(L);
- IF DIGIT(CH) THEN <<U:=CONS(CH,U); J:=J+1>> ELSE
- IF CH EQ '!. THEN <<APPEAR!.:=T ; J:=0 >> ELSE
- IF CH EQ '!/ THEN <<APPEAR!/:=T; V:=U; U:=NIL>> ELSE
- IF CH EQ '!- THEN SIGN:=-1 ELSE
- IF CH EQ 'E OR CH EQ 'D OR CH EQ 'B
- OR CH EQ '!e OR CH EQ '!d OR CH EQ '!b THEN GO TO JUMP;
- ENDL: IF L:=CDR(L) THEN GOTO LOOP ELSE GOTO MAKE;
- JUMP: WHILE L:=CDR(L) DO
- <<IF DIGIT( CH:=CAR(L)) OR CH EQ '!-
- THEN V:=CONS(CH,V) >>;
- L:=REVERSE(V);
- IF CAR(L) EQ '!- THEN M:=-COMPRESS( CDR(L))
- ELSE M:= COMPRESS(L);
- MAKE: U:=REVERSE(U);
- V:=REVERSE(V);
- IF APPEAR!/ THEN RETURN CONV!:R2BF
- ( MAKE!:RATNUM( SIGN*COMPRESS(V) , COMPRESS(U)) ,
- (IF !:PREC!: THEN !:PREC!: ELSE 50) );
- IF APPEAR!. THEN J:=-J ELSE J:=0;
- IF SIGN=1 THEN U:=COMPRESS(U) ELSE U:=-COMPRESS(U);
- RETURN MAKE!:BF(U,J+M);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE PRINT!:BF(NMBR,TYPE); %******************$
- %==========================================================$
- % This function prints a number "n" in the print-type TYPE.$
- % NMBR is a BIG-FLOAT representation of "n". $
- % TYPE is either 'N, 'I, 'E, 'F, 'L, 'R, meaning as: $
- % TYPE='N ... the internal representation is printed. $
- % TYPE='I ... the integer part is printed. $
- % TYPE='E ... <mantissa in form *.***>E<exponent>. $
- % TYPE='F ... <integer part>.<decimal part>. $
- % TYPE='L ... in a list form readable by READ!:LNUM. $
- % TYPE='R ... printed as a rational number. $
- % **** The number is printed by being inserted a blank $
- % **** after each five characters. Therefore, you $
- % **** can not use the printed numbers as input data, $
- % **** except when they are printed in type 'L. $
- %==========================================================$
- IF NOT( TYPE EQ 'N OR TYPE EQ 'I OR TYPE EQ 'E OR
- TYPE EQ 'F OR TYPE EQ 'L OR TYPE EQ 'R)
- OR NOT( BFP!:(NMBR)) THEN BFLERRMSG 'PRINT!:BF ELSE
- BEGIN INTEGER J,K; SCALAR U,V;
- IF ZEROP!:(NMBR) THEN NMBR:=MAKE!:BF(0,0);
- IF TYPE EQ 'I THEN GOTO ITYPE ELSE
- IF TYPE EQ 'E THEN GOTO ETYPE ELSE
- IF TYPE EQ 'F THEN GOTO FTYPE ELSE
- IF TYPE EQ 'L THEN GOTO LTYPE ELSE
- IF TYPE EQ 'R THEN GOTO RTYPE;
- NTYPE: PRINT(NMBR);
- RETURN T;
- ITYPE: U:=EXPLODE( CONV!:BF2I(NMBR));
- J:=0;
- WHILE U DO <<PRIN2( CAR(U)); U:=CDR(U); J:=J+1;
- IF J=5 THEN <<PRIN2(" "); J:=0>> >>;
- TERPRI(); RETURN T;
- ETYPE: U:=EXPLODE( ABS( J:=MT!:(NMBR)));
- IF NULL( CDR(U)) THEN RPLACD(U , LIST(0));
- IF J>=0 THEN U:=CONS( CAR(U) , CONS('!. , CDR(U)))
- ELSE U:=CONS('!- , CONS( CAR(U) , CONS('!.,CDR(U))));
- J:=0;
- WHILE U DO <<PRIN2( CAR(U)); U:=CDR(U); J:=J+1;
- IF J=5 THEN <<PRIN2(" "); J:=0>> >>;
- IF J=0 THEN <<PRIN2("E " ); J:=2>> ELSE
- IF J=1 THEN <<PRIN2(" E " ); J:=4>> ELSE
- IF J=2 THEN <<PRIN2(" E "); J:=0>> ELSE
- IF J=3 THEN <<PRIN2(" E " ); J:=0>> ELSE
- IF J=4 THEN <<PRIN2(" E "); J:=2>>;
- U:=EXPLODE( K:=ORDER!:(NMBR));
- IF K>=0 THEN U:=CONS('!+,U);
- WHILE U DO <<PRIN2( CAR(U)); U:=CDR(U); J:=J+1;
- IF J=5 THEN <<PRIN2(" "); J:=0>> >>;
- TERPRI(); RETURN T;
- FTYPE: U:=EXPLODE( ABS( MT!:(NMBR)));
- IF (J:=EP!:(NMBR)) >= 0 THEN
- <<V:=NIL; WHILE (J:=J-1)>=0 DO V:=CONS(0,V);
- U:=NCONC(U,V) >> ELSE
- IF (J:=ORDER!:(NMBR)+1) > 0 THEN
- <<V:=U; WHILE (J:=J-1)>0 DO V:=CDR(V);
- RPLACD(V , CONS('!.,CDR(V))) >>
- ELSE <<WHILE (J:=J+1)<=0 DO U:=CONS(0,U);
- U:=CONS(0 , CONS('!.,U)) >>;
- IF MT!:(NMBR) < 0 THEN U:=CONS('!-,U);
- J:=0;
- WHILE U DO <<PRIN2( CAR(U)); U:=CDR(U); J:=J+1;
- IF J=5 THEN <<PRIN2(" "); J:=0>> >>;
- TERPRI(); RETURN T;
- LTYPE: PRIN2(" '(");
- PRIN2( ORDER!:(NMBR));
- PRIN2(" ");
- U:=EXPLODE( MT!:(NMBR));
- J:=0;
- WHILE U DO <<PRIN2( CAR(U)); U:=CDR(U); J:=J+1;
- IF J>=5 AND U AND NOT( CAR(U) EQ '!0)
- THEN <<PRIN2(" "); J:=J-5>> >>;
- PRIN2(")"); TERPRI(); RETURN T;
- RTYPE: PRINT!:RATNUM( CONV!:BF2R(NMBR));
- RETURN T;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE PRINT!:RATNUM(RNMBR); %******************$
- %======================================================$
- % This function prints a rational number "n". $
- % RNMBR is a rational number representation of "n". $
- % **** The number is printed by being inserted a blank $
- % **** after each five characters. So, you can $
- % **** not use the printed numbers as input data. $
- %======================================================$
- IF NOT( RATNUMP!:(RNMBR)) THEN BFLERRMSG 'PRINT!:RATNUM ELSE
- BEGIN INTEGER J; SCALAR U,V;
- U:=NUMR!:(RNMBR);
- V:=DENM!:(RNMBR);
- IF V<0 THEN <<U:=-U; V:=-V>>;
- J:=0;
- U:=EXPLODE(U);
- WHILE U DO <<PRIN2( CAR(U)); U:=CDR(U); J:=J+1;
- IF J=5 THEN <<PRIN2(" "); J:=0>> >>;
- IF J=0 THEN <<PRIN2("/ " ); J:=2>> ELSE
- IF J=1 THEN <<PRIN2(" / " ); J:=4>> ELSE
- IF J=2 THEN <<PRIN2(" / "); J:=0>> ELSE
- IF J=3 THEN <<PRIN2(" / " ); J:=0>> ELSE
- IF J=4 THEN <<PRIN2(" / "); J:=2>>;
- V:=EXPLODE(V);
- WHILE V DO <<PRIN2( CAR(V)); V:=CDR(V); J:=J+1;
- IF J=5 THEN <<PRIN2(" "); J:=0>> >>;
- TERPRI(); RETURN T;
- END$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 2-1. Arithmetic manipulation routines. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE ABS!:(NMBR); %***************************$
- %===============================================$
- % This function makes the absolute value of "n".$
- % N is a BIG-FLOAT representation of "n". $
- %===============================================$
- IF MT!:(NMBR) > 0 THEN NMBR
- ELSE MAKE!:BF( -MT!:(NMBR) , EP!:(NMBR))$
- %*************************************************************$
- SYMBOLIC PROCEDURE MINUS!:(NMBR); %*************************$
- %=============================================$
- % This function makes the minus number of "n".$
- % N is a BIG-FLOAT representation of "n". $
- %=============================================$
- MAKE!:BF( -MT!:(NMBR) , EP!:(NMBR))$
- %*************************************************************$
- SYMBOLIC PROCEDURE PLUS!:(N1,N2); %*************************$
- %==========================================================$
- % This function calculates the sum of "n1" and "n2". $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- BEGIN INTEGER E1,E2;
- IF (E1:=EP!:(N1)) = (E2:=EP!:(N2)) THEN RETURN
- MAKE!:BF( MT!:(N1)+MT!:(N2) , E1)
- ELSE IF E1>E2 THEN RETURN MAKE!:BF
- ( MT!:( INCPREC!:(N1,E1-E2))+MT!:(N2) , E2)
- ELSE RETURN MAKE!:BF
- ( MT!:(N1)+MT!:( INCPREC!:(N2,E2-E1)) , E1);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE DIFFERENCE!:(N1,N2); %*******************$
- %==========================================================$
- % This function calculates the difference of "n1" and "n2".$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- BEGIN INTEGER E1,E2;
- IF (E1:=EP!:(N1)) = (E2:=EP!:(N2)) THEN RETURN
- MAKE!:BF( MT!:(N1)-MT!:(N2) , E1)
- ELSE IF E1>E2 THEN RETURN MAKE!:BF
- ( MT!:( INCPREC!:(N1,E1-E2))-MT!:(N2) , E2)
- ELSE RETURN MAKE!:BF
- ( MT!:(N1)-MT!:( INCPREC!:(N2,E2-E1)) , E1);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE TIMES!:(N1,N2); %************************$
- %==========================================================$
- % This function calculates the product of "n1" and "n2". $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- MAKE!:BF( MT!:(N1)*MT!:(N2) , EP!:(N1)+EP!:(N2))$
- %*************************************************************$
- SYMBOLIC PROCEDURE DIVIDE!:(N1,N2,K); %*********************$
- %==========================================================$
- % This function calculates the quotient of "n1" and "n2", $
- % with the precision K, by rounding the ratio of "n1" $
- % and "n2" at the (K+1)th place. $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- % K is any positive integer. $
- %==========================================================$
- BEGIN
- N1:=CONV!:MT(N1 , K+PRECI!:(N2)+1);
- N1:=MAKE!:BF( MT!:(N1)/MT!:(N2) , EP!:(N1)-EP!:(N2));
- RETURN ROUND!:MT(N1,K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE EXPT!:(NMBR,K); %************************$
- %===============================================$
- % This function calculates the Kth power of "n".$
- % The result will become a long number if $
- % ABS(K) >> 1. $
- % NMBR is a BIG-FLOAT representation of "n". $
- % K is an integer (positive or negative). $
- % **** For calculating a power X**K, with non- $
- % **** integer K, please use TEXPT!:ANY. $
- %===============================================$
- IF K>=0 THEN
- MAKE!:BF( EXPT( MT!:(NMBR) , K) , EP!:(NMBR)*K)
- ELSE DIVIDE!:( MAKE!:BF(1,0) , EXPT!:(NMBR,-K) ,
- -PRECI!:(NMBR)*K)$
- %*************************************************************$
- SYMBOLIC PROCEDURE TPLUS!:(N1,N2); %************************$
- %==========================================================$
- % This function calculates the sum of "n1" and "n2" $
- % up to a precision specified by !:PREC!: or N1 or N2.$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2",$
- % otherwise they are converted to <BIG-FLOAT>'s. $
- %==========================================================$
- IF BFP!:( N1:=CONV!:A2BF(N1)) AND
- BFP!:( N2:=CONV!:A2BF(N2)) THEN ROUND!:MT
- ( PLUS!:(N1,N2) , (IF !:PREC!: THEN !:PREC!:
- ELSE MAX( PRECI!:(N1) , PRECI!:(N2))) )
- ELSE BFLERRMSG 'TPLUS!:$
- %*************************************************************$
- SYMBOLIC PROCEDURE TDIFFERENCE!:(N1,N2); %******************$
- %==========================================================$
- % This function calculates the difference of "n1" and "n2" $
- % up to a precision specified by !:PREC!: or N1 or N2.$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2",$
- % otherwise they are converted to <BIG-FLOAT>'s. $
- %==========================================================$
- IF BFP!:( N1:=CONV!:A2BF(N1)) AND
- BFP!:( N2:=CONV!:A2BF(N2)) THEN ROUND!:MT
- ( DIFFERENCE!:(N1,N2) , (IF !:PREC!: THEN !:PREC!:
- ELSE MAX( PRECI!:(N1) , PRECI!:(N2))) )
- ELSE BFLERRMSG 'TDIFFERENCE!:$
- %*************************************************************$
- SYMBOLIC PROCEDURE TTIMES!:(N1,N2); %***********************$
- %==========================================================$
- % This function calculates the product of "n1" and "n2" $
- % up to a precision specified by !:PREC!: or N1 or N2.$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2",$
- % otherwise they are converted to <BIG-FLOAT>'s. $
- %==========================================================$
- IF BFP!:( N1:=CONV!:A2BF(N1)) AND
- BFP!:( N2:=CONV!:A2BF(N2)) THEN ROUND!:MT
- ( TIMES!:(N1,N2) , (IF !:PREC!: THEN !:PREC!:
- ELSE MAX( PRECI!:(N1) , PRECI!:(N2))) )
- ELSE BFLERRMSG 'TTIMES!:$
- %*************************************************************$
- SYMBOLIC PROCEDURE TDIVIDE!:(N1,N2); %**********************$
- %==========================================================$
- % This function calculates the quotient of "n1" and "n2" $
- % up to a precision specified by !:PREC!: or N1 or N2.$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2",$
- % otherwise they are converted to <BIG-FLOAT>'s. $
- %==========================================================$
- IF BFP!:( N1:=CONV!:A2BF(N1)) AND
- BFP!:( N2:=CONV!:A2BF(N2)) THEN
- DIVIDE!:(N1 , N2 , (IF !:PREC!: THEN !:PREC!:
- ELSE MAX( PRECI!:(N1) , PRECI!:(N2))) )
- ELSE BFLERRMSG 'TDIVIDE!:$
- %*************************************************************$
- SYMBOLIC PROCEDURE TEXPT!:(NMBR,K); %***********************$
- %=====================================================$
- % This function calculates the Kth power of "n" up to $
- % the precision specified by !:PREC!: or NMBR. $
- % NMBR is a BIG-FLOAT representation of "n", $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is an integer (positive or negative). $
- % **** For calculating a power X**K, where K is not $
- % **** an integer, please use TEXPT!:ANY. $
- %=====================================================$
- IF BFP!:( NMBR:=CONV!:A2BF(NMBR)) AND FIXP(K) THEN
- IF K=0 THEN MAKE!:BF(1,0) ELSE
- IF K=1 THEN NMBR ELSE
- IF K<0 THEN TDIVIDE!:( MAKE!:BF(1,0) ,
- TEXPT!:(NMBR,-K) )
- ELSE TEXPT!:CAL(NMBR , K , (IF !:PREC!: THEN
- !:PREC!: ELSE PRECI!:(NMBR)) )
- ELSE BFLERRMSG 'TEXPT!:$
- SYMBOLIC PROCEDURE TEXPT!:CAL(NMBR,K,PREC);
- IF K=1 THEN NMBR ELSE
- BEGIN INTEGER K2; SCALAR U;
- U:=ROUND!:MT( TIMES!:(NMBR,NMBR) , PREC);
- IF K=2 THEN RETURN U ELSE
- IF (K-2*(K2:=K/2)) = 0 THEN RETURN
- TEXPT!:CAL(U,K2,PREC)
- ELSE RETURN ROUND!:MT
- ( TIMES!:(NMBR , TEXPT!:CAL(U,K2,PREC)) , PREC);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE QUOTIENT!:(N1,N2); %*********************$
- %==========================================================$
- % This function calculates the integer quotient of "n1" $
- % and "n2", just as the "QUOTIENT" for integers does. $
- % **** For calculating the quotient up to a necessary $
- % **** precision, please use DIVIDE!:. $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- BEGIN INTEGER E1,E2;
- IF (E1:=EP!:(N1)) = (E2:=EP!:(N2)) THEN RETURN
- MAKE!:BF( MT!:(N1)/MT!:(N2) , 0)
- ELSE IF E1>E2 THEN RETURN
- QUOTIENT!:( INCPREC!:(N1,E1-E2) , N2)
- ELSE RETURN
- QUOTIENT!:( N1 , INCPREC!:(N2,E2-E1));
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE REMAINDER!:(N1,N2); %********************$
- %==========================================================$
- % This function calculates the remainder of "n1" and "n2", $
- % just as the "REMAINDER" for integers does. $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- BEGIN INTEGER E1,E2;
- IF (E1:=EP!:(N1)) = (E2:=EP!:(N2)) THEN RETURN
- MAKE!:BF( REMAINDER( MT!:(N1) , MT!:(N2)) , E2)
- ELSE IF E1>E2 THEN RETURN
- REMAINDER!:( INCPREC!:(N1,E1-E2) , N2)
- ELSE RETURN
- REMAINDER!:( N1 , INCPREC!:(N2,E2-E1));
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE TEXPT!:ANY(X,Y); %***********************$
- %====================================================$
- % This function calculates the power x**y, where "x" $
- % and "y" are any numbers. The precision of $
- % the result is specified by !:PREC!: or X or Y.$
- % **** For a negative "x", this function returns $
- % **** -(-x)**y unless "y" is an integer. $
- % X is a BIG-FLOAT representation of "x", otherwise $
- % it is converted to a <BIG-FLOAT>. $
- % Y is either an integer, a floating-point number, $
- % or a BIG-FLOAT number, i.e., a BIG-FLOAT $
- % representation of "y". $
- %====================================================$
- IF FIXP(Y) THEN TEXPT!:(X,Y) ELSE
- IF INTEGERP!:(Y) THEN TEXPT!:(X , CONV!:BF2I(Y)) ELSE
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( BFP!:( Y:=CONV!:A2BF(Y))) THEN BFLERRMSG
- 'TEXPT!:ANY ELSE
- IF MINUSP!:(Y) THEN TDIVIDE!:( MAKE!:BF(1,0) ,
- TEXPT!:ANY(X , MINUS!:(Y)) ) ELSE
- BEGIN INTEGER N; SCALAR XP,YP;
- N:=(IF !:PREC!: THEN !:PREC!:
- ELSE MAX( PRECI!:(X) , PRECI!:(Y)) );
- IF MINUSP!:(X) THEN XP:=MINUS!:(X) ELSE XP:=X;
- IF INTEGERP!:( TIMES!:(Y , CONV!:I2BF(2))) THEN
- <<XP:=INCPREC!:(XP,1);
- YP:=TEXPT!:(XP , CONV!:BF2I(Y));
- YP:=TIMES!:(YP , SQRT!:(XP,N+1));
- YP:=ROUND!:MT(YP,N) >>
- ELSE
- <<YP:=TTIMES!:(Y , LOG!:(XP,N+1));
- YP:=EXP!:(YP,N) >>;
- RETURN (IF MINUSP!:(X) THEN MINUS!:(YP) ELSE YP);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE MAX!:(N1,N2); %**************************$
- %==========================================================$
- % This function returns the larger of "n1" and "n2". $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- IF GREATERP!:(N2,N1) THEN N2 ELSE N1$
- %*************************************************************$
- SYMBOLIC PROCEDURE MIN!:(N1,N2); %**************************$
- %==========================================================$
- % This function returns the smaller of "n1" and "n2". $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- IF LESSP!:(N2,N1) THEN N2 ELSE N1$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 2-2. Arithmetic predicates. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE GREATERP!:(N1,N2); %*********************$
- %==========================================================$
- % This function returns T if "n1" > "n2" else returns NIL. $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- BEGIN INTEGER E1,E2;
- IF (E1:=EP!:(N1)) = (E2:=EP!:(N2)) THEN
- RETURN (IF MT!:(N1) > MT!:(N2) THEN T ELSE NIL)
- ELSE IF E1>E2 THEN
- IF MT!:( INCPREC!:(N1,E1-E2)) > MT!:(N2)
- THEN RETURN T ELSE RETURN NIL
- ELSE IF MT!:(N1) > MT!:( INCPREC!:(N2,E2-E1))
- THEN RETURN T ELSE RETURN NIL;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE GEQ!:(N1,N2); %**************************$
- %==========================================================$
- % This function returns T if "n1" >= "n2" else returns NIL.$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- NOT( LESSP!:(N1,N2))$
- %*************************************************************$
- SYMBOLIC PROCEDURE EQUAL!:(N1,N2); %************************$
- %==========================================================$
- % This function returns T if "n1" = "n2" else returns NIL. $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- IF ZEROP!:( DIFFERENCE!:(N1,N2)) THEN T ELSE NIL$
- %*************************************************************$
- SYMBOLIC PROCEDURE LESSP!:(N1,N2); %************************$
- %==========================================================$
- % This function returns T if "n1" < "n2" else returns NIL. $
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- GREATERP!:(N2,N1)$
- %*************************************************************$
- SYMBOLIC PROCEDURE LEQ!:(N1,N2); %**************************$
- %==========================================================$
- % This function returns T if "n1" <= "n2" else returns NIL.$
- % N1 and N2 are BIG-FLOAT representations of "n1" and "n2".$
- %==========================================================$
- NOT( GREATERP!:(N1,N2))$
- %*************************************************************$
- SYMBOLIC PROCEDURE INTEGERP!:(X); %*************************$
- %===================================================$
- % This function returns T if X is a BIG-FLOAT $
- % representing an integer, else it returns NIL.$
- % X is any LISP entity. $
- %===================================================$
- IF BFP!:(X) THEN IF EP!:(X)>=0 OR
- EQUAL!:(X , CONV!:I2BF( CONV!:BF2I(X))) THEN T
- ELSE NIL
- ELSE NIL$
- %*************************************************************$
- SYMBOLIC PROCEDURE MINUSP!:(X); %***************************$
- %===================================================$
- % This function returns T if "x"<0 else returns NIL.$
- % X is any LISP entity. $
- %===================================================$
- IF BFP!:(X) AND MT!:(X) < 0 THEN T ELSE NIL$
- %*************************************************************$
- SYMBOLIC PROCEDURE ZEROP!:(X); %****************************$
- %===================================================$
- % This function returns T if "x"=0 else returns NIL.$
- % X is any LISP entity. $
- %===================================================$
- IF BFP!:(X) AND MT!:(X) = 0 THEN T ELSE NIL$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 3-1. Elementary CONSTANTS. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:PI(K); %*******************************$
- %====================================================$
- % This function calculates the value of the circular $
- % constant "PI", with the precision K, by $
- % using Machin's well known identity: $
- % PI = 16*atan(1/5) - 4*atan(1/239). $
- % Calculation is performed mainly on integers. $
- % K is a positive integer. $
- %====================================================$
- IF NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG '!:PI ELSE
- IF K<=20 THEN ROUND!:MT
- ( MAKE!:BF( 314159265358979323846 , -20) , K) ELSE
- BEGIN INTEGER K3,S,SS,M,N,X; SCALAR U;
- U:=GET!:CONST( '!:PI , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- SS:=N:=EXPT(10 , K3:=K+3)/5;
- X :=-5**2;
- M:=1;
- WHILE NOT(N=0) DO <<N:=N/X; SS:=SS+N/( M:=M+2) >>;
- S:=N:=EXPT(10,K3)/239;
- X:=-239**2;
- M:=1;
- WHILE NOT(N=0) DO <<N:=N/X; S:=S+N/( M:=M+2) >>;
- ANS: U:=ROUND!:MT( MAKE!:BF( 16*SS-4*S , -K3) , K);
- SAVE!:CONST( '!:PI , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:BIGPI(K); %****************************$
- %====================================================$
- % This function calculates the value of the circular $
- % constant "PI", with the precision K, by the $
- % arithmetic-geometric mean method. (See, $
- % R. Brent, JACM Vol.23, #2, pp.242-251(1976).) $
- % K is a positive integer. $
- % **** This function should be used only when you $
- % **** need "PI" of precision higher than 1000. $
- %====================================================$
- IF NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG '!:BIGPI ELSE
- BEGIN INTEGER K2,N; SCALAR DCUT,HALF,X,Y,U,V;
- U:=GET!:CONST( '!:PI , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- K2 :=K+2;
- HALF:=CONV!:S2BF("0.5");
- DCUT:=MAKE!:BF(10,-K2);
- X:=CONV!:I2BF( N:=1);
- Y:=DIVIDE!:(X , !:SQRT2(K2) , K2);
- U:=CONV!:S2BF("0.25");
-
- WHILE GREATERP!:( ABS!:(DIFFERENCE!:(X,Y)) , DCUT) DO
- <<V:=X;
- X:=TIMES!:( PLUS!:(X,Y) , HALF);
- Y:=SQRT!:( CUT!:EP( TIMES!:(Y,V) , -K2) , K2);
- V:=DIFFERENCE!:(X,V);
- V:=TIMES!:( TIMES!:(V,V) , CONV!:I2BF(N));
- U:=DIFFERENCE!:(U , CUT!:EP(V,-K2));
- N:=2*N>>;
- V:=CUT!:MT( EXPT!:( PLUS!:(X,Y) , 2) , K2);
- U:=DIVIDE!:(V , TIMES!:( CONV!:I2BF(4) , U) , K);
- SAVE!:CONST( '!:PI , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:E(K); %********************************$
- %=====================================================$
- % This function calculates the value of "e", the base $
- % of the natural logarithm, with the precision K,$
- % by summing the Taylor series for exp(x=1). $
- % Calculation is performed mainly on integers. $
- % K is a positive integer. $
- %=====================================================$
- IF NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG '!:E ELSE
- IF K<=20 THEN ROUND!:MT
- ( MAKE!:BF( 271828182845904523536 , -20) , K) ELSE
- BEGIN INTEGER K2,ANS,M,N; SCALAR U;
- U:=GET!:CONST( '!:E , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- K2:=K+2;
- M :=1;
- N :=EXPT(10,K2);
- ANS:=0;
- WHILE NOT(N=0) DO ANS:=ANS+( N:=N/( M:=M+1));
- ANS:=ANS+2*EXPT(10,K2);
- U:=ROUND!:MT( MAKE!:BF(ANS,-K2) , K);
- SAVE!:CONST( '!:E , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:E01(K); %******************************$
- %=====================================================$
- % This function calculates exp(0.1), the value of the $
- % exponential function at the point 0.1, with $
- % the precision K. $
- % K is a positive integer. $
- %=====================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:E01 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=EXP!:( CONV!:S2BF("0.1") , K);
- SAVE!:CONST( '!:E01 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:LOG2(K); %*****************************$
- %==============================================$
- % This function calculates log(2), the natural $
- % logarithm of 2, with the precision K. $
- % K is a positive integer. $
- %==============================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:LOG2 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=LOG!:( CONV!:I2BF(2) , K);
- SAVE!:CONST( '!:LOG2 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:LOG3(K); %*****************************$
- %==============================================$
- % This function calculates log(3), the natural $
- % logarithm of 3, with the precision K. $
- % K is a positive integer. $
- %==============================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:LOG3 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=LOG!:( CONV!:I2BF(3) , K);
- SAVE!:CONST( '!:LOG3 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:LOG5(K); %*****************************$
- %==============================================$
- % This function calculates log(5), the natural $
- % logarithm of 5, with the precision K. $
- % K is a positive integer. $
- %==============================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:LOG5 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=LOG!:( CONV!:I2BF(5) , K);
- SAVE!:CONST( '!:LOG5 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:LOG10(K); %****************************$
- %===============================================$
- % This function calculates log(10), the natural $
- % logarithm of 10, with the precision K. $
- % K is a positive integer. $
- %===============================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:LOG10 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=LOG!:( CONV!:I2BF(10) , K);
- SAVE!:CONST( '!:LOG10 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:LOGPI(K); %****************************$
- %===============================================$
- % This function calculates log(PI), the natural $
- % logarithm of "PI", with the precision K. $
- % K is a positive integer. $
- %===============================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:LOGPI , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=LOG!:( !:PI(K+2) , K);
- SAVE!:CONST( '!:LOGPI , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:SQRT2(K); %****************************$
- %===================================================$
- % This function calculates SQRT(2), the square root $
- % of 2, with the precision K. $
- % K is a positive integer. $
- %===================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:SQRT2 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=SQRT!:( CONV!:I2BF(2) , K);
- SAVE!:CONST( '!:SQRT2 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:SQRT3(K); %****************************$
- %===================================================$
- % This function calculates SQRT(3), the square root $
- % of 3, with the precision K. $
- % K is a positive integer. $
- %===================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:SQRT3 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=SQRT!:( CONV!:I2BF(3) , K);
- SAVE!:CONST( '!:SQRT3 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:SQRT5(K); %****************************$
- %===================================================$
- % This function calculates SQRT(5), the square root $
- % of 5, with the precision K. $
- % K is a positive integer. $
- %===================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:SQRT5 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=SQRT!:( CONV!:I2BF(5) , K);
- SAVE!:CONST( '!:SQRT5 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:SQRT10(K); %***************************$
- %====================================================$
- % This function calculates SQRT(10), the square root $
- % of 10, with the precision K. $
- % K is a positive integer. $
- %====================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:SQRT10 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=SQRT!:( CONV!:I2BF(10) , K);
- SAVE!:CONST( '!:SQRT10 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:SQRTPI(K); %***************************$
- %====================================================$
- % This function calculates SQRT(PI), the square root $
- % of "PI", with the precision K. $
- % K is a positive integer. $
- %====================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:SQRTPI , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=SQRT!:( !:PI(K+2) , K);
- SAVE!:CONST( '!:SQRTPI , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:SQRTE(K); %****************************$
- %===================================================$
- % This function calculates SQRT(e), the square root $
- % of "e", with the precision K. $
- % K is a positive integer. $
- %===================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:SQRTE , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=SQRT!:( !:E(K+2) , K);
- SAVE!:CONST( '!:SQRTE , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:CBRT2(K); %****************************$
- %=================================================$
- % This function calculates CBRT(2), the cube root $
- % of 2, with the precision K. $
- % K is a positive integer. $
- %=================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:CBRT2 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=CBRT!:( CONV!:I2BF(2) , K);
- SAVE!:CONST( '!:CBRT2 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:CBRT3(K); %****************************$
- %=================================================$
- % This function calculates CBRT(3), the cube root $
- % of 3, with the precision K. $
- % K is a positive integer. $
- %=================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:CBRT3 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=CBRT!:( CONV!:I2BF(3) , K);
- SAVE!:CONST( '!:CBRT3 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:CBRT5(K); %****************************$
- %=================================================$
- % This function calculates CBRT(5), the cube root $
- % of 5, with the precision K. $
- % K is a positive integer. $
- %=================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:CBRT5 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=CBRT!:( CONV!:I2BF(5) , K);
- SAVE!:CONST( '!:CBRT5 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:CBRT10(K); %***************************$
- %==================================================$
- % This function calculates CBRT(10), the cube root $
- % of 10, with the precision K. $
- % K is a positive integer. $
- %==================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:CBRT10 , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=CBRT!:( CONV!:I2BF(10) , K);
- SAVE!:CONST( '!:CBRT10 , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:CBRTPI(K); %***************************$
- %==================================================$
- % This function calculates CBRT(PI), the cube root $
- % of "PI", with the precision K. $
- % K is a positive integer. $
- %==================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:CBRTPI , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=CBRT!:( !:PI(K+2) , K);
- SAVE!:CONST( '!:CBRTPI , U); RETURN U;
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE !:CBRTE(K); %****************************$
- %=================================================$
- % This function calculates CBRT(e), the cube root $
- % of "e", with the precision K. $
- % K is a positive integer. $
- %=================================================$
- BEGIN SCALAR U;
- U:=GET!:CONST( '!:CBRTE , K);
- IF U = "NOT FOUND" THEN NIL ELSE RETURN U;
- U:=CBRT!:( !:E(K+2) , K);
- SAVE!:CONST( '!:CBRTE , U); RETURN U;
- END$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 3-2. Routines for saving CONSTANTS. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE GET!:CONST(CNST,K); %********************$
- %==================================================$
- % This function returns the value of constant CNST $
- % of the precision K, if it was calculated $
- % previously with, at least, the precision K, $
- % else it returns "NOT FOUND". $
- % CNST is the name of the constant (to be quoted). $
- % K is a positive integer. $
- %==================================================$
- IF ATOM(CNST) AND FIXP(K) AND K>0 THEN
- BEGIN SCALAR U;
- U:=GET(CNST , 'SAVE!:C);
- IF NULL(U) OR CAR(U)<K THEN RETURN "NOT FOUND"
- ELSE IF CAR(U)=K THEN RETURN CDR(U)
- ELSE RETURN ROUND!:MT(CDR(U),K);
- END
- ELSE BFLERRMSG 'GET!:CONST$
- %*************************************************************$
- SYMBOLIC PROCEDURE SAVE!:CONST(CNST,NMBR); %****************$
- %=================================================$
- % This function saves the value of constant CNST $
- % for the later use. $
- % CNST is the name of the constant (to be quoted).$
- % NMBR is a BIG-FLOAT representation of the value.$
- %=================================================$
- IF ATOM(CNST) AND BFP!:(NMBR) THEN
- PUT(CNST , 'SAVE!:C , CONS( PRECI!:(NMBR) , NMBR))
- ELSE BFLERRMSG 'SAVE!:CONST$
- %*************************************************************$
- SYMBOLIC PROCEDURE SET!:CONST(CNST,L); %********************$
- %=================================================$
- % This function sets the value of constant CNST. $
- % CNST is the name of the constant (to be quoted).$
- % L is a list of integers, which represents the $
- % value of the constant in the way described $
- % in the function READ!:LNUM. $
- %=================================================$
- SAVE!:CONST(CNST , READ!:LNUM(L))$
- %*************************************************************$
- SYMBOLIC$ %SETTING THE CONSTANTS ***************************$
- SET!:CONST( '!:PI , '( 0 3141 59265 35897 93238 46264
- 33832 79502 88419 71693 99375 105820 9749 44592 30781
- 64062 86208 99862 80348 25342 11706 79821 48086 51328
- 23066 47093 84460 95505 82231 72535 94081 28481 1174
- 5028410 2701 93852 11055 59644 62294 89549 30381 96442
- 88109 8) )$
- SET!:CONST( '!:E , '( 0 2718 28182 84590 45235 36028
- 74713 52662 49775 72470 93699 95957 49669 67627 72407
- 66303 53547 59457 13821 78525 16642 74274 66391 93200
- 30599 21817 41359 66290 43572 90033 42952 60595 63073
- 81323 28627 943490 7632 33829 88075 31952 510190 1157
- 38341 9) )$
- SET!:CONST( '!:E01 , '( 0 1105 17091 80756 47624 81170
- 78264 90246 66822 45471 94737 51871 87928 63289 44096
- 79667 47654 30298 91433 18970 74865 36329 2) )$
- SET!:CONST( '!:LOG2 , '(-1 6931 47180 55994 53094 17232
- 12145 81765 68075 50013 43602 55254 1206 800094 93393
- 62196 96947 15605 86332 69964 18687 54200 2) )$
- SET!:CONST( '!:LOG3 , '( 0 1098 61228 866810 9691 39524
- 52369 22525 70464 74905 57822 74945 17346 94333 63749
- 42932 18608 96687 36157 54813 73208 87879 7) )$
- SET!:CONST( '!:LOG5 , '( 0 1609 43791 2434100 374 60075
- 93332 26187 63952 56013 54268 51772 19126 47891 47417
- 898770 7657 764630 1338 78093 179610 7999 7) )$
- SET!:CONST( '!:LOG10 , '( 0 2302 58509 29940 456840 1799
- 14546 84364 20760 11014 88628 77297 60333 27900 96757
- 26096 77352 48023 599720 5089 59829 83419 7) )$
- SET!:CONST( '!:LOGPI , '( 0 1144 72988 5849400 174 14342
- 73513 53058 71164 72948 12915 31157 15136 23071 47213
- 77698 848260 7978 36232 70275 48970 77020 1) )$
- SET!:CONST( '!:SQRT2 , '( 0 1414 21356 23730 95048 80168
- 872420 96980 7856 96718 75376 94807 31766 79737 99073
- 24784 621070 38850 3875 34327 64157 27350 1) )$
- SET!:CONST( '!:SQRT3 , '( 0 17320 5080 75688 77293 52744
- 634150 5872 36694 28052 53810 38062 805580 6979 45193
- 301690 88000 3708 11461 86757 24857 56756 3) )$
- SET!:CONST( '!:SQRT5 , '( 0 22360 6797 74997 89696 40917
- 36687 31276 235440 6183 59611 52572 42708 97245 4105
- 209256 37804 89941 441440 8378 78227 49695 1) )$
- SET!:CONST( '!:SQRT10, '( 0 3162 277660 1683 79331 99889
- 35444 32718 53371 95551 39325 21682 685750 4852 79259
- 44386 39238 22134 424810 8379 30029 51873 47))$
- SET!:CONST( '!:SQRTPI, '( 0 1772 453850 9055 16027 29816
- 74833 41145 18279 75494 56122 38712 821380 7789 85291
- 12845 91032 18137 49506 56738 54466 54162 3) )$
- SET!:CONST( '!:SQRTE , '( 0 1648 721270 7001 28146 8486
- 507878 14163 57165 3776100 710 14801 15750 79311 64066
- 10211 94215 60863 27765 20056 36664 30028 7) )$
- SET!:CONST( '!:CBRT2 , '( 0 1259 92104 98948 73164 7672
- 106072 78228 350570 2514 64701 5079800 819 75112 15529
- 96765 13959 48372 93965 62436 25509 41543 1) )$
- SET!:CONST( '!:CBRT3 , '( 0 1442 249570 30740 8382 32163
- 83107 80109 58839 18692 53499 35057 75464 16194 54168
- 75968 29997 33985 47554 79705 64525 66868 4) )$
- SET!:CONST( '!:CBRT5 , '( 0 1709 97594 66766 96989 35310
- 88725 43860 10986 80551 105430 5492 43828 61707 44429
- 592050 4173 21625 71870 10020 18900 220450 ) )$
- SET!:CONST( '!:CBRT10, '( 0 2154 4346900 318 83721 75929
- 35665 19350 49525 93449 42192 10858 24892 35506 34641
- 11066 48340 80018 544150 3543 24327 61012 6) )$
- SET!:CONST( '!:CBRTPI, '( 0 1464 59188 75615 232630 2014
- 25272 63790 39173 85968 55627 93717 43572 55937 13839
- 36497 98286 26614 56820 67820 353820 89750 ) )$
- SET!:CONST( '!:CBRTE , '( 0 1395 61242 50860 89528 62812
- 531960 2586 83759 79065 15199 40698 26175 167060 3173
- 90156 45951 84696 97888 17295 83022 41352 1) )$
- %*************************************************************$
- %*************************************************************$
- %** **$
- %** 4-1. Elementary FUNCTIONS. **$
- %** **$
- %*************************************************************$
- %*************************************************************$
- %*************************************************************$
- SYMBOLIC PROCEDURE SQRT!:(X,K); %***************************$
- %===================================================$
- % This function calculates SQRT(x), the square root $
- % of "x", with the precision K, by Newton's $
- % iteration method. $
- % X is a BIG-FLOAT representation of "x", x >= 0, $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %===================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR MINUSP!:(X) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'SQRT!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(0) ELSE
- BEGIN INTEGER K2,NCUT,NFIG; SCALAR DCUT,HALF,DY,Y,Y0,U;
- K2 :=K+2;
- NCUT:=K2-(ORDER!:(X)+1)/2;
- HALF:=CONV!:S2BF("0.5");
- DCUT:=MAKE!:BF(10,-NCUT);
- DY :=MAKE!:BF(20,-NCUT);
- Y0:=CONV!:MT(X,2);
- IF REMAINDER( EP!:(Y0) , 2) = 0 THEN
- Y0:=MAKE!:BF( 3+2*MT!:(Y0)/25 , EP!:(Y0)/2)
- ELSE Y0:=MAKE!:BF( 10+2*MT!:(Y0)/9 , (EP!:(Y0)-1)/2);
- NFIG:=1;
- WHILE NFIG<K2 OR GREATERP!:( ABS!:(DY) , DCUT) DO
- <<IF (NFIG:=2*NFIG) > K2 THEN NFIG:=K2;
- U :=DIVIDE!:(X,Y0,NFIG);
- Y :=TIMES!:( PLUS!:(Y0,U) , HALF);
- DY:=DIFFERENCE!:(Y,Y0);
- Y0:=Y>>;
- RETURN ROUND!:MT(Y,K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE CBRT!:(X,K); %***************************$
- %===================================================$
- % This function calculates CBRT(x), the cube root $
- % of "x", with the precision K, by Newton's $
- % iteration method. $
- % X is a BIG-FLOAT representation of any real "x", $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %===================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'CBRT!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(0) ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( CBRT!:( MINUS!:(X) , K)) ELSE
- BEGIN INTEGER K2,NCUT,NFIG,J; SCALAR DCUT,THRE,DY,Y,U;
- K2 :=K+2;
- NCUT:=K2-(ORDER!:(X)+2)/3;
- THRE:=CONV!:I2BF(3);
- DCUT:=MAKE!:BF(10,-NCUT);
- DY :=MAKE!:BF(20,-NCUT);
- Y:=CONV!:MT(X,3);
- IF (J:=REMAINDER( EP!:(Y) , 3)) = 0 THEN
- Y:=MAKE!:BF( 5 + MT!:(Y)/167 , EP!:(Y)/3) ELSE
- IF J=1 OR J=-2 THEN
- Y:=MAKE!:BF( 10+ MT!:(Y)/75 , (EP!:(Y)-1)/3)
- ELSE Y:=MAKE!:BF( 22+2*MT!:(Y)/75 , (EP!:(Y)-2)/3);
- NFIG:=1;
- WHILE NFIG<K2 OR GREATERP!:( ABS!:(DY) , DCUT) DO
- <<IF (NFIG:=2*NFIG) > K2 THEN NFIG:=K2;
- U :=CUT!:MT( TIMES!:(Y,Y) , NFIG);
- U :=DIVIDE!:(X , U , NFIG);
- J :=ORDER!:( U:=DIFFERENCE!:(U,Y))+NCUT-K2;
- DY:=DIVIDE!:(U , THRE , MAX(1,NFIG+J));
- Y :=PLUS!:(Y,DY) >>;
- RETURN ROUND!:MT(Y,K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE EXP!:(X,K); %****************************$
- %=================================================$
- % This function calculates exp(x), the value of $
- % the exponential function at the point "x", $
- % with the precision K, by summing terms of $
- % the Taylor series for exp(z), 0 < z < 1. $
- % X is a BIG-FLOAT representation of any real "x",$
- % otherwise it is converted to a <BIG-FLOAT>.$
- % K is a positive integer. $
- %=================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'EXP!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(1) ELSE
- BEGIN INTEGER K2,M; SCALAR ONE,Q,R,Y,YQ,YR,SAVE!:P;
- K2 :=K+2;
- ONE:=CONV!:I2BF(1);
- Q:=CONV!:I2BF( M:=CONV!:BF2I( Y:=ABS!:(X)));
- R:=DIFFERENCE!:(Y,Q);
- IF ZEROP!:(Q) THEN YQ:=ONE
- ELSE << SAVE!:P:=!:PREC!:; !:PREC!::=K2;
- YQ:=TEXPT!:( !:E(K2) , M);
- !:PREC!::=SAVE!:P>>;
- IF ZEROP!:(R) THEN YR:=ONE ELSE
- BEGIN INTEGER J,N; SCALAR DCUT,FCTRIAL,RI,TM;
-
- DCUT:=MAKE!:BF(10,-K2);
- YR:=RI:=TM:=ONE;
-
- M:=1;
- J:=0;
- WHILE GREATERP!:(TM,DCUT) DO
- <<FCTRIAL:=CONV!:I2BF( M:=M*( J:=J+1));
- RI:=CUT!:EP( TIMES!:(RI,R) , -K2);
- N :=MAX(1 , K2-ORDER!:(FCTRIAL)+ORDER!:(RI));
- TM:=DIVIDE!:(RI,FCTRIAL,N);
- YR:=PLUS!:(YR,TM); IF REMAINDER(J,10)=0 THEN
- YR:=CUT!:EP(YR,-K2) >>;
- END;
- Y:=CUT!:MT( TIMES!:(YQ,YR) , K+1);
- RETURN (IF MINUSP!:(X) THEN DIVIDE!:(ONE,Y,K)
- ELSE ROUND!:LAST(Y) );
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE LOG!:(X,K); %****************************$
- %===================================================$
- % This function calculates log(x), the value of the $
- % logarithmic function at the point "x", with $
- % the precision K, by summing terms of the $
- % Taylor series for log(1+z), 0 < z < 0.10518. $
- % X is a BIG-FLOAT representation of "x", x > 0, $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %===================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- MINUSP!:(X) OR ZEROP!:(X) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'LOG!: ELSE
- IF EQUAL!:(X , CONV!:I2BF(1)) THEN CONV!:I2BF(0) ELSE
- BEGIN INTEGER K2,M; SCALAR EE,ES,ONE,SIGN,L,Y,Z,SAVE!:P;
- K2 :=K+2;
- ONE:=CONV!:I2BF(1);
- EE :=!:E(K2);
- ES :=!:E01(K2);
- IF GREATERP!:(X,ONE) THEN <<SIGN:=ONE; Y:=X>>
- ELSE <<SIGN:=MINUS!:(ONE); Y:=DIVIDE!:(ONE,X,K2) >>;
- IF LESSP!:(Y,EE) THEN <<M:=0; Z:=Y>>
- ELSE <<IF (M:=(ORDER!:(Y)*23)/10) = 0 THEN Z:=Y
- ELSE << SAVE!:P:=!:PREC!:; !:PREC!::=K2;
- Z:=DIVIDE!:(Y , TEXPT!:(EE,M) , K2);
- !:PREC!::=SAVE!:P>>;
- WHILE GREATERP!:(Z,EE) DO
- <<M:=M+1; Z:=DIVIDE!:(Z,EE,K2) >> >>;
- L:=CONV!:I2BF(M);
- Y:=CONV!:S2BF("0.1");
- WHILE GREATERP!:(Z,ES) DO
- <<L:=PLUS!:(L,Y); Z:=DIVIDE!:(Z,ES,K2) >>;
- Z:=DIFFERENCE!:(Z,ONE);
- BEGIN INTEGER N; SCALAR DCUT,TM,ZI;
- Y:=TM:=ZI:=Z;
- Z:=MINUS!:(Z);
- DCUT:=MAKE!:BF(10,-K2);
- M:=1;
- WHILE GREATERP!:( ABS!:(TM) , DCUT) DO
- <<ZI:=CUT!:EP( TIMES!:(ZI,Z) , -K2);
- N :=MAX(1 , K2+ORDER!:(ZI));
- TM:=DIVIDE!:(ZI , CONV!:I2BF( M:=M+1) , N);
- Y :=PLUS!:(Y,TM); IF REMAINDER(M,10)=0 THEN
- Y:=CUT!:EP(Y,-K2) >>;
- END;
- Y:=PLUS!:(Y,L);
- RETURN ROUND!:MT( TIMES!:(SIGN,Y) , K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE LN!:(X,K); %*****************************$
- %=================================================$
- % This function calculates log(x), the value of $
- % the logarithmic function at the point "x", $
- % with the precision K, by solving $
- % x = exp(y) by Newton's method. $
- % X is a BIG-FLOAT representation of "x", x > 0, $
- % otherwise it is converted to a <BIG-FLOAT>.$
- % K is a positive integer. $
- %=================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- MINUSP!:(X) OR ZEROP!:(X) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'LN!: ELSE
- IF EQUAL!:(X , CONV!:I2BF(1)) THEN CONV!:I2BF(0) ELSE
- BEGIN INTEGER K2,M; SCALAR EE,ONE,SIGN,Y,Z,SAVE!:P;
- K2 :=K+2;
- ONE:=CONV!:I2BF(1);
- EE :=!:E(K2+2);
- IF GREATERP!:(X,ONE) THEN <<SIGN:=ONE; Y:=X>>
- ELSE <<SIGN:=MINUS!:(ONE); Y:=DIVIDE!:(ONE,X,K2) >>;
- IF LESSP!:(Y,EE) THEN <<M:=0; Z:=Y>>
- ELSE <<IF (M:=(ORDER!:(Y)*23)/10) = 0 THEN Z:=Y
- ELSE << SAVE!:P:=!:PREC!:; !:PREC!::=K2;
- Z:=DIVIDE!:(Y , TEXPT!:(EE,M) , K2);
- !:PREC!::=SAVE!:P>>;
- WHILE GREATERP!:(Z,EE) DO
- <<M:=M+1; Z:=DIVIDE!:(Z,EE,K2) >> >>;
- BEGIN INTEGER NFIG,N; SCALAR DCUT,DX,DY,X0;
-
- DCUT:=MAKE!:BF(10,-K2);
- DY :=MAKE!:BF(20,-K2);
- Y:=DIVIDE!:( DIFFERENCE!:(Z,ONE) ,
- CONV!:S2BF("1.72") , 2);
-
- NFIG:=1;
- WHILE NFIG<K2 OR GREATERP!:( ABS!:(DY) , DCUT) DO
- <<IF (NFIG:=2*NFIG) > K2 THEN NFIG:=K2;
- X0:=EXP!:(Y,NFIG);
- DX:=DIFFERENCE!:(Z,X0);
- N :=MAX(1 , NFIG+ORDER!:(DX));
- DY:=DIVIDE!:(DX,X0,N);
- Y :=PLUS!:(Y,DY) >>;
- END;
- Y:=PLUS!:( CONV!:I2BF(M) , Y);
- RETURN ROUND!:MT( TIMES!:(SIGN,Y) , K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE SIN!:(X,K); %****************************$
- %=================================================$
- % This function calculates sin(x), the value of $
- % the sine function at the point "x", with $
- % the precision K, by summing terms of the $
- % Taylor series for sin(z), 0 < z < PI/4. $
- % X is a BIG-FLOAT representation of any rael "x",$
- % otherwise it is converted to a <BIG-FLOAT>.$
- % K is a positive integer. $
- %=================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'SIN!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(0) ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( SIN!:( MINUS!:(X) , K)) ELSE
- BEGIN INTEGER K2,M; SCALAR PI4,SIGN,Q,R,Y;
- K2 :=K+2;
- M :=PRECI!:(X);
- PI4:=TIMES!:( !:PI(K2+M) , CONV!:S2BF("0.25"));
- IF LESSP!:(X,PI4) THEN <<M:=0; R:=X>>
- ELSE <<M:=CONV!:BF2I( Q:=QUOTIENT!:(X,PI4));
- R:=DIFFERENCE!:(X , TIMES!:(Q,PI4)) >>;
- SIGN:=CONV!:I2BF(1);
- IF M>=8 THEN M:=REMAINDER(M,8);
- IF M>=4 THEN <<SIGN:=MINUS!:(SIGN); M:=M-4>>;
- IF M=0 THEN GOTO SN ELSE IF M=1 THEN GOTO M1 ELSE
- IF M=2 THEN GOTO M2 ELSE GOTO M3;
- M1: R:=CUT!:MT( DIFFERENCE!:(PI4,R) , K2);
- RETURN TIMES!:(SIGN , COS!:(R,K));
- M2: R:=CUT!:MT(R,K2);
- RETURN TIMES!:(SIGN , COS!:(R,K));
- M3: R:=CUT!:MT( DIFFERENCE!:(PI4,R) , K2);
- SN: BEGIN INTEGER J,N,NCUT; SCALAR DCUT,FCTRIAL,RI,TM;
-
- NCUT:=K2-MIN(0 , ORDER!:(R)+1);
- DCUT:=MAKE!:BF(10,-NCUT);
- Y:=RI:=TM:=R;
- R:=MINUS!:( CUT!:EP( TIMES!:(R,R) , -NCUT));
-
- M:=J:=1;
- WHILE GREATERP!:( ABS!:(TM) , DCUT) DO
- <<J:=J+2;
- FCTRIAL:=CONV!:I2BF( M:=M*J*(J-1));
- RI:=CUT!:EP( TIMES!:(RI,R) , -NCUT);
- N :=MAX(1 , K2-ORDER!:(FCTRIAL)+ORDER!:(RI));
- TM:=DIVIDE!:(RI,FCTRIAL,N);
- Y :=PLUS!:(Y,TM); IF REMAINDER(J,20)=0 THEN
- Y:=CUT!:EP(Y,-NCUT) >>;
- END;
- RETURN ROUND!:MT( TIMES!:(SIGN,Y) , K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE COS!:(X,K); %****************************$
- %=================================================$
- % This function calculates cos(x), the value of $
- % the cosine function at the point "x", with $
- % the precision K, by summing terms of the $
- % Taylor series for cos(z), 0 < z < PI/4. $
- % X is a BIG-FLOAT representation of any real "x",$
- % otherwise it is converted to a <BIG-FLOAT>.$
- % K is a positive integer. $
- %=================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'COS!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(1) ELSE
- IF MINUSP!:(X) THEN COS!:( MINUS!:(X) , K) ELSE
- BEGIN INTEGER K2,M; SCALAR PI4,SIGN,Q,R,Y;
- K2 :=K+2;
- M :=PRECI!:(X);
- PI4:=TIMES!:( !:PI(K2+M) , CONV!:S2BF("0.25"));
- IF LESSP!:(X,PI4) THEN <<M:=0; R:=X>>
- ELSE <<M:=CONV!:BF2I( Q:=QUOTIENT!:(X,PI4));
- R:=DIFFERENCE!:(X , TIMES!:(Q,PI4)) >>;
- SIGN:=CONV!:I2BF(1);
- IF M>=8 THEN M:=REMAINDER(M,8);
- IF M>=4 THEN <<SIGN:=MINUS!:(SIGN); M:=M-4>>;
- IF M>=2 THEN SIGN:=MINUS!:(SIGN);
- IF M=0 THEN GOTO CS ELSE IF M=1 THEN GOTO M1 ELSE
- IF M=2 THEN GOTO M2 ELSE GOTO M3;
- M1: R:=CUT!:MT( DIFFERENCE!:(PI4,R) , K2);
- RETURN TIMES!:(SIGN , SIN!:(R,K));
- M2: R:=CUT!:MT(R,K2);
- RETURN TIMES!:(SIGN , SIN!:(R,K));
- M3: R:=CUT!:MT( DIFFERENCE!:(PI4,R) , K2);
- CS: BEGIN INTEGER J,N; SCALAR DCUT,FCTRIAL,RI,TM;
-
- DCUT:=MAKE!:BF(10,-K2);
- Y:=RI:=TM:=CONV!:I2BF(1);
- R:=MINUS!:( CUT!:EP( TIMES!:(R,R) , -K2));
-
- M:=1;
- J:=0;
- WHILE GREATERP!:( ABS!:(TM) , DCUT) DO
- <<J:=J+2;
- FCTRIAL:=CONV!:I2BF( M:=M*J*(J-1));
- RI:=CUT!:EP( TIMES!:(RI,R) , -K2);
- N :=MAX(1 , K2-ORDER!:(FCTRIAL)+ORDER!:(RI));
- TM:=DIVIDE!:(RI,FCTRIAL,N);
- Y :=PLUS!:(Y,TM); IF REMAINDER(J,20)=0 THEN
- Y:=CUT!:EP(Y,-K2) >>;
- END;
- RETURN ROUND!:MT( TIMES!:(SIGN,Y) , K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE TAN!:(X,K); %****************************$
- %=================================================$
- % This function calculates tan(x), the value of $
- % the tangent function at the point "x", $
- % with the precision K, by calculating $
- % sin(x) or cos(x) = sin(PI/2-x). $
- % X is a BIG-FLOAT representation of any real "x",$
- % otherwise it is converted to a <BIG-FLOAT>.$
- % K is a positive integer. $
- %=================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'TAN!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(0) ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( TAN!:( MINUS!:(X) , K)) ELSE
- BEGIN INTEGER K2,M; SCALAR ONE,PI4,SIGN,Q,R;
- K2 :=K+2;
- ONE:=CONV!:I2BF(1);
- M :=PRECI!:(X);
- PI4:=TIMES!:( !:PI(K2+M) , CONV!:S2BF("0.25"));
- IF LESSP!:(X,PI4) THEN <<M:=0; R:=X>>
- ELSE <<M:=CONV!:BF2I( Q:=QUOTIENT!:(X,PI4));
- R:=DIFFERENCE!:(X , TIMES!:(Q,PI4)) >>;
- IF M>=4 THEN M:=REMAINDER(M,4);
- IF M>=2 THEN SIGN:=MINUS!:(ONE) ELSE SIGN:=ONE;
- IF M=1 OR M=3 THEN R:=DIFFERENCE!:(PI4,R);
- R:=CUT!:MT(R,K2);
- IF M=0 OR M=3 THEN GOTO M03 ELSE GOTO M12;
- M03: R:=SIN!:(R,K2);
- Q:=DIFFERENCE!:(ONE , TIMES!:(R,R));
- Q:=SQRT!:( CUT!:MT(Q,K2) , K2);
- RETURN TIMES!:(SIGN , DIVIDE!:(R,Q,K));
- M12: R:=SIN!:(R,K2);
- Q:=DIFFERENCE!:(ONE , TIMES!:(R,R));
- Q:=SQRT!:( CUT!:MT(Q,K2) , K2);
- RETURN TIMES!:(SIGN , DIVIDE!:(Q,R,K));
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE ASIN!:(X,K); %***************************$
- %==================================================$
- % This function calculates asin(x), the value of $
- % the arcsine function at the point "x", $
- % with the precision K, by calculating $
- % atan(x/SQRT(1-x**2)) by ATAN!:. $
- % The answer is in the range [-PI/2 , PI/2]. $
- % X is a BIG-FLOAT representation of "x", IxI <= 1,$
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %==================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- GREATERP!:( ABS!:(X) , CONV!:I2BF(1)) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'ASIN!: ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( ASIN!:( MINUS!:(X) , K)) ELSE
- BEGIN INTEGER K2; SCALAR ONE,Y;
- K2 :=K+2;
- ONE:=CONV!:I2BF(1);
- IF LESSP!:( DIFFERENCE!:(ONE,X) , MAKE!:BF(10,-K2))
- THEN RETURN ROUND!:MT
- ( TIMES!:( !:PI(K+1) , CONV!:S2BF("0.5")) , K);
- Y:=CUT!:MT( DIFFERENCE!:(ONE , TIMES!:(X,X)) , K2);
- Y:=DIVIDE!:(X , SQRT!:(Y,K2) , K2);
- RETURN ATAN!:(Y,K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE ACOS!:(X,K); %***************************$
- %==================================================$
- % This function calculates acos(x), the value of $
- % the arccosine function at the point "x", $
- % with the precision K, by calculating $
- % atan(SQRT(1-x**2)/x) if x > 0 or $
- % atan(SQRT(1-x**2)/x) + PI if x < 0. $
- % The answer is in the range [0 , PI]. $
- % X is a BIG-FLOAT representation of "x", IxI <= 1,$
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %==================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- GREATERP!:( ABS!:(X) , CONV!:I2BF(1)) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'ACOS!: ELSE
- BEGIN INTEGER K2; SCALAR Y;
- K2:=K+2;
- IF LESSP!:( ABS!:(X) , MAKE!:BF(50,-K2))
- THEN RETURN ROUND!:MT
- ( TIMES!:( !:PI(K+1) , CONV!:S2BF("0.5")) , K);
- Y:=DIFFERENCE!:( CONV!:I2BF(1) , TIMES!:(X,X));
- Y:=CUT!:MT(Y,K2);
- Y:=DIVIDE!:( SQRT!:(Y,K2) , ABS!:(X) , K2);
- RETURN (IF MINUSP!:(X) THEN ROUND!:MT
- ( DIFFERENCE!:( !:PI(K+1) , ATAN!:(Y,K)) , K)
- ELSE ATAN!:(Y,K) );
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE ATAN!:(X,K); %***************************$
- %====================================================$
- % This function calculates atan(x), the value of the $
- % arctangent function at the point "x", with $
- % the precision K, by summing terms of the $
- % Taylor series for atan(z) if 0 < z < 0.42. $
- % Otherwise the following identities are used: $
- % atan(x) = PI/2 - atan(1/x) if 1 < x and $
- % atan(x) = 2*atan(x/(1+SQRT(1+x**2))) $
- % if 0.42 <= x <= 1. $
- % The answer is in the range [-PI/2 , PI/2]. $
- % X is a BIG-FLOAT representation of any real "x", $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %====================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'ATAN!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(0) ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( ATAN!:( MINUS!:(X) , K)) ELSE
- BEGIN INTEGER K2; SCALAR ONE,PI4,Y,Z;
- K2 :=K+2;
- ONE:=CONV!:I2BF(1);
- PI4:=TIMES!:( !:PI(K2) , CONV!:S2BF("0.25"));
- IF EQUAL!:(X,ONE) THEN RETURN ROUND!:MT(PI4,K);
- IF GREATERP!:(X,ONE) THEN RETURN ROUND!:MT
- ( DIFFERENCE!:( PLUS!:(PI4,PI4) ,
- ATAN!:( DIVIDE!:(ONE,X,K2) , K+1)) , K);
- IF LESSP!:(X , CONV!:S2BF("0.42")) THEN GOTO AT;
- Y:=PLUS!:(ONE , CUT!:MT( TIMES!:(X,X) , K2));
- Y:=PLUS!:(ONE , SQRT!:(Y,K2));
- Y:=ATAN!:( DIVIDE!:(X,Y,K2) , K+1);
- RETURN ROUND!:MT( TIMES!:(Y , CONV!:I2BF(2)) , K);
- AT: BEGIN INTEGER M,N,NCUT; SCALAR DCUT,TM,ZI;
- NCUT:=K2-MIN(0 , ORDER!:(X)+1);
- Y:=TM:=ZI:=X;
- Z:=MINUS!:( CUT!:EP( TIMES!:(X,X) , -NCUT));
- DCUT:=MAKE!:BF(10,-NCUT);
- M:=1;
- WHILE GREATERP!:( ABS!:(TM) , DCUT) DO
- <<ZI:=CUT!:EP( TIMES!:(ZI,Z) , -NCUT);
- N :=MAX(1 , K2+ORDER!:(ZI));
- TM:=DIVIDE!:(ZI , CONV!:I2BF( M:=M+2) , N);
- Y :=PLUS!:(Y,TM); IF REMAINDER(M,20)=0 THEN
- Y:=CUT!:EP(Y,-NCUT) >>;
- END;
- RETURN ROUND!:MT(Y,K)
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE ARCSIN!:(X,K); %*************************$
- %==================================================$
- % This function calculates arcsin(x), the value of $
- % the arcsine function at the point "x", with $
- % the precision K, by solving $
- % x = sin(y) if 0 < x <= 0.72, or $
- % SQRT(1-x**2) = sin(y) if 0.72 < x, $
- % by Newton's iteration method. $
- % The answer is in the range [-PI/2 , PI/2]. $
- % X is a BIG-FLOAT representation of "x", IxI <= 1,$
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %==================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- GREATERP!:( ABS!:(X) , CONV!:I2BF(1)) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'ARCSIN!: ELSE
- IF ZEROP!:(X) THEN CONV!:I2BF(0) ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( ARCSIN!:( MINUS!:(X) , K)) ELSE
- BEGIN INTEGER K2; SCALAR DCUT,ONE,PI2,Y;
- K2 :=K+2;
- DCUT:=MAKE!:BF(10 , -K2+ORDER!:(X)+1);
- ONE :=CONV!:I2BF(1);
- PI2 :=TIMES!:( !:PI(K2+2) , CONV!:S2BF("0.5"));
- IF LESSP!:( DIFFERENCE!:(ONE,X) , DCUT) THEN
- RETURN ROUND!:MT(PI2,K);
- IF GREATERP!:(X , CONV!:S2BF("0.72")) THEN GOTO AC
- ELSE GOTO AS;
- AC: Y:=CUT!:MT( DIFFERENCE!:(ONE , TIMES!:(X,X)) , K2);
- Y:=ARCSIN!:( SQRT!:(Y,K2) , K);
- RETURN ROUND!:MT( DIFFERENCE!:(PI2,Y) , K);
- AS: BEGIN INTEGER NFIG,N; SCALAR CX,DX,DY,X0;
- DY:=ONE;
- Y :=X;
- NFIG:=1;
- WHILE NFIG<K2 OR GREATERP!:( ABS!:(DY) , DCUT) DO
- <<IF (NFIG:=2*NFIG) > K2 THEN NFIG:=K2;
- X0:=SIN!:(Y,NFIG);
- CX:=DIFFERENCE!:(ONE , TIMES!:(X0,X0));
- CX:=CUT!:MT(CX,NFIG);
- CX:=SQRT!:(CX,NFIG);
- DX:=DIFFERENCE!:(X,X0);
- N :=MAX(1 , NFIG+ORDER!:(DX));
- DY:=DIVIDE!:(DX,CX,N);
- Y :=PLUS!:(Y,DY) >>;
- END;
- RETURN ROUND!:MT(Y,K);
- END$
- %*************************************************************$
- SYMBOLIC PROCEDURE ARCCOS!:(X,K); %*************************$
- %====================================================$
- % This function calculates arccos(x), the value of $
- % the arccosine function at the point "x", with $
- % the precision K, by calculating $
- % arcsin(SQRT(1-x**2)) if x > 0.72 and $
- % PI/2 - arcsin(x) otherwise by ARCSIN!:. $
- % The answer is in the range [0 , PI]. $
- % X is a BIG-FLOAT representation of "x", IxI <= 1, $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %====================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- GREATERP!:( ABS!:(X) , CONV!:I2BF(1)) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'ARCCOS!: ELSE
- IF LEQ!:(X , CONV!:S2BF("0.72")) THEN
- ROUND!:MT( DIFFERENCE!:
- ( TIMES!:( !:PI(K+1) , CONV!:S2BF("0.5")) ,
- ARCSIN!:(X,K) ) , K)
- ELSE ARCSIN!:( SQRT!:( CUT!:MT
- ( DIFFERENCE!:( CONV!:I2BF(1) , TIMES!:(X,X)) ,
- K+2) , K+2) , K)$
- %*************************************************************$
- SYMBOLIC PROCEDURE ARCTAN!:(X,K); %*************************$
- %==================================================$
- % This function calculates arctan(x), the value of $
- % the arctangent function at the point "x", $
- % with the precision K, by calculating $
- % arcsin(x/SQRT(1+x**2)) by ARCSIN!: $
- % The answer is in the range [-PI/2 , PI/2]. $
- % X is a BIG-FLOAT representation of any real "x", $
- % otherwise it is converted to a <BIG-FLOAT>. $
- % K is a positive integer. $
- %==================================================$
- IF NOT( BFP!:( X:=CONV!:A2BF(X))) OR
- NOT( FIXP(K)) OR K<=0 THEN BFLERRMSG 'ARCTAN!: ELSE
- IF MINUSP!:(X) THEN
- MINUS!:( ARCTAN!:( MINUS!:(X) , K))
- ELSE ARCSIN!:( DIVIDE!:(X , SQRT!:( CUT!:MT
- ( PLUS!:( CONV!:I2BF(1) , TIMES!:(X,X)) ,
- K+2) , K+2) , K+2) , K)$
- END;
|