From decwrl!elroy.jpl.nasa.gov!usc!cs.utexas.edu!uunet!allbery Sat Mar 10 15:35:48 PST 1990 Article 1398 of comp.sources.misc: Path: decwrl!elroy.jpl.nasa.gov!usc!cs.utexas.edu!uunet!allbery From: hp@relay.EU.net@vmars.UUCP Newsgroups: comp.sources.misc Subject: v11i023: RPL Part 3 of 3 Message-ID: <80862@uunet.UU.NET> Date: 10 Mar 90 19:57:58 GMT Sender: allbery@uunet.UU.NET Lines: 1366 Approved: allbery@uunet.UU.NET (Brandon S. Allbery - comp.sources.misc) Posting-number: Volume 11, Issue 23 Submitted-by: hp@relay.EU.net@vmars.UUCP Archive-name: rpl/part03 #! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'Readme' <<'END_OF_FILE' A Few Tips for installing: X Turbo-C: X Use the Project file rpl.prj. X Model must be large. X Append .\port to your include directories X (You will need hjp.h) X Disable the warnings: X Suspicious pointer conversion. X Possibly incorrect assignment. X (Because there will be many of them). X I used a profiler posted last summer to X comp.sources.misc. If you don't have it, remove X prof_start () from rpl.c and lprof.obj from the X project file. X UNIX: X The Makefile should work. In the directory port X are a few files our system (ULTRIX 2.1) didn't X have in it's C library. If your C library has them, X you may get errors from the linker. X I wrote a summary of the troubles I had porting the X Program to ULTRIX into the file porting.tips. It X may be helpful to you. X Other systems: X You are on your own, sorry. X Bug Reports, wishes, questions, ideas, ... please send to: X hp@honey.tuwien.ac.at X...!mcsun!tuvie!asupa!honey!hp hp@vmars.uucp END_OF_FILE if test 939 -ne `wc -c <'Readme'`; then echo shar: \"'Readme'\" unpacked with wrong size! fi # end of 'Readme' fi if test ! -d 'port' ; then echo shar: Creating directory \"'port'\" mkdir 'port' fi if test -f 'port/cabs.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/cabs.c'\" else echo shar: Extracting \"'port/cabs.c'\" \(6712 characters\) sed "s/^X//" >'port/cabs.c' <<'END_OF_FILE' X/* X * Copyright (c) 1985 Regents of the University of California. X * All rights reserved. X * X * Redistribution and use in source and binary forms are permitted X * provided that the above copyright notice and this paragraph are X * duplicated in all such forms and that any documentation, X * advertising materials, and other materials related to such X * distribution and use acknowledge that the software was developed X * by the University of California, Berkeley. The name of the X * University may not be used to endorse or promote products derived X * from this software without specific prior written permission. X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED X * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. X * X * All recipients should regard themselves as participants in an ongoing X * research project and hence should feel obligated to report their X * experiences (good or bad) with these elementary function codes, using X * the sendbug(8) program, to the authors. X */ X X#ifndef lint static char sccsid[] = "@(#)cabs.c 5.3 (Berkeley) 6/30/88"; X#endif /* not lint */ X X/* HYPOT(X,Y) X * RETURN THE SQUARE ROOT OF X^2 + Y^2 WHERE Z=X+iY X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) X * CODED IN C BY K.C. NG, 11/28/84; X * REVISED BY K.C. NG, 7/12/85. X * X * Required system supported functions : X * copysign(x,y) X * finite(x) X * scalb(x,N) X * sqrt(x) X * X * Method : X * 1. replace x by |x| and y by |y|, and swap x and X * y if y > x (hence x is never smaller than y). X * 2. Hypot(x,y) is computed by: X * Case I, x/y > 2 X * X * y X * hypot = x + ----------------------------- X * 2 X * sqrt ( 1 + [x/y] ) + x/y X * X * Case II, x/y <= 2 X * y X * hypot = x + -------------------------------------------------- X * 2 X * [x/y] - 2 X * (sqrt(2)+1) + (x-y)/y + ----------------------------- X * 2 X * sqrt ( 1 + [x/y] ) + sqrt(2) X * X * X * X * Special cases: X * hypot(x,y) is INF if x or y is +INF or -INF; else X * hypot(x,y) is NAN if x or y is NAN. X * X * Accuracy: X * hypot(x,y) returns the sqrt(x^2+y^2) with error less than 1 ulps (units X * in the last place). See Kahan's "Interval Arithmetic Options in the X * Proposed IEEE Floating Point Arithmetic Standard", Interval Mathematics X * 1980, Edited by Karl L.E. Nickel, pp 99-128. (A faster but less accurate X * code follows in comments.) In a test run with 500,000 random arguments X * on a VAX, the maximum observed error was .959 ulps. X * X * Constants: X * The hexadecimal values are the intended ones for the following constants. X * The decimal values may be used, provided that the compiler will convert X * from decimal to binary accurately enough to produce the hexadecimal values X * shown. X */ X X#if defined(vax)||defined(tahoe) /* VAX D format */ X#ifdef vax X#define _0x(A,B) 0x/**/A/**/B X#else /* vax */ X#define _0x(A,B) 0x/**/B/**/A X#endif /* vax */ X/* static double */ X/* r2p1hi = 2.4142135623730950345E0 , Hex 2^ 2 * .9A827999FCEF32 */ X/* r2p1lo = 1.4349369327986523769E-17 , Hex 2^-55 * .84597D89B3754B */ X/* sqrt2 = 1.4142135623730950622E0 ; Hex 2^ 1 * .B504F333F9DE65 */ static long r2p1hix[] = { _0x(8279,411a), _0x(ef32,99fc)}; static long r2p1lox[] = { _0x(597d,2484), _0x(754b,89b3)}; static long sqrt2x[] = { _0x(04f3,40b5), _0x(de65,33f9)}; X#define r2p1hi (*(double*)r2p1hix) X#define r2p1lo (*(double*)r2p1lox) X#define sqrt2 (*(double*)sqrt2x) X#else /* defined(vax)||defined(tahoe) */ static double r2p1hi = 2.4142135623730949234E0 , /*Hex 2^1 * 1.3504F333F9DE6 */ r2p1lo = 1.2537167179050217666E-16 , /*Hex 2^-53 * 1.21165F626CDD5 */ sqrt2 = 1.4142135623730951455E0 ; /*Hex 2^ 0 * 1.6A09E667F3BCD */ X#endif /* defined(vax)||defined(tahoe) */ X double hypot(x,y) double x, y; X{ X static double zero=0, one=1, X small=1.0E-18; /* fl(1+small)==1 */ X static ibig=30; /* fl(1+2**(2*ibig))==1 */ X double copysign(),scalb(),logb(),sqrt(),t,r; X int finite(), exp; X X if(finite(x)) X if(finite(y)) X { X x=copysign(x,one); X y=copysign(y,one); X if(y > x) X { t=x; x=y; y=t; } X if(x == zero) return(zero); X if(y == zero) return(x); X exp= logb(x); X if(exp-(int)logb(y) > ibig ) X /* raise inexact flag and return |x| */ X { one+small; return(x); } X X /* start computing sqrt(x^2 + y^2) */ X r=x-y; X if(r>y) { /* x/y > 2 */ X r=x/y; X r=r+sqrt(one+r*r); } X else { /* 1 <= x/y <= 2 */ X r/=y; t=r*(r+2.0); X r+=t/(sqrt2+sqrt(2.0+t)); X r+=r2p1lo; r+=r2p1hi; } X X r=y/r; X return(x+r); X X } X X else if(y==y) /* y is +-INF */ X return(copysign(y,one)); X else X return(y); /* y is NaN and x is finite */ X X else if(x==x) /* x is +-INF */ X return (copysign(x,one)); X else if(finite(y)) X return(x); /* x is NaN, y is finite */ X#if !defined(vax)&&!defined(tahoe) X else if(y!=y) return(y); /* x and y is NaN */ X#endif /* !defined(vax)&&!defined(tahoe) */ X else return(copysign(y,one)); /* y is INF */ X} X X/* CABS(Z) X * RETURN THE ABSOLUTE VALUE OF THE COMPLEX NUMBER Z = X + iY X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) X * CODED IN C BY K.C. NG, 11/28/84. X * REVISED BY K.C. NG, 7/12/85. X * X * Required kernel function : X * hypot(x,y) X * X * Method : X * cabs(z) = hypot(x,y) . X */ X double cabs(z) struct { double x, y;} z; X{ X return hypot(z.x,z.y); X} X double z_abs(z) struct { double x,y;} *z; X{ X return hypot(z->x,z->y); X} X X/* A faster but less accurate version of cabs(x,y) */ X#if 0 double hypot(x,y) double x, y; X{ X static double zero=0, one=1; X small=1.0E-18; /* fl(1+small)==1 */ X static ibig=30; /* fl(1+2**(2*ibig))==1 */ X double copysign(),scalb(),logb(),sqrt(),temp; X int finite(), exp; X X if(finite(x)) X if(finite(y)) X { X x=copysign(x,one); X y=copysign(y,one); X if(y > x) X { temp=x; x=y; y=temp; } X if(x == zero) return(zero); X if(y == zero) return(x); X exp= logb(x); X x=scalb(x,-exp); X if(exp-(int)logb(y) > ibig ) X /* raise inexact flag and return |x| */ X { one+small; return(scalb(x,exp)); } X else y=scalb(y,-exp); X return(scalb(sqrt(x*x+y*y),exp)); X } X X else if(y==y) /* y is +-INF */ X return(copysign(y,one)); X else X return(y); /* y is NaN and x is finite */ X X else if(x==x) /* x is +-INF */ X return (copysign(x,one)); X else if(finite(y)) X return(x); /* x is NaN, y is finite */ X else if(y!=y) return(y); /* x and y is NaN */ X else return(copysign(y,one)); /* y is INF */ X} X#endif END_OF_FILE if test 6712 -ne `wc -c <'port/cabs.c'`; then echo shar: \"'port/cabs.c'\" unpacked with wrong size! fi # end of 'port/cabs.c' fi if test -f 'port/hjp.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/hjp.h'\" else echo shar: Extracting \"'port/hjp.h'\" \(652 characters\) sed "s/^X//" >'port/hjp.h' <<'END_OF_FILE' X/**** X hjp.h X X Functions, Macros and Types often used X X****/ X/**** X Modification history: X X 1.0 88-07-06 Not really the first release which was lost some time ago. hjp X 1.1 89-06-13 ushort added; X version changed to Version X eprintf added hjp X 1.2 90-03-06 #ifndef _TYPES_ added for UNIX-compatibility. X X****/ X X#ifndef HJP X X#define HJP X X#define schar static char X typedef unsigned char uchar; /* 8 bit */ X X#ifndef _TYPES_ typedef unsigned short ushort; /* 16 bit */ typedef unsigned int uint; /* 16 bit */ X#endif typedef unsigned long ulong; /* 32 bit */ X X#define Version(a) schar VER[] = a; X void prof_start (char * filename); X X#endif END_OF_FILE if test 652 -ne `wc -c <'port/hjp.h'`; then echo shar: \"'port/hjp.h'\" unpacked with wrong size! fi # end of 'port/hjp.h' fi if test -f 'port/itoa.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/itoa.c'\" else echo shar: Extracting \"'port/itoa.c'\" \(943 characters\) sed "s/^X//" >'port/itoa.c' <<'END_OF_FILE' X/* X itoa.c X X replacement for itoa and ultoa functions X (nonexistent on UNIX) X*/ X static char digits [] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"; X char * itoa (int value, char * string, int radix) X{ X char buffer [33]; X char * s = string, * b = buffer; X X if (radix < 2 || radix > 36) return 0; X if (value < 0) { X * s ++ = '-'; X value = - value; X } X X while (value) { X * b ++ = digits [value % radix]; X value /= radix; X } X if (b == buffer) { /* value was zero */ X *s ++ = '0'; X } else { X while (--b >= buffer) { X * s ++ = * b; X } X } X *s = 0; X X return string; X} X char * ultoa (unsigned long value, char * string, int radix) X{ X char buffer [33]; X char * s = string, * b = buffer; X X if (radix < 2 || radix > 36) return 0; X X while (value) { X * b ++ = digits [value % radix]; X value /= radix; X } X if (b == buffer) { /* value was zero */ X *s ++ = '0'; X } else { X while (--b >= buffer) { X * s ++ = * b; X } X } X *s = 0; X X return string; X} X END_OF_FILE if test 943 -ne `wc -c <'port/itoa.c'`; then echo shar: \"'port/itoa.c'\" unpacked with wrong size! fi # end of 'port/itoa.c' fi if test -f 'port/mem.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/mem.h'\" else echo shar: Extracting \"'port/mem.h'\" \(67 characters\) sed "s/^X//" >'port/mem.h' <<'END_OF_FILE' X/* X mem.h X X*/ X X#define memmove(dst, src, cnt) bcopy(src, dst, cnt) END_OF_FILE if test 67 -ne `wc -c <'port/mem.h'`; then echo shar: \"'port/mem.h'\" unpacked with wrong size! fi # end of 'port/mem.h' fi if test -f 'port/port.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/port.h'\" else echo shar: Extracting \"'port/port.h'\" \(533 characters\) sed "s/^X//" >'port/port.h' <<'END_OF_FILE' X X/* X port.h X X this file contains types and constants which may have to be X changed for other compilers or machines. X X*/ X X#ifndef I_port X#define I_port X X#if defined (HUGE_VAL) && defined (__TURBOC__) X#undef HUGE_VAL X#define HUGE_VAL 1.7976931348623157e+308 X#endif X X#ifndef HUGE_VAL X#define HUGE_VAL 1.7976931348623157e+308 X#endif X X#ifndef TINY_VAL X#define TINY_VAL 2.2250738585072014e-308 X#endif X X#ifdef unix X struct complex {double x, y;}; extern double cabs (struct complex z); X X#endif X X#ifdef mips X#define host_mips X#endif X X#endif END_OF_FILE if test 533 -ne `wc -c <'port/port.h'`; then echo shar: \"'port/port.h'\" unpacked with wrong size! fi # end of 'port/port.h' fi if test -f 'port/process.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/process.h'\" else echo shar: Extracting \"'port/process.h'\" \(25 characters\) sed "s/^X//" >'port/process.h' <<'END_OF_FILE' X/* X process.h X X dummy X*/ END_OF_FILE if test 25 -ne `wc -c <'port/process.h'`; then echo shar: \"'port/process.h'\" unpacked with wrong size! fi # end of 'port/process.h' fi if test -f 'port/stdlib.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/stdlib.h'\" else echo shar: Extracting \"'port/stdlib.h'\" \(100 characters\) sed "s/^X//" >'port/stdlib.h' <<'END_OF_FILE' X/* X stdlib.h X X*/ X X#define max(x,y) ((x) > (y) ? (x) : (y)) X#define min(x,y) ((x) < (y) ? (x) : (y)) END_OF_FILE if test 100 -ne `wc -c <'port/stdlib.h'`; then echo shar: \"'port/stdlib.h'\" unpacked with wrong size! fi # end of 'port/stdlib.h' fi if test -f 'port/strerror.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/strerror.c'\" else echo shar: Extracting \"'port/strerror.c'\" \(316 characters\) sed "s/^X//" >'port/strerror.c' <<'END_OF_FILE' X/* X strerror.c X X replacement for Turbo-Cs strerror function X (nonexistent on UNIX systems) X*/ X extern int sys_nerr; extern char * sys_errlist []; X char * strerror (int errno) X{ X char unknown [] = "Unknown error"; X X if (errno < 0 || errno >= sys_nerr) { X return unknown; X } else { X return sys_errlist [errno]; X } X} END_OF_FILE if test 316 -ne `wc -c <'port/strerror.c'`; then echo shar: \"'port/strerror.c'\" unpacked with wrong size! fi # end of 'port/strerror.c' fi if test -f 'port/strtoul.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/strtoul.c'\" else echo shar: Extracting \"'port/strtoul.c'\" \(411 characters\) sed "s/^X//" >'port/strtoul.c' <<'END_OF_FILE' X/* X strtoul.c X X replacement for Turbo-Cs strtoul function. X (missing on UNIX) X X*/ X X#include X unsigned long strtoul (char * s, char ** endptr, int radix) X{ X unsigned long ul = 0; X int limit = radix > 9 ? radix - 10 + 'A' : radix + '0'; X int c; X X while ((c = toupper (*s)) && isalnum (c) && c < limit) { X ul = ul * radix + (c > '9' ? c + 10 - 'A' : c - '0'); X s ++; X } X * endptr = s; X return ul; X} END_OF_FILE if test 411 -ne `wc -c <'port/strtoul.c'`; then echo shar: \"'port/strtoul.c'\" unpacked with wrong size! fi # end of 'port/strtoul.c' fi if test -f 'port/support.c' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'port/support.c'\" else echo shar: Extracting \"'port/support.c'\" \(17016 characters\) sed "s/^X//" >'port/support.c' <<'END_OF_FILE' X/* X * Copyright (c) 1985 Regents of the University of California. X * All rights reserved. X * X * Redistribution and use in source and binary forms are permitted X * provided that the above copyright notice and this paragraph are X * duplicated in all such forms and that any documentation, X * advertising materials, and other materials related to such X * distribution and use acknowledge that the software was developed X * by the University of California, Berkeley. The name of the X * University may not be used to endorse or promote products derived X * from this software without specific prior written permission. X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED X * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE. X * X * All recipients should regard themselves as participants in an ongoing X * research project and hence should feel obligated to report their X * experiences (good or bad) with these elementary function codes, using X * the sendbug(8) program, to the authors. X */ X X#ifndef lint static char sccsid[] = "@(#)support.c 5.3 (Berkeley) 6/30/88"; X#endif /* not lint */ X X/* X * Some IEEE standard 754 recommended functions and remainder and sqrt for X * supporting the C elementary functions. X ****************************************************************************** X * WARNING: X * These codes are developed (in double) to support the C elementary X * functions temporarily. They are not universal, and some of them are very X * slow (in particular, drem and sqrt is extremely inefficient). Each X * computer system should have its implementation of these functions using X * its own assembler. X ****************************************************************************** X * X * IEEE 754 required operations: X * drem(x,p) X * returns x REM y = x - [x/y]*y , where [x/y] is the integer X * nearest x/y; in half way case, choose the even one. X * sqrt(x) X * returns the square root of x correctly rounded according to X * the rounding mod. X * X * IEEE 754 recommended functions: X * (a) copysign(x,y) X * returns x with the sign of y. X * (b) scalb(x,N) X * returns x * (2**N), for integer values N. X * (c) logb(x) X * returns the unbiased exponent of x, a signed integer in X * double precision, except that logb(0) is -INF, logb(INF) X * is +INF, and logb(NAN) is that NAN. X * (d) finite(x) X * returns the value TRUE if -INF < x < +INF and returns X * FALSE otherwise. X * X * X * CODED IN C BY K.C. NG, 11/25/84; X * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85. X */ X X X#if defined(vax)||defined(tahoe) /* VAX D format */ X#include X static unsigned short msign=0x7fff , mexp =0x7f80 ; X static short prep1=57, gap=7, bias=129 ; X static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ; X#else /* defined(vax)||defined(tahoe) */ X static unsigned short msign=0x7fff, mexp =0x7ff0 ; X static short prep1=54, gap=4, bias=1023 ; X static double novf=1.7E308, nunf=3.0E-308,zero=0.0; X#endif /* defined(vax)||defined(tahoe) */ X double scalb(x,N) double x; int N; X{ X int k; X double scalb(); X X#ifdef national X unsigned short *px=(unsigned short *) &x + 3; X#else /* national */ X unsigned short *px=(unsigned short *) &x; X#endif /* national */ X X if( x == zero ) return(x); X X#if defined(vax)||defined(tahoe) X if( (k= *px & mexp ) != ~msign ) { X if (N < -260) X return(nunf*nunf); X else if (N > 260) { X extern double infnan(),copysign(); X return(copysign(infnan(ERANGE),x)); X } X#else /* defined(vax)||defined(tahoe) */ X if( (k= *px & mexp ) != mexp ) { X if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf); X if( k == 0 ) { X x *= scalb(1.0,(int)prep1); N -= prep1; return(scalb(x,N));} X#endif /* defined(vax)||defined(tahoe) */ X X if((k = (k>>gap)+ N) > 0 ) X if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k< -prep1 ) X /* gradual underflow */ X {*px=(*px&~mexp)|(short)(1<>gap)-bias); X#else /* defined(vax)||defined(tahoe) */ X if( (k= *px & mexp ) != mexp ) X if ( k != 0 ) X return ( (k>>gap) - bias ); X else if( x != zero) X return ( -1022.0 ); X else X return(-(1.0/zero)); X else if(x != x) X return(x); X else X {*px &= msign; return(x);} X#endif /* defined(vax)||defined(tahoe) */ X} X finite(x) double x; X{ X#if defined(vax)||defined(tahoe) X return(1); X#else /* defined(vax)||defined(tahoe) */ X#ifdef national X return( (*((short *) &x+3 ) & mexp ) != mexp ); X#else /* national */ X return( (*((short *) &x ) & mexp ) != mexp ); X#endif /* national */ X#endif /* defined(vax)||defined(tahoe) */ X} X double drem(x,p) double x,p; X{ X short sign; X double hp,dp,tmp,drem(),scalb(); X unsigned short k; X#ifdef national X unsigned short X *px=(unsigned short *) &x +3, X *pp=(unsigned short *) &p +3, X *pd=(unsigned short *) &dp +3, X *pt=(unsigned short *) &tmp+3; X#else /* national */ X unsigned short X *px=(unsigned short *) &x , X *pp=(unsigned short *) &p , X *pd=(unsigned short *) &dp , X *pt=(unsigned short *) &tmp; X#endif /* national */ X X *pp &= msign ; X X#if defined(vax)||defined(tahoe) X if( ( *px & mexp ) == ~msign ) /* is x a reserved operand? */ X#else /* defined(vax)||defined(tahoe) */ X if( ( *px & mexp ) == mexp ) X#endif /* defined(vax)||defined(tahoe) */ X return (x-p)-(x-p); /* create nan if x is inf */ X if (p == zero) { X#if defined(vax)||defined(tahoe) X extern double infnan(); X return(infnan(EDOM)); X#else /* defined(vax)||defined(tahoe) */ X return zero/zero; X#endif /* defined(vax)||defined(tahoe) */ X } X X#if defined(vax)||defined(tahoe) X if( ( *pp & mexp ) == ~msign ) /* is p a reserved operand? */ X#else /* defined(vax)||defined(tahoe) */ X if( ( *pp & mexp ) == mexp ) X#endif /* defined(vax)||defined(tahoe) */ X { if (p != p) return p; else return x;} X X else if ( ((*pp & mexp)>>gap) <= 1 ) X /* subnormal p, or almost subnormal p */ X { double b; b=scalb(1.0,(int)prep1); X p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);} X else if ( p >= novf/2) X { p /= 2 ; x /= 2; return(drem(x,p)*2);} X else X { X dp=p+p; hp=p/2; X sign= *px & ~msign ; X *px &= msign ; X while ( x > dp ) X { X k=(*px & mexp) - (*pd & mexp) ; X tmp = dp ; X *pt += k ; X X#if defined(vax)||defined(tahoe) X if( x < tmp ) *pt -= 128 ; X#else /* defined(vax)||defined(tahoe) */ X if( x < tmp ) *pt -= 16 ; X#endif /* defined(vax)||defined(tahoe) */ X X x -= tmp ; X } X if ( x > hp ) X { x -= p ; if ( x >= hp ) x -= p ; } X X#if defined(vax)||defined(tahoe) X if (x) X#endif /* defined(vax)||defined(tahoe) */ X *px ^= sign; X return( x); X X } X} double sqrt(x) double x; X{ X double q,s,b,r; X double logb(),scalb(); X double t,zero=0.0; X int m,n,i,finite(); X#if defined(vax)||defined(tahoe) X int k=54; X#else /* defined(vax)||defined(tahoe) */ X int k=51; X#endif /* defined(vax)||defined(tahoe) */ X X /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ X if(x!=x||x==zero) return(x); X X /* sqrt(negative) is invalid */ X if(x1.0) t=1; /* b>1 : Round-to-(+INF) */ X if(t>=0) q+=r; } /* else: Round-to-nearest */ X else { X s *= 2; x *= 4; X t = (x-s)-1; X b=1.0+3*r/4; if(b==1.0) goto end; X b=1.0+r/4; if(b>1.0) t=1; X if(t>=0) q+=r; } X end: return(scalb(q,n)); X} X X#if 0 X/* DREM(X,Y) X * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) X * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS) X * INTENDED FOR ASSEMBLY LANGUAGE X * CODED IN C BY K.C. NG, 3/23/85, 4/8/85. X * X * Warning: this code should not get compiled in unless ALL of X * the following machine-dependent routines are supplied. X * X * Required machine dependent functions (not on a VAX): X * swapINX(i): save inexact flag and reset it to "i" X * swapENI(e): save inexact enable and reset it to "e" X */ X double drem(x,y) double x,y; X{ X X#ifdef national /* order of words in floating point number */ X static n0=3,n1=2,n2=1,n3=0; X#else /* VAX, SUN, ZILOG, TAHOE */ X static n0=0,n1=1,n2=2,n3=3; X#endif X X static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390; X static double zero=0.0; X double hy,y1,t,t1; X short k; X long n; X int i,e; X unsigned short xexp,yexp, *px =(unsigned short *) &x , X nx,nf, *py =(unsigned short *) &y , X sign, *pt =(unsigned short *) &t , X *pt1 =(unsigned short *) &t1 ; X X xexp = px[n0] & mexp ; /* exponent of x */ X yexp = py[n0] & mexp ; /* exponent of y */ X sign = px[n0] &0x8000; /* sign of x */ X X/* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */ X if(x!=x) return(x); if(y!=y) return(y); /* x or y is NaN */ X if( xexp == mexp ) return(zero/zero); /* x is INF */ X if(y==zero) return(y/y); X X/* save the inexact flag and inexact enable in i and e respectively X * and reset them to zero X */ X i=swapINX(0); e=swapENI(0); X X/* subnormal number */ X nx=0; X if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;} X X/* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */ X if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;} X X nf=nx; X py[n0] &= 0x7fff; X px[n0] &= 0x7fff; X X/* mask off the least significant 27 bits of y */ X t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t; X X/* LOOP: argument reduction on x whenever x > y */ loop: X while ( x > y ) X { X t=y; X t1=y1; X xexp=px[n0]&mexp; /* exponent of x */ X k=xexp-yexp-m25; X if(k>0) /* if x/y >= 2**26, scale up y so that x/y < 2**26 */ X {pt[n0]+=k;pt1[n0]+=k;} X n=x/t; x=(x-n*t1)-n*(t-t1); X } X /* end while (x > y) */ X X if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;} X X/* final adjustment */ X X hy=y/2.0; X if(x>hy||((x==hy)&&n%2==1)) x-=y; X px[n0] ^= sign; X if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;} X X/* restore inexact flag and inexact enable */ X swapINX(i); swapENI(e); X X return(x); X} X#endif X X#if 0 X/* SQRT X * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT X * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE X * CODED IN C BY K.C. NG, 3/22/85. X * X * Warning: this code should not get compiled in unless ALL of X * the following machine-dependent routines are supplied. X * X * Required machine dependent functions: X * swapINX(i) ...return the status of INEXACT flag and reset it to "i" X * swapRM(r) ...return the current Rounding Mode and reset it to "r" X * swapENI(e) ...return the status of inexact enable and reset it to "e" X * addc(t) ...perform t=t+1 regarding t as a 64 bit unsigned integer X * subc(t) ...perform t=t-1 regarding t as a 64 bit unsigned integer X */ X static unsigned long table[] = { X0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740, X58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478, X21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, }; X double newsqrt(x) double x; X{ X double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */ X long mx,scalx,mexp=0x7ff00000; X int i,j,r,e,swapINX(),swapRM(),swapENI(); X unsigned long *py=(unsigned long *) &y , X *pt=(unsigned long *) &t , X *px=(unsigned long *) &x ; X#ifdef national /* ordering of word in a floating point number */ X int n0=1, n1=0; X#else X int n0=0, n1=1; X#endif X/* Rounding Mode: RN ...round-to-nearest X * RZ ...round-towards 0 X * RP ...round-towards +INF X * RM ...round-towards -INF X */ X int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070 X * and a National 32081 & 16081 X */ X X/* exceptions */ X if(x!=x||x==0.0) return(x); /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */ X if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */ X if((mx=px[n0]&mexp)==mexp) return(x); /* sqrt(+INF) is +INF */ X X/* save, reset, initialize */ X e=swapENI(0); /* ...save and reset the inexact enable */ X i=swapINX(0); /* ...save INEXACT flag */ X r=swapRM(RN); /* ...save and reset the Rounding Mode to RN */ X scalx=0; X X/* subnormal number, scale up x to x*2**54 */ X if(mx==0) {x *= b54 ; scalx-=0x01b00000;} X X/* scale x to avoid intermediate over/underflow: X * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */ X if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;} X if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;} X X/* magic initial approximation to almost 8 sig. bits */ X py[n0]=(px[n0]>>1)+0x1ff80000; X py[n0]=py[n0]-table[(py[n0]>>15)&31]; X X/* Heron's rule once with correction to improve y to almost 18 sig. bits */ X t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; X X/* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */ X t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; X t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; X X/* twiddle last bit to force y correctly rounded */ X swapRM(RZ); /* ...set Rounding Mode to round-toward-zero */ X swapINX(0); /* ...clear INEXACT flag */ X swapENI(e); /* ...restore inexact enable status */ X t=x/y; /* ...chopped quotient, possibly inexact */ X j=swapINX(i); /* ...read and restore inexact flag */ X if(j==0) { if(t==y) goto end; else t=subc(t); } /* ...t=t-ulp */ X b54+0.1; /* ..trigger inexact flag, sqrt(x) is inexact */ X if(r==RN) t=addc(t); /* ...t=t+ulp */ X else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */ X y=y+t; /* ...chopped sum */ X py[n0]=py[n0]-0x00100000; /* ...correctly rounded sqrt(x) */ end: py[n0]=py[n0]+scalx; /* ...scale back y */ X swapRM(r); /* ...restore Rounding Mode */ X return(y); X} X#endif END_OF_FILE if test 17016 -ne `wc -c <'port/support.c'`; then echo shar: \"'port/support.c'\" unpacked with wrong size! fi # end of 'port/support.c' fi if test ! -d 'samples' ; then echo shar: Creating directory \"'samples'\" mkdir 'samples' fi if test -f 'samples/a.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/a.rpl'\" else echo shar: Extracting \"'samples/a.rpl'\" \(38 characters\) sed "s/^X//" >'samples/a.rpl' <<'END_OF_FILE' X<< DUP 0 > << 1 - A >> IFT >> 'A' STO END_OF_FILE if test 38 -ne `wc -c <'samples/a.rpl'`; then echo shar: \"'samples/a.rpl'\" unpacked with wrong size! fi # end of 'samples/a.rpl' fi if test -f 'samples/ack.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/ack.rpl'\" else echo shar: Extracting \"'samples/ack.rpl'\" \(386 characters\) sed "s/^X//" >'samples/ack.rpl' <<'END_OF_FILE' X<< X; X; ACK -- compute Ackermann function X; X -> a b ; two arguments X << IF a 0 == ; if first is zero ... X THEN X b 1 + ; it is easy X ELSE X IF b 0 == ; if second is zero ... X THEN X a 1 - 1 ACK ; it is a little more complicated. X ELSE ; if neither is zero ... X a 1 - a b 1 - ACK ACK ; it gets highly recursive. X ENDIF X ENDIF X >> X>> 'ACK STO END_OF_FILE if test 386 -ne `wc -c <'samples/ack.rpl'`; then echo shar: \"'samples/ack.rpl'\" unpacked with wrong size! fi # end of 'samples/ack.rpl' fi if test -f 'samples/ack0.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/ack0.rpl'\" else echo shar: Extracting \"'samples/ack0.rpl'\" \(166 characters\) sed "s/^X//" >'samples/ack0.rpl' <<'END_OF_FILE' X<< X -> a b X << IF a 0 == X THEN X b 1 + X ELSE X IF b 0 == X THEN X a 1 - 1 ACK X ELSE X a 1 - a b 1 - ACK ACK X ENDIF X ENDIF X >> X>> 'ACK STO END_OF_FILE if test 166 -ne `wc -c <'samples/ack0.rpl'`; then echo shar: \"'samples/ack0.rpl'\" unpacked with wrong size! fi # end of 'samples/ack0.rpl' fi if test -f 'samples/bench1.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/bench1.rpl'\" else echo shar: Extracting \"'samples/bench1.rpl'\" \(54 characters\) sed "s/^X//" >'samples/bench1.rpl' <<'END_OF_FILE' X<< TIME 1 10000 FOR a NEXT TIME - NEG >> 'BENCH1' STO END_OF_FILE if test 54 -ne `wc -c <'samples/bench1.rpl'`; then echo shar: \"'samples/bench1.rpl'\" unpacked with wrong size! fi # end of 'samples/bench1.rpl' fi if test -f 'samples/bench2.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/bench2.rpl'\" else echo shar: Extracting \"'samples/bench2.rpl'\" \(163 characters\) sed "s/^X//" >'samples/bench2.rpl' <<'END_OF_FILE' X<< TIME X 15 FIB DROP X TIME - NEG X>> 'BENCH2' STO X X<< X IF DUP 2 > X THEN X DUP 1 - FIB SWAP 2 - FIB + X ELSE X DROP 1 X ENDIF X>> 'FIB' STO END_OF_FILE if test 163 -ne `wc -c <'samples/bench2.rpl'`; then echo shar: \"'samples/bench2.rpl'\" unpacked with wrong size! fi # end of 'samples/bench2.rpl' fi if test -f 'samples/edit.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/edit.rpl'\" else echo shar: Extracting \"'samples/edit.rpl'\" \(230 characters\) sed "s/^X//" >'samples/edit.rpl' <<'END_OF_FILE' X<< X; X; EDIT -- edit object on level 1 using vi X; X "scratch" PURGE ; remove existing scratchfile X "scratch" SAVE ; save object to it X "vi scratch" SYSTEM ; edit it X "scratch" LOAD ; end get it back to level 1 X>> 'EDIT' STO END_OF_FILE if test 230 -ne `wc -c <'samples/edit.rpl'`; then echo shar: \"'samples/edit.rpl'\" unpacked with wrong size! fi # end of 'samples/edit.rpl' fi if test -f 'samples/hyp.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/hyp.rpl'\" else echo shar: Extracting \"'samples/hyp.rpl'\" \(34 characters\) sed "s/^X//" >'samples/hyp.rpl' <<'END_OF_FILE' X<< SQ SWAP SQ + SQRT >> 'HYP' STO END_OF_FILE if test 34 -ne `wc -c <'samples/hyp.rpl'`; then echo shar: \"'samples/hyp.rpl'\" unpacked with wrong size! fi # end of 'samples/hyp.rpl' fi if test -f 'samples/sin.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/sin.rpl'\" else echo shar: Extracting \"'samples/sin.rpl'\" \(58 characters\) sed "s/^X//" >'samples/sin.rpl' <<'END_OF_FILE' X"sintest.rpl" LOAD TIME 'a STO SINTEST CLEAR TIME a - OFF END_OF_FILE if test 58 -ne `wc -c <'samples/sin.rpl'`; then echo shar: \"'samples/sin.rpl'\" unpacked with wrong size! fi # end of 'samples/sin.rpl' fi if test -f 'samples/sintest.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/sintest.rpl'\" else echo shar: Extracting \"'samples/sintest.rpl'\" \(133 characters\) sed "s/^X//" >'samples/sintest.rpl' <<'END_OF_FILE' X<< X "sintest.out" PURGE X -6 6 FOR x X -6 6 FOR y X x y R->C SIN RE X "sintest.out" SAVE X NEXT X NEXT X>> 'SINTEST' STO END_OF_FILE if test 133 -ne `wc -c <'samples/sintest.rpl'`; then echo shar: \"'samples/sintest.rpl'\" unpacked with wrong size! fi # end of 'samples/sintest.rpl' fi if test -f 'samples/test.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/test.rpl'\" else echo shar: Extracting \"'samples/test.rpl'\" \(108 characters\) sed "s/^X//" >'samples/test.rpl' <<'END_OF_FILE' X"TRON X"BENCH1.RPL" LOAD BENCH1 X'BENCH1' PURGE X"BENCH2.RPL" LOAD BENCH2 X'BENCH2' PURGE X'FIB' PURGE CLEAR OFF END_OF_FILE if test 108 -ne `wc -c <'samples/test.rpl'`; then echo shar: \"'samples/test.rpl'\" unpacked with wrong size! fi # end of 'samples/test.rpl' fi if test -f 'samples/test1.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/test1.rpl'\" else echo shar: Extracting \"'samples/test1.rpl'\" \(29 characters\) sed "s/^X//" >'samples/test1.rpl' <<'END_OF_FILE' X"BENCH1.RPL" LOAD BENCH1 OFF END_OF_FILE if test 29 -ne `wc -c <'samples/test1.rpl'`; then echo shar: \"'samples/test1.rpl'\" unpacked with wrong size! fi # end of 'samples/test1.rpl' fi if test -f 'samples/test2.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/test2.rpl'\" else echo shar: Extracting \"'samples/test2.rpl'\" \(56 characters\) sed "s/^X//" >'samples/test2.rpl' <<'END_OF_FILE' X"BENCH2.RPL" LOAD BENCH2 X'BENCH2' PURGE X'FIB' PURGE OFF END_OF_FILE if test 56 -ne `wc -c <'samples/test2.rpl'`; then echo shar: \"'samples/test2.rpl'\" unpacked with wrong size! fi # end of 'samples/test2.rpl' fi if test -f 'samples/tf.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/tf.rpl'\" else echo shar: Extracting \"'samples/tf.rpl'\" \(49 characters\) sed "s/^X//" >'samples/tf.rpl' <<'END_OF_FILE' X<< IF THEN "TRUE" ELSE "FALSE" ENDIF >> 'TF' STO END_OF_FILE if test 49 -ne `wc -c <'samples/tf.rpl'`; then echo shar: \"'samples/tf.rpl'\" unpacked with wrong size! fi # end of 'samples/tf.rpl' fi if test -f 'samples/wuerfel.rpl' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'samples/wuerfel.rpl'\" else echo shar: Extracting \"'samples/wuerfel.rpl'\" \(98 characters\) sed "s/^X//" >'samples/wuerfel.rpl' <<'END_OF_FILE' X<< X 0 'S STO X 1 1000 FOR I X I 6 / 5 6 / I 1 - ^ * X S + DUP 'S STO X PRINT X NEXT X>> 'WUERFEL STO END_OF_FILE if test 98 -ne `wc -c <'samples/wuerfel.rpl'`; then echo shar: \"'samples/wuerfel.rpl'\" unpacked with wrong size! fi # end of 'samples/wuerfel.rpl' fi echo shar: End of shell archive. exit 0