Check-ins Using Download Hex

Artifact 0fa5c5ceb32143b27842386b2281be259d6ba71917d6923fc8481f9529bd389b:

%. MATHLIB.RED - Some useful mathematical functions for PSL
%
% Most of these routines not very heavily tested. 
% Contributions from Galway, Griss, Irish, Morrison, and others.
%
%  MATHLIB.RED, 16-Dec-82 21:56:52, Edit by GALWAY
%   Various fixes and enhancements too numerous for me to remember.
%   Includes fixes in SQRT function, modifications of RANDOM and other
%   functions to bring them more in line with Common Lisp, addition of MOD
%   and FLOOR.
%  <PSL.UTIL>MATHLIB.RED.13, 13-Sep-82 08:49:52, Edit by BENSON
%  Bug in EXP, changed 2**N to 2.0**N
%  <PSL.UTIL>MATHLIB.RED.12,  2-Sep-82 09:22:19, Edit by BENSON
%  Changed all calls in REDERR to calls on STDERROR
%  <PSL.UTIL>MATHLIB.RED.2, 17-Jan-82 15:48:21, Edit by GRISS
%  changed for PSL

% Should these names be changed so that they all begin with an F or some
% other distinguishing mark?  Are they in conflict with anything?  Or should
% we wait until we have packages?

% Consider using Sasaki's BigFloat package -- it has all this and more, to
% arbitrary precision.  The only drawback is speed.

%***************** Constants declared as NewNam's ****************************

% We can't use these long ones in Lisp1.6 'cause the reader craps out (and
% it would truncate instead of round, anyway).  These are here for reference
% for implementation on other machines.
% put('NumberPi,'NewNam,3.14159265358979324);
% put('NumberPi!/2,'NewNam,1.57079632679489662);
% put('NumberPi!/4,'NewNam,0.785398163397448310);

BothTimes <<
put('Number2Pi,'NewNam,6.2831853);
put('NumberPi,'NewNam,3.1415927);
put('NumberPi!/2,'NewNam,1.5707963);
put('NumberPi!/4,'NewNam,0.78539816);
put('Number3Pi!/4,'NewNam,2.3561945);
put('Number!-2Pi,'Newnam,-6.2831853);
put('Number!-Pi,'NewNam,-3.1415927);
put('Number!-Pi!/2,'NewNam,-1.5707963);
put('Number!-Pi!/4,'NewNam,-0.78539816);

put('SqrtTolerance,'NewNam,0.0000001);
put('NumberE, 'NewNam, 2.718281828);
put('NumberInverseE, 'NewNam, 0.36787944);     % 1/e
put('NaturalLog2,'NewNam,0.69314718);
put('NaturalLog10,'NewNam,2.3025851);
put('TrigPrecisionLimit,'NewNam,80);

>>;
%********************* Basic functions ***************************************

lisp procedure mod(M,N);
% Return M modulo N.  Unlike remainder function--it returns positive result
% in range 0..N-1, even if M is negative.  (Needs more work for case of
% negative N.)
begin scalar result;
    result := remainder(M,N);
    if result >= 0 then
        return result;
    % else
    return
        N + result;
end;

lisp procedure Floor X;
% Returns the largest integer less than or equal to X.  (I.e. the "greatest
% integer" function.)
if fixp X then
  X
else begin scalar N;
  N := fix X;
  % Note the trickiness to compensate for fact that (unlike APL's "FLOOR"
  % function) FIX truncates towards zero.
  return if X = float N then N else if X>=0 then N else N-1;
end;

lisp procedure Ceiling X;
% Returns the smallest integer greater than or equal to X.
if fixp X then
  X
else begin scalar N;
  N := fix X;
  % Note the trickiness to compensate for fact that (unlike APL's "FLOOR"
  % function) FIX truncates towards zero.
  return if X = float N then N else if X>0 then N+1 else N;
end;

lisp procedure Round X;
% Rounds to the closest integer.
% Kind of sloppy -- it's biased when the digit causing rounding is a five,
% it's a bit weird with negative arguments, round(-2.5)= -2.
if fixp X then
  X
else 
  floor(X+0.5);

%***************** Trigonometric Functions ***********************************

% Trig functions are all in radians.  The following few functions may be used
% to convert to/from degrees, or degrees/minutes/seconds.

lisp procedure DegreesToRadians x;
x*0.017453292; % 2*pi/360

lisp procedure RadiansToDegrees x;
  x*57.29578;    % 360/(2*pi)

lisp procedure RadiansToDMS x;
% Converts radians to a list of degrees, minutes, and seconds (rounded, not
% truncated, to the nearest integer).
begin scalar Degs,Mins;
  x := RadiansToDegrees x;
  Degs := fix x;
  x := 60*(x-Degs);
  Mins := fix x;
  return list(Degs,Mins, Round(60*(x-Mins)))
end;

lisp procedure DMStoRadians(Degs,Mins,Sex);
% Converts degrees, minutes, seconds to radians.
% DegreesToRadians(Degs+Mins/60.0+Sex/3600.0)
DegreesToRadians(Degs+Mins*0.016666667+Sex*0.00027777778);

lisp procedure sin x;
% Accurate to about 6 decimal places, so long as the argument is 
% of commensurate precision.  This will, of course, NOT be true for
% large arguments, since they will be coming in with small precision.
begin scalar neg;
  if minusp x then <<
    neg := T;
    x := - x >>;
  if x > TrigPrecisionLimit then
    LPriM "Possible loss of precision in computation of SIN";
  if x > NumberPi then
    x := x-Number2Pi*fix((x+NumberPi)/Number2Pi);
  if minusp x then <<
    neg := not neg;
    x :=  -x >>;
  if x > NumberPi!/2 then
    x := NumberPi-x;
  return if neg then -ScaledSine x else ScaledSine x
end;

lisp procedure ScaledSine x;
% assumes its argument is scaled to between 0 and pi/2.
begin scalar xsqrd;
  xsqrd := x*x;
  return x*(1+xsqrd*(-0.16666667+xsqrd*(0.0083333315+xsqrd*(-0.0001984090+
              xsqrd*(0.0000027526-xsqrd*0.0000000239)))))
end;

lisp procedure cos x;
% Accurate to about 6 decimal places, so long as the argument is 
% of commensurate precision.  This will, of course, NOT be true for
% large arguments, since they will be coming in with small precision.
<< if minusp x then
     x := - x;
   if x > TrigPrecisionLimit then
     LPriM "Possible loss of precision in computation of COS";
   if x > NumberPi then
     x := x-Number2Pi*fix((x+NumberPi)/Number2Pi);
   if minusp x then
     x := - x;
   if x > NumberPi!/2 then
     -ScaledCosine(NumberPi-x)
   else
     ScaledCosine x >>;

lisp procedure ScaledCosine x;
% Expects its argument to be between 0 and pi/2.
begin scalar xsqrd;
  xsqrd := x*x;
  return 1+xsqrd*(-0.5+xsqrd*(0.041666642+xsqrd*(-0.0013888397+
              xsqrd*(0.0000247609-xsqrd*0.0000002605))))
end;

lisp procedure tan x;
% Accurate to about 6 decimal places, so long as the argument is 
% of commensurate precision.  This will, of course, NOT be true for
% large arguments, since they will be coming in with small precision.
begin scalar neg;
  if minusp x then <<
    neg := T;
    x := - x >>;
  if x > TrigPrecisionLimit then
    LPriM "Possible loss of precision in computation of TAN";
  if x > NumberPi!/2 then
    x := x-NumberPi*fix((x+NumberPi!/2)/NumberPi);
  if minusp x then <<
    neg := not neg;
    x := - x >>;
  if x < NumberPi!/4 then
    x := ScaledTangent x
  else
    x := ScaledCotangent(-(x-numberpi!/2));
  return if neg then -x else x
end;

lisp procedure cot x;
% Accurate to about 6 decimal places, so long as the argument is 
% of commensurate precision.  This will, of course, NOT be true for
% large arguments, since they will be coming in with small precision.
begin scalar neg;
  if minusp x then <<
    neg := T;
    x := - x >>;
  if x > NumberPi!/2 then
    x := x-NumberPi*fix((x+NumberPi!/2)/NumberPi);
  if x > TrigPrecisionLimit then
    LPriM "Possible loss of precision in computation of COT";
  if minusp x then <<
    neg := not neg;
    x := - x >>;
  if x < NumberPi!/4 then
    x := ScaledCotangent x
  else
    x := ScaledTangent(-(x-numberpi!/2));
  return if neg then -x else x
end;

lisp procedure ScaledTangent x;
% Expects its argument to be between 0 and pi/4.
begin scalar xsqrd;
  xsqrd := x*x;
  return x*(1.0+xsqrd*(0.3333314+xsqrd*(0.1333924+xsqrd*(0.05337406 +
           xsqrd*(0.024565089+xsqrd*(0.002900525+xsqrd*0.0095168091))))))
end;

lisp procedure ScaledCotangent x;
% Expects its argument to be between 0 and pi/4.
begin scalar xsqrd;
  xsqrd := x*x;
  return (1.0-xsqrd*(0.33333334+xsqrd*(0.022222029+xsqrd*(0.0021177168 +
           xsqrd*(0.0002078504+xsqrd*0.0000262619)))))/x
end;

lisp procedure sec x;
1.0/cos x;

lisp procedure csc x;
1.0/sin x;

lisp procedure sinD x;
sin DegreesToRadians x;

lisp procedure cosD x;
cos DegreesToRadians x;

lisp procedure tanD x;
tan DegreesToRadians x;

lisp procedure cotD x;
cot DegreesToRadians x;

lisp procedure secD x;
sec DegreesToRadians x;

lisp procedure cscD x;
csc DegreesToRadians x;

lisp procedure asin x;
begin scalar neg;
  if minusp x then <<
    neg := T;
    x := -x >>;
  if x > 1.0 then
    stderror list("Argument to ASIN too large:",x);
  return if neg then CheckedArcCosine x - NumberPi!/2 
		else NumberPi!/2 - CheckedArcCosine x
end;

lisp procedure acos x;
begin scalar neg;
  if minusp x then <<
    neg := T;
    x := -x >>;
  if x > 1.0 then
    stderror list("Argument to ACOS too large:",x);
  return if neg then NumberPi - CheckedArcCosine x
		else CheckedArcCosine x
end;

lisp procedure CheckedArcCosine x;
% Return cosine of a "checked number", assumes its argument is in the range
% 0 <= x <= 1.
sqrt(1.0-x)*(1.5707963+x*(-0.2145988+x*(0.088978987+x*(-0.050174305+
        x*(0.030891881+x*(-0.017088126+x*(0.0066700901-x*(0.0012624911))))))));

lisp procedure atan x;
if minusp x then
  if x < -1.0 then
    Number!-Pi!/2 + CheckedArcTangent(-1.0/x)
  else
    -CheckedArcTangent(-x)
else
  if x > 1.0 then
    NumberPi!/2 - CheckedArcTangent(1.0/x)
  else
    CheckedArcTangent x;

lisp procedure acot x;
if minusp x then
  if x < -1.0 then
    -CheckedArcTangent(-1.0/x)
  else
    Number!-Pi!/2 + CheckedArcTangent(-x)
else
  if x > 1.0 then
   CheckedArcTangent(1.0/x)
  else
    NumberPi!/2 - CheckedArcTangent x;

lisp procedure CheckedArcTangent x;
begin scalar xsqrd;
  xsqrd := x*x;
  return x*(1+xsqrd*(-0.33333145+xsqrd*(0.19993551+xsqrd*(-0.14208899+
             xsqrd*(0.10656264+xsqrd*(-0.07528964+xsqrd*(0.042909614+
	     xsqrd*(-0.016165737+xsqrd*0.0028662257))))))))
end;

lisp procedure asec x;
acos(1.0/x);

lisp procedure acsc x;
asin(1.0/x);

lisp procedure asinD x;
RadiansToDegrees asin x;

lisp procedure acosD x;
RadiansToDegrees acos x;

lisp procedure atanD x;
RadiansToDegrees atan x;

lisp procedure acotD x;
RadiansToDegrees acot x;

lisp procedure asecD x;
RadiansToDegrees asec x;

lisp procedure acscD x;
RadiansToDegrees acsc x;

%****************** Roots and such *******************************************

lisp procedure sqrt N;
% Simple Newton-Raphson floating point square root calculator.
% Not waranted against truncation errors, etc.
begin integer answer,scale;
  N:=FLOAT N;
  if N < 0.0 then stderror list("SQRT given negative argument:",N);
  if zerop N then
    return N;
  % Scale argument to within 1e-10 to 1e+10;
  scale := 0;
  while N > 1.0E10 do
  <<
    scale := scale + 1;
    N := N * 1.0E-10 >>;
  while N < 1.0E-10 do
  <<
    scale := scale - 1;
    N := N * 1.0E10 >>;
  answer := if N>2.0 then (N+1)/2
         else if N<0.5 then 2/(N+1)
         else N;

  % Here's the heart of the algorithm.
  while abs(answer**2/N - 1.0) > SqrtTolerance do
    answer := 0.5*(answer+N/answer);
  return answer * 10.0**(5*scale)
end;

%******************** Logs and Exponentials **********************************

lisp procedure exp x;
% Returns the exponential (ie, e**x) of its floatnum argument as
% a flonum. The argument is scaled to
% the interval -ln  2 to  0, and a  Taylor series  expansion
% used (formula 4.2.45 on page 71 of Abramowitz and  Stegun,
% "Handbook of Mathematical  Functions").
begin scalar N;
  N := ceiling(x / NaturalLog2);
  x := N * NaturalLog2 - x;
  return 2.0**N * (1.0+x*(-0.9999999995+x*(0.4999999206+x*(-0.1666653019+
        x*(0.0416573475+x*(-0.0083013598+x*(0.0013298820+
        x*(-0.0001413161))))))))
end;


lisp procedure log x;
% See Abramowitz and Stegun, page 69.

 if x <= 0.0 then
   stderror list("LOG given non-positive argument:",x)
 else if x < 1.0 then
   -log(1.0/x)
 else
 % Find natural log of x > 1;
 begin scalar nextx, ipart;      % ipart is the "integer part" of the
                                 % logarithm.
   ipart := 0;

   % Keep multiplying by 1/e until x is small enough, may want to be more
   % "efficient" if we ever use really big numbers.
   while (nextx := NumberInverseE * x) > 1.0 do
   <<
       x := nextx;
       ipart := ipart + 1;
   >>;

   return
       ipart +
       if x < 2.0 then
         CheckedLogarithm x
       else
         2.0 * CheckedLogarithm(sqrt(x));
 end;
 
lisp procedure CheckedLogarithm x;
% Should have 1 <= x <= 2.  (i.e. x = 1+y  0 <= y <= 1)
<< x := x-1.0;
    x*(0.99999642+x*(-0.49987412+x*(0.33179903+x*(-0.24073381+x*(0.16765407+
         x*(-0.09532939+x*(0.036088494-x*0.0064535442))))))) >>;

lisp procedure log2 x;
log x / NaturalLog2;

lisp procedure log10 x;
log x / NaturalLog10;

%********************* Random Number Generator *******************************

% The declarations below  constitute a linear,  congruential
% random number generator (see  Knuth, "The Art of  Computer
% Programming: Volume 2: Seminumerical Algorithms", pp9-24).
% With the given  constants it  has a period  of 392931  and
% potency  6.    To   have  deterministic   behaviour,   set
% RANDOMSEED.
%
% Constants are:        6   2
%    modulus: 392931 = 3 * 7 * 11
%    multiplier: 232 = 3 * 7 * 11 + 1
%    increment: 65537 is prime
%
% Would benefit from being recoded in SysLisp, when full word integers should
% be used with "automatic" modular arithmetic (see Knuth).  Perhaps we should
% have a longer period version?
% By E. Benson, W. Galway and M. Griss

fluid '(RandomSeed RandomModulus);

RandomModulus := 392931;
RandomSeed := remainder(time(),RandomModulus);

lisp procedure next!-random!-number;
% Returns a pseudo-random number between 0 and RandomModulus-1 (inclusive).
RandomSeed := remainder(232*RandomSeed + 65537, RandomModulus);

lisp procedure Random(N);
% Return a pseudo-random number uniformly selected from the range 0..N-1.
% NOTE that this used to be called RandomMod(N).  Needs to be made more
% compatible with Common LISP's random?
    fix( (float(N) * next!-random!-number()) / RandomModulus);

procedure FACTORIAL N;   % Simple factorial
 Begin scalar M;
    M:=1;
    for i:=1:N do M:=M*I;
    Return M;
 end;


% Some functions from ALPHA_1 users

lisp procedure Atan2D( Y, X );
    RadiansToDegrees Atan2( Y, X );

lisp procedure Atan2( Y, X );
<<
    X := float X; Y := Float Y;

    if X = 0.0 then			% Y axis.
	if  Y >= 0.0  then  NumberPI!/2  else  NumberPi + NumberPI!/2

    else if X >= 0.0 and Y >= 0.0 then	% First quadrant.
	Atan( Y / X )

    else if X < 0.0 and Y >= 0.0 then	% Second quadrant.
	NumberPI - Atan( Y / -X )

    else if X < 0.0 and Y < 0.0 then	% Third quadrant.
	NumberPI + Atan( Y / X )

    else				% Fourth quadrant.
	Number2Pi - Atan( -Y / X )
>>;

lisp procedure TransferSign( S, Val );
% Transfers the sign of S to Val by returning abs(Val) if S >= 0,
% otherwise -abs(Val).
    if S >= 0 then abs(Val) else -abs(Val);

lisp procedure DMStoDegrees(Degs,Mins,Sex);
% Converts degrees, minutes, seconds to degrees
% Degs+Mins/60.0+Sex/3600.0
    Degs+Mins*0.016666667+Sex*0.00027777778;

lisp procedure DegreesToDMS x;
% Converts degrees to a list of degrees, minutes, and seconds (all integers,
% rounded, not truncated).
begin scalar Degs,Mins;
  Degs := fix x;
  x := 60*(x-Degs);
  Mins := fix x;
  return list(Degs,Mins, round(60*(x-Mins)))
end;

end;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]