koreth@ssyx.ucsc.edu (Steven Grimm) (06/04/89)
Submitted-by: dbrooks@osf.org (David Brooks)
Posting-number: Volume 2, Issue 41
Archive-name: fplib20/part01
This is revision 2.0 of FPLIB, a full floating-point library for use
with Sozobon C on the Atari ST. It fixes some minor bugs in revision
1.0, improves the accuracy and speed of all floating point operations,
and includes better documentation.
#!/bin/sh
# shar: Shell Archiver (v1.22)
#
# This is part 1 of a multipart archive
# do not concatenate these parts, unpack them in order with /bin/sh
#
# Run the following text with /bin/sh to create:
# ATOF.C
# DNOTES
# EXP.C
# FABS.C
# FMOD.C
# FPADD.S
# FPCMP.S
# FPDIV.S
# FPLIB.H
# FPMUL.S
# FP_PRT.C
# ITRIG.C
# LIBM.uue
# LOG.C
# MATH.H
# POW.C
# README
# RNOTES
# SPECS
# SQRT.C
# TRIG.C
#
if test -r s2_seq_.tmp
then echo "Must unpack archives in sequence!"
next=`cat s2_seq_.tmp`; echo "Please unpack part $next next"
exit 1; fi
sed 's/^X//' << 'SHAR_EOF' > ATOF.C &&
X/* ATOF function for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X/* Hack warning. Hack warning. Hack warning... */
X/* Somehow I feel scanf should call atof, not the other way around. */
X/* Guys??? */
X
Xfloat atof(s)
Xchar *s;
X{ float res;
X
X sscanf(s, "%g", &res);
X return res;
X}
X
X_scandum()
X{
X _scanf(); /* Force floating version to load */
X}
SHAR_EOF
chmod 0600 ATOF.C || echo "restore of ATOF.C fails"
sed 's/^X//' << 'SHAR_EOF' > DNOTES &&
X Design notes for FPLIB Release 2.0
X ==================================
X
X
XNUMBER REPRESENTATION
X---------------------
X
XThe floating point format used is Motorola's "Fast Floating Point":
X
X -----------------------------------------
X | 24 bits mantissa |S| 7 bits exp |
X -----------------------------------------
X
XThe mantissa has an implicit binary point just off to the left. The
Xexponent is binary, stored excess-64. In other words, the actual
Xmathematical exponent has 64 added, and the resulting 7-bit number is
Xunsigned.
X
XNegative numbers are stored as sign-and-magnitude.
X
XAll floating-point operations expect normalized numbers, and generate
Xnormalized numbers if the inputs are normalized. A number is
Xnormalized if the most significant bit is 1 (except for 0.0, which is
Xstored as 32 bits of zero, never negative). Also, a nonzero mantissa
Xwith a biased exponent of 0 is not normalized, letting us test the low
Xbyte or whole long as a test for 0.0. This all means that the most
Xsignificant bit is actually redundant, and could be left off, giving
Xone extra bit of precision. Oh well.
X
XSome of these routines make use of the assumption that all arguments
Xare normalized, and may generate slightly bogus results if that isn't
Xtrue. Release 1.0 assumed that a biased exponent of zero could have
Xa nonzero mantissa; that has been fixed.
X
XSome examples: 1.0 is hex 80000041 , or 0.5 * 2 ** 1 (binary 0.1,
Xexponent 1 + bias of 64). -0.75 is hex C00000C0 (note the sign bit is
Xinverted, but no other bit is).
X
X
XACCURACY AND LIMITS
X-------------------
X
XThe accuracy is in the 24-bit mantissa, which gives between 7.22 and
X7.52 decimal digits of relative error. Decimal constants are given to
X8 digits, and the compiler is trusted to convert them right (it
Xgenerally does a reasonable job).
X
XThe range is in the exponent. The smallest nonzero and largest
Xpositive normalized numbers are ~5.4e-20 and ~9.2e18.
X
XRounding of ambiguous numbers uses round-to-even (see D. E. Knuth's
XSeminumerical Algorithms for a discussion on this). This means that
Xif the trailing portion is exactly one half of the least significant
Xbit, the lsb is forced to be the nearer even number. We often keep
Xseveral insignificant bits to be sure this is done correctly.
X
X
XCODING TRICKS
X-------------
X
XIt's quite normal in a library of this type to resort to bit-twiddling
Xin the internal representation. To help this, FPLIB.H defines a C
Xunion, thus:
X
Xunion {
X unsigned long ul;
X float f;
X char sc[4];
X };
X
XThe number can be looked at as a floating or integer long (for handling
Xbits in the mantissa), the sign is the sign of sc[3], and the exponent
Xis the rest of sc[3]. Other shortcuts are commented in the inividual
Xroutines.
X
X
XMATHEMATICAL INFORMATION
X------------------------
X
XWell, this is mostly what I remember from my undergraduate days. Many
Xof these functions have a precise mathematical formula or algorithm
Xthat will give the result, but is too slow for practical use. Instead,
Xwe look for an approximation that is quick to compute, and gives us
Xjust enough precision and no more. After all, why compute to 100
Xdigits accuracy when you can only represent less than 8? Some of these
Xapproximations look very strange, but they all have sound mathematical
Xbases. The name Chebyshev springs to mind.
X
XDedicated mathematicians have struggled for decades to bring you these
Xformulae and, just as important, a predictable accuracy of better than
X7.5 decimal digits. Because of this, we can execute a known number of
Xoperations instead of loop-and-test (as you would in, for example, a
Xsimple Newton's iteration for square root).
X
XThe approximations are generally polynomials, usually avoid division, and
Xare valid only over a small range of the input parameter. We use simple
Xmathematical identities to extend the range. For example, the cyclic
Xnature of the trigonometric functions means we only need calculate in the
Xrange 0..pi/2. For another example, we can do much of the square root
Xcalculation by the simple expedient of dividing the binary exponent by 2.
X
X
XRANDOM LIBRARY INFORMATION
X--------------------------
X
XRelease 1.0 of Sozobon C had a few problems in the floating support:
X
X- Two negative floating numbers compared wrong (returned > when it
X should be <). I corrected that. However, the other routines don't
X assume the correction has been installed in the library (which it
X wasn't, at first...).
X
X- The library was not optimally ordered, and had to be scanned twice.
X I straightened that out.
X
X- Add and subtract were truncating. Rounding is generally a better idea.
X
Xatof() calls sscanf(). Ick.
X
X-------------------------------------------------------------------------------
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
XMarch 1989
X
XDavid W. Brooks
X36 Elm Street
XMedfield, MA 02052
X
XBROOKS@CSSS-A.PRIME.COM
X{uunet,mit-eddie}!csss-a.prime.com!brooks
SHAR_EOF
chmod 0600 DNOTES || echo "restore of DNOTES fails"
sed 's/^X//' << 'SHAR_EOF' > EXP.C &&
X/* EXP and hyperbolic trig functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <fplib.h>
X#include <errno.h>
X
X/* EXP: */
X/* */
X/* Uses the synthesis (2**n) * (2**[m/8]) * b**g */
X/* where b = 2**(1/8) and 8*n+m+g = arg/ln(b). b**g has a cubic */
X/* approximation: source and accuracy unknown. */
X/* */
X/* Beware a bug in the standard release: can't compare two negative */
X/* floating numbers. */
X
Xfloat exp(a) register float a;
X{ fstruct res;
X register int aint;
X
X static fstruct huge = {HUGE_AS_INT};
X static float powtab[] = {1.0, 1.09050773, 1.1892071, 1.2968396,
X 1.41421356, 1.5422108, 1.68179283, 1.8340081};
X
X if (a > 43.5)
X { errno = ERANGE;
X return huge.f;
X }
X
X if ((a + 43.5) < 0.0)
X return 0.0; /* Underflow */
X
X if ((aint = (int)(a *= 11.5415603)) < 0) /* 8/ln(2) */
X --aint; /* Correct mathematically */
X a = (float)aint - a; /* -(frac part) */
X
X res.f = 1.0 + a * (-0.0866439378 + a * (0.003750577 +
X a * -0.11321783e-3));
X res.sc[3] += aint >> 3; /* Mult by 2**n */
X return res.f * powtab[aint&7]; /* Mult by 2**[m/8] */
X}
X
X/* HYPERBOLIC FUNCTIONS: virtually free */
X
Xfloat sinh(a) float a;
X{ register float ea = exp(a);
X
X return 0.5 * (ea - (1.0 / ea));
X}
X
Xfloat cosh(a) float a;
X{ register float ea = exp(a);
X
X return 0.5 * (ea + (1.0 / ea));
X}
X
Xfloat tanh(a) float a;
X{ register float e2a;
X
X e2a = exp(a);
X e2a = e2a * e2a; /* exp-squared-a */
X return (e2a - 1.0) / (e2a + 1.0);
X}
X
SHAR_EOF
chmod 0600 EXP.C || echo "restore of EXP.C fails"
sed 's/^X//' << 'SHAR_EOF' > FABS.C &&
X/* FABS, LDEXP, and FREXP functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X/* Mustn't allow the standard declaration of fabs()... */
X#define MATH_H
X
X#include <fplib.h>
X#include <errno.h>
X
X/* FABS: directly manipulate the sign bit, punning the type. Actually */
X/* the argument and return are both float. */
X
Xunsigned long fabs(a)
Xunsigned long a;
X{ return a & 0xFFFFFF7FL; }
X
X/* LDEXP: again direct manipulation. */
X
Xfloat ldexp(a, n)
Xfstruct a;
Xint n;
X{
X/* "long" to handle int overflow/underflow... */
X register long exp = (long)(a.sc[3] & EXP_MASK) + n;
X
X if (exp <= 0L)
X return 0L; /* Underflow */
X
X if (exp > 0x7FL)
X { errno = ERANGE; /* Overflow */
X a.ul |= 0xFFFFFF7FL; /* Preserve sign */
X return a.f;
X }
X
X a.sc[3] = (a.sc[3] & 0x80) | exp;
X return a.f;
X}
X
X/* FREXP: split into convenient mantissa and exponent */
X
Xfloat frexp(a, ip)
Xfstruct a;
Xint *ip;
X{ if (a.sc[3] == 0)
X *ip = 0;
X else
X { *ip = (a.sc[3] & 0x7f) - BIAS;
X a.sc[3] = (a.sc[3] & 0x80) | BIAS; /* preserve sign,
X force exponent */
X }
X
X return a.f;
X}
SHAR_EOF
chmod 0600 FABS.C || echo "restore of FABS.C fails"
sed 's/^X//' << 'SHAR_EOF' > FMOD.C &&
X/* CEIL, FLOOR, MODF, and FMOD functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <fplib.h>
X
X/* This macro returns true if the number is so large that is has no */
X/* representable fraction. */
X
X#define ISLGINT(a) (((a).sc[3] & EXP_MASK) > (BIAS + MANT_BITS))
X
X/* CEIL: smallest integral value not less than a. */
X
Xfloat ceil(a)
Xfstruct a;
X{ register float aint;
X
X if (ISLGINT(a))
X return a.f;
X
X if ((aint = (float)((long)a.f)) != a.f && a.sc[3] >= 0)
X ++aint;
X
X return aint;
X}
X
X/* FLOOR: largest integral value not greater than a. */
X
Xfloat floor(a)
Xfstruct a;
X{ register float aint;
X
X if (ISLGINT(a))
X return a.f;
X
X if ((aint = (float)((long)a.f)) != a.f && a.sc[3] < 0)
X --aint;
X
X return aint;
X}
X
X/* Just for once, we luck out: negative numbers behave as needed... */
X
Xfloat modf(a, ip)
Xfstruct a;
Xfloat *ip;
X{ long ipart;
X
X *ip = (ISLGINT(a))?a.f:(float)(long)a.f;
X
X return a.f - *ip;
X}
X
Xfloat fmod(n, d)
Xfloat n, d;
X{ fstruct quot;
X long iquot;
X
X quot.f = n / d;
X if (ISLGINT(quot))
X return 0.0;
X
X iquot = (long)quot.f;
X return n - (float)iquot * d;
X}
SHAR_EOF
chmod 0600 FMOD.C || echo "restore of FMOD.C fails"
sed 's/^X//' << 'SHAR_EOF' > FPADD.S &&
X* Copyright (c) 1988 by Sozobon, Limited. Author: Johann Ruegg
X*
X* Permission is granted to anyone to use this software for any purpose
X* on any computer system, and to redistribute it freely, with the
X* following restrictions:
X* 1) No charge may be made other than reasonable charges for reproduction.
X* 2) Modified versions must be clearly marked as such.
X* 3) The authors are not responsible for any harmful consequences
X* of using this software, even if they result from defects in it.
X*
X*
X* MODIFICATIONS:
X*
X* D.W.Brooks Mar 89 Implemented rounding, inlined one-call subroutines,
X* and several other speedups
X*
X* Rounding in the ambiguous case (tail byte = 0x80) uses round-to-even, as
X* Knuth insists (TAOCP3SNA).
X*
X* fpsub
X*
X .globl fpsub
X .globl _fpsub
Xfpsub:
X_fpsub:
X eor.b #$80,11(sp) * negate B, fall through (can produce -0.0)
X*
X* fpadd
X*
X .globl fpadd
X .globl _fpadd
Xfpadd:
X_fpadd:
X move.l d3,a1 * save d3
X move.l d4,a2 * and d4
X move.l 4(sp),d0 * First arg
X move.l 8(sp),d1 * Second arg
X move.l d0,d4 * Keeps the result's sign
X moveq.l #$7f,d3
X
X move.b d0,d2 * Compare signs only
X eor.b d1,d2
X bpl same_sign
X* different signs
X and.b d3,d0 * Force positive
X and.b d3,d1
X
X move.l d0,d2 * note d2 is copy of d0
X cmp.b d0,d1 * Compare magnitudes
X bne sk1
X cmp.l d0,d1
Xsk1:
X ble norev1
X move.l d1,d0 * swap d0/d1...d0 was saved in d2...
X move.l d2,d1
X move.b d0,d2 * ...now d2.b is copy of new d0.b
X not.b d4 * Invert sign
Xnorev1:
X* Begin unsigned subtract. Here we have (d0-d1), d0 >= d1, both positive.
X* d2.b has result exp. Subtract 0.0 works without needing a special case.
X
X move.b d0,d3
X sub.b d1,d3 * diff of exps (in low byte)
X
X clr.b d0 * Extract mantissas
X clr.b d1
X
X cmp.b #25,d3
X bgt out * Nonsignificant addend: d0 is correct mantissa
X lsr.l d3,d1 * Shift uses d3 mod 64
X sub.l d1,d0
X
X beq out0 * Zero: all done!
X bmi outsub * Already normalized
X* normalize loop
Xnloop:
X subq.b #1,d2
X ble underfl
X add.l d0,d0
X bpl nloop
X
Xoutsub:
X add.l #$80,d0 * Force a rounding step. This can't overflow.
X * (rum-ti-rum...can it?)
X*
X* Common exit point for checking rounding after adding 0x80. If the tail
X* was not 0x80, it's already the correct result. Otherwise, the tail is now
X* 0, and we implement the round-to-even.
X*
Xoutround:
X tst.b d0
X bne out
X and.w #$FE00,d0
X*
X* Common exit point. d0 has the mantissa, d2 the exponent and d4 the sign.
X*
Xout:
X move.b d2,d0
X and.b #$80,d4
X or.b d4,d0 * fix sign
X*
X* Come here when d0 has the result
X*
Xout0:
X move.l a1,d3 * restore
X move.l a2,d4
X rts
X
X* underflow
Xunderfl:
X moveq.l #0,d0
X bra out0
X*
X* same signs: do add.
X*
Xsame_sign:
X and.b d3,d0 * Force both positive
X and.b d3,d1
X
X move.l d0,d2
X cmp.b d0,d1 * Compare magnitudes
X bne sk2
X cmp.l d0,d1
Xsk2:
X ble norev2
X move.l d1,d0 * Swap d0/d1...
X move.l d2,d1
X move.b d0,d2 * ...and d2.b is a copy of d0.b
Xnorev2:
X* Begin unsigned add: (d0+d1), d0 >= d1
X sub.b d1,d0
X move.b d0,d3 * diff of exps
X
X clr.b d0 * Extract mantissas
X clr.b d1
X
X cmp.b #25,d3
X bge out * Nonsignificant addend: d0 has result mantissa
X lsr.l d3,d1
X add.l d1,d0
X bcc around
X roxr.l #1,d0 * Mantissa overflowed the word
X addq.b #1,d2
X bmi oflo * Exponent too hot to handle
X
Xaround: * Round the fraction.
X add.l #$80,d0 * This rounds unless d0.b was 0x80
X bcc outround * No overflow, but do final rounding
X roxr.l #1,d0 * Actually loads 0x800000nn
X addq.b #1,d2 * Compensate
X bpl out * Ok if no exponent overflow
X
Xoflo:
X moveq.l #$ffffffff,d0 * Handle exponent overflow, maintaining sign
X moveq.l #$7f,d2
X bra out
X
X .end
SHAR_EOF
chmod 0600 FPADD.S || echo "restore of FPADD.S fails"
sed 's/^X//' << 'SHAR_EOF' > FPCMP.S &&
X* Copyright (c) 1988 by Sozobon, Limited. Author: Johann Ruegg
X*
X* Permission is granted to anyone to use this software for any purpose
X* on any computer system, and to redistribute it freely, with the
X* following restrictions:
X* 1) No charge may be made other than reasonable charges for reproduction.
X* 2) Modified versions must be clearly marked as such.
X* 3) The authors are not responsible for any harmful consequences
X* of using this software, even if they result from defects in it.
X*
X* MODIFICATIONS:
X*
X* D.W.Brooks 3/5/89 Fixed fpcmp (when both -ve) and fpneg (when zero)
X*
X
X .globl fpcmp
X .globl _fpcmp
Xfpcmp:
X_fpcmp:
X move.l 4(sp),d0
X move.l 8(sp),d1
X* Start DWB
X move.b d0,d2 * See if both signs are negative
X and.b d1,d2
X bpl L0
X move.l d0,d2 * They are: switch arguments so the rest works
X move.l d1,d0
X move.l d2,d1
XL0:
X* End DWB
X cmp.b d1,d0
X bne L1
X cmp.l d1,d0
XL1:
X rts
X
X .globl fpneg
X .globl _fpneg
Xfpneg:
X_fpneg:
X move.l 4(sp),d0
X
X* DWB removed: tst.b d0 (unnecessary to test for 0.0)
X
X beq L2
X eor.b #$80,d0
XL2:
X rts
SHAR_EOF
chmod 0600 FPCMP.S || echo "restore of FPCMP.S fails"
sed 's/^X//' << 'SHAR_EOF' > FPDIV.S &&
X* FPDIV.S -- floating point divide
X* Copyright = David Brooks, 1989 All Rights Reserved
X*
X* This version works by a normal bitwise shift/subtract loop, with separate
X* handling of the exponent and a final normalization. 26 bits are computed,
X* to allow one normalization and one rounding step. A certain amount of
X* black magic and handwaving are involved, but please note that at no time
X* do any bits leave the confines of the computer.
X*
X* Round-to-even is used for the ambiguous rounding case.
X*
X* float fpdiv(num, denom); (_fpdiv for explicit call; fpdiv for compiled /)
X
X .globl fpdiv
X .globl _fpdiv
X
Xfpdiv:
X_fpdiv:
X move.l d3,a1 * We don't need no wimpy stack frame
X moveq.l #0,d0
X move.l 4(sp),d1
X beq exit * Numerator 0 - all done
X move.l 8(sp),d2
X beq div0 * Divide by 0 - overflow
X clr.b d1 * Work with mantissas
X clr.b d2
X move.l #$2000000,d3 * Position of msb of 26-bit field
Xcmpbit:
X cmp.l d2,d1 * Compare against new divisor
X blo nobit
X add.l d3,d0 * Add in this bit
X sub.l d2,d1 * and adjust
Xnobit:
X lsr.l #1,d2
X lsr.l #1,d3
X bne cmpbit * Done 26 bits?
X
X lsl.l #6,d0 * Reposition
X tst.l d1 * If there was a remainder...
X beq doexp
X addq.l #2,d0 * record a memory of it. The 2 bit
X * .won't be lost by either of the
X * ..subsequent normalizations
Xdoexp:
X moveq.l #$7F,d1 * Calculate new exponent
X move.l d1,d2
X and.b 7(sp),d1
X and.b 11(sp),d2
X sub.w d2,d1
X add.w #$41,d1 * Adjust
X tst.l d0 * Check for already normalized or zero
X bmi normok
Xnormloop:
X subq.w #1,d1 * Do one normalize step (there can't
X add.l d0,d0 * .be more)
Xnormok:
X add.l #$80,d0 * Round up in most cases
X bcc round1
X roxr.l #1,d0 * The rounding caused overflow, sigh
X addq.w #1,d1
Xround1:
X tst.b d0 * See if trailer was exactly 0x80
X bne ckexp
X and.w #$FE00,d0 * It was: round to even
X
Xckexp: * Check exponent for sanity
X tst.w d1
X ble underflow
X cmp.w #$7F,d1
X bgt overflow
Xsetexp:
X move.b d1,d0 * Set exponent
X
X move.b 7(sp),d1 * Get signs
X move.b 11(sp),d2
X eor.b d2,d1 * Form new sign
X and.b #$80,d1 * Extract it
X or.b d1,d0
Xexit:
X move.l a1,d3
X rts
Xdiv0:
Xoverflow:
X moveq.l #$FFFFFFFF,d0
X moveq.l #$7F,d1
X bra setexp
X
Xunderflow:
X moveq.l #0,d0
X bra exit
X
X .end
SHAR_EOF
chmod 0600 FPDIV.S || echo "restore of FPDIV.S fails"
sed 's/^X//' << 'SHAR_EOF' > FPLIB.H &&
X/* Common definitions for FPLIB */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#ifndef MATH_H
X#include <math.h> /* Now they're all defined */
X#endif
X
X/* Some facts about the representation */
X
X#define BIAS 0x40
X#define MANT_BITS 24
X#define EXP_MASK 0x7F
X#define HUGE_AS_INT 0xFFFFFF7F
X
X/* Redefine a 32-bit float number. Can initialize with unsigned long. */
X
Xtypedef union {
X unsigned long ul;
X float f;
X char sc[4]; /* Assuming char signed */
X } fstruct;
SHAR_EOF
chmod 0600 FPLIB.H || echo "restore of FPLIB.H fails"
sed 's/^X//' << 'SHAR_EOF' > FPMUL.S &&
X* Copyright (c) 1988 by Sozobon, Limited. Author: Johann Ruegg
X*
X* Permission is granted to anyone to use this software for any purpose
X* on any computer system, and to redistribute it freely, with the
X* following restrictions:
X* 1) No charge may be made other than reasonable charges for reproduction.
X* 2) Modified versions must be clearly marked as such.
X* 3) The authors are not responsible for any harmful consequences
X* of using this software, even if they result from defects in it.
X*
X* MODIFICATION:
X*
X* D.W.Brooks Apr 89 Tightened up some coding, and looked after some
X* trailing bits for the sake of textbook rounding.
X*
X* fpmul
X*
X .globl _fpmult
X .globl _fpmul
X .globl fpmul
Xfpmul:
X_fpmul:
X_fpmult:
X move.l d3,a0 * save scratch registers
X move.l d4,a1
X move.l d5,a2
X move.l 4(sp),d0
X move.l 8(sp),d1
X
X move.l #$7F,d2 * Calculate the presumed exponent
X move.l d2,d3
X and.w d0,d2
X beq ret * Mult by 0: d0 is correct
X and.w d1,d3
X beq ret0 * Mult by 0: set d0
X add.w d3,d2
X sub.w #$40,d2 * Remove the extra bias
X
X clr.b d0 * Calculate the absolute mantissa product
X clr.b d1
X move.w d0,d3 * Get least significant bits of result
X mulu.w d1,d3
X swap d3 * d3 range now 0...fe01
X
X move.w d0,d4
X move.w d1,d5
X swap d0
X swap d1
X mulu d0,d5 * calculate inner portion
X mulu d1,d4
X add.l d3,d4 * no carry since d3 <= fe01 && d4 <= fffe0001
X add.l d4,d5
X bcc nocar1
X move.w d5,d3 * Save trailing bits
X move.w #1,d5
X bra t1
Xnocar1:
X move.w d5,d3
X clr.w d5
Xt1:
X swap d5
X
X mulu d1,d0 * calculate most significant part
X add.l d5,d0
X bcc nocar2
X roxr.l #1,d0
X addq.w #1,d2 * Normalize down
Xnocar2:
X bmi norm
X add.l d0,d0 * only need at most 1 shift: started norm AB
X subq.w #1,d2
Xnorm:
X add.l #$80,d0 * Round uppppp....
X bcc nocar3
X roxr.l #1,d0
X addq.w #1,d2
Xnocar3:
X tst.b d0 * Check if trailer was exactly 0x80 (now 0x00)
X bne rebuild
X tst.w d3 * Check back with spare bits
X bne rebuild
X and.w #$FE00,d0 * Round to even
Xrebuild:
X tst.w d2 * Reconstruct the number. First test overflow
X ble underflow * Could be mult by 0
X cmp.w #$7F,d2
X bgt overflow
Xexpsign:
X move.b d2,d0 * Stuff exponent in
X move.b 7(sp),d1 * Calculate sign
X move.b 11(sp),d2
X eor.b d1,d2
X and.b #$80,d2
X or.b d2,d0
Xret:
X move.l a2,d5
X move.l a1,d4
X move.l a0,d3
X rts
X
Xret0:
Xunderflow:
X move.l #0,d0 * Underflow, return 0
X bra ret
X
Xoverflow:
X move.l #$7F,d2 * Overflow, return +/- huge
X move.l #$FFFFFFFF,d0
X bra expsign
X
X .end
SHAR_EOF
chmod 0600 FPMUL.S || echo "restore of FPMUL.S fails"
sed 's/^X//' << 'SHAR_EOF' > FP_PRT.C &&
X/* Modified: */
X/* DWB Apr 89 Handle unnormalized numbers */
X/* Fix exponent if .9... overflows */
X
X#define MAXEC 40
X
Xstatic char *f_buf;
Xstatic int f_upper;
X
Xstatic char ec_buf[MAXEC+10];
Xstatic int ec_sign, ec_exp;
X
Xfp_print(x, fmtc, prec, ptmp)
Xfloat x;
Xchar *ptmp;
X{ register long lx, m;
X register int e;
X
X/* These routines have trouble with unnormalized numbers. */
X
X lx = *(long*)&x;
X if ((m = lx & 0xffffff00L) == 0 || (e = lx & 0x7f) == 0)
X x = 0.0;
X else if (m > 0L)
X { do
X { if (--e == 0)
X { x = 0.0;
X goto uflo;
X }
X m <<= 1;
X } while (m > 0L);
X *(long*)&x = m | e | (lx & 0x80L);
Xuflo:
X }
X
X f_buf = ptmp;
X f_upper = 0;
X
X switch (fmtc) {
X case 'E':
X f_upper = 1;
X case 'e':
X e_print(x, prec);
X break;
X case 'F':
X case 'f':
X f_print(x, prec);
X break;
X case 'G':
X f_upper = 1;
X case 'g':
X g_print(x, prec);
X break;
X }
X}
X
Xstatic
Xe_print(x, prec)
Xfloat x;
X{
X int nsig;
X register char *p;
X
X if (prec < 0)
X nsig = 7;
X else
X nsig = prec+1;
X
X ec_pr(x, nsig, 0);
X
X p = f_buf;
X if (ec_sign)
X *p++ = '-';
X *p++ = ec_buf[0];
X *p++ = '.';
X if (nsig > 1)
X strcpy(p, &ec_buf[1]);
X p += strlen(p);
X *p++ = f_upper ? 'E' : 'e';
X ec_exp--;
X if (ec_exp < 0) {
X *p++ = '-';
X ec_exp = -ec_exp;
X }
X if (ec_exp < 10)
X *p++ = '0';
X sprintf(p, "%d", ec_exp);
X}
X
Xstatic
Xf_print(x, prec)
Xfloat x;
X{
X int nsig, nz, i;
X register char *p;
X
X if (prec < 0)
X nsig = 6;
X else
X nsig = prec;
X
X ec_pr(x, -nsig, 0);
X
X p = f_buf;
X if (ec_sign)
X *p++ = '-';
X if (ec_exp < 1) {
X *p++ = '0';
X } else {
X strncpy(p, ec_buf, ec_exp);
X p += ec_exp;
X }
X if (prec != 0 || nsig)
X *p++ = '.';
X if (nsig == 0) {
X *p = 0;
X return;
X }
X
X if (ec_exp < 0) {
X nz = -ec_exp;
X if (nz > nsig)
X nz = nsig;
X for (i=0; i<nz; i++)
X *p++ = '0';
X nsig -= nz;
X if (nsig > 0) {
X strncpy(p, ec_buf, nsig);
X p += nsig;
X }
X *p = 0;
X } else {
X strcpy(p, &ec_buf[ec_exp]);
X }
X}
X
Xstatic
Xg_print(x, nsig)
Xfloat x;
X{
X int prec;
X
X if (nsig < 0)
X nsig = 6;
X if (nsig < 1)
X nsig = 1;
X
X ec_pr(x, 1, 1);
X
X if (ec_exp < -3 || ec_exp > nsig)
X e_print(x, nsig-1);
X else {
X prec = nsig - ec_exp;
X f_print(x, prec);
X }
X}
X
X/*
X * given x, ndig
X * if ndig is > 0, indicates number of significant digits
X * else -ndig is number of digits we want to the right of dec. pt.
X * return the following:
X * appropriate number of digits of significance in ec_buf
X * ec_sign true if x was negative
X * ec_exp indicates the decimal point relative to leftmost digit
X */
Xstatic
Xec_pr(x, ndig, trunc)
Xfloat x;
X{
X int isneg;
X int nhave;
X long part;
X int exp, newexp;
X int ndigcpy = ndig;
X float rem;
X char tbuf[20];
X
X /* ndig must be >= 1 and <= MAXEC */
X if (x < 0.0) {
X isneg = 1;
X x = -x;
X } else
X isneg = 0;
X
X /* get some digits */
X somedig(x, &part, &rem, &exp);
X
X sprintf(ec_buf, "%ld", part);
X nhave = strlen(ec_buf);
X exp = nhave + exp;
X
X if (ndig <= 0) {
X ndig = -ndig;
X ndig += exp;
X }
X
X if (ndig < 1)
X ndig = 1;
X else if (ndig > MAXEC)
X ndig = MAXEC;
X
X /* get some more digits */
X while (nhave < ndig+1 && nhave < MAXEC) {
X if (rem == 0.0) {
X while (nhave < ndig+1)
X ec_buf[nhave++] = '0';
X ec_buf[nhave] = 0;
X break;
X }
X
X x = rem;
X somedig(x, &part, &rem, &newexp);
X
X sprintf(tbuf, "%ld", part);
X newexp = strlen(tbuf) + newexp;
X while (newexp++)
X ec_buf[nhave++] = '0';
X strcpy(&ec_buf[nhave], tbuf);
X nhave = strlen(ec_buf);
X }
X
X if (fround(ndig, trunc))
X { if (ndigcpy <= 0)
X strcpy(ec_buf+ndig, "0");
X exp++;
X }
X
X ec_sign = isneg;
X ec_exp = exp;
X}
X
Xstatic
Xfround(n, trunc)
X{
X char *p;
X int oflo = 0;
X
X p = &ec_buf[n];
X if (*p >= '5' && !trunc) {
X p--;
X while (p >= ec_buf) {
X *p += 1;
X if (*p < '9')
X goto done;
X *p = '0';
X p--;
X }
X ec_buf[0] = '1';
X oflo = 1;
X }
Xdone:
X ec_buf[n] = 0;
X return oflo;
X}
X
Xstatic
Xsomedig(x, lp, remp, expp)
Xfloat x;
Xlong *lp;
Xfloat *remp;
Xint *expp;
X{
X int bexp, dexp;
X long ipart;
X float rem;
X
X bexp = fgetexp(x);
X dexp = 0;
X
X while (bexp > 31) {
X x *= 1E-3;
X dexp += 3;
X bexp = fgetexp(x);
X }
X while (bexp < 10) {
X x *= 1E3;
X dexp -= 3;
X if (dexp < -24) {
X ipart = 0;
X dexp = 0;
X rem = 0.0;
X goto iszero;
X }
X bexp = fgetexp(x);
X }
X fsplit(x, &ipart, &rem);
Xiszero:
X *lp = ipart;
X *remp = rem;
X *expp = dexp;
X}
X
Xstatic
Xfgetexp(x)
Xfloat x;
X{
X char *p;
X int i;
X
X p = (char *)&x;
X i = p[3] & 0x7f;
X i -= 0x40;
X return i;
X}
X
Xstatic
Xfsplit(x, vp, rp)
Xfloat x, *rp;
Xlong *vp;
X{
X long ival;
X float rem;
X int bexp, neg;
X
X ival = *(long *)&x;
X neg = ival & 0x80;
X ival &= ~0xff;
X
X bexp = fgetexp(x);
X
X if (bexp < 1) {
X ival = 0;
X rem = x;
X } else {
X ival = (unsigned long)ival >> (32-bexp);
X if (neg)
X ival = -ival;
X rem = x - (float)ival;
X }
X
X *vp = ival;
X *rp = rem;
X}
SHAR_EOF
chmod 0600 FP_PRT.C || echo "restore of FP_PRT.C fails"
sed 's/^X//' << 'SHAR_EOF' > ITRIG.C &&
X/* ASIN, ACOS, ATAN, and ATAN2 functions for Sozobon C. */
X/* Copyright = David Brooks, 1989 All Rights Reserved */
X
X#include <fplib.h>
X#include <errno.h>
X
X/* ARCTAN: */
X/* */
X/* The method is from Hart et al: "Computer Approximations". It uses a */
X/* divide-and-conquer algorithm. The argument range is divided into */
X/* five using the partition points tan({1,3,5,7}*pi/16), a polynomial */
X/* valid over +/- pi/16 is calculated, and other magic is used to */
X/* reposition the result. Precision is >8 places. */
X
Xfloat atan(a) fstruct a;
X{ fstruct absval;
X register float tx0, tsq, temp;
X register int part;
X register char sign;
X
X static float adj[] = {0.0, 0.414213562, 1.0, 2.414213562};
X static float atof[] = {0.0, 0.39269908, 0.78539816,
X 1.178097245, 1.57079633};
X
X sign = a.sc[3] & 0x80;
X a.sc[3] &= 0x7f; /* get fabs(a) */
X tx0 = a.f;
X
X/* Figure out the partition */
X
X part = tx0<0.66817864?(tx0<0.19891237?0:1):(tx0<1.49660576?2:
X tx0<5.0273395?3:4);
X if (part == 4)
X tx0 = -1.0 / tx0;
X else if (part != 0)
X { temp = adj[part];
X tx0 = (tx0 - temp) / (tx0 * temp + 1.0);
X }
X
X/* Here's the calculation */
X
X tsq = tx0 * tx0 + 1.67784279;
X a.f = (0.93833093 / tsq + 0.44075154) * tx0 + atof[part];
X
X if (a.sc[3] != 0) /* Negate if negative */
X a.sc[3] ^= sign;
X return a.f;
X}
X
X/* ARCSIN and ARCCOS use the standard identities. There's some less */
X/* optimal code here becausethe released Sozobon C can't properly */
X/* compare two negative floating numbers. */
X
Xfloat asin(a) float a;
X{
X if (a > 1.0 || (a+1.0) < 0.0)
X { errno = EDOM;
X return 0.0;
X }
X
X return atan(a / sqrt(1.0 - a*a));
X}
X
Xfloat acos(a) float a;
X{ float temp;
X
X if (a > 1.0 || (a+1.0) < 0.0)
X { errno = EDOM;
X return 0.0;
X }
X
X temp = atan(sqrt(1.0 - a*a) / a);
X if (a >= 0.0)
X return temp;
X else
X return M_PI + temp;
X}
X
X/* ATAN2: */
X/* */
X/* Computes atan(quotient), and returns that for positive cos, else */
X/* extends the range depending on the sin. */
X
Xfloat atan2(s, c) fstruct s, c;
X{ register float r;
X
X if (c.sc[3] == 0) /* Infinite argument */
X return (s.sc[3]<0)?-M_PI_2:M_PI_2;
X
X r = atan(s.f / c.f);
X
X if (c.sc[3] >= 0)
X return r; /* Range -pi/2..pi/2 */
X if (s.sc[3] >= 0)
X return M_PI + r;
X return -M_PI + r;
X}
SHAR_EOF
chmod 0600 ITRIG.C || echo "restore of ITRIG.C fails"
sed 's/^X//' << 'SHAR_EOF' > LIBM.uue &&
Xtable
X !"#$%&'()*+,-./0123456789:;<=>?
X@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_
Xbegin 644 LIBM.ARC
XM&@A,24)-+D$ \@ "&H@ $R0 ,,2Z6"H!:94 ,_\K0D9/FC(LW !( " !@z
XM0$( D P$L 4 &* ,!HJ"@ (@--$!4XL?(/'I)S& ! < %@ 1) P/@ )0Ry
XM%\L%#%AB:=E&P@Q^_X!!) ' )-O # X\>*DS@L>HO(! O*")0(GN89">=("x
XMR#]^3P%5#%+U*U:M3[X$$R8SP)F*S.P * O6A=FLD AL[?J/WM-D#[*1Y5&Jw
XM2S9@3T?-H-7C::8,?'(]]53"$HVR],[FW?H"@&:]3SIC_LS9,U[0HNUF/EV:v
XM=.BYJEUWE@W[+EJ^\)38G7<& XP!,0$(.&, $M _^8I(4!<P?'D(.S" P/ u
XMGS^ZGTT8'_09"U>O^$ YP5Q@%L48<?2+8N/MFN^_*C[<\<>.^LG8,R-K/ATt
XM:%4 5]UWQ5IMG:$ "#R(!4 02S7U'P)_!<;@?9V%49T-8(BS'VP 1B@8;1;Zs
XM@X$+@#P8H@,D?M9@'5_,D88;0XD"40%?C/'&'#&6%P C7] 1!HQYR0B)><Y\r
XM488<<KB!$ "(S&@&'&W4P49Y34)2P)-CM%'4D%5>"4<89)!!Y8Q?/&D&'6],q
XMV669<+"!IAE,.@G'''6(,::53Y*1AEQ<0F3 >4!.J2@ P!*J*&!SHCHD 4 p
XMP $ BQ8*Z:2-#IIHHR1,>B@ E4IJ::,L:#JII99N^BD F48*:*>K7NKJJ:1*o
XMBJFHK=) JZNIBMJHK8AN.FNOG + JZBGYHKHKK<62ZRGDJHZ:@' 6FDJLZWRn
XM<&NC1"Q[;**MQNJK4=<*JVVBW@;+1+*RDIMNJ-$V&VRIHZ9[ZKGM.@MKNO2.m
XM.^VV5&B:QD %'5010PXI%-%$%2&@0$4:532 (>5E"LD!($$B$DDFH:3232\%l
XM1Y-P+/W14@LP(O //D^MTL O5-G5WGU=$+C06P' \-0RK[#B WNNQ<S6S LQk
XM (&U +QD P"VR$0 G " 4E^"%07S'GA/63,,(V2YS#.%IMWF%3V_E0,$$L! j
XM0]Z0$(! %01/AZ-*"S]@]A0\-LP!A'M<X^TU@+JU= 8!-K44'8!');6B$D.\i
XM]0\=88]=]MD$J-T"%"6EMV!9\&R-UFQ<8Z[Y9J]-?3)U_W@P]$+R^<7#+Y h
XM8$S6)_<LLUL(G"Z ?-(\!<HW]PSQ.6@^%UB[M0/(!\OI!,B'RH:69PV@:U[,g
XM;H8*_>75_(-Z@Q[\S [\] \G0Q6UXH/8=P[U@MD30&+YFP-@83^Z@.'+2%]=f
XM;Q7TTE-OK?7G.Q]@^[+[F5NZ!Q3PY45\3'&*5<AW/_,I:$( /$V*N/8^$\#'e
XM!9!C'^BX,#OB1.<K8.!$A+0A!/L,!0R*H!\_^O8 ,V"@;PM@PX%X )AL2"V!d
XM"*IA$,#@%Y8P0(.H<1\ ^($+[\"'A6TPP ?CPPP8]F2$);0+/SX#!EJ,4!@Fc
XMS L8A$0 :DP@&3]XH"8$D0<A# 4PU_@!8+ !!$LL8QQ!T.$7PN C( U 2#0*b
XM@XOL*"3SV&*.-L+1D/HH "7,L8XRJ$@?!R !+&GI3E[2$Y\&T*4G16E*?<+3a
XME\($R2?1R4Z9I-&1DK2D-<TA#G*@PYW^I*Y-N5)1L&05KM#UK7+]:ESELA>Uz
XMWK6N<<G2EN("5J.LY2QA+H26O*QE+XU)3%T%TYG-W%>\I%FO<#'$6;_<);*8y
XMB<Q39<N8[)IF*U^E350A4YK9;)6QQ"G-<MWK6^M\I2S3&:QH4LJ<U214H_)5x
XMS5,-:UO>#%>_L F <!:3G/ \9S(7>BI[JE.AC (7L-J9KG_>,Z#&'.@]#2K.w
XM=\XRG]QB:+KLF<N0RA*@[@(F1U]94GC!BY[#="8PF]FIFI;'IA'-Z0!L"H!3v
XMIE)@"B'80PZ&$0$<+2,;&9(Z(+( (X1$A2\8 S-@("BF. #*UI)2XP& ,2Qu
XM@0'>,R!1,B(.2\@$"+\)CELJ@HH.R.0'+9&"[WA B%^P(VX\, 0J0C&84- #t
XM'[\ 8FLZU\ (>JU^4A3L:[((&OA,P2[_6$]A09<:VX ./K\1!Q*F (1H0':)s
XMI(.!6P$ A):03@LM0J4J!TFF42JI(I6$DI0ZN4DQA3)/>UHE2]TUT9"VJU$,r
XM@.@NI751>1GW6Q0(UZ.JV:CD&G.Y'2UG/%T%78HFE+F.TA<NAPNMBL#A#7<0q
XMV$.$:C")8"0 Z&!84@?@"J9NZ6)?.4E*A$"$J;Z0)0]P@T\*&#ZC($4I"80Ap
XM&Z(( 9DQSCM_>$ /C-!M%BA@P@(ZU#$ XRSF0<'( . R!;PGQ].5GV 2)\1o
XMO;)"(IR! ;KYQS^(8YR9J!BT %C<%[Y[AQQ9R4A(>FV<;FPF-*F)3$]RTQO@n
XM%-LL;8F29&+#&]Z2E]A>$I)&PL.6"% E5K8TI5@&KG9WVRKGLC-8U0W6=%VZm
XMTP++-*7=5:DS+5JM+8>4N)8"0!FD+-Z@-F2HYJU( 4*%5(Y (&(0.0 JGDH2l
XML S@*>T%@!$4&ST!GH$"WA-*7D"0@Z%<=7S[(P"B*[+H#P\6+=MSRUC!P(RGk
XMX$(5.R""8KO2H0][YVA>%< 4?&<5) PCJ@P.\7U8#:'^/246<NA"+YY2#V;Tj
XM =/:<<;MC +Q7*N?<[N&F4Y]+]IBXX?OS%'9R&;5ID<0&QD,UMY(K<VQFY%i
XM/OJ JE5"!(P+\D)!+:M?_VAC[M"=9D7J!E"(_.!N>*]GWH2U;!!5E,!\(R!$h
XM/+B@*Q K\+V0V'X-7ZROI9B^%?UC2']XX"(J$8P@8.(&[PA"*?QAC2#4@@#Rg
XM"$(QSH""(%R#"/L(@CJ \;HH3QF/+7H1&LK3QP#8H$8WVCEKAQ0 //3H1T*_f
XMHXYTX<@CK\FUI903F&R+9$T^.91L.E.:=DRC(+^)ZYZL$VTEJ5N3CO.6!XUHe
XM.!MUS7!YU)_A&C-OB_LMDI;SF[<Z:*OX&2R-*K-;QQ5IFXW)9E<Y-%B%;VBXd
XM+ HO+CNS[02UNW412OE&<0&9C?+[/0_OSG("\U27QV[A8]K/X9I^\.*$Z3-1c
XM^BW-MRKTQ<5IF74J>YXJV2!+LG/!_(2PAHA#O0,#Q;@#[0="]Z5ZFO9T6CJXb
XMW^]-N-)YP?3)4Z[84+\% -AXRLGAT;+GW:?1!:H(+UH@D[P<_WPEG+CW#UN7a
XMB/,%95JC^*X?GECY'Q9E4!M ^E6QB6U0A0?4( D?$#<NDSZE\3N?-FWU5EFKz
XM<4(Q]D,\ &P5\ -:Y1* , Y T#HR$0AD\VQYT6"7]36C55H<Q@/;UWW5UEA4y
XM8X($P'VK1F*DPP',PP/>L 750(#Z5AUN (($8'%Y\06W9V,T<GLQ .*I"."x
XMT'10!G6P)2>?1%M3-W:Y=5NRA4E5UW5M\G55]G>NXF5I5V;0Q79N5TY>&%V w
XMAUS=)"_=-4]GAT]F&"P<U7FULEUUYTN5UX:'!RJ8!P!Y6%![R'G+1'F-YUO!v
XM@G>H9RIKV(9Q&'C*HDM<6%*D=T\ 8 9M\ 9D %3'M'M$Y3"R 'P+47R0L 5u
XM U_\4($=(Q-_ LR@05O@ "#<Q4)Q$#K]X&Z-A0*%W$/YF@:X$04 '%3=!KPt
XM 5JNH$*F"!.HJ(H P(JN>#^QN$ ?QH.W^(L.UD&[>!-LT(L3)W#!*$5@,(R%s
XM5HS!D8JK^ 8>\(H@ ",, * <(X%]@)L\V'C\XRUF!?> 0;00#_PUX[*!Q\Kr
XMH6*GJ&'(R(H&@$ .8G]Y\7XND&"_PX,'^&%@@ PU4@9I8(4X9P9*]@9R<(1Yq
XMH09?4(ED "=#UT5EXI$\IR,<D'4^QG5L(F1$)B=&!H6<1(5/2(5D1X575W56p
XMMDO41(C M$YZB$ZK]V4OA8=ON)-75GH/M2WA!$R)AV6LYRIX=Y1OII/CY'C8o
XMLH<^Z8=(Z2I\)Y6#2%"Y(H=PZ$RN!Y608@9A( 9S@(GDQ7OG10O -R2PP%08n
XM9C$;\HH5=G$#\ ?X9@7\,0#@B(ID(P,^A 3! V!H!U*, AN4%4 <&D)9 L6m
XM]@=O8 &1-B,+AW$ \HI@P V!.1-D PB#T $DR)G$L"%3Y0#LR !" 6OF%4<l
XM8XP:IH&DM9HP B?"0B <%:D:5K04"9I*4A*=V-L0 9S5A0 P$424"9R<)P:k
XM210X1DI-&!% .8AD]H6]E9UT%R^-$H:P=)U?9I14"9[86953N2]T-&1L>6?Ej
XMA3!#@@MQ.0 L,20$$&AU28I(X ;RYA\?MD1]61$$UQ0E<'V'-&0E>6,M,@8_i
XM0@9UP$4@T")SH*!NT)((&J$_0J$Y^2UM=RI>9IWI$F8'1!!N0 =F(#"E8!%5h
XMD"V;" &X 6>* #PP%2@!%^R<!(>P $F ",0@ (P0@'TI003\!8"P ,,, $ g
XM4 ),HPI"4 1"X#*C!0+_D \]P =ZJ0]"H)_\,*5ZZ0^CI0#_< ]+X [_< M"f
XM8 1", 12L A% 0$3(&Z3!@1" 00 '!,1PX$*=L2ILP( /6H@##, S2 0Ve
XM$ !*@*7$T0-#\ _^ ;*@#AF $Y 0!0F@]4-5HP$*55A*C_H ]@ MX.@'(d
XMR -N@HP < 9J@(P@< 9$@(PJ< 9&@(PK< 8X@(PM< 8J@(PN< 8:L)H/PJ8Oc
XM 5=.T#I<\ 120 1I2CHR,*4+ 4(QY@4Z8"T!0#I6T %3,#HQY@0VP &%^A5Gb
XM@ ":NJ@QY@.:RJDQ9@-3L B=BHQ>,"5(@XQ%<'V$@(Q'< 9\$)!F0 FC!0%?a
XMVH_W\*MST0!*8 1GP 4@T(O21@#>P \ H !^S=)P X64:0 @ =F\$,-\ (-z
XMX!H_8!?^L+&;VH_YX*_]B _^^@,%@+$=!D0)"P 6P$:DHP-% 0,< :@T )4y
XM\ ^R\+ 6@;$>ZP^^"@AP51:R\!DKRP!:.B2+R@G(* 9M4!WN@(QX "/^8 ZTx
XM*09D4P[KN <(T(/0P!&6P!'(P!%BLA"J,#!E2[8< ;'",;8"H+9MRQ&7)QQNw
XMZ[;")QQG^[9RRQ'DD+8<00Q\F[<+00/RP0@P0 !'@0-XL"(\ *W4T0_F\ -]v
XMBK$-*P(%ZX$(J[ ,*[ (H+/)4QWJH+.W4QWB ,40#9@H VC)0!?*A^<,%HFu
XM<P_R(0F+NQ"6&J7*:@";JJPF Z[^( ?R00<@P 14P @:,&@!$ 3ML!!" %-t
XM>C)/BJG5X0@J0 6<"P;JT+C:\ 4S)J(DBG,U0@=Y =EH))2@B86HI)S,!!Ls
XM< <9N2;;2Z+%N29/LKTOLB9"9B&9E*'CY*'B*2AIEF:>\K^OU%M:UEM?:9[Er
XM&2W=Y;]Z)Y1.Z< )C,!YUU$ /'GDN;]4:<'CB68#/,%=%IX2W"H@>GINYHABq
XM!L%IV,#?LI0!;,)B25 =T2@Q+!S!,L,V7,,X+,,Y3,,ZW,,\_,,W[,-!#,0Xp
XM['@N/'>$!\)B><!PUL+[*\#?F<'F22H5_&4]):$EBA"$<*($D*)N61$28"M]o
XM!B#S.0"8$&B,8'Q(<!(?@*,EPZ, X*-'H 048 8*$*>'VQXXA(XOL #!FHQ/n
XM@&%*0 1+(@Y7JJ@,D%Q(NA#&(!?/JF+*&@ GHZQ>:J5&$,F*BLGV<+0# !;[m
XM,P"MDY__( >RD6FB' :N$000&P!" P J8,>'[ ^'; ]2P BUO*84X*;DAJ=Tl
XM*A/# 0*Z0*7_8 8"L*1X2@%ZRJ<+.PS%$*@Z$ ]4 2C@PQ38,N,@,NZK#:\k
XM'!S,$33'K "E!0*AR@9F0*H#RZIGH 3'C M < 8B,,A?I01'0%4PP ,J)A_Aj
XMH:-]_ ,]**Q/ ,'(!_FP*M^W#K>(<CZZ0]_@P7MA0&U3,=PX@_<(#^Q/,NDi
XM@P?5# :F@(Q;T :@ +5NH FTN057F[5;ZP1="RE.D# JP-(N#2 M#=,OK3X5h
XM80!G;!'G8A$[H=,500 [+0!0P!'](M-$+0!(\-(($-,$@&%$G=0O?0 *2]0$g
XM,-1.+=5"0#I(<,SXC,V0H\URVLLJ\1;^( ](0 P(,A300 A>/:?! 0&*DS9\f
XM7-" C&&XT!9L8 WRP0RS&ZW5@0CR00B'K ]E?=9YD=9K#=8 0]E?=9YD=9Ke
XM#=80SE@ X08 @Q<+G$ "@.4)#@"X^ (P 8;& #^%L(! ($"-%X (!$ 40(d
XM$ @0@$DD<(8 #QX3E;CYI\_,1"T])"J@*60FOX@ 2M@,P" G $T"1I$J)"Ac
XM0X@2 [3!*"2(NW]A8+B@B018-&!"(" A!@*F2&B$2)I$J1( @C, !B#04 !b
XM$ @?@! ,C.,DUPA"7!YH@*($B)G$.SB%]B, 25%;/:CMSCP&30P"DPIH@3(a
XM&3"'SW#1RA78"P O!OP-#.6)5:PB&DA"!"V0";P%#R9<V/#AV+*!T:H]F5+Bz
XM6PQT[>+5RQ?K:I&#"X.A1!>,.IG__,'MITHI!A50&!0,\(^?8Q #B "ZK<1(y
XM8A6#I!RYV$^1CD%@[%4'\,\'4!AFV*"4;DWU!M5$;=@"%!YG4 4%F<X8!0"x
SHAR_EOF
echo "End of part 1, continue with part 2"
echo "2" > s2_seq_.tmp
exit 0