koreth@ssyx.ucsc.edu (Steven Grimm) (06/04/89)
Submitted-by: dbrooks@osf.org (David Brooks)
Posting-number: Volume 2, Issue 42
Archive-name: fplib20/part02
#!/bin/sh
# this is part 2 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
# file LIBM.uue continued
#
CurArch=2
if test ! -r s2_seq_.tmp
then echo "Please unpack part 1 first!"
exit 1; fi
( read Scheck
if test "$Scheck" != $CurArch
then echo "Please unpack part $Scheck next!"
exit 1;
else exit 0; fi
) < s2_seq_.tmp || exit 1
sed 's/^X//' << 'SHAR_EOF' >> LIBM.uue
XM/&V5QD7_L&$4!#1=Y(\PV)F'D3/B ;!331@1 U0+5P#35@ ] %&>3?YX M0*w
XM-*KBH21&(:6/AR),V"%&&&PXI#\.]%2>AP9X*,!O9A$@7$G$M1700+DQQ=M3v
XM0%%D$0#\:,211R")-!Q;*[7TTD1*.J$D/TI*80< 0KCR#Q9U_D.%DE"$%A<Fu
XMF6W6&5P ./+?5 '( 50.XZFA)!.D,9$.$%U]%!A)I#DQ#J5> 0,#!+!$!5< t
XM,(1ZU*@6F#H H>Z8RE)<OIA: *&2F&K &2Y$.FFE90X !*9;:<KI5SU \8\3s
XM VKIU$,8FI'IIKR&- @)"D!!5S[J(*==@#L \T/6T'QPE94C'MG:N#^(T6Zr
XM_#Q'@#>,2=""C$AX!UY[9RB@@IU(R'?1/ICL^P\28-#S@ED#*& N%>ZV)@,!q
XMDD@$ Q+!U#:OGAB6B[#"6V'1\!,/1PS Q!4'<C$64LP$!88L'"RMPJH!)I+#p
XM$$M,<6VW8<C BBV^&*.Z'NK#LXM1_<R/A^@ Y8(9"AA%'A,>&BK1"F<8L"+3o
XM1@WP#]08X0%E<&E1B28 +"%7UUT@F%H&C<STN+6'2E2GC@Y;X:.D/4)L)8>]n
XM?19Q!@X_$-"<'.X.=C$77;8A ;Y^*&$G%SU'-:@#V-D#ESXT"%'!1?JXL.T9m
XM FA 5I:[+7N@5 8(&ID%GP. >>6@:^ZOZR8$/GCA3QR.+P2._P,YT3%Y9@H8l
XMQ)#F1<'&@V'($E?)(04C#% @!<V :#)\V[3$ST 752]?0N.:<"-DMI+W[T<k
XM2E 5SZ2D$:/J0"<(86I 9P!JD0"M"0^&-J0GQ\(#:@"M<0&@9_)8RT4D,%6j
XMY"$+(+A!!%5H! (5^ \& H$-#""-/%)600T PSHR!X8=/&]R$WD9_0 @RHTi
XMYR%9K"D &"+#U\X2MK44AR-P\8<DDE4Z W6)#=:X"#[4(41O & 1TPB <@&h
XM@ (T\8F0.( ! ##%*E+QBE.$Q "RN$4L3L2+5A2 %\5(1B8N\0 0",:D8C$g
XM)4+B 4[48AP/$("0)""+":BC%O58@ D0 $+8$ #0!"7+Y@!#E\@PQOB0HJ0f
XM%,"0B'3#'!@9&"! \@MSL$,8V-!(D0#B"Y@<0QC<8!-1.!*3=)##&."0!P @e
XMXI1?<$,=VC"'5+X2$H_\PACHD <XE,&5IZRE*M$@!V#BD@UMJ ,; G/+1QXRd
XME*.\I0'0J$5J#L":V$RC-M4H+6YFTYO;[*8XJRG.;.*2F@7@YCD=B4XT.K&.c
XMV70B"<+I1!9HDYSR#"<YR6G.:X:3FTZDP3_OZ<]O#G2?!=5G0L'I3@ (E*'Cb
XM7&A$U9E0?@X4HOTT*$01BM%R>M2@'/VH1RUJT7@>Y8EB=&)*3ZI2E+JTI3!Ea
XMJ4Q72M.7SM2F-8UI3F^J4YRZ5*,B#>A%)ZI-H78TI"55*$4S2M2&/G2H))7Hz
XM4@D*3J,>M*+^1&HZB;K/K29UJDB5:E$=.E2F@M2?3H1)1^,H47R*]*AB=2(1y
XME I5L0*4K' %9U0;JM:R<K6B6^VH7L4:4J!.]:SJ]&I6LZI8IJZSL']U+&/-x
XMV5AV1M:C5MTH5O-*57KB]:I@72QBF_K6RX85J.^DJV8ERUK3;M:S\.SL8 W;w
XMV:2F=K:N':U;F5!;T8*VMW<-J5Q_&]J&"A< <\5M<<7)UKLF5[5TW>ME6TO2v
XMRDJWG9C][&JCZ]O'MK6[OMUN6KF[7:!"%K*T3>]TH>O6OIKWM<H=Z6NK.]KEu
XM+O>ZI?TJ<65;5]PF=9VZ_:Y9G_C4[:*7O5G5Q#/A( <ZN&"1=41B D)B@ #8t
XM(BX/0$5<P*"!N!R@F%IT@+00QHB0*, )5H@+$L[! 0!8Q2TF@ !Z@,&+?V@1s
XM "B(GP-FW*L_Z U5-$"R2"2P2,P!T<"8P06JE% FC%(*;RFQ=,Y9Y22>0(r
XM9T""VLY U.9X0P>,-49$L.$;P !T[P@A/J "X .&!<;@G#4= P&!"V^<TNq
XMB/-$:%!G;-P9S@@( Q3JK(P9,$:/ R &)!AZ"\V&0S>4 (A%OD//JC.# P8p
XMLB,1@"$UC(,0DG8#.I!0#! 0 !"%$,$ A-$K0 B$"]& @ @2$EH#@[")&)o
XM098) %T8Y0#9$<(9TN4/0,M9 $QHC0F83 E,'D .&B)!1C3 @OD0UH&F#8 n
XM7% B\OCC#1!P6369(C,"M(;<@6E(5ZX!#&=+BR52@ !,BF !8$SAV3IPMQ9Um
XMP 8):+L%,,B!M'1 !"#, !@#9T K%3!P-A! VS#X@<#W+>Z$H5LD\#*1";Z0l
XMYC7+&P!E ,KGLP *0#A9_[ D!!NK8]<YV:9 9"01 P0[#.0QA\%UXNQ)T('k
XM93-;W]"6-K45[N@!\#LNL("XKAG@F,\]G#&>FXD_S( (4Q@XQVO ]#Y38: j
XM#[S@1B,-/V A]3<8H =2+X^2]$$:?9"]/&Q0 ,0Y^"-VB#T24M_6&SQ [(H;i
XMX.+O8HP";KYNL$P #([P^KZCI44#D 0)XDI(N5N#(3A(?.!^!WS&!S]QHQ<>h
XM0S) NT&RXR%YW-H?2V?#V6=ND&[O>O5?9, /8!* V4_$V)!_@L+_T8^!MR&#g
XMG=<!V0WR!@J0? I T'F> PV ?M@BZR1?"+Y1WO?ER]D?; #.&ZM#=$SH(+Bf
XM)@"@/]:%8 A#5!.QA]-PO7QW-61>;FF./9KCCN;@8^<"H$+&PUVWS"/, .Z2e
XM<0?Q?^[W!&&W%?P #3.A?C)B#TOG!A) <F!W%QFT%?80#;GA>@P <S@1%09!d
XM<B9'0?S00&R@!26" G G!>%G+N[P,5=@?N@' (( @D @@B2(' C(0?S >%M$c
XM$A $C# /V+'@XY7$A!P9VW& '(&)@TQ*"2 '?+@!AH0?B1'A*?A+AQD#S-Pb
XM%?X@+3BPA?]@#P.7=2TP$TB& ,U!#_1G?SL7 "2P?^9R?PCS)#.A#0'(& S0a
XM''98;@VA@/] #RCD/AQ$#X."@^6A@T3H@T (!MR@AT,(#);2>(\'!1_C=P5Hz
XM-!X""R=(?!AA"D:18PS@(9R@)"F'$1S <FOW9&@(B;U2A$, #2V0!.5!$.7!y
XM ! #70!D]8$ I@!AXP!7;"#PS@"MA6'F4@(#X1;[0HC+>8 VR @_''BMBFx
XM%A#0=EFWC#\X,L"X), P 8P1 ]BV?C^B#+=F,-:'= V!<F]2(A^0'6^ @<Sw
XM"!+P!C<P?I.7.W>! %;@$\;&!^DH(RD'#26R$]C'!0>C#W&A!/;(&OCH%E?@v
XM$^+Q#^@ =RR0)W;C$YHF$G\X("!@)_B@%#)F)_2@%,87C%EG;!CPCZ3'"6G(u
XMAM87A,; <K'6 D@PBSX1!.AP1"B!!(#0%G_P,_I 4A(VRW%?J0=>4XD\;"t
XM#R" @#0F+8" 0@)@)S=6'L:F"BII#]T6AFS@ 7G"#V7H%O_PD;G1D3;)$6:)s
XM#] '$R"P$ O8%2* @."0"BT0!.6!'?20+TU9'NIQ,?Q@+NU2;BX " LI$A>#r
XM#R(' T'!!F2PF(U9 HT) )?$8&G@!J:$2ZBD2JQD3(\D3&Q0!@713)A$F6Y q
XM!V8@FL+D!JN4!\UT2&/0!ANA1:T)!VY0!G AFXYT2,FT3+B)F8?$!G3P!C8Qp
XMFW-0!V(@+=+46L$55_)E7(KE1$C44-%)6OA5G? %7")U6N1U5="I4-T97\%5o
XM68(U5M898-BEG$\TG=_96>O957<UG0W% /E1-(#4>MI4.MYG\R51.])4/KIn
XM5NKY17?58K15GM3978&5G>&5GF,%GQ'5GN!U6*1EG_S9H-M)7^MU7A'*GA7:m
XM4/,T8%75H0Q*H?(Y5?DIHOED4@Y*6>CY6 GJ72 J+?H)H1(:HP"*G1J*7<VUl
XMG_;D7Y.UGP666T)Z6R/Z1 2ZGP%:GSBZH>8Y3M8EGMF9H#OJ5CUZH=<Y8/]Yk
XMG?J)GA)*I!8UHP;16Q"ZI9:EG6:*5@*ZGT=Z7$I*70*670XJ75\JHO'I64EZj
XMH4_JIEA*IQGZ6MUIH/A%IO_%I#):H0#6IRQ*GKXU7O?UIN<)18L* .[%52\:i
XMJ3 JI(,:7]K9J#:*H._56E.*II<*5I6*76?JHQ*:H]@Y7"8* +QE3@[P3',@h
XM2F[P8'%! %4P5Y!081=V%+*P81TV$3D&"0M@!"96CG#& 7HET,&"'%AE_CPg
XM9!@P!06! 85X,,S #,W0 W'X,8?)B' & >XRF-UZCZ<!F-XJ(_S@. "@/G%Af
XM",WZK&.I:PI *)PP9-M@#OW0CMI#C%K$<!30!AY@+N8XENDZE@3@KPG3KF" e
XM#?1Z!EHP9-FZK1N99Q0 C.VJ!( @$ 1;K@SYK=X0A^AZCWZI:V^!9^7A+GZYd
XME^VB9G4P9* 0%T3PAV! "D-6"0*0#V) L_Q!!I<TJZ-TF<X$!\ IG)RIF\J$c
XMG+D)!V% !HW9FP5P2+5YFP.0G$3E1"7*H9=J6U?;K@A6IZW:IAQZI*!:II:Zb
XM57%4JOLIMFYUJ'.ZM6_J5FP[7V@ZMZ.5HA,*15)J64 &!Z\)#TMW$AM%)KFZ\;&!P 25,[N[B+G>-97MWDy
XMFF_@!G9@JZ5 -KE*N+TJ *^$N,$: !U&K(U[8E9 $^R* &R@ D8!"-G1E&[1x
XME_EH;!I 0?Y ,:]V,?JPET?ILJ2H?2BF8N> 3BV?$RS$64& !B0=38 $Z Qw
XM:0;P:5,P!/]# #1F8TB$!,, "/C!#;?VEP-@ LO'8U!Y!BJ@ (, E:D1!C[Pv
XM?NK:._R 067AK%ID A]PC0BH'F @#QFL 1Q\E\'(!K+& W]@P/S#LF! #>^[u
XM%_&+%\M7P#=V&]Z+ $^I17_P$$ !!&W0;(-0OZYJ9OGKLD !!F] '#P#V#@t
XM#I@@$4 VRY$$30%>"P"NC!B)!4M*5T2H=4!\%YQB*A Y!D!L&Y3)<IE&\<s
XMG'4 '.L#+YKMO:EH)^ZI.!96B&J7X&:782\H#L:R*BJJ*FZH8F\G_@5R(><r
XM7H=$!FE0O(LDN P8;M:N!.!!,"JM,3:8HXK D.P$9-[NG! I1[!FX@! $@q
XM! )@ CP@1EHD"X)0!E(R")(@".*@".*P"&8@#^V "$H0"*!#!8 @!W] @$@p
XM#*=Q ,1P&@L@"4(@#1(1!!J[!@0P!4$ #8 0B=-"!@0@#@LA!=H,9 3,!UVHo
XM!$'P!B(@"Q+Q!VY@ 1 0 !(0S11 S;7PS##A:@%@ @FP9E7,S&! #AN1PI5\n
XMR;?KFW!@R7."QQ2FJ<RI7[_[6Y'L61;-7W^<T1AU2$U+!H$;89NLO'$1 +@0m
XMRK<<$HQK8@JCN@"P *9, D1@NI)+N3EF!W\@/;,@ &I !<"P:@,P3[\JNK( l
XM"&\@:YTK -)C! \' )80RW1B83"1 6[@ >:0"@N!S""P!@@P!0+P!B SFJ k
XM#^(,"$J@SBK1A1 @ ,CPSX! *- B=6!PEM#T$M#$,- $4M $=M=ISKN=*Sj
XM$ %@ 70RU;90U6S #%E=&V1@ .8, %(@ &O@ &=-!K 0V9.M!K)0Q73P!R)Wi
XM2,5YG%H$NZ)MG+R+2R#MM U- 8 22&=VJY]3PD:HQ=]V_N%VX]<7WZ%V^.Yh
XMT8-L5]N$M&PPTDG$R;QZTH:@T@,0L\2J B8& D,@ D1 D5@TZF, )]M @) g
XM#)Z!!\80!&? !]0P!)2@Q4(PU38@%\& !$/P$#H0 $A *4&@#," #,& #8.@f
XM#81 !@A@ P6@ [0'!@0 X$)0!$A0!,!0,85 SIHM!-P,SE,@!)==SN<L!&E-e
XM=4HP!&:PSNTL!&)-"_-\U6V=SVFTSS'=T\0 UP)0P2B0 ": 2-N/3PV55,d
XIQL1-!Q+MT+N9VE$+!SU>VA.-T<*MR+U]Y,!]Y,95Y-#5T5>JY($!&@#Pc
X b
Xend
SHAR_EOF
chmod 0600 LIBM.uue || echo "restore of LIBM.uue fails"
sed 's/^X//' << 'SHAR_EOF' > LOG.C &&
X/* LOG functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <fplib.h>
X#include <errno.h>
X
X/* LOG: */
X/* */
X/* This yanks out the exponent, performs a polynomial valid over +/- */
X/* sqrt(2), and adds the scaled exponent back in. Domain error on */
X/* arguments <= 0. */
X
Xfloat log(a) fstruct a;
X{ register float t1, t2, t3;
X static fstruct nhuge = {0xFFFFFFFFL}; /* Largest negative */
X
X if (a.f <= 0.0)
X { errno = EDOM; /* Impossible really */
X return nhuge.f;
X }
X
X if (a.f < M_SQRT2 && a.f > M_SQRT1_2)
X { t3 = 0.0; /* arg between sqrt(2) and sqrt(2)/2 */
X t2 = (t1 = a.f - 1.0) + 2.0;
X }
X else
X { t3 = (float)(((int)a.sc[3] << 1) - 0x81) * 0.34657359;
X /* ...(log2(a) * 2 + 1) * ln(2)/2 */
X a.sc[3] = BIAS; /* Scale back to range */
X t2 = a.f + M_SQRT1_2;
X t1 = a.f - M_SQRT1_2;
X }
X
X t1 /= t2;
X t2 = t1 * t1;
X
X return t3 + t1 * (2.0000008 + t2 * (0.66644078 + t2 * 0.41517739));
X}
X
X/* LOG10 is easy. */
X
Xfloat log10(a) float a;
X{
X return M_LOG10E * log(a);
X}
SHAR_EOF
chmod 0600 LOG.C || echo "restore of LOG.C fails"
sed 's/^X//' << 'SHAR_EOF' > MATH.H &&
X/*
X * MATH.H Math function declarations
X */
X
X#ifndef MATH_H
X#define MATH_H
X
Xextern float sin(), cos(), tan();
Xextern float asin(), acos(), atan(), atan2();
Xextern float exp(), sinh(), cosh(), tanh();
Xextern float log(), log10(), pow();
Xextern float sqrt();
Xextern float ceil(), floor(), fabs();
Xextern float ldexp(), frexp(), modf(), fmod();
Xextern float atof();
X
X/* Some useful constants, generally to 8-9 digits */
X
X#define M_E 2.71828183
X#define M_LOG2E 1.44269504
X#define M_LOG10E 0.43429448
X#define M_LN2 0.69314718
X#define M_LN10 2.30258509
X#define M_PI 3.14159265
X#define M_PI_2 1.570796327
X#define M_PI_4 0.785398163
X#define M_1_PI 0.318309886
X#define M_2_PI 0.636619772
X#define M_2_SQRTPI 1.128379167
X#define M_SQRT2 1.41421356
X#define M_SQRT1_2 0.70710678
X/* This is what the compiler should turn into 0xFFFFFF7F... */
X#define MAXFLOAT ((float)9.2233715e+18)
X#define HUGE MAXFLOAT
X
X#define _ABS(x) ((x) < 0 ? -(x) : (x))
X
X#endif MATH_H
SHAR_EOF
chmod 0600 MATH.H || echo "restore of MATH.H fails"
sed 's/^X//' << 'SHAR_EOF' > POW.C &&
X/* POW function for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X/* */
X/* This uses exp and log in the obvious way. A domain error occurs if */
X/* a=0 and p<=0, or a<0 and p is not integral. */
X
X#include <fplib.h>
X#include <errno.h>
X
Xfloat pow(a, p)
Xfstruct a, p;
X{ fstruct r;
X register long pint;
X register int sign = 0;
X
X if (a.sc[3] == 0 && p.sc[3] <= 0)
X { errno = EDOM;
X return 0.0;
X }
X
X if (a.sc[3] < 0) /* ...if a.f < 0 */
X { pint = (long)p.f;
X if ((float)pint != p.f) /* ... if nonintegral */
X errno = EDOM;
X sign = pint & 1; /* ... carry on regardless */
X a.sc[3] &= EXP_MASK;
X }
X
X r.f = exp(log(a.f) * p.f);
X if (sign && (r.sc[3] != 0))
X r.sc[3] |= 0x80;
X return r.f;
X}
SHAR_EOF
chmod 0600 POW.C || echo "restore of POW.C fails"
sed 's/^X//' << 'SHAR_EOF' > README &&
X FLOATING POINT LIBRARY FOR SOZOBON C
X
X Rev 2.0 -- prepared April 30, 1989
X
XThe files contained in this archive are a set of floating point
Xfunctions for Sozobon C. We provide:
X
X sqrt exp floor
X sin sinh modf
X cos cosh fmod
X tan tanh fabs
X atan log ldexp
X atan2 log10 frexp
X asin pow atof
X acos ceil math.h
X
Xand modified versions of:
X
X fpcmp fpadd fpmul
X fpdiv fp_prt
X
XThese are available in source. There is also a new LIBM.A that has the
Xnew routines included (and the other routines re-arranged to remove the
Xneed to scan the library twice).
X
XThe complete archive also contains design notes, release notes for
Xrelease 2.0, and a brief set of specifications.
X
XOnly LIBM.A and MATH.H are needed for normal use of the library.
XLIBM.A should go in your LIB directory, and MATH.H in your INCLUDE
Xdirectory.
X
XThese routines may be used without restriction, but I retain the copyright
Xon them, and the Only True Upgrades will be those I make. If you have any
Ximprovements, please send them to me.
X
XDavid W. Brooks
X36 Elm Street
XMedfield, MA 02052, USA
X
Xdbrooks@osf.org
Xuunet!osf.org!brooks
SHAR_EOF
chmod 0600 README || echo "restore of README fails"
sed 's/^X//' << 'SHAR_EOF' > RNOTES &&
X Release notes for FPLIB Release 2.0
X ===================================
X
XThis release contains fixes to FPLIB release 1.0, and fixes to the
Xstandard Sozobon release 1.0 floating support.
X
XIn the standard support:
X
X- The original floating point add/subtract, multiply, and divide have
X been rewritten for better speed and accuracy.
X
X- Bugs in the floating version of printf have been fixed.
X
XThe remaining information is of interest only if you have been using
XFPLIB release 1.0. The following changes have been made:
X
X- MATH.H has been expanded to provide symbolic definitions of some
X common constants.
X
X- A few minor bugs have been fixed:
X
X * Some routines could return -0.0; the sign has been corrected.
X
X * Some routines could return a number with a biased exponent of
X zero; they now return true 0.0.
X
X * Out of range results are now signaled by ERANGE, not EDOM.
X
X * frexp() returned an incorrect value on overflow, and ldexp()
X didn't handle overflow or underflow correctly.
SHAR_EOF
chmod 0600 RNOTES || echo "restore of RNOTES fails"
sed 's/^X//' << 'SHAR_EOF' > SPECS &&
XFPLIB implements the set of routines described in K&R, Second Edition,
Xpp 250-251. In the following descriptions, x and y are of type float,
Xand n is of type int. All functions (including ceil and floor) return
Xa value of type float.
X
XRemember to declare any routine you use, either explicitly or by
Xincluding MATH.H. If you don't, the compiler will assume the function
Xreturns a value of type int, and will perform an unwanted conversion.
X
XA domain error occurs if an argument is outside the range for which
Xthe function is defined (e.g. sqrt(-1.0)). In this case, errno is
Xset to EDOM and the result is usually 0.0. A range error occurs if
Xthe result cannot be represented as a float (e.g. exp(100.0)). In
Xthis case, errno is set to ERANGE and the result is HUGE_VAL. If the
Xfunction succeeds, errno is not zeroed.
X
XIn all trigonometric functions, angles are expressed in radians. Pi
Xradians equal 180 degrees.
X
XAll routines work in all three resolutions :-)
X
Xsin(x) sine of x.
Xcos(x) cosine of x.
Xtan(x) tangent of x.
Xasin(x) inverse sine of x in the range -pi/2 to pi/2.
Xacos(x) inverse cosine of x in the range 0 to pi.
Xatan(x) inverse tangent of x in the range -pi/2 to pi/2.
Xatan2(x,y) inverse tangent of x/y in the range -pi to pi. The
X angle is chosen so that x has the sign of the sine and
X y has the sign of the cosine.
Xsinh(x) hyperbolic sine of x.
Xcosh(x) hyperbolic cosine of x.
Xtanh(x) hyperbolic tangent of x.
Xexp(x) e to the power of x.
Xlog(x) logarithm, base e, of x: x > 0.
Xlog10(x) logarithm, base 10, of x: x > 0.
Xpow(x,y) x to the yth power. A domain error occurs if x=0 and
X y<=0, or if x<0 and y is not an integer.
Xsqrt(x) square root of x: x>=0.
Xceil(x) smallest whole number not less than x.
Xfloor(x) largest whole number not greater than x.
Xfabs(x) absolute value of x.
Xldexp(x,n) x, times 2 to the nth power.
Xfrexp(x,pn) splits x into a fraction between 0.5 and 0.999..., which
X is the returned value, and a power of 2, which is stored
X in the int pointed to by pn. If x is 0, both parts of
X the result are 0.
Xmodf(x,py) splits x into integral and fractional parts, both with
X the same sign as x. The fraction is returned and the
X integral part is stored in the float pointed to by py.
Xfmod(x,y) remainder of x/y, with the same sign as y.
Xatof(ps) converts the string ps into a float.
X
XThe file MATH.H also defines several mathematical constants.
SHAR_EOF
chmod 0600 SPECS || echo "restore of SPECS fails"
sed 's/^X//' << 'SHAR_EOF' > SQRT.C &&
X/* SQRT function for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X/* */
X/* The method is from Hart et al: "Computer Approximations". We reduce */
X/* the range to [0.25..1.0], then apply the polynomial: */
X/* 0.25927688 + 1.0520212n - 0.31632214n^2, giving 2.3 digits */
X/* approximation. Then we do two Newton's iterations (giving ~9 digits */
X/* precision). */
X/* */
X/* A negative argument gives a return value of 0.0 and sets EDOM. */
X
X#include <fplib.h>
X#include <errno.h>
X
X/* Sadly we can't use register for objects of type fstruct. */
X
Xfloat sqrt(n)
Xfstruct n;
X{ register int exp; /* Separate unbiased exponent. */
X fstruct s;
X
X if (n.sc[3] == 0) /* Check for 0.0 */
X return 0.0;
X
X if ((exp = n.sc[3]) < 0) /* Test argument's sign. */
X { errno = EDOM;
X return 0.0;
X }
X
X exp -= BIAS; /* Get unbiased exponent. */
X if (exp & 1)
X { n.sc[3] = BIAS-1; /* Convert to 0.25..0.5. */
X ++exp;
X }
X else
X n.sc[3] = BIAS; /* Or convert to 0.5..1. */
X
X s.f = 0.25927688 + n.f * (1.0520212 + n.f * -0.31632214);
X
X s.f = s.f + n.f / s.f; /* One Newton. */
X s.sc[3] -= 1; /* (here's the divide by 2). */
X s.f = s.f + n.f / s.f; /* The other Newton. */
X s.sc[3] += ((unsigned int)exp >> 1) - 1;
X /* Divide by 2 and insert
X half the original exponent. */
X return s.f;
X}
SHAR_EOF
chmod 0600 SQRT.C || echo "restore of SQRT.C fails"
sed 's/^X//' << 'SHAR_EOF' > TRIG.C &&
X/* SIN, COS, and TAN functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <fplib.h>
X#include <errno.h>
X
X/* SIN:
X/* */
X/* The method is from Hart et al: "Computer Approximations". It reduces */
X/* the range to 0..pi/2, scaled to 0..1, and applies: */
X/* 1.57079629t - 0.64596336t^3 + 0.0796884805t^5 - 0.0046722279t^7 */
X/* + 0.00015082056t^9 (precision of 8.5 digits) and */
X/* finally fixes the sign. */
X/* */
X/* MODIFICATIONS: */
X/* DWB 03/15/89 Don't try to negate x if it's zero. */
X/* DWB 04/30/89 Fix cos(0.0) */
X
Xfloat sin(a) fstruct a;
X{ register float t, ipart, tsq;
X fstruct x;
X register unsigned long tint;
X register char sign;
X
X sign = a.sc[3] & 0x80;
X a.sc[3] ^= sign; /* sin(negative x) = -sin(-x) */
X
X/* If the argument is huge, the result is rather arbitrary -- and, */
X/* besides, the following code will break. */
X
X if (a.sc[3] >= (BIAS + MANT_BITS))
X { errno = ERANGE;
X return 0.0;
X }
X
X t = a.f * M_2_PI; /* Express in quarter-turns */
X
X/* Get integer part and fraction. Optimise the common case, in the
X first quadrant */
X
X if (t <= 1.0)
X tint = 0L;
X else
X { tint = t; /* tint = whole quarterturns */
X ipart = (float)tint;
X t = (tint & 1)?ipart - t + 1.0:t - ipart; /* Get fraction */
X }
X tsq = t * t;
X x.f = t * (1.57079629 + tsq * (-0.64596336 + tsq * (0.0796884805
X + tsq * (-0.46722279e-2 + tsq * 0.15082056e-3))));
X
X if (x.sc[3] != 0) /* Don't negate zero! */
X { if ((int)tint & 2) /* Quadrants 3 and 4 */
X x.sc[3] |= 0x80; /* Set negative */
X if (sign != 0)
X x.sc[3] ^= 0x80; /* Negate */
X }
X
X return x.f;
X}
X
X/* COS: */
X/* */
X/* Although approximations are available, this will do fine. I'm not */
X/* proud of the hack to fix the only rational argument with a rational */
X/* result, but then pride is the first deadly sin... */
X
Xfloat cos(a) float a;
X{
X return a==0.0?1.0:sin(M_PI_2 - a);
X}
X
X/* TAN: */
X/* (v-e-r-y s-l-o-w-l-y) */
X
Xfloat tan(a) float a;
X{
X return sin(a) / sin(M_PI_2 - a);
X}
SHAR_EOF
chmod 0600 TRIG.C || echo "restore of TRIG.C fails"
rm -f s2_seq_.tmp
echo "You have unpacked the last part"
exit 0