home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Library Routines / Text Examples / Miscellaneous Routines < prev    next >
Encoding:
Text File  |  1992-10-06  |  11.4 KB  |  342 lines  |  [TEXT/NCII]
  1. ####################################
  2. #  Miscellaneous Routines
  3. #
  4. # These routines don't use any xCOD's, just NCII's built-in
  5. #  programming language.  
  6. #
  7. #   Truncate(x) = Int(x*(1+1e-19) - .5) # Returns integer <= x.;
  8. #   FractionalPart(x) = x - Truncate(x)# Returns fractional part of x>0.;
  9. #   Round(x, DecimalDigits) =  int(x*10^DecimalDigits)/10^DecimalDigits # rounds off x.;
  10. #   Mod(a,b) = a - b*Truncate(a/b) # Returns a mod b.;
  11. #   function IsEven(n) # Returns 1 if n is even, 0 if odd.
  12. # 
  13. #   Length(V) = sqrt(V•V) # length of a vector. ;
  14. #  AngleBetween(V,W) = acos( V•W /length(V) /length(W) ) # angle from vector V to W, in radians.;
  15. #   function Cross3D(X,Y)  # Returns the 3D cross product of two 3D vectors x,y.
  16. #   Hermitian(theMat) = Transpose(theMat†);
  17. #   Commutator(a,b) = a•b - b•a;
  18. #   function Trace(theMatrix) # Returns the sum of the diagonal elements of a matrix.
  19. #
  20. #   Differentiate(f,x,dx) = [ f(x+dx/2) - f(x-dx/2) ]/dx  # Returns df/dx at x[1…N]
  21. # 
  22. #   Mean(V) = sum(V)/size(V) # Returns the mean of an array V;
  23. #   Variance(V) = Mean( | V-Mean(V) |^2 ) # Returns the sigma^2 of an array V;
  24. #   StandardDeviation(V) = sqrt(Mean(|V|^2) - |Mean(V)|^2) #Returns standard deviation of an array V;
  25. #   Gaussian(x,x0,sigma)=1/sqrt(2πsigma^2) exp(-0.5{(x-x0)/sigma}^2);
  26. # 
  27. #   Log(x) = ln(x)/ln(10) # log base 10 of x.;
  28. #    function LogGrid(min, max, N) # Makes an N point logarithmic grid.
  29. #
  30. #   Sigma(x) = (1 + Sign(x+1e-16))/2 # is 1 for x=>0, 0 for x<0.  
  31. #   aTan2(s,c) = if c>=0 then atan(s/c) else sign(s+1e-18)*π+atan(s/c)# is -π<theta<=π from s=sin(theta), c=cos(theta);
  32. #   FixRoundOff(x) = (x+1)-1 # Returns x rounded off to about 1e-19.;
  33. #  program BeepIfComplex(z) # Beeps if z is complex ;
  34. #  program ConvertToComplex(x) # Converts x to a complex quantity.
  35. #   SetNegativeToZero(A) =  (Sign(A+1e-16) + 1)/2 * A     #Zeros the elements of an array A which are < 0.
  36. # 
  37. # 
  38. #  This text file explains and gives examples of
  39. #  some of the routines in the file All Library Routines, which
  40. #  should be Opened before trying any of these examples.
  41. ####################################
  42.  
  43.  
  44. ##############
  45. #   Integer Arithmetic
  46.  
  47. # Truncate(x) truncates the fractional part of a number.
  48. # Int (a built-in) rounds to the nearest even integer, 
  49. # while 1e-19 is near the smallest SANE number eps such that
  50. # (1+eps)-eps <> 0.
  51.   Truncate(x) = Int(x*(1+1e-19) - .5) # Returns integer <= x.;
  52.  
  53.   FractionalPart(x) = x - Truncate(x)# Returns fractional part of x>0.;
  54.  
  55.   Round(x, DecimalDigits) =  int(x*10^DecimalDigits)/10^DecimalDigits # rounds off x.;
  56.  
  57. # Examples:
  58. Truncate(2.9872)
  59.   2 
  60. FractionalPart(2.9872)
  61.   0.9872 
  62. Round(2.9872, 3)
  63.   2.987 
  64. Round(2.9872, 2)
  65.   2.99 
  66.  
  67.  
  68. # A simpler way to implement Trunc is to use the fact
  69. # that array indeces are truncated. Therefore, if you know
  70. # the range of integers of interest, just define an integer
  71. # array in that range and subscript it with the value you want
  72. # to truncate.  For example, if you only use integers 1…100,
  73.   Trunc[1…100] = 1…100;
  74. # Then Trunc(1.9) = Trunc[1] = 1,  Trunc(2.3) = 2, etc.
  75.  
  76.  
  77. # Mod(a,b) is "a mod b", i.e. the remainder of a/b.
  78.   Mod(a,b) = a - b*Truncate(a/b) # Returns a mod b.;
  79.  
  80.  
  81. # Also using Int, we can test for even/odd.
  82. # IsEven returns 1 if n is even, 0 if odd. 
  83.   function IsEven(n) # Returns 1 if n is even, 0 if odd.
  84. .   begin
  85. .   if n=int(n+0.5) then
  86. .     IsEven=1
  87. .   else
  88. .     IsEven=0
  89. .   end
  90. .   
  91.  
  92.  
  93. ####################################
  94. #   Matrix  and Vector
  95. # Also see the xCOD files for Inverse, Determinant, and Eigens., 
  96. # which have many more matrix routines, and the sample file
  97. # Vetor and Matrix.
  98.  
  99. # Find the length of a one-dimensional vector. ("•" is opt-8, or use the Shortcuts menu.)
  100.  Length(V) = sqrt(V•V) # length of a vector. ;
  101.  
  102. # The angle between two vectors.
  103.  AngleBetween(V,W) = acos( V•W /length(V) /length(W) ) # angle from vector V to W, in radians.;
  104.  
  105. # The three-dimensional cross product.
  106. function Cross3D(X,Y)  # Returns the 3D cross product of two 3D vectors x,y.
  107. .   var Z
  108. .   #  Input:   X,Y[1…3] 3-D vectors
  109. .   #  Output: Cross3D = X cross Y
  110. .   begin
  111. .      if size(x) <> 3 then
  112. .        Cross3D = " ERROR: X must be a 3D vector for Cross3D"
  113. .      else if size(y) <> 3 then
  114. .         Cross3D = " ERROR: Y must be a 3D vector for Cross3D"
  115. .      else
  116. .      begin
  117. .        z[1] = x[2] y[3] - x[3] y[2];
  118. .        z[2] = x[3] y[1] - x[1] y[3];
  119. .        z[3] = x[1] y[2] - x[2] y[1];
  120. .        Cross3D = z
  121. .      end
  122. .    end
  123. .    
  124.  
  125. # The Hermitian is the conjugate transpose of a matrix.
  126.   Hermitian(theMat) = Transpose(theMat†);
  127.  
  128. # This should only be applied to square matrices.
  129.   Commutator(a,b) = a•b - b•a;
  130.  
  131.  function Trace(theMatrix) # Returns the sum of the diagonal elements of a matrix.
  132. .  var d,Answer
  133. .  begin
  134. .      d = dimension(theMatrix)
  135. .      Answer = 0
  136. .      if size(d)<>2 then 
  137. .         Answer = " ERROR: argument to Trace must be a matrix."
  138. .      else if d[1]<>d[2] then
  139. .          Answer = " ERROR: matrix for Trace must be square."
  140. .      else
  141. .         for j=1…d[1] do Answer = Answer + theMatrix[j,j]
  142. .     Trace = Answer
  143. .  end
  144. . 
  145.  
  146.  
  147. ####################################
  148. #   Calculus
  149.  
  150. # All of these are crude but fairly simple, and use only
  151. # NCII's built in programming language. 
  152. #
  153. #  Since the xCOD based routines in the Discrete Calculus file
  154. #  do most of these same task quicker and more accurately, 
  155. #  all but the first of these have been commented out, leaving
  156. #  the rest only for examples of the NCII programming language.
  157.  
  158.  
  159. # The numerical derivative of a function f(x), returned
  160. #  at the values of the array x. dx is any "small" number.
  161.   Differentiate(f,x,dx) = [ f(x+dx/2) - f(x-dx/2) ]/dx  # Returns df/dx at x[1…N]
  162. # For example, 
  163.   f(x) = sin(x);  x =0…2π@.01;  dsin = DifferentiateFunc(f,x,1e-6)
  164.  
  165. # # This one returns an array. x should be monotomic with y=f(x).
  166. # function DifferentiateArray(y,x)
  167. # .   var y0, y1, n0, n, DA,DA0,DA1, ans
  168. # .   begin
  169. # .     n = size(y);  n0 = FirstIndex(y); nN=n0+n-1;
  170. # .     y0 = y[n0…(nN-1)];  y1 = y[(n0+1)…nN];
  171. # .     DA = (y1-y0)/(x[2]-x[1]);  # derivs at in-between points.
  172. # .     DA0 = DA[1…nN-2]; DA1=DA[2…nN-1];    
  173. # .     ans = y;
  174. # .     ans[(n0+1)…(nN-1)] = (DA0+DA1)/2;  # interpolate middle pts
  175. # .     ans[n0] = (3 DA0[1] - DA0[2])/2;        # extrapolate start
  176. # .     ans[nN] = (3 DA1[nN-2] - DA1[nN-3])/2;  # extrapolate end.
  177. # .     DifferentiateArray = ans
  178. # .   end
  179. # .   
  180.  
  181. # # This assumes that y and x are arrays with y=f(x) and that
  182. # # x is monotomic. Simple trapezoidal integration is done.
  183. # function IntegrateArray(y,x) = sum(y) *(x[2]-x[1]);
  184. # .   var n0, nN, ans
  185. # .   begin
  186. # .     n0 = FirstIndex(y);
  187. # .     nN = n0 + Size(y)-1;
  188. # .     ans = sum(y) - { y[n0]+y[nN] }/2
  189. # .     ans = ans*{ x[2] - x[1] }
  190. # .     IntegrateArray = ans
  191. # .   end
  192. # .   
  193.  
  194. # function Integrate(F,from,to);
  195. # .   var x,dx,y
  196. # .   begin
  197. # .     dx = (to-from)/300;   # we'll do the integral with 300 points.
  198. # .     x = from…to@dx;
  199. # .     y = F(x);
  200. # .     IntegrateFunc = IntegrateArray(y,x)
  201. # .   end
  202. # .   
  203.  
  204. # # Finally an indefinite integral. Both y and x are arrays,
  205. # # y=f(x) and x is monotomic.
  206. # # Because of the FOR loop, this one is slow.
  207. # function IndefiniteIntegral(y,x)
  208. # .   var n0,j, ans
  209. # .   begin
  210. # .     n0 = FirstIndex(x)
  211. # .     ans = y # save space for answer
  212. # .     ans[n0] = 0  #  integral with same limits
  213. # .     for j=(n0+1) to n0+size(x)-1 do
  214. # .       ans[j] = IntegrateArray(y[n0…j],x[n0…j])
  215. # .     IndefiniteIntegral = ans
  216. # .   end
  217. # .   
  218. # # Example: 
  219.  #  x = 0…2π@.01; y = sin(x);  IntegralY = IndefiniteIntegral(y,x)
  220. # # Then plotting [x,y] ; [x,IntegralY] will show 
  221. # #  sin(x) and (cos(x)-1).   (The -1 happens because 
  222. # #  the integration is only to within a constant which makes the 
  223. # #  integral 0 at the smallest value of x.)
  224.  
  225.  
  226.  
  227. ####################################
  228. #   Statistics
  229. # Also see the files for Sorting and Mean, as well as
  230. # Special Functions (which has the Error Function and Chi Square probability).
  231.  
  232. # Standard measures of an array of values.
  233.   Mean(V) = sum(V)/size(V) # Returns the mean of an array V;
  234.   Variance(V) = Mean( | V-Mean(V) |^2 ) # Returns the sigma^2 of an array V;
  235.   StandardDeviation(V) = sqrt(Mean(|V|^2) - |Mean(V)|^2) # Returns standard deviation of an array V;
  236.  
  237. # A "bell-curve" with mean x0 and Std. Dev. sigma.
  238.   Gaussian(x,x0,sigma)=1/sqrt(2πsigma^2) exp(-0.5{(x-x0)/sigma}^2);
  239. #
  240. #    Example:
  241. dx=0.01;x = -4…4@dx; y = Gaussian(x,0,1);    # Now plot [x,y] .
  242. #
  243. sum(y) dx        # normalization of gaussian
  244.   0.99993799 
  245. #
  246. sum(y x^2) dx
  247.   0.9988873     # variance of x.
  248.  
  249.  
  250. ####################################
  251. #   Log Utilities
  252.  
  253. # In general, log A base B is ln(a)/ln(b). 
  254. # ("e" is not defined by default, but "e=exp(1)" will define it.)
  255. Log(x) = ln(x)/ln(10) # log base 10 of x.;
  256.  
  257. # This is useful for making an even grid when doing
  258. # log plots.
  259. #
  260.  function LogGrid(min, max, N) # Makes an N point logarithmic grid.
  261. .   var x,dx, logMax, logMin
  262. .    # Input:    
  263. .    #              min = lowest grid value
  264. .    #              max = highest grid value
  265. .    #              N = number of points
  266. .    # Output:  
  267. .    #             LogGrid = [min, min*Z, min*Z*Z, ..., max]
  268. .    #                    (evenly spaced points in a log plot).
  269. .   begin
  270. .     logMin = log(min); logMax = log(max);
  271. .     dx = (logMax - logMin)/(N-1);
  272. .     x = 10^( logMin…logMax @ dx)
  273. .     LogGrid = x
  274. .   end
  275. .   
  276.  
  277.  
  278. ####################################
  279. #   Other
  280.  
  281. # Sigma is sometimes called the "step" function.
  282.   Sigma(x) = (1 + Sign(x+1e-16))/2 # is 1 for x=>0, 0 for x<0. ;
  283. # for example, 
  284. Table( -2:2, Sigma(-2:2) )
  285. # 
  286. # - - - - - - - - - - - - - - - 
  287.    -2                0
  288.    -1                0
  289.    0                  1
  290.    1                  1
  291.    2                  1
  292.  
  293.  
  294. # aTan2 is takes 2 arguments, sin(theta), cos(theta), and returns -π< theta < π.
  295. # This uses an "if" as part of a one line function. A full-fledged function might be
  296. # clearer, but not as much fun.
  297. aTan2(s,c) = if c>=0 then atan(s/c) else sign(s+1e-18)*π+atan(s/c)# is -π<theta<=π from s=sin(theta), c=cos(theta);
  298.  
  299.  
  300.  
  301. # Return 1 if an object exists, 0 if it doesn't.
  302. function Exists(Obj) # Return 1 if Obj exists, 0 if not.
  303. .   begin
  304. .     if Obj=Obj then Exists=1 else Exists=0
  305. .   end
  306. .   
  307.  
  308. # Sometimes useful for setting small numbers due to
  309. # round off errors to zero.  This works for numbers,
  310. # arrays, or matrices. However, if x is complex
  311. # this will only affect the real part; (i+x)-i will do the
  312. # same for the complex part. (And will make it complex
  313. # if it isn't already.)
  314.   FixRoundOff(x) = (x+1)-1 # Returns x rounded off to about 1e-19.;
  315.  
  316. # Usually real to complex conversions are handled 
  317. # automatically; however, it is possible to tell how and object
  318. # is stored.  
  319.  program BeepIfComplex(z) # Beeps if z is complex ;
  320. .   begin
  321. .      if Size(z) <> SizeAsIfReal(z) then beep
  322. .   end
  323. .   
  324.  
  325. # This makes x complex, which means it will 
  326. # take up twice as much space.  Usually this conversion
  327. # is done automatically when needed.  Under some circumstances
  328. # (such as "if x=1") it will be converted back to real 
  329. # automatically if the complex part is zero.
  330.  program ConvertToComplex(x) # Converts x to a complex quantity.
  331. .   begin
  332. .     x  = (x+i)-i
  333. .   end
  334. .   
  335.  
  336. # The Sign function may be used to do the work of an IF on arrays.  
  337. #  Sign(x) is 1 for x>0, 0 for x=0, and -1 for x<0.  
  338. #  On array, it does Sign on each element of the array to produce a new array.
  339. #  This example sets all the elements of an array which are less than zero to zero.
  340.  SetNegativeToZero(A) =  (Sign(A+1e-16) + 1)/2 * A     #Zeros the elements of an array A which are < 0.
  341.  
  342.